Simulating Brain Tumor Heterogeneity with a Multiscale Agent-Based Model: Linking Molecular Signatures, Phenotypes and Expansion Rate

We have extended our previously developed 3D multi-scale agent-based brain tumor model to simulate cancer heterogeneity and to analyze its impact across the scales of interest. While our algorithm continues to employ an epidermal growth factor receptor (EGFR) gene-protein interaction network to determine the cells' phenotype, it now adds an explicit treatment of tumor cell adhesion related to the model's biochemical microenvironment. We simulate a simplified tumor progression pathway that leads to the emergence of five distinct glioma cell clones with different EGFR density and cell 'search precisions'. The in silico results show that microscopic tumor heterogeneity can impact the tumor system's multicellular growth patterns. Our findings further confirm that EGFR density results in the more aggressive clonal populations switching earlier from proliferation-dominated to a more migratory phenotype. Moreover, analyzing the dynamic molecular profile that triggers the phenotypic switch between proliferation and migration, our in silico oncogenomics data display spatial and temporal diversity in documenting the regional impact of tumorigenesis, and thus support the added value of multi-site and repeated assessments in vitro and in vivo. Potential implications from this in silico work for experimental and computational studies are discussed.


INTRODUCTION
We present the extension of our previously developed three-dimensional (3D) agent-based, multi-scale brain tumor model to incorporate a simplified glioma progression pathway in an effort to study the impact of clonal heterogeneity on tumor growth dynamics. For the micromacroscopic environment, we expanded the initial setup [1][2][3] by simulating the interaction between the chemoattractant concentration and the degree of cell adhesion. Like before [3], a set of chemoattractants such as glucose, oxygen tension and transforming growth factor alpha ( α TGF ) diffuses [4-6] along a gradient that depends both on the initial conditions as well as on the cells' consumption and secretion ( α TGF ). Also, the glioma cells' considerably high migration rate in vivo, and even more so, in vitro [7] is impacted by the α TGF concentration to properly reflect available experimental data [8][9][10][11]. Specifically, we modeled the chemotactically attracted cells as follows: once the concentration of the so called homotype attractant [12], in this case α TGF , substantially exceeds that of the heterotype attractant [13] such as glucose, the homotype attractant ( α TGF ) plays a much more important role for directing cell attraction and vice versa. For the molecular environment, each cell is equipped with an epidermal growth factor receptor (EGFR) gene-protein interaction pathway [1-3] that is implicitly integrated with a simplified cell cycle module [14]. The pathway molecule 'activated phospholipase Cγ' (PLCγ) is critical in determining the cell phenotype which represents experimental results reported in part by [15] for the case of breast cancer. However, other than in our previous works where heterogeneity arose only from phenotypic alterations of an otherwise mono-clonal tumor, here, we have integrated a simplified cancer progression model (based on [16]) to generate a heterogeneous brain tumor that eventually consists of five distinct clonal populations that emerged through mutational events. These clones differ in their EGFR density in addition to being able to adjust their phenotypic behavior readily to the micro-environmental conditions. To properly reflect available experimental data, the subclones with higher EGF receptor density are more aggressive [17], and operate with an inefficient metabolism [18,19]; yet, they achieve a higher chemotactic search precision to move faster along the environment's least resistance, most permission and highest attraction paths [20]. We will discuss the results and their potential implications, and introduce plans for future work.

