Current trends in mathematical modeling of tumor–microenvironment interactions: a survey of tools and applications

Introduction

In its simplest description, a tumor is comprised of an expanding population of transformed cells that is supported by a surrounding microenvironment termed the tumor stroma. Cancer begins with a single normal cell that undergoes initiation, expands in cell number through tumor promotion and progresses to an invasive stage. Throughout tumor progression, tumor cell behavior is altered due to an accumulation of both genetic and epigenetic changes that influence a variety of cellular functions, including proliferation, migration, invasion, metabolism and survival. At the same time, the tumor cells are engaged in a complex interplay with the surrounding tumor microenvironment (Figure 1). Major stroma cell types include: fibroblasts that provide a provisional matrix and source of growth factors; endothelial cells that facilitate new vasculature to meet nutritional and oxidative requirements; and immune cells that have competing antitumorigenic surveillance roles as well as tumor-promoting pro-growth, pro-angiogenic and pro-tumor invasive roles. Coordinating the entire process are biochemical cues (i.e. growth factors, chemokines, metabolites) and biophysical properties (i.e. matrix tension, compression, stromal fiber alignment) that regulate cell behavior and impact tumor progression. In order to study this complex process, cancer is dissected and considered at different scales of complexity that include: minimalist systems used to investigate the dynamics of protein–protein interactions; monolayer cell system models that focus on a single cell-type behavior; three-dimensional (3D) cell or organotypic studies that examine interactions of multiple cell types and incorporate tissue structure in a complex matrix environment; and whole animal tumor studies. Integrating this information into a common description is as complicated as the disease itself. Computational formulations of cancer are providing cancer researchers with invaluable tools that are capable of integrating this complexity into organizing principles. These mathematical models (1) enable one to integrate measures made at different levels to generate simple-to-multifaceted computational descriptions of the cancer disease process; and (2) provide tools that of themselves can be used to generate new hypotheses regarding tumor initiation, progression and treatment.

OPEN IN VIEWER

Figure 1

A scheme of tumor microenvironment components. A normal duct composed from one layer of epithelial cells and surrounded by the basement membrane (BM) is filled with a mass of tumor cells occupying the lumen. A small cohort of invasive tumor cells is breaking through the BM. Tumor microenvironment consists of: various fibers that are aligned differently depending on whether they are located near the invading tumor cells or near the intact BM; multiple stromal cells, such as fibroblasts, macrophages, lymphocytes, endothelial and lymphatic cells; and the interstitial fluid. Tumor cells are exposed to gradients of various chemical cues

The first papers laying the basis for the future mathematical models of tumor growth were published at the beginning of the last century by Hill, ¹ who studied the role of diffusion of dissolved substances, such as oxygen and lactic acid, on the physiological processes in solid tissues, and by Mayneord, ² who investigated the effects of different distributions of actively dividing cells. Mayneord demonstrated the change in growth dynamics from exponential to linear when the growth of the entire tumor volume is gradually restricted to an outer rim of tissue of decreasing thickness. This phenomenon of diffusion limited tumor growth that leads to the emergence of a layered tumor structure composed of a necrotic core, a ring of quiescent cells and a rim of proliferating cells has been reproduced by virtually every mathematical model, whether continuous or individual cell based. The class of continuous models is well suited to represent tumor growth on a tissue scale where local properties of individual cells can be averaged and are assumed to change gradually (for a review see references ^3–8 and references therein). Using these models one can investigate how the whole tumor mass evolves and interacts with the surrounding environment via diffusive (such as oxygen, nutrients, waste products) and/or mechanical (such as internal and external tissue stress and compression) factors. These models can trace changes in the overall morphology of the growing tumor tissue and can also distinguish between a limited number of cell subpopulations modeled as densities of different cell types sharing common space; however, they are more suited to represent localized averages of cell properties than specific changes in individual cells. Therefore, a class of individual-cell-based models has been introduced, in which each cell is treated as a separate entity with individually defined cell life processes, cell–cell interactions and changes in cell phenotype or genotype (for a review see references ^6–9 and references therein). This class includes models in which individual cells are represented by either single points (defined on or off lattice) or as deformable bodies composed from either several lattice sides or intracellular elements not following the lattice geometry (off lattice), or as elastic bodies filled with fluid (see Figure 2 for examples). Since in these models the properties of individual cells can be defined independently, and the fate of each cell can be traced separately, these models are easier to parameterize and validate using experimental data. Thus, individual-cell-based models are used to investigate the emergent behaviors of smaller cell populations. These types of models are less useful for modeling large cell populations that correspond to typical in vivo tumor volumes, as the computational requirement to follow each cell within the tumor mass would be quite large. A comparison of relative strengths and weaknesses of the various models is presented in Figure 3. Generally, the kind of modeling approach that should be used to address a particular biological question depends on the level of detail that needs to be included in the computational model.

OPEN IN VIEWER

Figure 2

A schematic representation of different mathematical models of multicellular systems. (a) In a continuous model one considers densities of various types of cells, such as growing, viable and necrotic, and cell properties are averaged over a mass of cells of a given type; (b) in a cellular automata model each cell is represented as a separate entity with individually regulated cell properties that occupies one site of a regular square lattice; (c) in an off-lattice center-based particle model cells are represented as individual non-overlapping points; (d) in a cellular Potts model every cell is composed of several lattice sides allowing for representation of cell shape and its deformation; and (e) in the IBCell model every cell can acquire variable shapes and is modeled as an elastic body immersed in the viscous incompressible fluid

