A MULTIPLE TIME-SCALE COMPUTATIONAL MODEL OF A TUMOR AND ITS MICRO ENVIRONMENT

. Experimental evidence suggests that a tumor’s environment may be critical to designing successful therapeutic protocols: Modeling interactions between a tumor and its environment could improve our understanding of tumor growth and inform approaches to treatment. This paper describes an eﬃcient, ﬂexible, hybrid cellular automaton-based implementation of numerical solutions to multiple time-scale reaction-diﬀusion equations, applied to a model of tumor proliferation. The growth and maintenance of cells in our simulation depend on the rate of cellular energy (ATP) metabolized from nearby nutrients such as glucose and oxygen. Nutrient consumption rates are functions of local pH as well as local concentrations of oxygen and other fuels. The diﬀusion of these nutrients is modeled using a novel variation of random-walk techniques. Furthermore, we detail the eﬀects of three boundary update rules on simulations, describing their eﬀects on computational eﬃciency and biological realism. Qualitative and quantitative results from simulations pro- vide insight on how tumor growth is aﬀected by various environmental changes such as micro-vessel density or lower pH, both of high interest in current cancer research.

Abstract. Experimental evidence suggests that a tumor's environment may be critical to designing successful therapeutic protocols: Modeling interactions between a tumor and its environment could improve our understanding of tumor growth and inform approaches to treatment. This paper describes an efficient, flexible, hybrid cellular automaton-based implementation of numerical solutions to multiple time-scale reaction-diffusion equations, applied to a model of tumor proliferation. The growth and maintenance of cells in our simulation depend on the rate of cellular energy (ATP) metabolized from nearby nutrients such as glucose and oxygen. Nutrient consumption rates are functions of local pH as well as local concentrations of oxygen and other fuels. The diffusion of these nutrients is modeled using a novel variation of randomwalk techniques. Furthermore, we detail the effects of three boundary update rules on simulations, describing their effects on computational efficiency and biological realism. Qualitative and quantitative results from simulations provide insight on how tumor growth is affected by various environmental changes such as micro-vessel density or lower pH, both of high interest in current cancer research.
1. Introduction. Fast, biologically realistic hybrid cellular automaton-based tumor simulations could benefit cancer treatment by allowing researchers to simulate proposed treatment strategies in order to determine the most effective treatment. In thin, effectively two-dimensional slices of tissue, tumor development can be efficiently, realistically simulated using a hybrid cellular automaton (HCA), a grid whose elements evolve based on a system of differential equations. This paper describes the development and implementation of an HCA-based mathematical model of tumor development that includes blood flow, nutrient diffusion, and cellular consumption of nutrients.
Time and space are both discretized in this HCA model; time passes in small, fixed-length steps, and each grid element corresponds to a small region of space, represented by local state data such as nutrient concentrations and cell populations. All grid elements are concurrently updated based on their local state and those of neighboring elements at the previous time step. That is, S x,y (t + 1) = f (S x,y (t), S x−1,y (t), S x+1,y (t), S x,y−1 (t), S x,y+1 (t)) (1) where S x,y (t) is the state of element (x, y) at time t. Cellular automata (CAs), originally described in [74], have a grid with infinite length in both dimensions; HCAs differ from CAs only in that some modeled quantities (e.g., nutrient concentrations) are continuously rather than discretely described, but the spatial representations of HCAs are based on grids exactly as fully discrete CAs are. In these CAs or CA-based models, all grid elements have the same number of neighbors, so they are all updated using the same set of rules. Practical implementations of CAs have finite grid sizes, however, so edge and corner elements have fewer neighbors than interior elements. Special boundary rules must therefore be devised to govern the updating of these elements. In this paper we discuss the implementation of a HCA model with three types of boundary rules.

