A Review of Cell-Based Computational Modeling in Cancer Biology

Cancer biology involves complex, dynamic interactions between cancer cells and their tissue microenvironments. Single-cell effects are critical drivers of clinical progression. Chemical and mechanical communication between tumor and stromal cells can co-opt normal physiologic processes to promote growth and invasion. Cancer cell heterogeneity increases cancer’s ability to test strategies to adapt to microenvironmental stresses. Hypoxia and treatment can select for cancer stem cells and drive invasion and resistance. Cell-based computational models (also known as discrete models, agent-based models, or individual-based models) simulate individual cells as they interact in virtual tissues, which allows us to explore how single-cell behaviors lead to the dynamics we observe and work to control in cancer systems. In this review, we introduce the broad range of techniques available for cell-based computational modeling. The approaches can range from highly detailed models of just a few cells and their morphologies to millions of simpler cells in three-dimensional tissues. Modeling individual cells allows us to directly translate biologic observations into simulation rules. In many cases, individual cell agents include molecular-scale models. Most models also simulate the transport of oxygen, drugs, and growth factors, which allow us to link cancer development to microenvironmental conditions. We illustrate these methods with examples drawn from cancer hypoxia, angiogenesis, invasion, stem cells, and immunosurveillance. An ecosystem of interoperable cell-based simulation tools is emerging at a time when cloud computing resources make software easier to access and supercomputing resources make large-scale simulation studies possible. As the field develops, we anticipate that high-throughput simulation studies will allow us to rapidly explore the space of biologic possibilities, prescreen new therapeutic strategies, and even re-engineer tumor and stromal cells to bring cancer systems under control.


INTRODUCTION
Cancer is a complex systems problem that involves interactions between cancer cells and their tissue microenvironments. [1][2][3] Therapeutic approaches that narrowly focus on cancer cells frequently lead to disappointing outcomes, including resistance, tissue invasion, and treatment failure. Such failures are partly due to the unexpected behaviors that emerge from the dynamical systems of cancer tissues. Therapies act as selective pressures, even while cancer cells use increased genetic variability to broadly sample survival strategies and adapt. 3,4 Chronic hypoxia, another selective pressure, leads to metabolic changes, selection for cancer stem cells that resist treatment, invasion, and angiogenesis. [4][5][6] Tumor cells communicate biochemically and biomechanically with stromal cells, which allows them to co-opt normal physiologic processes. [1][2][3]7,8 Mathematical models can serve as "virtual laboratories" with fully controlled conditions where scientists and clinicians can investigate the emergent clinical behaviors that result from basic cell hypotheses and can evaluate new therapeutic strategies. 1,9 This review surveys cell-based methods for simulating cancer. Also known as discrete models, agent-based models, or individual-based models, cell-based models simulate individual cell behaviors within tissue environments. These models have several advantages. Each cell agent can track a fully independent state with individual parameters that reflect heterogeneity in cancer. Modelers can directly implement cell rules that reflect observations of single-cell behavior and cell-cell interactions, which allow us to translate biologic hypotheses to mathematical rules quickly; run simulation experiments that explore the emergent behaviors of these hypotheses; and compare against new data to confirm, reject, or iteratively improve the underlying hypotheses. 1,9,10 along a rigid grid and off-lattice models that have no such restriction. Figure 1 classifies most cell-based modeling approaches. Table 1 lists major open source modeling packages.