OPEN IN VIEWER

Figure 3

A comparison of relative strengths and weaknesses of various types of mathematical models discussed in this Minireview

In this Minireview, we discuss certain components of tumor microenvironment that in our opinion play a crucial role in tumor initiation, progression and treatment. In particular, we will focus on the topological structure of the host tissue, and the physical, chemical, fluid and cellular components of the stroma. We limit our discussion to mathematical models that focus on a cell as the main modeling unit. We present the current state of cell-focused mathematical modeling of different tumor–stroma interactions, and discuss what modeling approaches could be undertaken to complement the in vitro and in vivo experiments. Because of the limited size of this Minireview, we will refer mostly to the papers that have been published in the last five years. Moreover, since mathematical models dealing with several aspects of tumor development, such as avascular growth, angiogenesis and tumor cell metabolism, have been recently surveyed in-depth, we will only briefly write about them here and refer our readers to corresponding reviews.

Host tissue structure

One of the fundamental features of tumors is that the origin of epithelial tumors begins with one initiated cell. It is technically challenging, experimentally, yet more feasible mathematically to follow this one individual cell and its subsequent progeny through all the avenues of genetic and phenotypic changes through the early stages of tumor promotion and later through progression to the invasive cancer. To model computationally the initiation and early stages of carcinogenesis and cell escape from tissue homeostasis, one needs to consider the geometry of the host tissue. For example, when modeling adenocarcinomas, the tumors arising from glandular epithelial cells, one may include a structure of the normal duct that consists of a single layer of epithelial cells enclosing a hollow lumen and surrounded by the basement membrane, a specialized type of the matrix that separates the gland from the surrounding connective tissue. These transformed epithelial cells must acquire specific properties allowing them to overcome a tightly controlled epithelial organization in order to disrupt the epithelial structure, fill the lumen, break through the basement membrane and invade the surrounding tissue to form secondary metastatic tumors.

The formation of different patterns of ductal carcinoma in situ (DCIS), and mechanisms leading to the transition from DCIS to an invasive tumor, have been investigated by several authors. Rejniak and Dillon have simulated the formation of four distinct patterns of DCIS (micropapillary, cribriform, tufting and solid) using two single-cell-based models in which epithelial and tumor cells were represented as elastic fully deformable bodies filled with viscous incompressible fluid (see Figure 2), comparable in concept to water-filled balloons, that were capable of interacting with one another via a set of cell membrane receptors. The cells could also interact with the basement membrane modeled as a stiff elastic band enclosing the epithelial duct with a fluid-filled lumen (Figure 4a). By choosing different proliferative conditions for the growing tumor cells, they classified the final patterns with respect to gradual dedifferentiation of tumor cells corresponding to epithelial-like, intermediate and tumor-like cell characteristics. ¹⁰ Dillon et al. ¹¹ and Rejniak ¹² have extended the models presented in reference ¹⁰ by introducing the background nutrient concentration that influences cell growth and/or death. Emergence of different patterns of DCIS has also been investigated by Shumate and El-Shenawee ¹³ using the hexagonal cellular automata model, in which each cell occupies one hexagonal grid site and adheres to at most six neighbors (‘honeycomb’ like). By setting specific rules of cell division regulated by contact inhibition and the expression level of the major cell junctional protein, E-cadherin, the authors were able to simulate several DCIS patterns commonly seen in clinical research. Bankhead and co-workers have proposed a model based on the cellular automata simulation in which cells of different types (luminal, myoepithelial, mammary stem cells) and the basement membrane structure are represented as points on the square lattice and can act according to predefined rules on cell proliferation, death, movement and other epigenetic changes (Figure 4b). This model examined the role of different progenitor cell hierarchies in the evolution of DCIS and its aggressiveness, and showed that the structure of progenitor hierarchy drastically effects genetic heterogeneity of DCIS and its ability to sustain aggressive cancerous growth. ¹⁴ Franks et al. ¹⁵ represented the tumor as a growing cell mass and investigated interactions between the expansion forces created by tumor cell proliferation and the stresses that develop in the basement membrane. In this model, matrix degradation decreased the basement membrane's resistance to dilation and correlated with the tumor cells' ability to breach the duct wall and invade into the surrounding stroma. Therapeutic effects of inhibitors of matrix metalloproteinases (MMPi) designed to inhibit enzymes involved in the degradation of basement membranes and extracellular matrices in the vicinity of growing tumors have been investigated by Ribba et al. ¹⁶ The authors proposed a multiscale model coupling the continuous representation of both the healthy tissue and the tumor mass with a discrete cell cycle including a restriction point at which the tumor cells can change from a proliferative to quiescent state depending on their local overpopulation. This model offered an explanation for why MMPi treatment in clinical trials was not as promising as the results obtained in animal studies. Computer simulations showed that as cancer cells become insensitive to antigrowth signals, the efficacy of MMPi treatment decreases, and is not able to prevent the basement membrane breach by the invading tumor cells (Figure 4c).

OPEN IN VIEWER

Figure 4