2.
Background. There is a long history of mathematical models of tumor growth that leads to this work. Early models consisted of systems of differential equations, such as [22], in which differential equations model effects of immune response to tumor growth, and [54], in which delays are introduced into the model. Previous models with differential equation-based representations of different cell populations, however, have ignored spatial characteristics such as blood vessel location and tissue heterogeneity. Furthermore, clinical insights can be gained from models that incorporate and integrate factors from a variety of biological systems that influence tumor development and treatment; models such as [3,52], for example, incorporate elements ranging from angiogenic factors to cytotoxic immune cell dynamics and pH gradients. In particular, evidence suggests that local conditions such as nutrient availability, oxygen concentration, and local pH have key roles in carcinogenesis and tumor behavior, affecting tumor growth, malignancy, and response to treatment [72], and our model includes the effects of these factors on cell growth.
In this paper, we present a hybrid cellular automaton model of early solid tumor growth from an energy budget perspective, where nutrient consumption and ATP production are explicitly represented as dependent on oxygen concentration, available fuel concentration, and pH. The approach models tumor development over both time and space, and it incorporates the effects of variations in blood supply as well as the diffusion of oxygen and other nutrients on the morphology and aggressiveness of the tumor. Simulations of the model are compared with experimental data showing concentrations of hypoxic cells and tumor morphologies consistent with those observed in both in vitro and in vivo studies ( [58], [38]).
Several approaches to the spatial modeling of tumors in the micro-environment have appeared in the literature. One approach uses partial differential equations and/or level-set methods to describe the evolution of the tumor boundary [11,16]; others describe mechanical models of individual cells linked together [55], and others bridge the gap between the continuum and discrete models, such as [59,44,57,56] and [75]. Pertinent to the model presented in this paper is [10], which describes a mathematical model of the effect of cellular adhesion on tumor growth. Highly adhesive tumor cells are likely to result in a compact, circular tumor, whereas less adhesive cells may produce a more diffuse, branching tumor structure.
2.1. CA models of tumor growth. In [1], a cellular automaton-based model of tumor growth, the effects of differential blood flow and vessel adaptation are modeled. In this model, blood is treated as a non-Newtonian fluid which flows through a network of blood vessels. For example, red blood cells are distributed unevenly between the two daughter vessels in a network bifurcation; the daughter vessel with higher flow velocity receives a greater number of red blood cells. Furthermore, the blood vessel network adapts in order to best provide nutrients to the tissue. For example, a decrease in red blood cell flow rate may trigger an increase in blood vessel radius. Zero-flux boundary conditions are used in this simulation.
Another cellular automaton-based tumor simulation, described in [18,45], includes the immune response. Two types of immune cells are modeled: natural killer (NK) cells and cytotoxic T lymphocytes (CTLs). NK cells wander the tissue area until encountering a tumor cell, which they then attack. When a tumor cell is lysed in this manner, additional CTLs are recruited to the same site, where they will attack other tumor cells in the immediate neighborhood. Immune cells in tumor-free locations die or randomly move to neighboring locations.
IMMSIM [53] is a cellular automaton-based simulation of immune system response. The IMMSIM model simulates discrete events such as the destruction of an individual antigen by an immune cell, rather than taking a continuous, differential equation-based approach. As in [1,18,45], biological components such as antigens and antibodies diffuse among grid elements. If two components are present in the same grid element, an interaction event such as antigen-antibody binding may occur. Each type of interaction has specific rules by which the probabilities of final outcomes are determined; the end result is then decided stochastically.
ParImm [4], a further development of IMMSIM, is another cellular automatonbased immune system simulator. Like IMMSIM, ParImm models cellular and molecular components of the immune system. In ParImm, cells are implemented as stochastic finite state machines, with states corresponding to possible cellular conditions, such as normal, resting, infected, or dead. Cellular interactions may result in a change in cellular state. For example, an interaction between an immune cell and an infected cell might might change the infected cell's state to "dead." Several recent automata models allow cellular behavior to depend on local chemical concentrations. Dormann et al. reproduced the layered structure of avascular V-79 tumors using cellular automata governed by chemical signaling [24], while Patel incorporated H + ion production rates and randomly placed blood vessels to explore optimal micro-vessel densities and pH levels for tumor invasion [51]. In addition, Alarcon employed vascular dynamics in a cellular automaton model to explore the role of blood flow and oxygen distribution [1].
Smallbone et al. [63] use a hybrid cellular-automaton model to explore the interactions between the tumor and its environment, and the effect that these have on carcinogenesis and invasiveness. This model builds on these earlier models by incorporating the interdependency between tumor cell metabolism and local chemical concentrations. In particular, it takes into consideration the delivery of oxygen and other fuels via micro-vessels, the constriction of blood vessels due to pressure from the surrounding cells, cell-cell adhesion, the synthesis of nutrients for the production of ATP, the budgeting of cellular energy towards growth and maintenance processes, and the production of acidic byproducts during glycolysis. In this model, cell proliferation and death is a deterministic function of available energy, in contrast to CA models in the literature which use probabilistic rules [27,45,18,63].
Several non-biological cellular automaton-based simulation methods [31,43,30] also exist. While these simulations vary widely in their applications, all CA-based simulation methods require boundary rules. It is likely that the boundary rules discussed in this paper could be applied to other biological or non-biological gridbased methods.
2.2. Tumor metabolism. While tumor cells consume a variety of fuels including glutamine, palmitate, oleate and others, oxygen and glucose are the most significant contributors to cellular energy [33]. The model discussed by Patel et al. in [51] considers tumor metabolism as characteristically glycolytic, where cells rely on anaerobic pathways for the majority of ATP production, regardless of local oxygen availability. This concept of a glycolytic phenotype was first proposed in 1929 by Warburg after he found that mouse ascites cancer cells produced high amounts of lactate even with abundant oxygen available, evidence of aerobic glycolysis [76]. Currently, this concept underlies many approaches to basic cancer research as well as cancer prognosis and treatment.
However, recent experimental evidence challenges the perception of a glycolytic phenotype [76]. Guppy et al. [33] conclude that under aerobic conditions, the metabolism of MCF-7 tumor cells is not different from that of many normal cells. If tumor cells were characteristically glycolytic, they would rely on glucose for the majority of their ATP production, even with normal oxygen availability. On the other hand, because tumor cells often consume a markedly high rate of glucose, PET (positron emission tomography) has become a very popular technique for detecting and measuring a wide variety of tumors, supporting the concept of a glycolytic phenotype [49]. Gatenby and Gillies [29] have proposed a plausible Darwinian scenario in which the glycolytic phenotype evolves as a response to intermittent hypoxia in pre-malignant populations.
Some argue that the oxygen and glucose consumption rates associated with tumors are more likely due to a heightened "Pasteur effect," where low oxygen conditions force tumor cells to metabolize glucose less efficiently. In fact, the role of tumor oxygenation has received attention recently for a variety of reasons, especially after some researchers have found correlations between low oxygenation and metastatic progression [67], treatment resistance, and patient survival. Oxygen is delivered to tissue via an efficient network of microvasculature, but with increasing size, tumors often compromise blood flow which causes a decrease in oxygen perfusion, leading to hypoxic regions distributed heterogeneously within the tumor, [72], [46], [37], [23]. Our model extends previous work by simulating the progression of hypoxia, its influence on local cell metabolism, and the consequences for solid tumor growth. Our assumptions on the differences in cell metabolism between normal cells and tumor cells differ somewhat from those made in the model of Smallbone et al. [63], in that both tumor cells and normal cells can metabolize either aerobocially or glycolytically. This situation corresponds to what is described as the "second phase of growth" in the evolutionary scenario described in [63]. It supports the hypothesis set forth in [29] by showing that the cell metabolism is sensitive to very small changes in oxygen and hydrogen ion concentrations, and that a small number  of malignant cells can modify the micro-environment enough to allow the formation of a significant tumor mass.

