Biology Background

The development of a primary solid tumor starts from a single cell that proliferates in an inappropriate manner, dividing repeatedly to form a cluster of tumor cells. Nutrient and waste diffusion lengths limit the size of such avascular tumor spheroids to about 1 mm [1]. The central region of the spheroid becomes necrotic, with a surrounding layer of cells whose hypoxia triggers a cascade of hypoxia-inducible factor-1 (HIF-1) and vascular-endothelial-growth-factor (VEGF)-mediated signaling events that initiate tumor vascularization by promoting growth and extension (angiogenesis) of nearby blood vessels [2]. This general sequence occurs in many types of tumors. Vascularized tumors are able to grow to a much larger size than spheroids and are more likely to spread and metastasize using blood vessels as pathways for invasion [3].

Both fetal and adult angiogenesis is primarily a response to hypoxia [2], [4]–[7]. In adults, angiogenesis plays key roles during tissue repair and remodeling, e.g. during wound healing and expansion of tissues during the female reproductive cycle.

The level of HIF-1α DNA-binding activity in the nucleus varies exponentially with oxygen tension over the physiological range in mammalian cells [8]. Cells exposed to hypoxic conditions accumulate activated HIF-1α in their nuclei in less than 2 minutes [9]. HIF-1α changes the expression levels of numerous hypoxia-dependent genes including those responsible for oxygen transport, vascular regulation, angiogenesis, glucose uptake, and glycolysis (reviewed in [10]). HIF activation also upregulates key angiogenesis regulatory signaling molecules including VEGF-A and its receptors VEGFR-1 and VEGFR-2 [2]. VEGF-A has diffusive, ECM-bound, and semi-ECM-bound isoforms which differ in weight and heparin-binding affinity [11]. Tumors secreting different VEGF-A isoforms induce the formation of morphologically different neo-vascular beds [12].

Endothelial cells form two distinct types which respond differently to VEGF-A. Non-dividing tip cells form filopodia and migrate towards sources of VEGF-A, while non-migrating stalk cells proliferate but do not form filopodia. [13]–[15]. Two cell types are functionally necessary because if every endothelial cell were a tip cell, the vessels would disintegrate, while uniform division of stalk cells would fail to form vessels in the correct pattern [13], [14].

Cell-adhesion plays a crucial role in the formation and stabilization of nascent blood vessels (see [16], [17] for reviews). Formation of cell-cell adhesion junctions via cadherins like VE-cadherin inhibits the chemotaxis response of endothelial cells to VEGF-A at endothelial cell-endothelial cell boundaries (contact-inhibited chemotaxis) and increases the stability of those boundaries [18]. The growth rate of cultured endothelial cells decreases as the area of VE-cadherin junctions increases (contact inhibition of growth)[19].

Normal new blood vessels and tumor-induced blood vessels differ greatly in morphology and function. Normal new vessels recruit pericytes and vascular smooth-muscle cells to sites adjacent to the endothelial cells to stabilize the vessel [20]. Tumor-induced blood vessels often lack a hierarchical arrangement and have irregular diameters, high tortuosity, random branching, and defective endothelial barrier function [20]–[22].

Computational Background

Because of the range of scales involved in cancer biology, cancer simulations employ a wide range of techniques depending on their biological focus. Standard partial-differential-equation (PDE) continuum models include scales down to the grid size used to solve the equations. Continuum multi-grid and adaptive-mesh techniques can cover scales from µm to 10 cm. Hybrid models, which use cellular automata, agent-based or multi-cell techniques to represent individual cells and continuum PDEs to represent diffusion of molecules, can capture more detail than continuum models spanning scales from microns to millimeters (for comprehensive reviews of mathematical models of tumor-growth and angiogenesis see [23]–[27] and references therein).