PREVIOUS WORKS
Earlier tumor models focused primarily either on the micro-macroscopic or on the molecular level to simulate tumor growth and migration. For instance, for the micro-macroscopic scale, Mansury et al. [21][22][23] proposed a two-dimensional (2D) agent-based model in which the spatio-temporal expansion of malignant brain tumor cells is guided by environmental heterogeneities (mechanical confinement, toxic metabolites and nutrient sources) in order to gain more insight into the systemic effect of such cellular chemotactic search precision modulations. For the molecular scale, Schoeberl et al. [24] built up a computational model to describe the dynamics of the MAPK cascade activated by EGF receptors and Araujo et al. [25] developed an EGFR signaling network model that not only shows the protein dynamics in the EGFR pathway but also predicts the therapeutic effect of inhibitors on the molecular network. In an effort to develop a cell cycle module, Alarcon et al. [14] extended the work of Tyson and Novak [26] to model the cancer cell cycle process under hypoxic conditions. Additionally, to simulate mutation among different clones, Sole and Deisboeck [16] revised the previous quasispecies model of Swetina and Schuster [27] to research the impact of genetic instability on tumor growth. While these models undoubtedly are useful at their particular scale of interest, the quest of studying the tumor as a biological system that crosses scales requires a multiscale modeling approach. Therefore, in [28] and [1, 2] our group developed a 2D multi-scale agent-based tumor model to integrate both, the molecular and micro-macroscopic environments. Specifically, in [1] we put forward an EGFR gene-protein interaction pathway in which activated PLCγ is employed to determine the cell phenotype.
Subsequently, in [2] we then compared the biological behavior of glioma cell clones with varying EGFR density. Building on this work, in [3] we recently introduced a 3D multi-scale agent-based brain tumor model which explicitly includes a cell cycle description. However, this multi-scale model has a number of shortcomings that include simulating a monoclonal, and thus a genetically stable tumor, and it also suffers from a rather simplistic description of the tissue microenvironment. In an effort to model these features more realistically the subsequent section will describe the setup of a much more sophisticated 3D multi-scale brain tumor model that includes explicitly the phenotype of cell adhesion as it relates to chemoattractant [10] and incorporates specifically the emergence of heterogeneous tumor cell clones.

MATHEMATICAL MODEL
As before, this multi-scale agent-based brain tumor model includes macro-microscopic and molecular environments that interact with each other. The following sections describe these environments and their relationships in detail.

Clonal heterogeneity
Here, a slice of virtual brain tissue is modeled with a 100*100*100 3D rectangular lattice as depicted in Figure 1. A continually replenished nutrient source, representing a blood vessel, is located in Cube 4.

Figure 1
To yield tumor heterogeneity, an element of genetic instability is introduced that is based on the mutation model presented by [16]. Our 'progression pathway' here incorporates five glioma cell clones (A, B, C, D and E) that emerge sequentially (Figure 2), and that are distinct in their EGFR receptor density. Based on experimental data, in our setup, cell clones with higher EGFR receptor density are more aggressive [17,29,30] and operate with an inefficient metabolism [18,19]; yet, they achieve a higher chemotactic search precision to move faster along the environment's least resistance, most permission and highest attraction paths [20].

Figure 2
In the beginning, five hundred clone A cells are initialized as close as possible to the center of the lattice as each cell occupies one lattice point. As a first iteration, the model implements a linear progression pathway in assuming that clone A can only mutate to B, clone B to C, C to D, and D to E. The respective mutation rates depend on the cell's proliferation period and are listed in Table 1.  [13]. Our setup places a blood vessel in the center of Cube 4 to supply glucose, α TGF and oxygen for the developing brain tumor. Initially, these chemoattractants are dispersed by normal distribution described in Eq. 1-3: In  [3,14], o σ is the parameter controlling the dispersion of the oxygen tension level, and ijk d is the L-infinity [21] distance to the center of Cube 4. During the simulation, at each time step the cell will take up glucose to maintain its metabolism and secrete α TGF [1] triggering both paracrine and autocrine biological behavior that is described by Eq. 4-5: In Eq. 4-5, t represents the time step, T S is the α TGF secretion rate [32] and n r is the cell's glucose uptake coefficient [33] for each clone. Because [19] argues that the metabolism of cancer cells is approximately eight times larger than the metabolism of normal cells and since [33] discovered that for U87 human glioma cells the glucose uptake rate is 0.77 pmol/h per cell, we assign each clone (A, B, C, D and E) a different glucose uptake rate related to its EGFR density, which is listed in the Table 2.
and o D stands for the diffusion coefficient of oxygen [4].

