Theoretical Biology and Medical Modelling

Background: The epidermal growth factor receptor (EGFR) is frequently overexpressed in many cancers, including non-small cell lung cancer (NSCLC). In silico modeling is considered to be an increasingly promising tool to add useful insights into the dynamics of the EGFR signal transduction pathway. However, most of the previous modeling work focused on the molecular or the cellular level only , neglecting the crucial feedback between these scales as well as the interaction with the heterogeneous biochemical microenvironment. Results: We developed a multiscale model for investigating expansion dynamics of NSCLC within a two-dimensional in silico microenvironment. At the molecular level, a specific EGFR-ERK intracellular signal transduction pathway was implemented. Dynamical alterations of these molecules were used to trigger phenotypic changes at the cellular level. Examining the relationship between extrinsic ligand concentrations, intrinsic molecular profiles and microscopic patterns, the results confirmed that increasing the amount of available growth factor leads to a spatially more aggressive cancer system. Moreover, for the cell closest to nutrient abundance, a phase-transition emerges where a minimal increase in extrinsic ligand abolishes the proliferative phenotype altogether. Conclusion: Our in silico results indicate that in NSCLC, in the presence of a strong extrinsic chemotactic stimulus (and depending on the cell


Background
Non-small cell lung cancer (NSCLC) remains at the top of the list of cancer-related deaths in the United States [1]. The epidermal growth factor receptor (EGFR) is frequently overexpressed in NSCLC [2,3]. Binding of epidermal growth factor (EGF) or transforming growth factor alpha (TGFα) to the extracellular domain of EGFR produces a from the spread of primary tumor cells into the surrounding tissue [9], quantitative measurements of the relationship between the level of the growth factors and the resulting tumor expansion is crucial -all the more so, since EGFR has emerged as an attractive therapeutic target for patients with advanced NSCLC [10].
A number of EGFR-related intracellular signal transduction pathways have been studied [11][12][13][14][15][16], including NSCLC [17], and corresponding computational models at the molecular-level have been developed. These quantitative works mainly focused on signal-response relationships between the binding of EGF to EGFR and the activation of downstream proteins in the signaling cascade. With these in silico approaches, experimentally testable hypotheses can be made on signaling events controlling divergent cellular responses such as cell proliferation, differentiation, or apoptosis [18,19]. However, most signaling works did not yet consider the cellular level (see [20,21] for a review), and, conversely, only a few recent EGF/EGFR-mediated cellular-level models have started to incorporate a simple molecular level in studying e.g., cell migration in breast cancer [22], cell proliferation [23], and autocrine receptor-ligand dynamics [24,25]. We argue that a more detailed understanding of a complex cancer system requires integrating both molecular-and cellular-level works to properly examine multicellular dynamics. To our knowledge, to date, no multiscale model of NSCLC has been developed or published.
Our group has been developing multiscale models to investigate highly malignant brain tumors as complex dynamic and self-organizing biosystems. Since this NSCLC model builds on these works, we will briefly review some milestones. First, an agent-based model for studying the spatio-temporal expansion of virtual glioma cells in a two-dimensional (2D) environment was built and the relationship between rapid growth and extensive tissue infiltration was investigated [26,27]. This 'micro-macro' framework was then extended 'top-down' by incorporating an EGFR molecular interaction network [28] so that molecular dynamics at the protein level could be related to multi-cellular tumor growth patterns [29]. Most recently, an explicit cell cycle description was implemented to study in more detail brain tumor growth dynamics in a three-dimensional (3D) context [30]. These previous works have provided a computational paradigm in which biological processes have been successfully simulated from the molecular scale up to the cellular level and beyond. This progress led us to test the platform's applicability to and flexibility for other cancer types as well.
In this paper, we have therefore extended these previous modeling works to the case of NSCLC. Necessary modifi-cations include at the molecular level the implementation of a NSCLC-specific EGFR-ERK signal transduction pathway. A novel, data-driven switch that is operated by two key molecules, i.e. phospholipase Cγ (PLCγ) and extracellular signal-regulated kinase (ERK), processes the phenotypic decision at the cellular level. The aim of this in silico work is to provide insights into the externally triggered molecular-level dynamics that govern phenotypic changes and thus impact multicellular patterns in NSCLC. In the following sections, we will first show the detailed design of the model before we present and then discuss the simulation results.

Molecular Signaling Pathway
The kinetic model of the implemented NSCLC-specific molecular signaling pathway, which consists of 20 molecules, is shown in Fig. 1. These proteins, including both receptor (EGFR) and non-receptor kinases (e.g., PLCγ and protein kinase C (PKC) [31,32], Raf, mitogen-activated protein kinase kinase (MEK), and ERK [33][34][35]), have been experimentally or clinically proven to play an important role in NSCLC tumorigenesis. Although in reality these molecules fulfil their functions by interacting with a multitude of other molecular species from many distinct pathways [36,37], we choose to start with these proteins not only because of their significance in the case of NSCLC but also since most of their kinetic parameters can be found in the literature. Also, it is reasonable to reduce the number of involved molecules as a starting point for modeling [38]. Amongst these proteins, both PLCγ and ERK are of particular interest for determining the cell's phenotypic changes as we will detail below.
Kinetic equations are written in terms of concentrations and the reaction rates are functions of concentrations. The association and dissociation steps are characterized by first-order and second-order rate constants, respectively. We note that, although in reality chemical reactions of second or higher order are two-step processes, they are usually treated as a one-step process in mathematical modeling [39]. Our model is based on a total of 20 ordinary differential equations (ODEs) and uses exactly the same modeling techniques as other pathway analysis studies (see [11,12] for detailed definitions). For simplicity, the ODEs for different molecules were calculated by Eq. (1): where X i represents one of these 20 molecular pathway components. In Eq. (1), the change in concentration of molecule X i is the result of the reaction rates producing Kinetic model of the NSCLC-specific EGFR signaling pathway Figure 1 Kinetic model of the NSCLC-specific EGFR signaling pathway. The arrows represent the reactions specified in Tables  1 and 2.

Micro-Environment
The 2D virtual micro-environment is made up of a discrete lattice consisting of a grid with 200 × 200 points (Fig.  2). We use p(i, j) to express each point in the lattice, where i and j indicate the integer location in Euclidean terms. One single, distant nutrient source (simulating a crosssectional blood vessel) is located at p(150, 150). To start with, a number of M × N cells (in other words, an M-by-N matrix) are initialized in the center of the lattice (and this number can be set to meet different simulation purposes). Each grid point can be occupied by one cell only or remain empty at a time.
Three external chemical cues are employed in the model: EGF, glucose and oxygen tension. As we have done in previous studies [29,30], the nutrient source carries the highest value of these three diffusive cues, which implicates that it is the most attractive location for the chemotactically acting tumor cells. Then, by means of normal distribution, each grid point of the lattice is assigned a concentration profile of these three cues. The levels of these distributions are weighted by the distance, d ij , of a given grid point from the nutrient source. The distributions of these three cues are described by the following equations: Moreover, the three chemotactic cues continue to diffuse over the lattice throughout the entire process of a simulation with a fixed rate, using the following equation: where M represents one of the three external cues, and t represents a time step. The coefficients in Eqs. (2)(3)(4)(5) are listed in Table 3 (see also [30] for more details). It is evident then that the closer a given location is to the nutrient source, the higher the levels of the three cues will be at this grid point. Glucose will be continuously taken up by cells to support their metabolism. Only the nutrient source, p(150, 150), is replenished at each time step while all other grid points are not. In addition, cells take up both their own EGF and that secreted by adjoining cells in our model, because cancer cells act in both autocrine and paracrine manner in consuming EGF [40,41].
Each cell encompasses a self-maintained molecular interaction network (shown in Fig. 1) and the simulation system records the molecular composite profile at every time step to determine the cell's phenotype for the next step. In , , ,... . t

Reaction number Equation Kinetic parameter Reference
k 9 ·X 9 -k -9 ·X 10 k 9 = 1 k -9 = 0.03 [11] [12] between time steps, the chemical environment is being updated, including EGF and glucose concentration as well as oxygen tension (according to Eq. (5)). When the first cell reaches the nutrient source the simulation run is terminated.

Cellular Phenotype Decision
Four tumor cell phenotypes are considered in the model: proliferation, migration, quiescence and death. Cell death is triggered when the on site glucose concentration drops below 8 mM [42]. A cell turns quiescent when the on site glucose concentration is between 8 mM and 16 mM, when it does not meet conditions for migration or proliferation (see below), or when it cannot find an empty location to migrate to or proliferate into.
The most important two phenotypic traits for spatio-temporal expansion, i.e. migration and proliferation, are decided by evaluating the dynamics of the following critical intracellular molecules. (1) PLCγ is known to be involved in directing cell movement in response to EGF [43][44][45]; PLCγ dynamics are accelerated during migration Two-dimensional virtual micro-environment    [46]. Therefore, in our model, the rate of change of PLCγ (ROC PLC ) decides if a cell proceeds to migration or not. That is, if ROC PLC exceeds a certain set threshold, T PLC , the cell has the potential to migrate. (2) Similarly, the rate of change of ERK (ROC ERK ) decides if a cell proceeds with proliferation. ERK has been found experimentally to have a strong influence on cell proliferation [33,47,48], and transient activation of ERK with EGF leads to cell replication [49,50]. If a cell decides to migrate or proliferate, it will search for an appropriate location to move to or for its offspring to reside in. Candidate locations are those grid points surrounding the cell. Implementing a cell surface receptor-mediated chemotactic evaluation, the most appropriate location is detected by using a 'search-precision' mechanism [27] according to: where T ij represents the perceived attractiveness of location p(i, j), L ij represents the result of an evaluation function for location p(i, j) (see [27] for the definition of L ij ), and ε~ N(μ, σ 2 ) is an error term following a normal distribution with mean μ and variance σ 2 . ψ ∈ [0,1] denotes the search-precision parameter that for a given run is held constant for all cells. Briefly, for a given cell at a certain location, when ψ = 0 the cell performs a pure random walk, whereas when ψ = 1 the cell always selects the location with the highest glucose concentration. Based on previous results [27], we set ψ = 0.7 because this value tends to lead to the highest average velocity of the tumor's spatial expansion.
It is worth noting that even if ROC PLC or ROC ERK exceed their corresponding thresholds, it does not necessarily have to lead to cell migration or proliferation. Rather, if nowhere else to go, the cell remains quiescent and continues to search for an empty location at the next time step.
Any cell in the process of changing its phenotype will fall into one of these four categories: (i) ROC PLC < T PLC and ROC ERK < T ERK ; (ii) ROC PLC > T PLC and ROC ERK < T ERK ; (iii) ROC PLC < T PLC and ROC ERK > T ERK ; and (iv) ROC PLC > T PLC and ROC ERK > T ERK . Figure 3 lists these conditions and their phenotypic consequences, respectively. Following the first three cell decisions is straightforward; first, if a cell experiences condition (i) no phenotypic change results as both ROC PLC and ROC ERK remain below their corresponding thresholds; however, if a cell faces condition (ii) the cell migrates because of ROC PLC exceeding its threshold while in the presence of (iii) the cell proliferates due to ROC ERK exceeding its threshold. However, for (iv), and in the absence of any specific experimental data, i.e. for the case that both ROC PLC and ROC ERK exceed their corresponding thresholds, we explored two hypotheses: 'rule A' yielding migration advantage (i.e., the cell decides to migrate) whereas 'rule B' resulting in a proliferation advantage (i.e., the cell decides to proliferate). For simplicity, decision rules for the first three conditions are referred to 'general rules', while rules A and B are referred to 'special rules' hereafter. In the following section, we will describe the corresponding simulation results.

Results
Our algorithm was implemented in C/C++. A total of 49 seed cells were initially set up in the center of the lattice, and these cells were arranged in a 7 × 7 square shape (i.e., M = 7 and N = 7, see Fig. 2 for the configuration of the seed cells). We defined cell IDs from 0 to 48 (left to right, bottom to top). To investigate cell expansion dynamics, we moni-  Through the distinct microenvironmental conditions they face, these corner cells exemplify the impact of location on single cell behavior, while they however still grasp the nature of the entire system. As described before, both rules A and B were tested for each different simulation condition. Figure 4 shows two simulation results for rules A and B, respectively. The simulations were conducted with a standard EGF concentration of 2.56 nM. Note that this concentration is derived from the literature [51,52] and has been rescaled to fit our model as a benchmark starting point for further simulations. In the upper panel of Fig.  4(a) for rule A, tumor cells first display on site proliferation prior to exhibiting extensive migratory behavior towards the nutrient source. However, for rule B (lower panel), cells remain stationary proliferative throughout, thereby increasing the tumor radius yet without substantial mobility-driven spatial expansion. The run time for the latter case (rule B) was considerably longer than for rule A. Based on the criterion chosen for terminating the run, i.e. the first cell reaching the nutrient source, this result is somewhat expected since rule A favors migration whereas rule B promotes proliferation. This is further supported by analysis of the evolution of the various phenotypes and the change of [total] cell numbers ( Fig. 4(b)). While both rules generate all three cell phenotypes (proliferation (dark blue), migration (red), and quiescence (green)), rule A (left panel) indeed appears to result in a cancer cell population that exhibits a larger migratory fraction than the one emerging through rule B (right panel) which, however, yields a larger portion of proliferative cells (light blue). It is thus not surprising that for rule B, the [total] cell population of the tumor system exceeds the one achieved through rule A by a factor of 10.

Influence of Decision Rules on Phenotypic Changes
To better understand the significance of each rule for the tumor system, we have investigated its influence on generating the intended phenotype. Figure 5 shows the weight of rule A on migration (a), and that of rule B on proliferation (b). (The results are taken from the two simulation runs reported in Fig. 4). In Fig. 5(a) has influence on proliferation only ( Fig. 5(b)) and it contributes more to inducing proliferations than the corresponding general rule does.
However, as documented in the linear least square fittings, the rate at which rule A causes an increase in migration exceeds by far the one by which rule B induces an increase in proliferation. This indicates that the influence of rule A on increasing migrations is more substantial than that of rule B on increasing proliferations. Being particularly interested in gaining insights into spatially aggressive tumors, we continue in the following with investigating the implications of rule A on microscopic and molecular level dynamics of the cancer system.  corner and p(110,110) is the upper right corner (5 runs, data not shown).

Discussion & Future Works
While using mathematical models to investigate the behavior of signaling networks is hardly new, understanding a complex biosystem, such as a tumor, by focusing on the analysis of its molecular or cellular level separately or exclusively is insufficient, particularly if it excludes the interaction with the surrounding tissue. Recent analyses of signaling pathways in mammalian systems have revealed that highly connected sub-cellular networks generate signals in a context dependent manner [53]. That is, biological processes take place in heterogeneous and highly structured environments [54] and such extrinsic conditions alone can induce the transformation of cells inde- pendent of genetic mutations as has been shown for the case of melanoma [55]. Taken together, modeling of cancer systems requires the analysis and use of signaling pathways in a simulated cancer environment (context) across different spatial-temporal scales.
Our group has been focusing on the development of such multiscale models for studying highly malignant brain tumors [27,29,30,56]. Here, on the basis of these previous works, we presented a 2D multiscale agent-based model to simulate NSCLC. Specifically, we monitored how, dependent on microenvironmental stimuli, molecular profiles dynamically change, and how they affect a single NSCLC cell's phenotype and, eventually, the resulting multicellular patterns.
Proceeding top-down in our analysis, we first evaluated the multicellular readout of molecular 'decision' rules A and B (versus general rules; Fig. 3). The patterns of a more stationary, concentrically growing cancer system (following rule B) are quite different from the rapid, chemotacticallyguided, spatial expansion that can be seen in the tumor regulated by rule A (Fig. 4(a)). Not surprisingly, the latter also operates with many more migratory albeit overall less Weight of decision rules on changing cell phenotypes [total] cells ( Fig. 4(b)). Furthermore, examining in more detail the influence of the two distinct rules on their respective phenotypic yield, we found that the impact of rule A on increasing cell migration is more substantial than rule B's influence on furthering proliferation (Fig. 5). This finding suggests that the migratory rule A can operate the cancer system through incrementally smaller changes (while the simulation system is more robust for rule B). Such sensitivity to migratory cues corresponds well with experimental data on the response of human breast cancer cells, which showed that a spatially successful expansive system reacts rather quickly to even miniscule changes in chemotactic directionality [57,58].
Continuing therefore with rule A, our effort was then geared to gain insights into tumor expansion dynamics not only with regards to extrinsic stimuli but also to cell geography, i.e. a cell's location relative to the replenished nutrient source. Most interestingly, we found a phase transition in the cancer cell closest to the nutrient source (i.e. Cell No 48, while none of the other three corner cells showed similar behavior). Specifically, for a tumor cell at this location, i.e., facing nutrient abundance, proliferation is completely abolished once the extrinsic EGF concentration exceeds a certain level. While this at first may seems rather unexpected, this finding however only confirms the experimentally sound notion that EGF stimulates the spatial expansion of a cancer system [5][6][7][8]. Moreover, with increasing EGF concentrations, the maxima of ROC PLC (Fig. 6) gradually occur earlier which seems to indicate that, under these conditions, the downstream signal is processed faster. Interestingly, such a 'no proliferation, just migration' behavior in the presence of chemo-attractant has indeed already been reported in several in vitro studies using a variety of cancer cell lines [59,60] as well as in non-cancerous human cells [61]. (While admittedly, for the reasons stated, rule B did not receive similar attention in our analysis), we nonetheless argue that, on the basis of our results and the experimental reports they seem to correspond with, rule A and thus a migratory decision prompted by a [ROC PLC > T PLC and ROC ERK > T ERK ] condition is a reasonable outcome for the signaling process taking place in NSCLC also in vitro and in vivo.
Changes at the molecular level for Cell No 48 with an increasing extrinsic EGF concentration (rule A) Figure 6 Changes at the molecular level for Cell No  However, moving the model closer to reality will require a multitude of adjustments, one of which is its ability to account for up-or down-regulation in key molecules as a result of tumorigenesis. As a first step, and since experimental data on over-expression of EGFR in a variety of cancer types, including NSCLC, are ample [62][63][64][65] we have begun to simulate the impact of an increasing number of receptors on the cancer system ( Fig. 7; simulations conducted with an EGFR concentration of 800 nM (per system)). Comparing this preliminary data with those reported in Fig. 6 (simulations conducted with an EGFR concentration of 80 nM (per system)), we find that an EGFR-overexpressing NSCLC tumor seems to operate with even more migration and does so earlier on. The result is a spatially even more aggressive cancer system, which seems to correspond well with the aforementioned experimental studies. And, intriguingly, while the phase transition itself is preserved, it however occurs already at a smaller EGF concentration, hence indicating that the increase in receptor density leads to an amplification of the downstream signal, which again corresponds well with experimental results in examining signaling activities gen-  [66]. Taken together, while preliminary, this finding demonstrates applicability and confirms flexibility of this multiscale platform, hence warrants its further expansion.
There are a number of research tracks that can and should be pursued in future works. First, it will be intriguing to see if, in the presence of a non-replenished nutrient source, the proliferative phenotype eventually can be recovered once extrinsic ligand concentrations fall beyond the phase-transition threshold. More generally, while most of the pathway's parameters, including rate constants and initial component concentrations were obtained from the experimental literature, this data naturally originated from a variety of often stationary experimental settings and different cell types. It therefore represents a less desirable and reliable input than time series data that come from one experimental setting only. Also, some parameters had to be estimated, much like in other well-established pathway models [11,12]. Taken together, future works will have to include not only proper experimental verification of the estimated parameters and evaluation of the simulation results but also, on the in silico side, techniques such as sensitivity analysis to help determine the effects of parameter uncertainties on model outcome [67] and to identify control points for experiment design [68]. While a pathway model cannot be a biological representation in every detail [38] we plan on adding, in incremental steps, other pathways of relevance for NSCLC such as e.g. PI3K/PTEN/AKT [69]. Moreover, simulating a more heterogeneous biochemical environment and implementing both cell-cell and cell-matrix interactions [70] are planned steps at the cellular level that should help representing the cancer system of interest in more detail.
Regardless, we believe that the current model already provides useful insights into NSCLC from a systematic view in terms of quantitatively understanding the relationship between extrinsic chemotactic stimuli, the underlying properties of signaling networks, and the cellular biological responses they trigger. Our results yield several experimentally testable hypotheses and thus further support the use of multiscale models in interdisciplinary cancer research. To our knowledge, this presents the first multiscale computational model of Non-Small Cell Lung Cancer and is thus potentially a significant first step towards realizing a fully validated in silico model for this devastating disease. TGFα: Transforming growth factor α.