A Review of Cell-Based Computational Modeling in Cancer Biology

1Cancer is a complex systems problem that involves interactions between cancer cells and their tissue microenvironments.- 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., Chronic hypoxia, another selective pressure, leads to metabolic changes, selection for cancer stem cells that resist treatment, invasion, and angiogenesis.- Tumor cells communicate biochemically and biomechanically with stromal cells, which allows them to co-opt normal physiologic processes.-,, 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.,

2This 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.,,

3Cell-based models represent individual cells with two main paradigms—lattice-based models that track cells 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.

4Lattice-based models can use regular structured meshes (eg, Cartesian [two- or three-dimensional [2D/3D], dodecahedral [3D]) or unstructured meshes. Structured meshes are simpler to implement, visualize, and combine with partial differential equation (PDE) solvers, but their structure can lead to grid biases. Unstructured meshes can avoid these issues but with greater complexity.

5We can further categorize lattice-based methods by their spatial resolution. In cellular automaton (CA) models, each lattice site can hold a single cell.- 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.- These methods usually update the lattice sites in a random order to reduce grid artifacts.,

6In lattice gas CA (LGCA) models, a single lattice site can contain multiple cells.,,, 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.,

7Some problems may require resolution of individual cell morphologies. Cellular Potts models (CPMs) use multiple lattice sites to represent each cell.,, 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.

8We 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.

9CBMs track each cell’s center of mass or volume, typically by using a single software agent per cell.-, 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.-, Most CBMs approximate cells as spheres; however, some approximate cells as deformable ellipsoids to better represent their morphologies.,

10CBMs can model cell morphology in greater detail by breaking cells into subcellular elements,: 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,; this allows modelers to incorporate heterogeneous details into individual clusters of cells but with greater computational efficiency than traditional CBMs.

11Vertex-based methods (eg, Fletcher et al) 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. 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. Level set methods have been applied to implicitly track the movement of cell boundaries, and VCell (see Connecting to Molecular Effects) recently added front-tracking capabilities., 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.

12Most cell-based models are hybrid discrete-continuum; they couple a discrete cell model to continuum models of the microenvironment.,, In general, these models use reaction-diffusion PDEs to simulate biotransport of oxygen, growth factors, and drugs. Ghaffarizadeh et al 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). In this framework, modelers write rules to relate individual cell phenotypes to local chemical substrate conditions.

13Many discrete models include systems of ordinary differential equations (ODEs) to model molecular processes in individual cells., VCell can simulate reacting flows of many proteins within a single detailed cell,, and many modeling packages (eg, Chaste, CompuCell3D, and EPISIM) 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 used small neural networks to simulate individual cell phenotypic "decisions" on the basis of microenvironmental inputs, whereas PhysiBoSS combines the Boolean network modeling approach of MaBoSS, with PhysiCell to simulate molecular processes in individual cells.

14We 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.

15Many 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 and Smallbone et al 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, extended earlier CA results by adapting IBCell (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 (Fig 2B).

16Macklin et al and Hyun and Macklin 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), 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 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 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 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).

17Tumor-induced angiogenesis allows lesions to grow to clinically detectable sizes. McDougall and colleagues, 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 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 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 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).

18Cai et al 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 mechanisms in vascular remodeling rather than on initial conditions. Wu et al extended an earlier hybrid discrete-continuum model (that was based on that of McDougall and colleagues,) to investigate the influence of interstitial fluid pressure, interstitial fluid flow, and lymphatic drainage on drug delivery in growing tumors. They found that elevated interstitial hydraulic conductivity and high interstitial fluid pressure limit the transvascular delivery of nutrients and therapeutics (Fig 3E).

19Models of cancer stem cells (CSCs) offer valuable insights into the driving forces of cancer biology. Fletcher and colleagues,, 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).

20Norton et al 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).

21Poleszczuk et al 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 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.

22Alfonso et al 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.

23Tektonidis et al 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.

24Hatzikirou et al 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 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.

25In related work, Böttger et al explored a 2D LGCA GoG model for GBM to provide a more quantitative parameter space analysis of a tumor’s invasive dynamics. Systematically varying model parameters for proliferation and motility led to counterintuitive results about invasion. Specifically, invasion speed depended on two competing processes: emptying space as a result of cell migration and filling space as a result of cell proliferation. In later work, Böttger et al used similar techniques to find that if cell motility decreases with increased cell density, then small tumors self-extinguish. On the other hand, increasing cell motility with increasing cell density leads to self-sustaining growth, similar to the Allee effect frequently observed in ecology.

26Kim et al 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.

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

28Reher et al 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.

29Kim and Othmer 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.

30Zhang et al 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).

31Anderson et al 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, Anderson et al showed that atypical behavior in a single cell can lead to eventual intraductal and stromal invasion (Fig 5C).

32In a sophisticated investigation of metastatic colonization, Araujo et al 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 osteoblast-osteoclast 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% inhibition closer to a theoretical maximum 100% inhibition could dramatically improve outcome.

33Individual interactions between tumor cells and immune cells are critical to immunosurveillance.,, Cell-based models are uniquely capable of examining these while considering the roles of stochasticity and heterogeneity.

34Kather et al developed a CA model to study complex interactions between tumor cells and immune cells. In their model, immune cells could randomly appear (an influx model), migrate toward tumor cells, and kill them. Immune cells were assumed to be exhausted after killing a maximum number of times and would induce tissue fibrosis that impaired cell migration. Tumor cells could proliferate, die, remain stationary, or migrate on the basis of the number of open neighboring lattice sites and their distance to the nearest open lattice site (a phenomenologic model of necrosis). Depending on the relative migration and proliferation parameters for tumor and immune cells, this model could exhibit diverse tumor-immune interactions, including successful immunosurveillance (eradicated), tumor encapsulation by fibrotic tissue (immune excluded), and ongoing immune responses (Fig 5D). Gong et al recently developed a 3D CA model to investigate the tumor-immune responses to programmed cell death-1 and programmed death-ligand 1 inhibition (Fig 5E).

35Ghaffarizadeh et al 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 tumor-released 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; this work showcases the potential for using supercomputers to explore large therapeutic design spaces in high throughput.

36As companions to this review, we maintain three online curated collections:

37Cloud-hosted examples of cell-based cancer models: https://www.scienceopen.com/collection/online-cancer-simulators

38Open source cell-based simulation frameworks: https://www.scienceopen.com/collection/open-source-agent-frameworks-biology

39Additional reviews of mathematical modeling in cancer: https://www.scienceopen.com/collection/math-modeling-in-cancer-reviews

40Cell-based methods can track single-cell traits and individual behaviors, which make them well suited for problems where single-cell effects are important, such as stem-cell hierarchies, heterogeneity, invasion, and tumor-immune interactions. They are ideal for hypothesis-driven computational experiments because biologic observations of single-cell behavior can be directly translated to agent rules.

41Lattice-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.

42Lattice-based methods are at risk for grid-based artifacts, and their biomechanical realism is limited; off-lattice models can readily incorporate biomechanics and off-lattice 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, off-lattice 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.

43This 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. 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, which would potentiate a community-driven ecosystem of interoperable tools that together exceed the sum of their parts. We are excited to imagine the new insights that are on the horizon.