Snapshots from simulations of various mathematical models of tumor growth. (a) A micropapillary pattern of ductal carcinoma in situ; IBCell model; reprinted from reference, ¹² with permission from Elsevier. (b) Growing solid DCIS with different cell subpopulations including epithelial, progenitors, apoptotic and malignant cells; cellular automata model; reprinted from reference, ¹⁴ with permission from Elsevier. (c) Invading epithelial cancer cells with basal membrane degradation; continuous model; reprinted from reference, ¹⁶ with permission from Elsevier. (d) Rearrangement of ECM fibers (bottom) by the migrating tumor cell density (top) with dominant cell orientation (middle); continuous model; reprinted from reference ³⁵ (Figure 6), with permission from Springer Science. (e) A capillary sprout morphology resulting from VEGF gradients and endothelial cell migration along the matrix fibers; cellular Potts model; reprinted from reference, ²² with permission from Elsevier. (f) Formation of finger-like cohorts of invasive tumor cells (bottom) in an non-homogeneous ECM (top); cellular automata model; reprinted from reference, ³⁷ with permission from Elsevier. (g) Clonal growth of tumor cells due to cell competition for space and nutrients, leading to the configuration of three different cell subpopulations: an inner core of necrotic cells, an outer rim of proliferating cells and a middle zone of quiescent cells between them (left) and gradient of nutrients concentration (right); IBCell model; reprinted from reference, ⁴⁵ with permission from Birkhauser-Verlag. (h) Tumor-induced angiogenesis and vascular tumor growth. Shown from left to right, top: tumor cell subpopulations, oxygen profile, variable intravascular radius and ECM density profile; bottom: mechanical pressure profile, intravascular pressure profile, hematocrit profile and TAF concentration; hybrid continuous model; reprinted from reference ²³ (Figures 7 and 8), with permission from Springer Science. (i) Creation and amplification of autologous morphogen gradients by matrix-binding properties of morphogen and subtle physiological flows (right) and with no interstitial flow present (left); continuous model; reprinted from reference, ⁵¹ with permission from Elsevier

Stromal cellular structure

The tumor microenvironment is comprised of a variety of cell types, such as endothelial cells, fibroblasts and immune cells, which constantly interact with malignant cells. These stromal cells provide a rich source of growth factors and cytokines, and play a crucial role in modifying the existing connective tissue and depositing new matrix. As the tumor cell population expands, at a certain point the existing vasculature can no longer meet the nutrient and metabolic requirements of the cells. This can lead to the activation of endothelial cells lining the nearby blood vessel, inducing endothelial cell proliferation and chemotactic migration towards the tumor. This results in the growth of new vessels from pre-existing vascular beds (the process of angiogenesis). Fibroblasts located in the surrounding connective tissue are activated to form myofibroblasts. These myofibroblasts are capable of synthesis, deposition and remodeling of the underlying tissue matrix. In addition, these cells produce many soluble paracrine growth factors that regulate tumor cell proliferation, morphology, survival and death.

Immune cells are readily attracted to the tumor site from the peripheral blood or lymphatic vasculature. They are responsible for detecting and destroying all entities that are foreign to the body, not necessarily tumor-related. However, tumor cells present antigens distinct from those of normal cells and the surrounding stroma that may alert the immune system and stimulate their surveillance action. The cells forming the adaptive immune response, such as cytotoxic T-lymphocytes (CD8+), are drawn to a location due to the presence of a specific antigen and induce apoptosis in cells presenting these ‘foreign’ antigens. On the other hand, the cells forming the innate immune system, such as natural killer (NK) cells or macrophages, infiltrate the tumor tissue and target tumor cells without the need to recognize a specific antigen before swinging into action. However, the role of tumor-associated macrophages in regulating tumor growth is complex, as they may inhibit tumor growth by producing antiangiogenic proteins, and kill tumor cells either directly by producing cytotoxic molecules or indirectly by inducing antitumor effects in NK cells or cytotoxic T-lymphocytes. On the other hand, macrophages and other immune cell types are able to produce an array of pro-tumor factors that can stimulate angiogenesis, promote tumor growth and facilitate metastasis. As with any organ environment, there are other supporting cell types that can play individualized roles and provide a rich source of growth factors and cytokines, such as neurons and adipocytes. Interactions between the growing tumors and various stromal cells have been of great interest to mathematical modelers; however, a larger body of research exists focused specifically on processes such as angiogenesis and immune surveillance. This research can define and establish biological roles for certain stromal cell types, such as endothelial and immune cells, whereas other types of stromal cells are often disregarded when mathematically modeling tumorigenesis due to the ongoing but not yet completed experimental investigations of the role these cell types play in tumorigenesis.