Spatial search process
Representing the biological mechanism of chemotaxis we have again implemented a search mechanism as put forward already in [21,22].
where ijk T stands for the perceived attractiveness of location (i,j,k), ijk L represents the correct, is an error term that is normally distributed with mean µ and variance 2 σ . The parameter n ψ is positive between zero and one, 1 0 ≤ ≤ n ψ , and represents the extent of the 'search precision' (with '0' exhibiting purely random walk behavior, and '1' being fully biased without any error in signal processing); the subscript of ψ represents the different clones. In accordance with experimental data [29] that report an acceleration of the tumor's spatio-temporal expansion rate with increasing EGFR density per cell, we argue that a cancer cell with high EGFR density should express a comparably high search precision (being lesser affected by the noise function). Thus, we assign a search precision value of 0.5 for clone A and B, and 0.7 for clone C, D and E, respectively.
The definitions of the Von Neumann and Moore neighborhoods [21] usually are restricted to 2D lattices. We therefore define 3D Von Neumann neighborhood, 3D Diagonal neighborhood, and 3D Moore neighborhood as detailed in Eq.10-12 with Here, (x 0 , y 0 , z 0 ) denotes the current location of the cell and (x ,y ,z) represents the coordinates of its neighborhood. We set the grid size to 10μm to reflect the (idealized) average diameter of a real tumor cell and each simulation step represents approximately 2.5 hours, similar to our previous study [3]. The maximum distance for a given tumor cell to migrate under a 3D Von Neumann neighborhood search is therefore 10μm/simulation step, ). Otherwise, the cell will select its next location from the more restrictive 3D Von ). We therefore incorporate in our expanded model here not only that the chemotactically-driven search precision of clones C, D, E is higher than that of clones A and B. Rather, in addition, we implement that a higher α TGF concentration can increase the cell's spatial "permission", thus reflecting its progressive loss of tumor cell-tumor cell adhesion as reported in [8][9][10][11], and as such implicitly adding an element of haptotaxis.
As discussed earlier, ijk L is the correct, non-erroneous evaluation of location, which includes homotype and heterotype chemoattractants. Eq. 13 models this biological behavior, i.e., if a cell resides in location (x 0 , y 0 , z 0 ) and the sum of its 3D Von Neumann neighborhood , then this cell is subjective to choose the homotype attractant, otherwise the cell is attracted by the heterotype cue.

Cell phenotype
As before, we simulate four different cellular phenotypes, i.e. migration, proliferation, quiescence and apoptosis -now, however, for each of the five different glioma cell clones.
We again implement the concept of 'dichotomy' for glioma cells, first put forward by [35] which argues on the basis of mounting experimental data that migration and proliferation are mutually exclusive behaviors that can happen in sequence, yet hardly at the same time in the same cell. Using breast cancer cells, Dittmar et al. [15]

Molecular environment
The sub-cellular setup is described in detail in [3]. Here, we therefore only briefly summarize this molecular environment, focus on its modifications, and note that the model's intracellular workflow is depicted in Figure 3.(a). Specifically, each cell is equipped with an EGFR geneprotein interaction network shown in Figure 3.(b), which is composed of 14 molecular species in addition to the five molecular species involved in our simplified cell cycle representation as listed in Table 3. These molecular species interact with each other according to mass balance equations, denoted by the general form: where X i is the mass of ith (i=1-14) molecules comprising the EGFR network and As an important extension of our previous works, a proliferation process has now been explicitly integrated into the model in that more aggressive cell clones are associated with a shorter cell doubling time.   To achieve a higher spatial resolution of analysis, we divide the macroscopic tumor into three microscopic regions, sectioning the distance between the tumor and the nutrient source such as shown in Figure 5.(a): region 1 is closest to the nutrient source, region 3 is farthest from the nutrient source, with region 2 being located in between.