Zheng et al. [28] have used an adaptive finite-element/level-set method to model solid tumor growth in combination with Anderson and Chaplain's hybrid model of angiogenesis [24]. Zheng's model treats tumor cells as a viscous fluid flowing through a porous medium obeying the Darcy-Stokes law. Zheng et al. have shown that both diffusional instability (competition of growth and surface tension) and co-option of the new anastomosed capillaries may be key glioma invasion mechanisms. Frieboes et al. [29] have used Zheng's level-set method in combination with Plank and Sleeman's hybrid continuum-discrete [30], lattice-free model of tumor angiogenesis to model the physiology and evolution of glioma neovasculature in 3D. Frieboes et al's model allowed them to correlate measurable tumor microenvironment parameters to cell phenotypes and potentially to tumor-scale growth and invasion. Cristini et al. [31] have also developed a continuum model of solid avascular tumors using a mixture model obeying Cahn-Hilliard-type equations. Cristini et al. found that taxis of tumor cells up gradients of nutrient produces fingering instabilities, especially when tumor cell proliferation is slow.

The effects of blood flow on vascular remodeling and tumor growth have been extensively studied by Owen et al. [32], Alarcon et al. [33], Bartha and Rieger [34], Welter et al. [35], McDougall et al. [36], Stephanou et al. [37], [38], Pries et al. [39]–[42] and Macklin et al. [26]. Other vasculature and transport-related effects remain to be studied [43], for example how tumor cells interact with endothelial cells and enter and exit the blood stream (intra/extravasation) and spread to other organs (metastasis) via blood vessels. The simple model we present in this initial paper neglects the crucial effects of detailed transport due to blood flow and the effects of flow-induced vascular remodeling. We discuss some of the resulting artifacts and missing behaviors below.

Because the Glazier-Graner-Hogeweg model (GGH, also known as the Cellular Potts Model, CPM) [44]–[46] handles cell-cell and cell-ECM adhesion and cell motility more naturally than many other modeling methods, GGH simulations may provide additional insight into tumor growth and the complex roles of angiogenesis. Poplawski et al. [47] developed a GGH simulation of avascular tumor based on Anderson's model [48] to study the effects of adhesion and nutrient transport on the morphology of avascular tumors. Poplawski et al. found that nutrient-deprived tumors are generally more invasive and that high tumor-ECM surface tension changes seaweed-like tumor-morphologies into dendritic morphologies. Rubenstein and Kaufman [49] have developed a GGH model of avascular tumors with explicit representation of two types of ECM fibers to study the effect of ECM on growth of glioma spheroids in vitro. Rubenstein and Kaufman showed that invasion is maximized at intermediate collagen concentrations, as occurs experimentally.

In our simulations, the vascular structure produces inhomogeneities in oxygen tension on scales larger than those appearing in continuum simulations depending on inhomogeneities in Extracellular Matrix (ECM) and smaller than those due to tissue structure. These inhomogeneities may affect tumor growth rates and morphology and the somatic evolution of metastatic potential within tumors. Thus, simulating vascular structure at the cell level is crucial to developing biomedically-meaningful tumor simulations.

In this paper we simulate 3D solid tumor growth and angiogenesis using the multi-cell GGH model. Our simulation omits many biological details, but provides a useful starting point for the construction of more realistic models. We focus on tumors where the vasculature remains peripheral to the growing tumor. Our major simplifications include: 1) We neglect the distinction between veins and arteries, anastomosis, and the possible presence of pericytes and smooth muscle cells. 2) Since oxygen-depleted areas of a tissue coincide with nutrient-depleted and energy-depleted areas [50], we assume that oxygen serves as the single limiting substrate field. 3) Cells become hypoxic or necrotic by simple thresholding depending on the local concentration of oxygen. 4) We neglect the basal metabolic consumption of oxygen by tumor cells. 5) We assume that the oxygen concentration in the capillaries is constant along the blood vessels, neglecting vessel diameter, blood flow rate and vessel collapse due to external pressure. 6) We assume that oxygen diffuses uniformly in the host tissue and tumor. 7) We caricature the results of the hypoxic signaling pathways as constant-rate secretion by hypoxic cells of a single long-diffusion-length isoform of VEGF-A like VEGF-A. 8) Since we do not model blood flow explicitly, we neglect its biomedically important effects on vascular remodeling and the maturation of nascent blood vessels. 9) Rather than model tip-cell selection explicitly, we distribute a certain number of inactive neovascular cells (terms in bold, e.g. neovascular, refer to specific simulation classes), which behave identically to vascular cells until the concentration of VEGF-A exceeds an activation threshold. 10) We adopt the vascular-patterning hypothesis that activated vascular endothelial cells secrete and chemotax towards a short-diffusion-length chemoattractant (for a discussion of possible chemoattractant candidates see [51], [52]). 11) We employ a simplified model of tight junctions between endothelial cells in the preexisting vasculature. 12) We do not represent the ECM and the cells in the surrounding normal tissue explicitly and neglect mechanisms related to ECM and normal tissue remodeling, like secretion of ECM proteins by tumor cells, secretion of ECM-degrading enzymes (matrix metalloproteinases), and lactic acid accumulation.