The process of tumor-induced angiogenesis has been investigated by many different mathematical models, including continuous and individual-cell-based ones (see recent reviews in references ^17–19 and references therein), shedding light on the mechanisms of new capillary sprouting, as well as interactions between endothelial and tumor cells. In general, the following scenario has been considered. Upon an increase in tumor size, the areas of nutrients deprivation (such as hypoxic regions) are formed, which induce tumor cells to secrete growth factors, e.g. vascular endothelial growth factor (VEGF) and basic fibroblast growth factor. These factors influence activation of endothelial cells lining the blood vessels and initiate the process of angiogenesis, thereby bringing essential nutrients to the tumor and facilitating its further vascular growth. The more recent computational models extended this scenario by including additional factors, such as the structure of the extracellular matrix (ECM), auxiliary molecular species and different mechanisms of capillary formation and growth. Milde et al. ²⁰ investigated the direct effects of the ECM and the soluble as well as matrix-bound growth factors on capillary growth. They constructed a 3D hybrid-particle model, in which (1) endothelial cells are represented as particles; (2) molecular species, such as fibronectin, MMPs and VEGF, are represented by their concentrations; and (3) the orientation of randomly placed bundles of collagen fibers that modulate the direction and velocity of endothelial cells are modeled as a vector field. A large-scale parameter study using computer simulations revealed that in the presence of matrix-bound VEGF, an increase in the number of branches is observed as a result of an increase in both the amount of distributed VEGF and cell–matrix adhesion. Merks et al. ²¹ investigated biophysical mechanisms behind the formation of new blood vessels and showed that contact-inhibited chemotaxis to an autocrine, secreted chemoattractant can explain aspects of both de novo and sprouting blood-vessel growth. The authors used the cellular Potts model (see Figure 2) to represent endothelial cells together with the continuous description of a chemoattractant secreted by endothelial cells that stimulated the movement and growth of other endothelial cells. Their results suggested that branching in aggregates of chemotacting endothelial cells could result from two separate effects of the same mechanism based entirely on experimentally observed cell-level properties, i.e. for low cell motilities the branching resembles a buckling instability in which the surface cells exert a force on the cluster's inner core pushing the cells to buckle outside. However, for larger cell motilities, the shallower chemoattractant gradients at protrusions make the endothelial cells there more likely to extend outward-directed pseudopods. Two other interesting models, by Bauer et al. ²² and Macklin et al. ²³ describing capillary sprouting in addition to other tumor-related phenomena, are described later in the text.

The interactions between tumor and immune cells were investigated by de Pillis and Radunskaya ²⁴ by using a continuous model and by Mallet and de Pillis ²⁵ by using a hybrid cellular automata model. In both cases the authors have studied the role of NK and cytotoxic T-lymphocytes cells in tumor surveillance. The continuous model focused on the long-duration impact of the immune system on tumor growth and indicated that to promote tumor regression, it may be necessary (although perhaps not sufficient) to focus on increasing CD8+ T-cell activity. In fact, the authors proposed that there may be a direct positive correlation between the patient-specific efficacy of the CD8+ T-cell response, as measured by cytotoxicity assays, and the likelihood of a patient favorably responding to immunotherapy treatments. To validate the model, the authors compared their results to previously published data from mouse and human studies, in particular, CD8+ T-tumor and NK-tumor lysis data from chromium release assays. In the individual-cell-based cellular automata model, the authors explored the impact of the host immune system on the morphology of the initial stages of solid tumor growth and the growth of undetectable metastases, and revealed that when the immune system is introduced to the model, tumor growth depends strongly on the combined effects of the innate and specific immune cells, the potency in killing tumor cells and the levels of tumor infiltration by the immune cells. The authors were able to identify parameter sets that correspond to strong immune systems of healthy individuals, capable of early tumor detection and weaker immune systems of individuals from whom tumors are able to grow to the point of metastasis. They have also determined the parameter regimes in which T-cell infiltration is effective in tumor destruction. Matzavinos et al. ²⁶ have investigated whether the tumor-infiltrating cytotoxic lymphocytes can control cancer dormancy, a dynamic equilibrium state between the tumor and immune cells, in a relatively small tumor without central necrosis. The authors applied a continuous model in which they considered a population of tumor cells, together with active and inactivated T-lymphocytes. The T-lymphocytes are assumed to migrate into the growing solid tumor and to interact with the tumor cells resulting in either tumor cell death or lymphocyte inactivation. This model identified critical parameters during which cancer cells are present in the tissue but are not clinically detectable for a long period of time; furthermore, the tumor cells can begin to grow progressively at a later time. Hence, this model could be used to estimate the time interval between the primary treatment of an immunogenic tumor and tumor recurrence.

New therapies that exploit the tendency of macrophages to infiltrate solid tumors have been investigated by Byrne and co-workers. ^27,28 They have developed continuous mathematical models of a growing avascular tumor spheroid. The tumor volume is composed of a mixture of tumor cells, macrophages and extracellular material, and tumor cell proliferation and death are regulated by nutrient diffusion. Macrophages are assumed to migrate through the tumor by a combination of random motion and chemotaxis. The chemoattractant is produced by the tumor cells under hypoxia. This modeling approach has been used to study the capability of specially engineered macrophages to induce tumor cell lysis, either directly by delivering cytotoxic factors under hypoxia or by releasing enzymes that activate an externally applied pro-drug that causes fatal damage to proliferating cells. The model showed that if the overall tumor burden (i.e. the number of tumor cells) increased when macrophages were present, the spheroid volume occupied by tumor cells declined. Moreover, the number of tumor cells in the hypoxic region becomes reduced as the tumor cells are displaced by the macrophages, thus making the tumor more responsive to standard chemotherapeutic agents that target proliferating cells. This may contribute to uncertainty about the use of macrophage infiltration as a prognostic indicator.

Stromal physical structure