Tumor heterogeneity:
In an attempt to quantify clonal heterogeneity we employ the socalled Shannon Index which is defined as where i p is the frequency of clone i within the tumor [37]. Figure 5.(b) displays this index for region 1, 2, and 3 at each time step. While the indices for all three regions generally increase (i.e., exhibit more clonal heterogeneity) region 1 experiences an unexpected drop around t = 100, which indicates the emergence of a more homogenous sub-section prior to a gain of clonal diversity similar to that of the other two geographic sub-regions.

Figure 5
Visualizing this type of emergent heterogeneity by orienting us on curve-crossing points seen in Figure 5.(b), we display in Figure 6.

Figure 6
Expansion rate per tumor region: To support this impression quantitatively, we measure the expansion rate for each sub-region by computing the average expansion radius for each region based on its center of mass. The results show that tumor region 1, i.e., closest to the nutrient source, expands much faster than regions 2 and 3 that are farther from the nutrient source (Figure 7). The lesser expansion rate of region 2 (versus region 3) is due to the fact that cells in region 2 find fewer empty locations to choose from.

Figure 7
Expansion rate per tumor clone: The next question to address is then how much of that regional specificity is due to its prevailing clonal population and as such we measure the expansion rate for clone A, B, C, D and E at each time step by computing the average expansion radius for each clone based on the center of each clone. From Figure 8 one can see that clone E (with the highest EGFR density per cell) exhibits the highest expansion rate increase. In contrast, although clone A and B (with lower EGFR density) originate first, and thus also start to expand much earlier than clone C, D and E, these two lower malignant clones keep a much smaller expansion rate for most of the simulation.
So far, the results support the notion that a higher density of the more aggressive clone E is responsible for the higher expansion rate increase in the tumor region that is closest to the nutrient source. This therefore warrants a more detailed look into the phenotypic dynamics of the cancerous clones in general, and of clone E in particular.

Figure 8
Phenotypic spectrum per tumor clone: Turning then to the microscopic analysis level by monitoring the phenotypic fractions of clone A, B, C, D and E at each time step (Figure 9. (a)-(e)) it is evident that clones with higher EGFR density are eventually comprised of a larger migratory fraction, while exhibiting a smaller proliferation and quiescence cell population. Specifically, clone E has the largest and most sustained migratory cell population, and, as the curve shows, the 'population' crossover between a more stationary, proliferationdominated to an expansive, migration-dominated clone happens earlier than in the other, less aggressive tumor clones. We also note that, qualitatively, the fluctuations in the curves increase from a fairly periodic pattern in the lower malignant clones A and B, to a much less regular, 'noisy' pattern in the higher malignant clones.