The Glazier-Graner-Hogeweg Model

Our simulation uses the Glazier-Graner-Hogeweg model (GGH, also known as the Cellular Potts Model, CPM), a multi-cell, lattice-based, stochastic model which describes biological cells and their interactions in terms of Effective Energies and constraints. GGH applications include modeling avascular tumor growth [47], [49], biofilms [55], chick limb growth [56], somitogenesis [57], blood flow and thrombus development [58] and angiogenesis [51], [59]–[61]. Each cell consists of a domain in a cell lattice of lattice sites, or voxels, at locations , sharing the same cell index, . Each cell has an associated cell type, . The Effective Energy governs changes in the cell configuration. The basic Effective Energy describes constant-volume cells interacting via differential adhesion:(1)where is the total number of lattice sites in cell , which is constrained to be close to the target volume (i.e., deviation of from increases the Effective Energy), and is the inverse compressibility of cells of type . is the contact energy per unit area between cells and . is the internal pressure in cell . To simulate the cytoskeletally-driven reorganization of cells and tissues, we model cell protrusions and retractions using a modified Metropolis dynamics. For each attempt of a cell to displace a neighbor, we select at random a cell boundary (source voxel) and a neighboring target voxel and calculate the change in the effective energy , if the source cell displaced the target cell at that voxel. If is negative, i.e., the change is energetically favorable, we make it. If is positive, we accept the change with probability , where describes the amplitude of cytoskeletal fluctuations. On a lattice with sites, displacement attempts represent our basic unit of time, one Monte Carlo Step (MCS). We describe secreted morphogens and oxygen macroscopically by concentration fields discretized at the resolution of the cell lattice. Fields evolve according to simple time-dependent partial-differential equations (PDEs). In addition to contact and volume energies, we include a chemotaxis term in the Effective Energy to model angiogenesis (see below).

We define a special, generalized cell representing the extracellular medium (ECM). The total volume and surface area of ECM are not constrained. ECM voxels can be both source voxels, e.g. during retraction of the trailing-edge of a cell, and target voxels, e.g. during formation of lamellipodia. Since ECM does not have cytoskeletal fluctuations, we use the amplitude of cytoskeletal fluctuations of the target or source cell to determine the acceptance of the displacement.

Our model contains three tumor-cell types: normal, hypoxic and necrotic (we use the designation tumor cell to refer to both normal and hypoxic cells). Normal tumor cells become hypoxic when the oxygen partial pressure () is less than 5 mmHg [22], [50], [62] and become necrotic when drops below 1 mmHg. Normal and hypoxic cells take up oxygen (see below) and proliferate at a rate which depends on the oxygen partial pressure according to a Michaelis-Menten form. To model tumor cell proliferation, we increase the cell's target volume according to:(2)where is the at the center-of-mass of the cell. We discuss the maximum growth rate, , in the next section. When a cell's volume reaches its doubling volume, we split the cell along the x–y plane. The two daughter cells receive equal target volumes of half of the target volume of the mother cell at mitosis. Necrotic cells lose volume at a constant rate:(3)