The ECM provides biological context and structural support for cells via a class of matrix receptors termed integrins that bind the ECM with some specificity, initiating signaling cascades that regulate various cell functions such as growth, survival and migration. Additionally, the mechanical tension of the ECM, itself, is informationally transmitted through this integrin engagement and receptor association with intracellular cytoskeleton elements. Furthermore, if one considers the porosity of the ECM, a result of both the density and covalent interactions of the ECM that resolves into a mesh-like structure, the matrix can serve as a footpath for cellular migration as well as a surmountable barrier. In order for transformed epithelial cells to become invasive, they need to actively migrate into the surrounding stromal tissue consisting of a meshwork of various cross-linked fibers, such as filaments of collagen, elastin, fibronectin and laminin, that together constitute the ECM. The organization and density of ECM fibers surrounding the normal epithelium differs from tumors. ²⁹ Moreover, the alignment and cross-linking of stromal fibers around tumors can be altered by both tumor-associated stromal fibroblasts ³⁰ and the invading tumor cells ²⁹ in such a way to facilitate tumor cell migration along the ECM fibers. ^31,32 Moreover, the elevated tissue rigidity (that can be modified in vitro by increasing the collagen concentration, the extent of gel cross-linking or the gel stiffness by applying exogenous force ³³ ) influences molecular signals that promote the malignant behavior of tumor cells. ³⁴