3.
Model development. The simulation represents a thin, effectively twodimensional slice of tissue with a cellular automaton data structure. Automaton grid elements represent small tissue areas and contain cell population and nutrient concentration data for those areas. To simulate the passage of time, all grid elements are updated in a two-step process: First, nutrient concentrations are updated by simulating diffusion and consumption by living cells. Second, cell populations are updated by simulating reproduction and death based on local nutrient concentrations in each grid element. In this implementation, cell growth, death and movement are determined by the amount of energy delivered via the blood vessels. This connects the "health" of the host directly to the growth of the tumor cells, and enables a connection between the macroscopic, holistic state of the larger organism with the local micro-environment of the tumor cells.
The diagram in Figure 1 indicates the model flow showing the inner diffusion loop inside the cell reproduction and movement loop. In this section we describe each of the modeling sub-steps in detail.
3.1. Model components: Nutrients, tissue and blood vessels. The model uses a multi-layered n × n grid to describe cell populations and chemical concentrations throughout a two-dimensional simulation space. Each grid element represents a physical volume of 175x175x40 µm, or 10 −6 ml. Thus the simulation space represents a thin layer of tissue 40 micrometers thick, a reasonable setting for studying the growth of some solid tumors, such as melanoma. Tumor cell and normal cell populations within each grid element are stored in n × n matrices, the matrix V n represents normal cells; the tumor cell population is stored in two matrices, divided into V t for viable cells, V nec for necrotic cells. Thus, V n (i, j) gives the number of normal cells in the grid element representing the ith row and the jth column of the discretized two-dimensional tissue. By modeling the way cell populations in each grid element interact with the local environment and neighboring cells, we are able to investigate early tumor growth in a qualitative and quantitative manner.
In vivo, blood vessels deliver nutrients (e.g., oxygen and other fuels) that cells consume, while these same cells produce waste (e.g., lactate) that is removed by the blood vessels. These two phenomena affect local extra-cellular chemical concentrations. Cell metabolism is dependent on the concentrations of nutrients and waste. In favorable conditions, cells can produce enough energy (ATP) from nearby nutrients to provide for cell maintenance and growth. In harsh conditions (e.g., in hypoxic regions), consumption rates are lowered and there is a significant decrease in doubling rate and an increase in cell death [14].
To simulate the micro-vascular network, a subset of grid elements are designated to contain blood vessels. Depending on the simulation desired, the distribution of blood vessels can be randomly generated, it can be prescribed according to some pattern, or it can be made to mimic an actual layer of tissue. In the case of a random distribution, the number of blood vessel elements is determined by experimental micro-vessel density values found in relevant literature, and might be one parameter of interest if the model is used to study the role of tumor angiogenesis, as was done in [51]. The size of the blood vessels can also be modified to be a portion of a grid element.
In a given simulation, all chemical concentrations generally begin at levels within normal ranges found in vivo. Small molecules will diffuse from the blood vessels at a rate dependent on the permeability of the vessel wall and on the gradient of chemical concentration. To store the local chemical concentrations within each grid element, three n × n matrices are used: the O matrix represents oxygen concentration [O 2 ]; G represents the concentration of glucose, representing the fuels available to the cells for maintenance and reproduction, and H represents the concentration of H + ions, determining the pH of the environment.
3.2. Nutrient diffusion: "Dirty Diffusion". The diffusion of oxygen, glucose and H + ions through the simulation space and the growth of cell populations are on two different time-scales: cell growth is on the order of days, but diffusion is on the order of seconds [51]. Thus, in our model, the diffusion process is modeled separately from cell proliferation, death and migration, allowing the simulation of processes on these two different time-scales. In this section, we describe the novel approximation methods used to speed up the simulation of the diffusion process, which occurs on the smaller time-scale and hence requires many more iterations than the cell-proliferation sub-routine.
Because small molecules diffuse on a time scale that is much faster than cell migration or proliferation, in many spatially dependent models of tumor growth that include the effects of nutrient concentration, the diffusion process is assumed to have reached its steady state so that the PDE can be replaced by an ODE. Alternatively, in a discretized version, the change in concentrations of nutrients is the solution of a large, sparse, system of linear equations (see Equation (5) ). This system is then typically solved using an iterative method, such as the SOR method [51] or an LU decomposition [45]. In other implementations (e.g., [15]), the diffusion step is modeled explicitly using a finite element method or, as in [16], numerical solutions are obtained using level set methods. In all of these cases, computational costs are large: simulations on reasonably large grids can take hours unless some simplifying assumptions are made, such as radial symmetry [11,39,15]. In our method we exploit the difference in time scales without assuming a steady state for the diffusing nutrients by using a local averaging technique. Briefly, this technique involves averaging concentrations over 3 × 3 non-overlapping neighborhoods, known as Margolus neighborhoods, [26], and is depicted graphically in Figure  2. The idea is to by-step the intermediate graphs in this figure, and jump immediately from the first one to the last one. One interpretation of this numerical method is that in one step it gives the expected value of the concentrations of molecules moving randomly in space.
The center grid elements of each neighborhood are updated at each sub-iteration; therefore, 9 sub-iterations are required to update all grid elements. If the Margolus Neighborhoods are chosen in a random order, spurious drift patterns are avoided, and the numerical scheme converges very quickly to solutions of the continuous diffusion equation, Equation (3), without giving up the explicit time dependency. Experimentally measured diffusion coefficients (see Table A), are used to estimate the amount of time represented by the local averaging step, which will vary for each molecule.
We call our method "Dirty Diffusion," denoting a coarse approximation to the continuous representation of diffusion. A comparison of the Dirty Diffusion method and the SOR method is shown in Figure 3. This figure depicts solutions to an equation of the form given in Equation (3), where concentrations are constant at certain grid elements representing blood vessels (darkest grid elements), and the consumption function, k, is zero. On the sides of the grid, the algorithm is modified depending on which boundary condition is chosen. See Section 3.7 for details on these boundary rules.
Details of the Dirty Diffusion method. We explain the Dirty Diffusion method in more detail by describing its implementation for the solution of one diffusion equation. Consider a concentration of small molecules such as oxygen, denoted by C(t, r), where t represents time and r represents the location in space. Suppose the molecules diffuse with constant diffusion coefficient, D, and are consumed by cells, represented by T (r), at a rate k(t, r, T (r)). Then the partial differential equation describing the evolution of the concentration is: which in two dimensions becomes (we omit the independent variables where the dependence is clear): Discretizing space into an n × n grid, and discretizing time at intervals of size ∆t gives the following discretized evolution rule for the concentration C i,j in the (i, j)th grid element: where the discretization of the partial derivative is: and ∆t is small relative to the diffusion rate, D. Note that the evolution rule depends on concentrations at the (i, j)th grid element and at the four neighboring elements. If the diffusion occurs at a much higher rate than consumption, then we might approximate the solution to Equation (3) by assuming that the nutrients first reach a steady state, and then are consumed. This is the approach taken by many modelers (e.g., in [51]). For the discretized equation, the steady state is a solution of n 2 equations of the form: where 1 ≤ i, j ≤ n. Thus, at the steady state of the diffusion portion of the equation, the concentration of nutrient in a given grid element is equal to the average of the concentrations in the four neighboring grid elements. We can therefore approach this steady state by successively replacing each grid element and its four neighbors by their average value: The execution of one averaging step, Equation (6), simulates the diffusion of the molecules in and out of the neighboring cells over a time interval ∆t which we can estimate as the expected time it takes for a molecule to diffuse across one grid element: If the rate of change of the cell population T is small relative to the rate of diffusion, then T (x, y, t) can be approximated by a constant in the consumption term of Equation (3). Thus, after one time step of length ∆t, where the time step is short compared to the cell proliferation rate, we have the following update rule: If ∆t is the length of the proliferation time step, then in between each update of the cell matrices, Equation (8) is applied ∆t/∆t times, or until an approximate steady state is reached, whichever comes first. See Figure 1 for a graphical representation of the program flow.
3.3. Nutrient diffusion: Boundary values at blood vessels, modeling vascular constriction and collapse. In our model, nutrients enter the tissue from simulated blood vessels and diffuse to locations with lower concentrations. Blood vessels are simulated by computing the effective nutrient concentration C x at grid elements x that contain blood vessels, based on the bloodstream nutrient concentration b and blood vessel permeability q. The nutrient flux from a blood vessel to a neighboring grid element y is defined as q(b − C y ), where C y is the nutrient concentration in element y. This expression is then combined with Equation (6) to produce the following formula for C x : For every element x that contains a blood vessel, the nutrient concentration in x is set to C x . The rule given by Equation (8) is then applied, causing nutrients to diffuse away from blood vessels into the tissue, or causing metabolic by-products to diffuse from the tissue into the vessel. Modifications can be made if the grid element is not entirely filled with blood vessel, the details of which are omitted here. We have developed Equation (8) which we use in the model to describe the diffusion and consumption/products of two nutrients and metabolic by-products. An important consideration in tumor growth is the collapse of vasculature inside the tumor due to constriction from the pressure of surrounding cells. This vascular collapse leads to two phenomena: angiogenesis [68,35] and metastasis due to hypoxic stress [12,8,42].
Constriction is modeled by setting the blood vessel permeability q as a function of t, the total number of tumor cells in the five-element neighborhood surrounding the vessel. The formula and a graph is given in Figure 4 Figure 5 shows the constriction of blood vessels in the region occupied by tumor cells.  Figure 4. Blood vessel permeability, q(t), as a function of t T , the cell population at which vessel constriction begins, and t max , the population at which the vessel collapses and blood flow is completely cut off.
3.4. Nutrient consumption model. Two types of nutrient molecules are explicitly modeled: oxygen and glucose. Metabolic waste is also modeled as free hydrogen ions created in the production of lactic acid. In the literature, cell metabolism and growth rates for EMT6/Ro tumor cells have been determined to be functions of local environmental factors such as glucose concentrations, oxygen availability, and pH [14]. Based on this experimental evidence, our nutrient consumption model has the following properties: • Glucose consumption increases with glucose concentration.
• Oxygen consumption increases with oxygen concentration.
• Up to a saturation threshold, glucose consumption increases as oxygen concentration decreases. • Oxygen consumption decreases as proton concentration increases.
This model is biologically motivated by two chemical processes used by cells to obtain energy: aerobic respiration and anaerobic respiration. Biological evidence for the model is presented in [14].
Recall that O denotes the current oxygen concentration (mM), G the current glucose concentration (mM), and H is the current ion concentration (mM), so that pH =-log10(H/1000). Casciari et al. [14] derived the following empirical equations for the rate of consumption of oxygen and glucose per cell (units are mol cell −1 s −1 ). We have discretized and modified slightly the original equations from Casciari's paper to give the following equations describing the oxygen and glucose consumed by one cell in a time ∆t.
The lower case letters are empirically determined constants whose values are listed in Table A. Based on experimental evidence [61,41], our model assumes that normal cells consume less glucose than tumor cells as oxygen levels decrease. Therefore, we introduce a parameter c ν into the first term of Equation (12) to model the cellspecific sensitivity of dG/dt to O 2 , where the subscript, ν, indicates the type of cell: normal cell, proliferating tumor cell, or necrotic tumor cell. The model also allows the maximum consumption rate, q, in Equation (13) to depend on cell type, again indicated by the subscript, ν. We have also modified the glucose consumption equation given in [14] under the assumption that the rate saturates as oxygen concentrations go to zero i.e. −dG/dt can never exceed a fixed constant, q ν . Assuming that the initial value of G is non-zero, the inclusion of a saturation term in Equation (12) insures that glucose levels will never reach zero. Therefore, we accept the empirically derived Equation (11) as it is given in [13] , without adding a saturation term.
Metabolic Byproducts and pH: Derivation of Equation (14). Low extra-cellular pH is a trademark of solid tumor growth. When glucose is metabolized, cells produce lactic acid as a byproduct, which lowers the pH within the tumor (although other contributing factors have been suggested as well [36]). The blood vessels regulate pH by removing acidic byproducts, but in solid tumors this phenomena is inhibited by a chaotic vasculature [60]. Other buffers also aid in maintaining a pH suitable for normal tissue. Guppy et al. [33] estimate that the ratio of lactate to glucose production is We note here that other approaches to the mathematical modeling of ATP production have been taken that involve the modeling of the metabolic reactions themselves, rather than starting with the empirically derived rates of consumption as we do here. See [73,6,7] for examples of this approach in models of tumor spheroids. We feel that our simpler approach is suitable in this context, since we are interested ultimately in the long-term macroscopic features of tumor growth, vascular constriction and invasion, and one of our goals is computational efficiency.
The empirical constants, represented by lower case letters in Equations (11), (12) and (13), were determined experimentally by Casciari et al. [14]. Specific values are given in appendix A.
3.5. Model of cell division and death. The cell's ability to withstand a harsh environment plays a large role in tumor progression and resistance to treatment. However, experimental evidence shows that different cell lines, when exposed to hypoxic conditions, vary in their ability to withstand energy deprivation rather than vary in anaerobic consumption rates [62]. Thus, it is important to incorporate energy concepts into our model where cells experience energy deprivation in regions of low O 2 and low pH while cells are able to undergo mitosis in more favorable conditions.
We use a generalized framework for cellular energy budgeting formulated by Kooijman [40] that has been applied to tumor-host energy interactions [71]. This approach uses three key assumptions: 1. Each cell requires a certain rate of energy uptake in order to maintain itself, prioritizing energy use for common cell functions. Let M denote this per-cell maintenance cost. 2. There is also a fixed energy cost for mitosis, implying that the growth rate of a cell population is proportional to the energy available for growth. Let g denote the per-cell cost of mitosis. 3. The maintenance and mitosis parameters, M and g, respectively, depend on cell type.
Thus, if A(t) is the rate of energy intake (mol. (grid element) −1 cycle −1 ), the growth rate for a given number of cells V at grid element (i, j) is where, in general, the energy required for maintenance and mitosis will depend on cell type. Thus, in what follows, we add a subscript to the parameters M and g to denote their dependence on cell type. Cells use ATP as a carrier of cellular energy storage, assembling ATP during metabolic processes and utilizing ATP for numerous cellular functions. Thus the energy intake rate, A(t), described above corresponds to ATP availability rates. If we assume that all of the available energy is used, cell growth and death can then be determined from Equation (15) by calculating the available ATP. To calculate the amount of ATP energy derived from nutrient consumption, we consider oxidative and glycolytic metabolism separately. First, we consider the aerobic synthesis of ATP. We can calculate the total oxidative ATP turnover from the rate of oxygen consumption and a P:O ratio of 2.36 [33]. We may also calculate the number of glucose molecules utilized in oxidative processes, G O , since the oxidization of one glucose molecule requires 6 molecules of O 2 . Once the amount of glucose metabolized aerobically has been calculated, we may deduce that the rest of the glucose was metabolized anaerobically. Glycolysis is a far more inefficient process, producing at most 2 ATP per glucose molecule. In the model, we assume that, on average, 4 ATP molecules are produced for every 3 glucose molecules consumed anaerobically, taking into account any metabolic inefficiencies that may be present [48]. Note that this is a conservative estimate, and this aspect of the model could be adjusted to mimic particular environments. See, for example, the discussion in [7] where it is assumed that 2 ATP are produced per glucose molecule. Letting AT P O and AT P G denote the number of available AT P molecules produced aerobically and anaerobically, respectively, we have the following equations: Note that the last equation is only valid when ∆G ≥ ∆G O , i.e., if some anaerobic glycolysis occurs; otherwise ∆AT P G = 0. Thus, when enough oxygen is present, there may be no ATP molecules produced anaerobically. In experiments such as those reported in [14], it is possible that other fuel substrates are present. Therefore, while Equations (11) -(13) might violate the constraint 6∆G > ∆O, this would correspond to the presence of other fuel substrates in the experimental set-up. In our model, however, we assume that the only nutrients present are glucose and lactate. A refinement of the present model might take other fuels into account. The total available ATP is the sum of Equations (16) and (17), which can then be used to calculate the change in the number of cells of type ν in grid element (i, j): We remark that the calculation of the ATP produced as a function of the nutrient consumption rates provides insight on the energy status of our simulated cell populations. The total amount of energy available to the host at any given time might, in principle, be a measure of the health of the organism, and could potentially be used for prognosis and/or treatment design.
See Appendix A for the values of the constants in all equations.
3.6. Model of cell movement. Each grid element represents a fixed spatial volume, so it can be filled with proliferating tumor cells, which can then invade neighboring grid elements. . However, cellular adhesion molecules, such as integrins and Tissue Factor (TF), play a large role in the mobility of tumor cells, affecting their interactions with neighboring cells and the surrounding extracellular matrix (ECM) [47]. When modeling the progression of a tumor through its different stages, it is necessary to consider the effects of cellular adhesion between cells and the ECM since the secretion or suppression of adhesion molecules is a possible mechanism for the onset of metastasis [2,34,28]. Mathematical models incorporating adhesion include [10] and [25], while [66] and [65] explore the relationship between cell-cell adhesion, pattern formation and disease prognosis. Turner [70,69], provides a method for estimating the diffusion coefficient for biological cells modeled as adhesive, deformable spheres by considering the "potential energy of interaction" between individual cells. Turner shows that the diffusion coefficient is proportional to the second derivative of the cells' energy density, E(n),which can be derived in terms of the biological parameters of individual cells (such as elasticity and adhesiveness). This provides a way to consider a population of discrete cells macroscopically. A polynomial of Laundau-Ginzburg form: An 2 2 + Bn 4 4 fits the solution of E(n) quite well near the equilibrium density, where n is the density of cells, A < 0 and B > 0. We modify Turner's formula for E(n) by setting it to zero below a minimum cell density, at which point adhesive forces have no effect, (see Figure 6).
Thus, below a specific minimal density (cells per grid element) n min , cells are assumed to make no physical contact, so E is zero. At densities greater than n min , cells begin to make contact and adhere, and E decreases to its minimum at the optimal cell density n eq . Above n eq , cell elasticity becomes more pronounced than adhesion, and E increases up to maximum possible density n max .
In our simulations, each new tumor cell is placed in its grid element of origin, or in one of the four neighboring elements, according to the following rules. Tumor cells are placed first in grid elements with populations less than n eq , with densely populated elements receiving cells before sparsely populated elements, in order to minimize total E over the local neighborhood. Once all elements in the five-element neighborhood containing the element of origin and the four adjacent elements are filled to n eq , the remaining cells are distributed randomly over the neighborhood until n max is reached in all grid elements. Any remaining new tumor cells are discarded due to overcrowding.
Adhesion between normal cells is assumed to be negligible compared with that of tumor cells, so no potential energy function is used for the distribution of normal cells. Newly generated normal cells are distributed randomly among the five-element neighborhood surrounding the grid element of origin.