Lattice-Based Methods
Lattice-based models can use regular structured meshes (eg, Cartesian 11 [two-or three-dimensional [2D/3D], dodecahedral [3D]) 12 or unstructured meshes. 13 Structured meshes are simpler to implement, visualize, and combine with partial differential equation (PDE) solvers, but their structure can lead to grid biases. 13 Unstructured meshes can avoid these issues 13 but with greater complexity.
We can further categorize lattice-based methods by their spatial resolution. In cellular automaton (CA) models, each lattice site can hold a single cell. [14][15][16][17] At each time step, each cell is updated with discrete lattice-based rules: remain, move to a neighboring lattice site, die (free a lattice site), or divide to place a daughter cell in a nearby site. [14][15][16][17] These methods usually update the lattice sites in a random order to reduce grid artifacts. 14,15 In lattice gas CA (LGCA) models, a single lattice site can contain multiple cells. 14,15,17,18 LGCA models track the number of cells that move through channels between individual lattice sites rather than the motion of each individual cell. They can simulate very large numbers of cells efficiently over long periods while also connecting to statistical mechanics theory; this facilitates analysis and provides a bridge to continuum methods that model cell densities or populations instead of single cells. 17,18 Some problems may require resolution of individual cell morphologies. Cellular Potts models (CPMs) use multiple lattice sites to represent each cell. 14,15,19 At each time step, CPMs visit each pixel (2D) or voxel (3D), test a random swap with a neighboring pixel/voxel, and accept or reject the swap (probabilistically) on the basis of whether it would reduce a global energy. Although CPMs can model cell morphologies and mechanics that cannot be incorporated in CA models, they are much more computationally intensive. Also, the calibration of Monte Carlo steps to physical time can be challenging. 20

Off-Lattice Methods
We can divide off-lattice models into center-based models (CBMs) that focus on cell volumes (or masses) and models that focus on cell boundaries. We can further classify these approaches by level of morphologic detail.
CBMs. CBMs track each cell's center of mass or volume, typically by using a single software agent per cell. [13][14][15]21 Some CBMs represent cells as points, whereas others explicitly model cell volumes. CBMs typically update the cells' positions by explicitly formulating the adhesive, repulsive, locomotive, and drag-like forces exchanged between cell centers. [13][14][15]21 Most CBMs approximate cells as spheres; however, some approximate cells as deformable ellipsoids to better represent their morphologies. 22,23 CBMs can model cell morphology in greater detail by breaking cells into subcellular elements 24,25 : Each cell is represented by multiple center-based agents that interact with adhesive and repulsive forces. These models better approximate cell biomechanics but at increased computational cost. Conversely, cells can be organized into clusters or functional units (eg, breast glands or colon crypts) that are simulated as agents that interact by mechanical forces or  other rule-based motions 26,27 ; this allows modelers to incorporate heterogeneous details into individual clusters of cells but with greater computational efficiency than traditional CBMs.

Off-Lattice
Boundary-tracking models. Vertex-based methods (eg, Fletcher et al 28 ) model cells as polygons (2D) or polyhedra (3D) and compute the forces that act on their vertices; they are particularly useful for modeling confluent tissues. 29 For greater spatial resolution, front-tracking methods, such as the immersed boundary method (IBM), solve PDEs for fluid flow inside and between cells and then advect boundary points along the cells' membranes in this flow. 30 Level set methods have been applied to implicitly track the movement of cell boundaries, 31 and VCell (see Connecting to Molecular Effects) recently added front-tracking capabilities. 32,33 These are among the most computationally intensive cell-based methods, but they are useful for coupling detailed cell mechanics to fluid and solid tissue mechanics.

Connecting to Molecular Effects
Most cell-based models are hybrid discrete-continuum; they couple a discrete cell model to continuum models of the microenvironment. 1,14,15 In general, these models use reaction-diffusion PDEs to simulate biotransport of oxygen, growth factors, and drugs. Ghaffarizadeh et al 34 developed BioFVM to solve diffusive transport of tens to hundreds of chemical substrates in 3D tissues; it is the underlying PDE solver for PhysiCell (a center-based simulation framework). 21 In this framework, modelers write rules to relate individual cell phenotypes to local chemical substrate conditions. 21 Many discrete models include systems of ordinary differential equations (ODEs) to model molecular processes in individual cells. 35,36 VCell can simulate reacting flows of many proteins within a single detailed cell, 32,33 and many modeling packages (eg, Chaste, 37 CompuCell3D, 38 and EPISIM 39 ) support systems biology markup language (SBML) to include systems of ODEs that simulate molecular effects in individual cells. Others use discrete models within individual agents: Gerlee and Anderson 40 used small neural networks to simulate individual cell phenotypic "decisions" on the basis of microenvironmental inputs, whereas PhysiBoSS 41 combines the Boolean network modeling approach of MaBoSS 42,43 with PhysiCell 21 to simulate molecular processes in individual cells.

EXAMPLES OF CELL-BASED MODELING IN CANCER BIOLOGY
We now explore a series of modeling themes that illustrate the use of cell-based modeling in cancer biology. Although we cannot comprehensively review all cell-based modeling in cancer (or even sample all major use cases for cell-based modeling), these themes are drawn from across the field to demonstrate scientific problems with significant cell-scale effects where cell-based models can yield new insights.

Hypoxia in Breast Cancer
Many groups have used cell-based models to investigate tumor growth in hypoxic tissues and more generally, the effect of diffusive transport limits. Gatenby et al 44 and Smallbone et al 45 used CAs to examine hypoxia-driven switching to invasive phenotypes in ductal carcinoma in situ (DCIS). They incorporated cellular metabolic adaptations to hypoxia, which allowed them to study early tumor invasion (Fig 2A). Anderson and colleagues 46,47 extended earlier CA results by adapting IBCell 30 (an IBM) to mouse mammary (EMT6/Ro) tumor cell proliferation in hypoxic tissues. As before, they found that hypoxic gradients could drive tissue invasion, but IBCell's improved modeling of cell adhesion and biomechanics predicted more rounded invasive tips 47 ( Fig 2B).
Macklin et al 50 and Hyun and Macklin 51 applied a CBM to study oxygen-driven proliferation and necrosis in solid-type DCIS with comedonecrosis. After calibrating to individual patient pathology data (tissue specimens immunostained for the Ki67 protein to detect cycling cells, cleaved caspase 3 to detect apoptosis, and annotated with viable rim sizes and cell density 50 ), they were able to simulate comedonecrosis and microcalcifications as emergent properties of the simulations along with realistic, constant rates of tumor advancement along the breast ducts. Ghaffarizadeh et al 21 refined the DCIS model and extended it to 3D as well as simulated the hypoxic interiors of hanging drop spheroids calibrated to match MCF-10A birth and death kinetics in culture (Figs 2C and D). As in early 3D work by Drasdo and Höhme 48 on EMT6/Ro cells (Fig 2E), they predicted a layered structure-an outer proliferative rim surrounding a quiescent perinecrotic region and an interior necrotic core. They were the first to predict networks of fluid-filled pores in the necrotic cores that emerge from the competing effects of necrotic cell shrinking and adhesion; these structures are observed in experimental models (Fig 2D  inset). Szymańska et al 49 used a CBM of EMT6 cells to simulate a growing tumor cord-a solid tumor that grows around a blood vessel. They predicted a similar three-layer structure but in reverse order-a proliferating core nearest the blood vessel, quiescent interior, and necrotic exterior ( Fig 2F).

Tumor-Induced Angiogenesis and Drug Delivery
Tumor-induced angiogenesis allows lesions to grow to clinically detectable sizes. 3 McDougall and colleagues 52,53 modeled sprouting angiogenesis with a CA model of vessel tip migration: Sprout tip agents followed chemotactic and haptotactic signals to migrate toward hypoxic tumor regions and left a trail of functional vessels. They incorporated a detailed vascular network flow model, including dynamic wall shear stress rules for vessel branching and anastomosis (vessel looping), and used this framework to explore therapeutic delivery from tumor-associated vasculatures ( Fig 3A). Bauer et al 54 used a CPM to simulate tumor-induced angiogenesis, by adding a detailed microenvironment, including extracellular matrix (ECM) and multiple vascular endothelial growth factor isoforms; they concluded that variations in the spatial distributions of proangiogenic factors greatly affect capillary morphology and that inhomogeneities in nonvascular tissue naturally lead to capillary anastomosis. Boas and Merks 55 used a CPM to investigate novel hypotheses on cell overtaking: Cell-cell biomechanical and chemical communication can cause endothelial cells in the stalk to assume the role of migrating tip cells by migrating to the front of an advancing vessel (Fig 3B). Shirinifard et al 56 used a CPM to investigate tumor growth with angiogenesis and showed that tumor size increases with increasing angiogenesis and that the tumors grow along the vasculature (Fig 3C).
Cai et al 57 used a CA model of tumor cells in a continuous ECM coupled with a discrete angiogenesis model that included flow effects and substrate perfusion from the vasculature ( Fig 3D). They showed that the final vessel configuration depends on emergent, dynamic feedback Cell-Based Computational Modeling in Cancer mechanisms in vascular remodeling rather than on initial conditions. Wu et al 59 extended an earlier hybrid discretecontinuum model 58 (that was based on that of McDougall and colleagues 52,53 ) to investigate the influence of interstitial fluid pressure, interstitial fluid flow, and lymphatic drainage on drug delivery in growing tumors. 59 They found that elevated interstitial hydraulic conductivity and high interstitial fluid pressure limit the transvascular delivery of nutrients and therapeutics (Fig 3E).

Cancer Stem Cells
Models of cancer stem cells (CSCs) offer valuable insights into the driving forces of cancer biology. Fletcher and colleagues 37,60,61 developed a 3D CBM of colonic crypts to explore the role of stem cells (in the bottom of the crypt) in colorectal carcinogenesis. Neighboring cells were connected by linear springs, and stem-cell division and differentiation were driven by Wnt gradients along the crypt axis. The geometry of the stem-cell hierarchy (proliferation at the crypt base, expansion and differentiation along the middle and top) created an overall base-to-top proliferative cell flux. This flux has an anticancer protective effect wherein it pushes any mutated cell and its progeny out of a crypt before they can spread throughout a crypt, unless the mutation occurs in a stem-cell niche (Fig 4A).
Norton et al 62 built a 3D CA model to examine the interaction between triple-negative breast cancer and stromal cells. Stem cells proliferated and differentiated into progenitor cells, and cancer cells exchanged chemical signals with fibroblasts and infiltrating macrophages. Among their results, they found that increasing the stromal effect on cancer cell proliferation decreased overall tumor size, whereas increasing the stromal effect on cancer cell migration increased tumor size (Fig 4B).
Poleszczuk et al 66 developed a 2D CA model of CSCs and nonstem cancer cells, which tracked four traits in each individual cell: migration rate, apoptosis, symmetric CSC division, and cancer cell proliferation potential. They found that increasing the cancer cell proliferation potential could reduce tumor growth because the increased cancer cell population competed with CSCs for space and inhibited CSC division. They also found that traits propagated radially from the centers of growing tumors; this has implications for biopsies of tumor heterogeneity (Fig 4C). Gao et al 63 used a CPM to investigate the role of glioma stem cells (GSCs) in glioblastoma growth and radiation therapy response (Fig 4D). They found that switching from asymmetric to symmetric division or fast GSC cycling was necessary to explain clinical observations of glioma repopulation after radiotherapy and that the expanded GSC fraction could reduce radiosensitivity.
Alfonso et al 67 also explored radiotherapy treatment paradigms with respect to a heterogeneous population of CSCs and cancer cells using a 3D CA model. They found that CSCs, which are typically more radioresistant, segregated to the center of the tumor across a range of proliferation and death parameters. This emergent phenomenon is due to the faster cycling time of the cancer cells compared with CSCs. When these cell arrangements were subjected to radiotherapy, they found that radiotherapy is more effective at tumor control when it is concentrated on the tumor center where CSCs are located rather than when it is spread homogeneously across the entire tumor.

The "Go or Grow" Hypothesis for Glioblastoma Multiforme
Tektonidis et al 68 examined the "go or grow" (GoG) hypothesis in gliomas, where tumor cells must make a "decision" between migration (go) or proliferation (grow). They modeled data from 3D spheroid cell cultures with a 2D LCGA model and attempted to recapitulate three experimental observations: nonidentical spreading rates of the invasive rim and central core, radially persistent and symmetric cell motion, and a highly proliferative central core compared with the remaining tumor. Tektonidis et al evaluated the emergent model behavior under a variety of cell phenotype rule sets to determine which rules were required to predict the three observations. They found that a proliferative-motility dichotomy (the GoG hypothesis), cell-cell repulsion, and density-dependent switching between the proliferative and motile states were required to match experimental observations. They concluded that disruption of the GoG mechanism to favor proliferation could limit the required tumor resection volume in surgical interventions.
Hatzikirou et al 64 investigated the GoG hypothesis in glioblastoma multiforme (GBM) using a 2D LGCA model. In their work, cells could divide, re-orient, migrate, or apoptose on the basis of local oxygenation conditions. They modeled the GoG hypothesis by switching hypoxic cells to a motile phenotype and reverting to a proliferative phenotype after escaping hypoxia. They found that increasing the cell bias toward proliferation increases overall tumor growth, but crossing a threshold could decrease overall tumor growth when motility is insufficient to open new space for cell division (Fig 4E). Comparable results were obtained by Gerlee and Nelander 69 using stochastic switching between the two phenotypes (migrate or proliferate). From their CA model, they derived a system of coupled PDEs to investigate further the relationships between cell-level parameters and tumor-scale dynamics.
In related work, Böttger et al 70  Kim et al 65 used a CBM to explore GBM using miR-451 as an intracellular detector of glucose, with an ODE model of miR-451 as the effector for selecting migration versus proliferation for glioma cells (Fig 4F). They found that cell migration depended not only on glucose, but also on mechanical spacing between cells. They also predicted that the placement of chemoattractants at the edges of a resected tumor could reduce GBM cell migration from the resection site. These GoG examples highlight the key role played in single-cell decisions in GBM and the potential for cell-based models that explore the clinical behaviors that emerge from single-cell effects. 72

Cancer Invasion and Epithelial-Mesenchymal Transition
Cancer invasion is essential to metastatic progression. 3,[72][73][74] Cancer cells acquire a motile phenotype to escape primary tumors and invade nearby tissues (as conceptually modeled in epithelial-to-mesenchymal transition [EMT] 75 ), invade and travel within blood and lymphatic vessels, 76,77 and finally colonize distant metastatic niches. 78 Because single-cell effects are critical, many cell-based models have investigated cancer invasion with or without explicit modeling of EMT.

Reher et al 79 developed a 2D
LGCA model to simulate the effects of a heterogeneous cell-cell adhesion in an epithelial layer, with decreased cell-cell adhesion representing one of the effects of EMT. Cell-cell adhesion was modeled by varying the initial and maximum number of adhesion receptors for each virtual cell. They found that increased adhesion heterogeneity as well as decoupling receptor number from environmental signals (cell-cell contact) lead to increased dissemination.
Kim and Othmer 80 developed a center-based model with a continuous ECM description to investigate the interactions of tumor cells, signal-secreting stromal cells, and the ECM. Stromal levels of fibroblast-secreted protein initiated EMT, which switches cells to an invasive, motile phenotype. Invasive cells also secreted a tumor-associated protease, which degrades the basement membrane and ECM. Degradation of the basement membrane results in more cell exposure to fibroblast-secreted protein, which results in more phenotypic switching; more invasion; more ECM degradation; and ultimately, collective cell invasion ( Fig 5A). These results suggest that inhibiting fibroblast secretions could affect invasive potential.
Zhang et al 81 used a 3D CA model to investigate tissue invasion and tumor cell heterogeneity in glioma. Using a subcellular signaling network, cells divide and move on the basis of external concentrations of nutrients and signaling molecules. They acquire oncogenic mutations upon division to produce new clones. Each successive clone is more proliferative and nutrient seeking, which increases overall invasive potential. By starting with a sphere of the least oncogenic cells in the center of the simulation, Zhang et al found that heterogeneity increased in all regions of the simulation followed by a decrease in heterogeneity in the regions that contain the nutrient source and eventual recovery of heterogeneity. The decrease was attributed to an outgrowth of the most oncogenic subclone, which potentially explains the asymmetric invasive growth seen in experimental and clinical reports (Fig 5B).
Anderson et al 82 used an IBM to investigate morphologic changes in breast acini as a result of variation in cell polarization and anoikis behaviors in response to cell contact with other cells and the basement membrane. In this, and additional work by Rejniak et al, 83 Anderson et al showed that atypical behavior in a single cell can lead to eventual intraductal and stromal invasion (Fig 5C).
In a sophisticated investigation of metastatic colonization, Araujo et al 86 developed a hybrid CA model of prostate cancer (PCa) bone metastasis. In their work, mesenchymal stromal cells differentiated into osteoblasts that created bone material, whereas osteoclasts degraded bone. In the absence of tumor cells, these cell populations coordinated through transforming growth factor-β and receptor activator of nuclear factor kappa beta ligand (RANKL) signaling to maintain healthy bone tissue. When PCa cell agents were introduced into the bone, tumor-secreted transforming growth factor-β promoted osteoblast activity, but osteoblastosteoclast feedbacks simultaneously degraded bone to release new growth factors that could drive additional PCa proliferation and invasion. Araujo et al studied the model to compare potential improvements of RANKL inhibitors and bisphosphonates (two standard-of-care treatments for bone metastases). They found that additional refinements to bisphosphonates would yield little clinical improvement over current treatments, whereas improvement of RANKL inhibitors from the current (model-estimated) 40%  80 (B) A three-dimensional (3D) cellular automaton (CA) model to study selection in heterogeneous brain cancers. Cells could mutate their signaling network parameters, which leads to more invasive clones. Adapted with permission from Zhang et al. 81 (C) An immersed boundary method of contact-based signaling and inhibition closer to a theoretical maximum 100% inhibition could dramatically improve outcome.

Tumor Immunosurveillance
Individual interactions between tumor cells and immune cells are critical to immunosurveillance. 3,7,87 Cell-based models are uniquely capable of examining these while considering the roles of stochasticity and heterogeneity.
Kather et al 84  Ghaffarizadeh et al 21 developed a 3D off-lattice model of an immune attack on a tumor with a heterogeneous oncoprotein (mutations increased both proliferation and immunogenicity). After simulating initial growth, they added simulated immune cells that chemotaxed toward tumorreleased immunostimulatory factors. Each immune cell tested for mechanical collision with cells, formed an (Hookean) adhesion, tested for immunogenicity, and stochastically attempted to induce tumor cell apoptosis. Attached immune cells eventually would succeed in killing the tumor cell and detach or continue the attempt before detaching and continuing their search for new targets. Their simulations showed initial tumor regression but eventual tumor regrowth when immune cells passed some tumor cells and formed large clumps near maxima of the immunostimulatory factor (Fig 5F). The authors have expanded their investigation to supercomputers to explore further the effect of stochastic migration on the overall efficacy of the immune response 10 ; this work showcases the potential for using supercomputers to explore large therapeutic design spaces in high throughput.

ADDITIONAL RESOURCES
As companions to this review, we maintain three online curated collections: Lattice-based methods are straightforward to implement and fast, which makes them useful for quickly testing new ideas. We show how lattice-based methods have yielded insights on cancer metabolism, cancer stem cells, angiogenesis, the GoG hypothesis, invasion, and cancer immunosurveillance. CA methods remain the most common method for modeling vascular networks, whereas CPMs have investigated the finer details of angiogenesis.
Lattice-based methods are at risk for grid-based artifacts, and their biomechanical realism is limited; off-lattice models can readily incorporate biomechanics and offlattice cell-cell interactions. As the hypoxia and immunosurveillance examples show, novel structures can arise from the interplay of stochasticity, transport limitations, mechanics, and single-cell characteristics. However, offlattice models often are computationally demanding, and they generally have many parameters to calibrate. Ultimately, no single method is best for all problems. Modelers should reproduce findings with multiple approaches to avoid algorithm-dependent biases. 88  85 (F) A 3D center-based model of immune responses to an immunostimulatory factor in a heterogeneous tumor (shaded by immunogenicity; yellow cells are most immunogenic). Immune cells (red) seek and adhere to cancer cells, test for immunogenicity, and induce apoptosis. The immune response failed after immune cells aggregated near a local maximum in the signaling factor, which allows the tumor to repopulate. Adapted with permission from Ghaffarizadeh et al. 21 This work was explored further with high-performance computing. 10 This is an exciting time to apply cell-based modeling to cancer. Increases in computational power are allowing larger simulation studies with greater sophistication, and high-throughput computing is enabling exploration of high-dimensional parameter spaces. 10 Open source platforms have lowered the barrier to entry for using sophisticated techniques (Table 1). In the future, we envision that cell-based modeling software will be increasingly user friendly, cloud hosted, and open to modular contributions from the community, 88 which would potentiate a community-driven ecosystem of interoperable tools that together exceed the sum of their parts. 1 We are excited to imagine the new insights that are on the horizon.