Different ECM alignments and their influence on two distinct cell motility modes (mesenchymal- and ameboid-like) have been investigated by Painter ³⁵ using a continuous model in which the orientation of stromal fibers is modeled as a direction field. The author showed that the processes of ECM contact guidance and remodeling can spatially organize cell populations to form chains of cells aligned along the matrix fibers (Figure 4d), and such patterns arise without additional environmental cues such as chemoattractants. By applying this model to tumor invasion, it has been shown that the interactions between invading cells and the ECM structure surrounding the tumor can have profound impact on the pattern and rate of cell infiltration. The effects of inhomogeneity in the ECM (modeled as anisotropic alignment of stromal fibers) on the formation of vascular sprouts during tumor-induced angiogenesis have been investigated by Bauer et al. ²² using a cellular Potts model. In this model, the configuration of stromal fibers is designed by distribution of thick fiber bundles (composed of nearby grid sites) at randomly chosen discrete orientations until the desired fibers density is achieved. The dynamics of growing endothelial cells is influenced by both mechanical (surface adhesion and matrix fibers' resistance to compression due to cell migration) and chemical (VEGF binding and uptake, tip cell activation by VEGF and matrix degradation) interactions with the ECM. Cells move in such a way to reduce the total energy of the system, i.e. endothelial cells will move to promote stronger over weaker adhesive bonds, and towards regions of higher chemical concentration (energy minimalization is the main principle of the cellular Potts model). This model showed that branching and anastomosis of a growing angiogenic vasculature emerge naturally (without any prescribed rule incorporated in the model) due to anisotropies in the stroma, such as variable matrix fiber alignment and the presence of other cells (Figure 4e). The influence of ECM cross-linking and fiber alignment on the formation and function of invadopodia (the slender, actin-rich protrusions extending into the ECM that are responsible for matrix degradation) was investigated by Enderling and collaborators in reference. ³⁶ The authors used a cellular automata approach to model both randomly distributed cross-linked ECM fibers and the cell invadopodia that can potentially grow in four orthogonal directions or remain stationary. The model combined with adequate experimentation revealed that cross-linking of the ECM inhibits the invadopodia-associated degradation and penetration of ECM. Moreover, it pointed out that the organization of collagen, as may occur in the stroma, forms gaps large enough so as to have little impact on invadopodia penetration/degradation. By contrast, dense ECM, such as occurs in basement membranes, is an effective obstacle to invadopodia penetration and degradation, particularly when cross-linked.

Other modelers have considered inhomogeneous properties of the stroma, such as its variable density, stiffness or exerted tension, in a more general way without directly including fine details of stromal cells and fiber alignment. Anderson et al. ³⁷ have discussed how the heterogeneity in the ECM density may act as a selective environmental pressure driving tumor invasion. They employed a cellular automata model, in which tumor development depends on the concentration of oxygen and production of matrix-degrading enzymes (MDEs), to compare several simulations of tumor growth that differs only in the structure of the surrounding stroma. Final morphologies obtained from these simulations revealed that the homogeneous environments promote the growth of tumor clusters with smooth round boundaries, whereas heterogeneous environments (in which the ECM density forms smooth gradients or changes randomly) drive growth of irregular, finger-like invasive cohorts of tumors cells (Figure 4f). Chaplain et al. ³⁸ have studied the potential role of a pathological cellular response to tissue compression. They showed how loss of contact inhibition to growth and of compression responsiveness can generate, by itself, clonal advantage that gives rise to the formation of hyperplasia and tumor growth. The authors employed a continuous model in which the space is occupied by a mixture of both normal and abnormal (pretumor) cells, the ECM produced by each type of cells, and the extracellular fluid in which the MDE can diffuse. Moreover, all cellular and extracellular components can influence each other's behavior by exerting stress due to cell growth and movement. The authors also showed that a stress-dependent production of ECM and MDE, in which either an increased production of ECM is not balanced by an increased release of MDEs or an excessive degradation of ECM is due to increased production of MDEs, can both lead to abnormal tissue development observed in many tumors.

Stromal chemical structure

The stroma surrounding normal and tumor tissues has a very complex chemical structure containing various diffusive or matrix-bound components. Most nutrients required for cell survival, such as oxygen or glucose, are transported through a vascular system and subsequently diffuse forming decreasing gradients within the tissue site. Many other proteins critical for cell life processes, such as growth factors, chemokines or MDEs, can be secreted by certain stromal cells and subsequently activated by tumor cells. As a result, these proteins form inhomogeneous patterns in the stroma. Moreover, all proteins that are not matrix bound can be transported either actively (via convection of the interstitial fluid) or passively (via diffusion) through the tissue. To maintain tissue homeostasis the balance between delivery and removal of each chemical factor, and its uptake and/or release by the cells is needed. However, in many cases the way of achieving such a balance may be quite complex since some (if not all) factors can simultaneously influence cell behavior and be altered by the very same cells. The complex non-linear relations between oxygen, glucose and extracellular pH levels in the cultured 3D multicellular spheroids have been investigated by Sutherland and co-workers, ^39,40 whereas the role of matrix proteinases in cancer development, invasion and suppression has been reviewed in references. ^41,42 The evolution of glycolytic and acid-resistant cell phenotypes in growing tumors as a critical step toward the emergence of invasive cancers has been discussed by Gatenby and Gilles. ^43,44

Mathematical models were used extensively to investigate relations between gradients of various metabolites and the emerging morphologies of developing tumors. The establishment of a three-layered tumor structure (consisting of proliferating, quiescent and necrotic cells) that is a consequence of nutrient depletion has been reproduced by virtually every kind of modeling approach, and has become a test problem for every newly developed mathematical model of solid tumor growth, whether continuous or single-cell-based (see reviews ^3–9 ). The typical situation presented in Figure 4g ⁴⁵ shows the gradient of oxygen with a highly depleted central area and a steady high concentration in the outside media (right panel) that coincide with locations of three cell sub-populations: an inner core of necrotic cells (cyan, left panel) that correlates with a low level of oxygen, a ring of viable quiescent cells (gray) and an outer rim of proliferating cells (red). Usually, the following terms are included in a mathematical reaction–diffusion model of metabolite kinetics: (i) a supply or production term, such as supply of oxygen or glucose from a capillary, or production of MDEs by the tumor cells; (ii) diffusion through the ECM or culture media as a means of an active transport that results in the creation of a metabolite gradient from its source towards the distant parts of the tissue; (iii) a decay term, naturally most metabolites have a tendency to dissipate over time; and (iv) a consumption or uptake term, such as consumption of oxygen or glucose by the viable cells, or depletion of enzymes used for ECM degradation. Once the cell metabolism is updated (uptake of oxygen and nutrients or production of waste products), cell phenotype may be changed depending on predefined conditions that relate cell metabolism with cell life processes, such as cell death (e.g. if the level of oxygen is low), proliferation (e.g. if the cell is not in contact inhibition) or migration (e.g. towards higher concentrations of chemoattractants).

Various mathematical approaches have been used to define cell response to multiple, and often contradictory, microenvironmental cues sensed by the cell. Silva et al. ⁴⁶ have investigated the hypothesis that increased systemic concentrations of pH buffers reduce intratumoral and peritumoral acidosis and, as a result, inhibit malignant growth. The authors used a hybrid 3D cellular automata to model tumor cells and reaction–diffusion equations to model several chemical species: glucose, oxygen, carbon dioxide, hydrogen ions, ATP, bicarbonate anions and a hypothetical buffer. All these chemical factors were coupled with each other (e.g. production of ATP and CO₂ depends on the levels of glucose and oxygen; production of lactate, hydrogen ions and ATP depends on the level of glucose) and together influenced cell fate (e.g. the cell death occurred when either the level of pH or the amount of ATP production fell below the prescribed thresholds; cell proliferation is proportional to ATP production rate provided there is an empty space in the cell vicinity to place a daughter cell). The simulations showed that increased concentration of serum bicarbonate can decrease intratumor and peritumor acidosis without altering blood pH, and that reduction in these pH gradients will also reduce tumor growth and invasion.

Anderson et al. ⁴⁷ have used three different mathematical models to investigate mechanisms underlying the emergence of finger-like cohorts of invasive tumor cells in relation to nutrient availability and cell metabolism. By reproducing similar cellular behavior over a range of parameters in three different mathematical models based on distinct mathematical formulations, they showed that the obtained results were model-independent and were a consequence of the underlying biological assumptions. In the Anderson cellular automata model each tumor cell is characterized by a predefined set of phenotypic traits, such as sensitivity to cell–cell adhesion, cell doubling time, rate of oxygen consumption, rate of production of MDEs and cell sensitivity to haptotactic cues that drive cell movement toward the area of higher matrix concentration. In general, these phenotypes are passed from a mother cell to its daughters; however, with a small probability, a daughter cell may acquire a different phenotype chosen randomly from a 100 predefined phenotypes. Tumor cell fate depends then on this phenotype and the concentration of oxygen and ECM density sensed by the cell from its immediate vicinity. The second model, developed by Gerlee, is a modification of the Anderson model in which each cell is equipped with a decision mechanism (a response network modeled as a feed-forward neural network) that determines cell behavior, such as cell proliferation, quiescence or apoptosis, based on cell ‘genotype’ and the microenvironmental cues, such as the number of neighbors and the concentration of oxygen. Cell ‘genotype’ consists of network parameters that are initially chosen by training the network to capture the essential behavior of real cancer cells, i.e. in the initial genotype the competition for space determines if the cell will proliferate or become quiescent, while the oxygen concentration influences the apoptotic response. When the cell divides the network parameters are copied to the daughter cell with a slight probability that some of them will be mutated, which implies that daughter cells may react differently to a similar microenvironment. Therefore, in the Gerlee model changes in cell phenotype are not random like in the Anderson model, but can be slightly modified when compared with the mother-cell phenotype. In the third model, developed by Rejniak, each cell is represented by a two-dimensional, deformable, fluid-filled membrane equipped with a set of cell membrane receptors. These receptors are capable of forming adhesive connections with other cells and with the ECM, sensing concentrations of nutrients in their immediate vicinity, and determining whether or not there is a space for cell expansion, i.e. if the membrane receptors are not engaged in either cell–cell or cell–ECM adhesion. Based on the configuration of different membrane receptors, the cell can undergo a receptor-driven change in its phenotype. If a cumulative concentration of nutrients sensed through all cell membrane receptors is below the prescribed threshold, the cell will enter into a necrotic state. However, when the concentration of nutrients is high enough to keep the cell viable, it will inspect the space available for its growth. If a certain number of cell receptors is engaged in cell–cell adhesion, thus resulting in cell contact inhibition, the cell remains quiescent; otherwise, it will acquire the fluid necessary for its growth and will expand towards the space not occupied by other cells. This receptor-driven cell dynamics allows one to relate cell life processes and changes in cell phenotype to the total or average concentrations of external factors sensed by the whole cell, or even to concentrations sensed by individual cell membrane receptors. All three models examined the impact of nutrient availability and cell metabolism (the intrinsic rate of nutrient consumption and cell resistance to starvation) on tumor morphology and aggressiveness, and showed that tumor fingering occurred when cell metabolism was high or when nutrient availability was low. Therefore, these simulations suggest that invasion is an outcome of the competition between the cells induced by a given microenvironment. This competition occurs only if there exists variation either within the cell population or/and within the microenvironment. Thus, invasion is an emergent property of malignant tumors at the tissue scale not at the individual cell scale.

Macklin et al. ²³ have developed a multiscale model of vascular solid tumor growth that couples a non-linear continuum model of solid tumor invasion (solved by using the ghost-cell/level-set methods), including cell–cell adhesion, cell–ECM adhesion, ECM degradation, tumor cell migration, proliferation and necrosis, with a hybrid model of tumor-induced angiogenesis that accounts for vascular remodeling due to wall sheer stress and mechanical stresses generated by the growing tumor and for blood flow through the vascular network. The invasion and angiogenesis models are coupled through the stromal chemical factors, such as the tumor angiogenic factors (TAF) that are released by tumor cells and influence the formation of a neo-vascular network, and the nutrients that extravasate from the vascular network, diffuse through the ECM and trigger further growth of the tumor. In addition, the vascular network and tumor expansion are also coupled via the ECM since both the tumor and the endothelial cells upregulate matrix-degrading proteolytic enzymes that cause localized degradation of the ECM, which in turn affects cell haptotactic migration. Through a series of computer simulations, the authors have traced shape instability of the growing tumor (with a layered structure of proliferating, quiescent and necrotic cells), the profiles of ECM density, mechanical pressure, oxygen and TAF concentrations in the ECM, as well as intravascular radii, pressure and hematocrit profiles (Figure 4h). Their results demonstrated that hydrostatic stress generated by tumor cell proliferation shuts down large portions of the vascular network affecting blood flow, subsequent network remodeling, delivery of nutrients to the tumor and tumor progression. In addition, ECM degradation by tumor cells had a dramatic effect on both the development of the vascular network and the growth response of the tumor. In particular, when the ECM degradation was significant, the newly formed vessels tended to encapsulate rather than penetrate the tumor and were thus less effective in delivering nutrients.

Stromal fluid structure

The interstitial space between cells is filled with interstitial fluid that consists of a water solvent containing amino acids, sugars, fatty acids, hormones and salts, as well as waste products from the cells. The interstitial water space in tumors (calculated as a fraction of all water space including cellular, interstitial and vascular) comprises a significant part of total tissue water (30–60% in carcinomas, 50–70% in skin, 20% in muscle ^48–50 ). It plays an important role in the exchange of oxygen, nutrients and waste products in normal tissues, as well as in the delivery of anticancer drugs in tumors. The transcapillary flow of water and other molecules is carried through the interstitial space due to the differences between pressures in capillaries and in the interstitium (interstitial fluid pressure gradient, IFP) and the drainage into the lymph vessels that prevents a build-up of tissue fluid. It has been noticed that the IFP is often increased in solid tumors, and may constitute a barrier for efficient drug delivery, since many anticancer drugs are transported through the interstitial space by convection rather than diffusion ⁵¹ ; however, the mechanisms leading to increased IFP levels in tumors are not fully understood. ^52,53

Swartz and co-workers have investigated theoretically ⁵⁴ and demonstrated experimentally ⁵⁵ the ability of slow interstitial flows to drive a significant asymmetry in pericellular protein concentrations to create transcellular chemical gradients in the direction of the flow. They have shown that when proteins are secreted in a matrix-binding form and transported by convective flow, autologous chemotactic gradients increasing in the direction of the flow can be cell-created (Figure 4i). Moreover, the combined influence of both factors (flow and matrix binding) has been amplified when compared with the effects of each factor alone. This mechanosensing mechanism allows the cell to gather information about the dynamic status of its microenvironment, and to create the autologous chemotactic gradient directing its migration into tissue-draining lymphatics. The effects of antiangiogenic therapy on the levels of the IFP and velocity have been mathematically modeled by Jain et al. ⁵⁶ This investigation revealed that upon the antiangiogenic treatment, interstitial convection within the tumor may increase, whereas fluid convection out of the tumor margin decreases. This can have a few of the following implications: increased drug convection within the tumor will improve treatment efficacy; decreased convection of growth factors from the tumor margin into the peritumor tissue may reduce angiogenesis and limit peritumor hyperplasia; and decreased convective flow out of the tumor will likely decrease dissemination of cancer cells to lymph nodes. In all cases the probability of lymphatic metastasis would be reduced. The impact of fibrillar assembly and tortuosity of the interstitium on the effective diffusion and convection of tracer molecules has been investigated by Ramanujan et al. ⁵⁷ Stromal cells residing in the interstitial space create obstacles for moving molecules that can be measured by the interstitial space tortouisity defined as a ratio between the effective path of a molecule's transport and a linear path length. Similarly, the assembly and cross-linking of ECM fibers can hinder molecular transport. The mathematical analysis revealed that collagen in physiological concentration presents a major barrier to molecular diffusion, especially for larger particles. Comparison of diffusion and convection data suggested that collagen may obstruct diffusion more than convection transport in tumors.

Discussion

We have discussed in this Minireview several structural and chemical components of the tumor microenvironment and highlighted the most prominent aspects of their roles in tumor progression. Among these are the various types of stromal cells, the interstitial fluid interpenetrating bundles of ECM fibers, the gradients of diffusive and matrix-bound metabolites and other chemical agents, and finally the structure of the host tissue in which these tumors are growing. We have also presented a survey of different mathematical and computational approaches that have been undertaken to model various aspects of tumor–microenvironment interactions. It is worth noting that every mathematical model is by its nature a simplification of the biological system it is assumed to represent, so we do not expect that one model will incorporate all elements of the stroma surrounding the tumor. Rather, the models presented here have focused on these parts of the tumor microenvironment that are crucial for investigating a given problem in tumor progression and/or treatment, in a manner similar to how biological experiments focus on the selected aspects of tumor growth and do not address in a single experiment all possible combinations of involved factors. We have presented different classes of mathematical models – continuous and individual-cell-based: those in which cell properties are averaged over large populations of cells, and those that model individual cells with individually regulated cell processes. The later models, however, can be subdivided into two classes depending reciprocally on the number of cells the model can handle and the included details of each individual cell structure. Thus, we presented examples of both classes of individual-cell-based models: those dealing with large cell populations with simplified cell geometry and those modeling small colonies of fully deformable cells with detailed representation of cell life processes, including cell–cell and cell–microenvironment interactions. We have summarized in Figure 3 the features of each class of model discussed in this Minireview.

In conclusion, the presented examples of mathematical models simulating tumor growth and treatment indicate that:

It is relatively easy to introduce the structure of the host tissue into a mathematical model (we discussed the cases of continuous, cellular automata, particle-based, cellular Potts and elastic cell-based models), easier than to design a laboratory experiment incorporating the same phenomena;

When working on the edge between experiments and computations one can experience certain difficulties that need to be overcome in order to design a successful integrated model. Some of these model weaknesses have to be addressed by the experimentalists, and others by the bio-mathematicians. All are challenging in their respective scientific fields:

Model parameterization: ideally, one would like to be able to parameterize the computational model with experimental measurements taken directly from the system that is being modeled. However, often certain parameters required by the model are technically difficult to measure, so the current practice is to combine the parameters collected either directly from the experiments or from the literature on similar experimental systems. In order to investigate how the departure from the chosen parameters values will influence the final computational results, one can perform a parameter sensitivity analysis (either analytically or by using simulations);

Other investigations of tumor–stroma interactions that can highly benefit from using the integrated approach and high-throughput simulations studies may address: the influence of stromal structure on breaking the natural tissue homeostasis that leads to tumor initiation; the relation between tumor physical and chemical microenvironment and antitumor drug penetration; the impact of tumor intrinsic heterogeneity on alterations in stromal composition that subsequently influences tumor development and treatment; and the effect of patient-specific tissue structure on antitumor treatment efficacy.

Mathematical models can be applied in addressing such complex features since by their nature they are suited to handle multiple variables with numerous parameters. It is relatively easy and inexpensive to simulate tumor growth and treatment in silico and to compare differences in simulation outcomes when such parameters are changed simultaneously and over a wide range of values, which is not often possible in the laboratory. However, such models can only be successful when they are built in close collaboration between experimental and computational laboratories, since truly integrated models need constant interactions: (1) first, to design a computational model based on sound biological assumptions; (2) to parameterize the model with experimental measurements; (3) to test experimentally model hypotheses; (4) to modify the computational model depending on experimental outcome; and (5) to repeat this integrated loop as necessary. This approach may provide a new insight on tumor progression leading to testable predictions and hypothesis. The integration of experimental and computational data that is supported by constantly improving techniques for experimental data acquisition, and constantly increasing computational power is in our view the way to utilize interdisciplinary research. This approach may advance our understanding of cancer and pave the way for novel antitumor interventions.

Author contributions: The authors have contributed equally to this review providing expertise in cancer biology (LJM) and mathematical modeling of cancer (KAR).