We remove zero-volume cells. Hypoxic cells secrete the field , which models the pro-angiogenic factor VEGF-A. Hypoxic cells stop secreting VEGF-A and become normal if is above 5 mmHg and become necrotic if drops below 1 mmHg. Hypoxic cells secrete at a constant normalized rate per MCS at each voxel in their volume. diffuses with diffusion constant and decays at a rate , so evolves according to:(4)where at voxels belong to hypoxic cells and elsewhere. Since cell motility is large, cells diffuse and rearrange fast enough to prevent artifacts due to their fixed cleavage plane [47].

Our model also contains two basic endothelial cell (EC) types: vascular and neovascular. We further distinguish two types of neovascular cells, inactive neovascular and active neovascular. Vascular cells build the preexisting mature vasculature. To represent tight junctions between endothelial cells in mature capillaries which maintain the integrity of blood vessels, linear springs connect the centers-of-mass of vascular cells. The springs obey the elastic constraint, , where is the equilibrium length of the connection and is the distance between the two neighbors. To model vascular rupture we set between the two neighbors when between the two neighbors is greater than . Inactive neovascular cells behave exactly like vascular cells. However, above a threshold value of inactive neovascular cells switch irreversibly into active neovascular cells. Active neovascular cells can proliferate and chemotax up gradients of VEGF-A. To account for the contact-inhibited growth of neovascular cells, when the common surface area with other vascular, inactive neovascular and active neovascular cells is less than a threshold, we increase the target volume of the active neovascular cells according to the Michaelis-Menten relation:(5)where is the concentration of VEGF-A at the center-of-mass of the cell and is a scaling constant describing the proportionality of the activation concentration to the concentration at which the growth rate is half its maximum. Active neovascular cells divide along the x–y plane when their volume reaches a specified doubling volume. After mitosis, both daughter cells are active neovascular and inherit half of their mother's target volume. Active neovascular cells at the tip of a sprout share two features with real tip-cells: 1) Compared to stalk cells, they have more free boundary which can respond to gradients of VEGF-A, dragging the rest of the sprout up the gradients. 2) Dragging reduces the contact area between active neovascular cells, promoting neovascular growth. Unlike real tip-cells, active neovascular cells at the tip of a sprout proliferate. Thus sprouts grow faster in steeper gradients of VEGF-A as long as the concentration of VEGF-A is well below the saturation concentration. We add a saturated Savil-Hogeweg-type chemotaxis term with contact inhibition to the basic of the Effective Energy to represent the net effect of preferential formation of pseudopods in response to the gradient of the chemoattractant field near the active neovascular cell's boundary [63]:(6)(7)(8)where is the degree of chemotactic response of the cell and s is a positive number which scales the VEGF-A concentration field relative to the neovascular activation threshold. .

Growth and chemotaxis of active neovascular cells up gradients of VEGF-A produce a dispersed growing population of neovascular cells rather than a connected capillary network [51]. To self-organize vascular and neovascular cells into a capillary-like structure we extended to 3D Merks' 2D model of angiogenesis in which endothelial cells self-organize into capillary-like networks in response to autocrine chemotaxis to a very short-diffusing chemoattractant [51]. We denote this chemoattractant field by . Vascular and neovascular cells secrete it at a constant rate at all their voxels, the chemoattractant degrades at a constant rate in the ECM, and diffuses at a constant rate everywhere:(9)where inside vascular and neovascular cells and elsewhere. Both vascular and neovascular cells in our model chemotax up gradients of . We include a linear version of eq. 8 [51] to model contact-inhibited chemotaxis:(10)(11)(12)

Autocrine chemotaxis produces a capillary-like structure, while the elastic constraint between vascular cells makes the preexisting vasculature more stable and less flexible.