Potential Energy Due to Adhesion Versus Cell Density in Grid Element
Cells don't touch at low density: no adhesive effect Cells just barely touch at this density ! adhesion begins to take effect Lowest energy: at higher densities than this, elasticity compensates for adhesion Figure 6. The potential energy of cell configurations is used to determine cell movement. Densities are given in cells per grid element. The potential energy function is zero for densities less than n min , at which density the cells are assumed to be "just touching". At higher densities, the potential function decreases, reaching a minimum at density n eq . For densities higher than n eq , the potential function increases steeply, reaching a maximum at density, n max . In this figure, n min = 50, n eq = 80, and n max = 100.   4.1. Nutrient consumption. The model parameters describing nutrient consumption and the production of lactate were initially calibrated by checking nutrient consumption rates in a variety of micro-environments to make sure that model predictions were consistent with experimental observation. Figure 10 shows oxygen and glucose consumption rates predicted by the model. In the model, tumor cells and normal cells consume oxygen at the same rate, but this rate does depend on local glucose, oxygen and hydrogen concentrations. Since normal pH levels are around 7.4 [64], it is the top left graph in each set that are the most significant to our current study. The upper set of graphs show glucose consumption rates as a function of oxygen concentration at four different levels of acidity. Consumption rates for tumor cells (solid lines) and normal cells (dashed lines) are drawn on the same axes for comparison. We see that the difference in consumption rates as most marked at pH levels close to normal, and when oxygen concentrations are low. In these graphs, the glucose concentration is fixed at 5.5 mM. The lower set of graphs shows oxygen consumption rates as a function of oxygen concentration at two different glucose concentrations and four different levels of acidity. Here we see that at lower glucose levels oxygen consumption rates are significantly higher than at normal glucose levels, consistent with [14]. As acidity increases, this difference diminishes.