DISCUSSION & CONCLUSIONS
In here, we present an expansion of our 3D multi-scale agent-based brain tumor model that has now been extended to (1) also incorporate an element of genetic instability to simulate glioma progression and to analyze its impact across multiple scales. Building on our previous works [1 -3], this model not only encompasses molecular, microscopic and macroscopic scales of observation, but now also integrates (2)  The tumor region adjacent to a blood vessel representing nutrient source (i.e., region 1) harbors early on the highest clonal heterogeneity (Figure 5.(b)), reflecting the spatio-temporal output of the implemented linear progression path. However, around time step t = 100 this geographic 'tip' region ( Figure 6.(a)-(c)) becomes increasingly homogeneous due to a temporary competitive advantage of the most aggressive clone E with regards to its chemosensing ability and high spatial expansion rate. During a period that is marked by variability across runs (increasing fluctuations in standard deviation), clone E prevails in its quest to move fastest into the direction of "(least resistance,) most permission and highest attraction" [20], and, in the process, impacting the spatio-temporal patterns of the entire tumor system.
However, the micro-environmental conditions closer to the core of the tumor, i.e. for cells in regions 2 and 3 are less selective for clone E and thus continuous tumor progression yields a sustained increase in the Shannon Index that eventually even exceeds the peak value for region 1. One can argue that, other than region 1, regions 2 and 3 allow for an element of temporary coexistence if not cooperation between the clones. Since regions 2 and 3 together add up to a larger volume than region 1 alone, this finding corresponds well with another recent claim of cancer cell cooperation fostering tumorigenesis [38], and also with the results reported by [37] on esophagus cancer which indicated that clonal heterogeneity is a greater risk factor for progression to cancer than the expansion of a homogeneous clone with a low Shannon Index value. This is confirmed in our own results as clone E's dominance of region 1 is relatively short-lived with the heterogeneity index soon returning to values comparable to the two proximal tumor regions. The question arises then as to what precisely are the effects of such clonal heterogeneity on tumor growth dynamics? First, due to their higher density of EGF receptors per cell, clones C, D and primarily E exhibit an increased expansion rate (Figure 8) over clones A and B. Correspondingly, Figure 9 confirms that these tumor clones with higher EGFR density switch earlier to an invasive phenotype-dominated cell population that operates with a larger fraction of migratory cells, yet with lesser proliferative and quiescent cells. Our in silico model shows that higher EGFR expression leads to faster expansion of tumor areas that harbor more such aggressive cell populations, in areas of nutrient abundance. This, in turn, will result in overall growth asymmetries (see region 1, clinics. In fact, applying a global genomic analysis Maher et al. [39] recently showed widespread differences between primary and secondary glioblastoma entities and reported that the latter group (which develops over multiple years from lower malignant astrocytomas [40]) is heterogeneous in its molecular pathogenesis. Our findings should therefore not be misinterpreted as questioning the value of such high-throughput data, rather, our in silico results argue that assessing the molecular profile repeatedly, and at different sites of the tumor, should add predictive power. Interestingly, we found that at time step 200 the percentage changes of some molecules (i.e., 3 X , 4 X , 5 X , 6 X and 10 X ) in region 2 and 3 is more substantial than in region 1 regardless of the switch per se. The reason for this is that the higher EGFR density in clones C-E results in higher concentrations of other downstream molecules as well. Hence, the relative molecular fluctuations should be lesser in clones C, D and E than in clones A and B. And, since region 1 is mostly comprised of cells derived from clones C-E, its molecular composite signature appears to be more preserved between the time points than that of the other two regions. However, over time many more aggressive cells originate from clones A and B also in the more proximal tumor regions, and we see therefore eventually a somewhat more homogeneously aggressive pattern where regions 2 and 3 start to resemble the distant region 1. As such, by monitoring the molecular signature dynamics per region, one can at least in silico visualize the gradual tumorigenic evolution of sub-regions of the tumor, in an effort to better understand how they may contribute to multicellular growth behavior.
Despite these considerable extensions, our modeling platform of course has still a number of shortcomings that can and should be addressed in future works. For example, we currently rely on a fixed mutation rate between each clone and mutations only occur linearly from A to E, without any non-linear 'shortcuts' and without assigning different weight to specific mutations. In subsequent research, we intend to implement a dynamical mutation rate that depends both on intrinsic and extrinsic or microenvironmental factors. Other extensions on the molecular level will rely on sensitivity analysis, or may include a more explicit treatment of vascular architecture. Nonetheless, this work significantly enhances our previously introduced platform [1-3] by including the element of genetic instability and integrating both cell cycle and cell doubling time into the cells' proliferation process. Given the fact that we begin to be able to validate some of the model's predictions with experimental studies, we firmly believe that this type of interdisciplinary approach holds significant promise to help improve our understanding of brain cancer progression and its impact on the tumor system's regional and global growth dynamics.     Depicted is the Shannon Index (y-axis) over time (x-axis) for these three distinct regions.
Shown is the result of 10 runs with standard deviation.