Vasculature transports oxygen to the host tissue and the tumor. We represent the oxygen partial pressure, , by a field . Since we initialize our simulations with a fully anastomosed preexisting vasculature and since not all neovascular cells form in closed loops, we assume that the oxygen partial pressure in the preexisting capillaries, , is higher than in the tumor-induced vasculature (see [64], [65], for more accurate blood-flow calculations). We assume that the oxygen field concentration neighboring a vessel changes proportionally to according to a solubility coefficient, . Available diffuses at a constant rate in the surrounding tissue and is consumed by tumor cells. Since oxygen consumption of human cells is almost constant until drops below 0.5–1 mmHg and the consumption rate is constant below this pressure, we use a piecewise-linear approximation of Michaelis-Menten kinetics to model oxygen consumption. Tumor cells take up oxygen at a rate proportional to with a maximum rate of . The cells near the vasculature consume oxygen at their maximum rate, but cells far from the vasculature have growth limited by the available oxygen and take up oxygen at a rate . To represent oxygen consumption by host cells in the normal tissue, which we do not model explicitly, we assume an oxygen consumption rate of saturating at a maximum rate of in the ECM. We include oxygen consumption by EC cells in the adjusted oxygen partial pressure. Thus oxygen evolves according to:(13)(14)where and inside vascular and neovascular cells and while and elsewhere. is the oxygen consumption rate for both normal and hypoxic cells. We have summarized the properties of the fields in Table 1 and cell types and their behaviors in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Diffusive molecules in the vascular tumor-growth simulation.

https://doi.org/10.1371/journal.pone.0007190.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Cell types in the simulations and their behaviors.

https://doi.org/10.1371/journal.pone.0007190.t002

Implementation Parameters and Initial Conditions

Our simulations use the open-source CompuCell3D simulation environment (http://www.compucell3d.org/) [66]. We ran our simulation on a lattice with periodic boundary conditions. One voxel is equivalent to 125 µm. The cell lattice represents a tissue with a volume of mm. Our average simulated tumor cell has a volume of about voxels or 3375 µm. We stored the cell-lattice configuration every 3 simulated hours and rendered each snapshot using the MATLAB volume-visualization functions (Our group is integrating post-rendering capabilities into CompuCell3D). Since rendering individual cells is computationally expensive, in Figures 1 and 4 and Supplemental Movies (Movie S1 and S2) we have only rendered boundaries between cells which differ in type. For demonstration purposes, we have rendered the boundaries of individual tumor cells in the simulation with angiogenesis on day 60 (Figure 2B) and in the simulation without angiogenesis on day 10 (Figure 2A).

Experimentally, tumor cells from lines such as U-87 human glioma [67]–[70] can move at a rate of about 0.35 µm/min. For typical parameter settings in our simulations tumor cells move at about 0.1 pixels/MCS. Equating the experimental and simulated mean cell speeds implies MCS h.

We start our simulation with a single normal tumor cell near the center of the cell lattice and a preexisting network of blood vessels (see Figure 1). We assume that 90 mmHg in the preexisting vasculature and 50 mmHg in the tumor-induced vasculature. We set the diffusion constant of oxygen in our simulations to 10 µm/s, about half its diffusion constant in water. Since hypoxic cells in our simulations become necrotic if drops below 1 mmHg, we set both and to 1. We set so the average of the tissue reaches an asymptotic value of 20 mmHg. We assume that the density of cells in the tumor is about 10 times the density of the cells distributed in the ECM (which we do not represent explicitly). We assume that both hypoxic and normal tumor cells in our model consume oxygen at a rate up to 3 times that of the cells in the surrounding tissue, thus  = . Higher oxygen consumption results in shorter oxygen penetration lengths and smaller solid tumors. We set 10 mmHg and choose so that the cell cycle of tumor cells at is 24 h. For these parameters our simulations produce solid tumors with maximum diameters of 200 µm.

Experimentally, the VEGF-A diffusion constant is about 10 µm/s in typical tissues and its decay rate is about 0.65 h [52]. We set the activation VEGF-A concentration  = 0.5, , and pick so that a neovascular cell that does not contact any vascular or neovacular cells has a minimum cell-cycle time of 24 h. Due to contact-inhibition of growth, neovascular cells incorporated in a tumor-induced vessel grow at a negligible rate.

We assume that integrins are down-regulated in tumor cells and that cell-cell adhesion via cadherins keeps the tumor solid, i.e. that the surface tension at the tumor-ECM interface is positive, (for definitions of surface tensions see [44], [45]). We set a positive surface tension between EC(vascular and neovascular) and tumor cells () which keeps the vasculature peripheral to the tumor. For the specific values of J and the other parameters, see the supplemental XML and python files (Code S1).