4.2.
Tumor growth and micro-environment. In this section we describe the results of some experiments with the model. Figure 9 shows the results of a simulation run for 365 days on a 300 by 300 grid with randomly distributed blood vessels beginning with a very small initial tumor population. Normal cell populations, shown on the left, are unaffected away from the tumor. Consistent with diffusion-limited growth, the tumor has developed a necrotic core and a proliferating rim. The center and right panels show that the tumor has become necrotic in the center.   Figure 11. Simulation showing the final tumor after 400 days (left) and total tumor population over time (right). In this case, tumor growth is diffusion limited, and the tumor will grow no further unless new blood vessels are formed.
of the tumor on the acidity of the surrounding tissue: far from the tumor, pH levels approach normal, but nearby and in the center of the tissue, the environment has become relatively acidic. Figure 11 illustrates that the inability of nutrients to penetrate the core of the tumor limits its ability to grow beyond a certain size, measured in terms of total population. In this simulation, adhesion parameters and the distribution of blood vessels are kept constant. If the tumor cells have the ability to encourage the growth of new vasculature (angiogenesis) or to reduce cell-cell adhesion, then the tumor has the ability to grow beyond this limiting population. These preliminary simulations suggest that the model adequately captures these features of tumor growth, so that mechanisms involved in the progression of the disease can then be explored.
All of our results reflect the assumption that tumor cells are able to metabolize glycolytically more readily than normal cells: glycolysis produces more lactate and is less efficient than aerobic metabolism, and so the O 2 concentrations and the pH of the tissue surrounding the tumor is decreased, and the tissue occupied by the tumor cells becomes hypoxic and more acidic. These changes in the micro-environment have consequences not only for the normal cells and the overall energy budget of the host, but can also affect the delivery and bio-distribution of potential treatments, a factor that we propose to explore with this model. We note that these results are consistent with observations from both in vitro and in vivo studies (e.g. [58], [38], [32]).

4.3.
The effect of glucose levels and adhesivity. Figure 12 shows that tumor growth is enhanced by increased glucose levels in the blood vessels. Additionally, a loss in adhesivity can enable the cells to migrate to more favorable regions of tissue, thereby escaping the diffusive limit. Figure 13 shows the available ATP at the end of the same simulations. The white contour, visible in the graph of ATP levels under "normal" conditions, encloses a cross-sectional area of tumor that is unable to sustain cell life. The contour is not visible in the other two graphs, since cells either have more glucose available (middle graph) or are able to migrate more freely towards areas of higher nutrient concentration (rightmost graph).   Figure 13. Graphs of the available ATP over the simulation grid at the end of the simulations shown in Figure 12. The contour, visible in the leftmost graph, encloses a cross-sectional area that is unable to sustain cell growth.
5. Discussion and future directions. We have developed a two-dimensional hybrid cellular automaton model of tumor growth in tissue that incorporates the effects of the cancer cells on the micro-environment. In particular, the model reflects a cancer cell's ability to metabolize glucose anaerobically and to survive in an acidic environment that can be hostile to normal cells. By calculating the energy available to cells locally, the model predicts that tumor cells have the advantage over normal cells and are able to invade surrounding tissue under certain conditions. The model has the ability to describe the effects of cell-cell adhesion on the invasiveness of the tumor, and it can also describe vascular constriction and collapse due to pressure from the growing tumor.
The goal of the model is to provide a computationally efficient framework for the investigation of the growth of specific tumors in tissue, of the mechanisms underlying cell metastasis, and of the effects of certain cancer treatments, including combination treatments made up of both chemotherapy and immunotherapy. The model is detailed enough to be calibrated to a specific tumor-cell type, once data describing cell metabolism, tissue vascularization, and tissue density is provided. Once the calibration is completed, the effects of certain treatments can be simulated. For example, the administration of a cytotoxic drug can be modeled by adding a state variable describing the delivery and diffusion of the drug through the tissue once the diffusion constant of the drug is known as well as its cytolytic potential. We have developed equations describing fractional cell kill by a chemotherapeutic drug in [19], [20] and [54], while we have modeled immunotherapy and combination therapy in [21], [22] and [17]. The equations describing the effects of the drug on the tumor cells and normal cells, as well as the interactions between immune cells and tumor cells can be used in this CA model by using them as local evolutionary rules. Recruitment and proliferation of immune cells must also be added to the current model. Some of the ideas from the authors' previous work in [18] describing a cellular automaton model can be used, although the probabilistic rules in that earlier model will be translated to deterministic laws, in line with this model's attempt to rely on tissue heterogeneity and the nonlinear nature of the interactions to generate the complicated and unpredictable behavior observed in actual tumor growth.
It is known that oxygen concentrations can be an important factor in the effectiveness of radiotherapy as well as in the progression of the disease by metastasis (e.g. [9], [12]). The model presented here could serve as a refinement of earlier models such as [5] and [50] that investigate the role of oxygen concentrations on radiotherapy and the delivery of chemotherapeutic agents. By appealing to the insights offered in these and other works, we plan to use this DEB-based model to connect the microscopic effects of therapy to overall host fitness.
The simulation could be extended to model tumor development in threedimensional tissue rather than an effectively two-dimensional tissue slice. This would require several modifications to simulation data structures and an adjustment to the Dirty Diffusion computation. Dirty Diffusion as described in this paper is a matrix multiplication-based operation; since matrix multiplication is defined only for two-dimensional matrices, running a single matrix multiplication as described in §B.0.1 would no longer be possible. Instead, the three dimensional n-by-n-by-n nutrient concentration matrix N would be treated as a stack of n two-dimensional matrices [N 0 , N 1 , . . . , N n−1 ]. Dirty diffusion would then be accomplished for each matrix N x by performing the following computation: Thus, addition of adjacent grid elements within a two-dimensional slice is accomplished using matrix multiplication as before, whereas addition of neighboring elements in adjacent slices is accomplished using matrix addition.
Appendix A. Parameter values.   Appendix B. Implementation of boundary rules. We have implemented three different boundary rules, illustrated in Figure 14. In general, simulations with periodic grid topology or zero-flux boundaries demonstrate biologically realistic tumor growth patterns. Simulations with fixed boundary state demonstrate biologically realistic growth only if the boundary's fixed nutrient concentrations are appropriately tuned. The three boundary rules can be implemented in similar ways and have similar computational costs. For completeness, we include here a discussion of the methods used to implement each type of boundary. It is certainly possible to perform simulations where the different boundary rules are combined should that be desirable. B.0.1. Periodic Grid Topology. In simulations with periodic grid topology, edge elements are assumed to be adjacent to the element at the opposite grid boundary; nutrient diffusion and cell movement between "neighboring" elements therefore occur across grid boundaries. Diffusion and cell movement in periodic grids is implemented using multiplication with two specific matrixes, A and B, shown in Figure 15. Figure 15 illustrates the simulation of one timestep of nutrient diffusion in an nby-n periodic grid. Grid-wide nutrient concentrations are stored in an n-by-n matrix N , with each entry containing the nutrient concentration in the corresponding grid element. To apply the diffusion rule in equation 8, the matrix multiplication A × N is carried out; as shown in Figure 15, each entry (r, c) in the resulting matrix selfPlusVertNeighbors contains the sum of the entries (r, c), ((r + 1) mod n, c), and ((r − 1) mod n, c), i.e., each entry contains its own nutrient concentration plus that of its vertical neighbors. Values for entries in which r = 0, the top row, are added to values in the second row and the bottom row; this multiplication therefore enforces vertical periodic grid topology. Figure 15. Implementation of periodic grid topology. N is a cellular automaton grid containing nutrient concentrations. The matrix multiplication A × N yields the matrix selfPlusVertNeighbors. The multiplication N × B yields horizontalNeighbors. These two matrices are then added, and the resulting matrix is divided by 5 to complete the diffusion calculation.
Similarly, the matrix multiplication N × B is then calculated to produce matrix horizontalNeighbors, which contains the sum of each grid element's horizontal neighbors while enforcing horizontal periodic grid topology. The two product matrices are then added, and their results are divided by 5 to produce the final result, in which each element contains the average nutrient concentration of it and its four adjacent elements, with element adjacency defined using periodic grid topology. The calculation of nutrient diffusion with periodic topology is summarized in equation 20: Where N t is the matrix of gridwide nutrient concentrations at timestep t.
B.0.2. Fixed Boundary State. The implementation of fixed boundary state is based on that of periodic grid topology. To simulate diffusion with fixed boundary state on an n-by-n grid, a list edgeVals of n precomputed nutrient concetrations is stored; this list specifies the fixed state of the grid boundary. equation 20 is then used to calculate matrix N t+1 from N t as if the grid had periodic topology. The resulting matrix N t+1 has invalid values in all edge elements, because the grid boundary state should be fixed; however, all other elements have valid nutrient concentrations. Finally, list edgeVals is copied into each of the four grid edges, restoring the fixed boundary state and completing the diffusion calculation. In order to produce a biologically realistic simulation, edgeVals must be accurately tuned. If the nutrient concentrations in edgeVals are high, the simulated tumor will grow outward toward the grid edges, which contain a plentiful, inexhaustible supply of nutrients. If the nutrient concentrations in edgeVals are low, grid-wide nutrient concentration levels may decrease below the threshold required for tumor survival.
Values in list edgeVals were set for an n-by-n grid using the following procedure. A tumor-free simulation was run on an n-by-n grid G with zero-flux boundaries for 500 timesteps, in order to allow grid-wide nutrient concentrations to reach a steady state in the absence of tumor cells. The nutrient concentrations in the edges of G were then used as the values for edgeVals in an n-by-n tumor simulation with fixed boundary state. B.0.3. Zero-Flux Boundary. In simulations with a zero-flux boundary, grid elements do not all have the same number of neighbors -edge and corner elements have fewer neighbors than central elements. Thus, the denominator of the diffusion rule in equation 8 is no longer the constant value 5; it is replaced by a value from the matrix numNeigh, illustrated in Figure 16: N t+1 (r, c) = N t (r, c) + (i,j)∈nbhd(r,c) N t (i, j) numNeigh(r, c) + 1 (21) Where N t (r.c) is the nutrient concentration in element (r, c) at time t, numNeigh(r, c) is element (r, c) of the numNeigh matrix, and nbhd(r, c) is the set of all elements (i, j) that are adjacent to element (r, c). For example, nbhd(0, 0) = {(1, 0), (0, 1)}. Because each element N t+1 (i, j) receives the average concentration in N t (i, j) and nbhd(i, j), the denominator in equation 21 must be equal to |nbhd(i, j)| + 1. With periodic grid topology, |nbhd(i, j)| = 4 for all grid elements; however, |nbhd(i, j)| varies with location in a grid with zero-flux boundaries, requiring the use of numNeigh. Matrix multiplication is used to calculate total neighborhood nutrient concentrations with zero-flux boundaries by using the matrices A and B , as shown in Figure 17.
B.0.4. Alternative Implementation. The previous sections describe an implementation of the three boundary rules based on matrix multiplication. Alternatively, a looping-based approach can be used, in which Equations (8) or (21) are directly applied to individual elements. This approach requires a separate computation for each grid element, so its computational complexity is in O(n 2 ) on an n-by-n grid.  . Implementation of zero-flux boundaries. N is a cellular automaton grid containing nutrient concentrations. The matrix multiplication A ×N yields the matrix selfPlusVertNeighbors. The multiplication N ×B yields horizontalNeighbors. These two matrices are then added, and the resulting matrix element-wise divided by matrix numNeigh to complete the diffusion calculation.
In contrast, the computational complexity of the matrix multiplication-based approach is in O(n 3 ); however, in practice, the matrix multiplication-based method performs much faster than the looping-based method.