Our model encompasses four biological scales (Figure ): the molecular scale, the cellular scale, the micro-environmental scale and the tissue scale. The molecular scale consists of the EGFR signaling and cell cycle pathways (Figure ). The cellular scale describes the phenotypic switch of the tumor cells which is determined by the molecular scale, the cell-cell and the cell-microenvironment interactions (Figure ). The micro-environmental scale provides a connection between the micro-vasculature (angiogenesis) and the tumor cells. Tumor cells consume glucose and oxygen nutrients. Transforming growth factor alpha (TGFα) secreted by tumor cells triggers the EGFR signaling pathway, which in turn determines cell’s phenotypic switch. Oxygen plays an important role in the cell cycle during the proliferation phase. The VEGF secreted by the tumor cells induces angiogenesis. Meanwhile, TGFα, glucose, oxygen and VEGF diffuse continuously in the tumor microenvironment. The tissue scale focuses on angiogenesis (Figure ). Blood vessel sprouts migrate and branch via tip endothelial cells' migration in response to VEGF chemotaxis and fibronectin haptotaxis. The new blood vessels supply tumor cells with glucose and oxygen to maintain their metabolism and further invasion.
Figure 1 Flow chart of multi-scale agent-based cancer modeling. The model encompasses four biological scales: the molecular scale, the cellular scale, the micro-environmental scale and the tissue scale. The molecular scale consists of the EGFR signaling and cell (more ...)
EGFR signaling pathway connected to cell cycle pathway with TKIs to block the EGFR signaling pathway. EGFR signaling pathway is connected to the cell cycle pathway, and TKIs block the EGFR signaling pathway.
Illustrations of phenotype switch "decision" of tumor cell. Please see details in the main text (Cellular scale: Phenotype switch of tumor cell as "agent").
Probability of migration of tip endothelial cell. The tip endothelial cell probabilistically moves up, down, right, or left, or stays at its current position.
We performed our simulations on a two-dimensional square lattice. The lattice size was set to L
200 representing a 4
5 mm length of a brain tissue slice. The lattice spacing is 20μm
, which is approximately the diameter of tumor cells. Initially a parent blood vessel along with 6 tip endothelial cells is located near the left boundary of the computational domain. Also a small cluster of active tumor cells were randomly distributed as shown in Figure . Each tumor cell is initialized with its age being between 0
24 hours randomly. The time step of the simulation is one hour. The details of the model parameters are in the Additional files 1
, and 5
of this paper.
Molecular scale: signaling pathway
The concentration of each component (Figure ) in the EGFR
] and cell-cycle pathways [21
] are described by a system of coupled ordinary differential equations (ODE
s) in the form
represents the production rate of Xi
_ the consumption rate. The parameters of ODE
s are in Additional files 1
, and 4
. Parameter sensitivity and model robustness analysis are presented in the results section. These methods enabled us to determine which parameters of the system are more sensitive as well as to examine if the system is robust enough to the small parameter perturbations.
Cellular scale: phenotype switch of tumor cells as "agents"
At every simulation step, each "agent" (i.e. a brain tumor cell) switches its phenotype according to the following rules:
1. First, each agent evaluates the concentration of glucose at its current location. If the concentration is greater than the cell active threshold, the agent becomes active and uses its EGFR signaling pathway to determine its phenotype. If the concentration of glucose is less than the dead threshold, the cell dies. If the concentration is between active and dead thresholds, the agent enters into a reversible quiescent state [10
2. Each active agent evaluates its migration potential (MP) by the following equation:
where d[PLCγ]/dt is the rate of change of the PLCγ concentration. If MP is greater than a threshold σPLCγ, the average rate of change of PLCγ concentration, the agent will choose the migration phenotype.
3. If MP is less than σPLCγ, the agent starts to proliferate. If the concentration of CDh1 is less than a threshold thr1 and the concentration of cycCDk is greater than the threshold thr2, the cell divides. After that, the cell chooses the most attractive free site
(detailed in equation (3)) in the neighborhood to deliver its offspring. If there is no empty neighborhood, the cell turns into a reversible quiescent state until free space becomes available.
4. Each agent chooses the "most attractive" location mentioned above according to the following probability:
is the glucose concentration at location j, Fj
is the fibronectin concentration at j
~N(0,1) is a normally distributed error term, the parameter ψ
(0, 1) represents the extent of the search precision, which is set to 0.7 [22
Microenvironmental scale: extracellular chemotaxis
Five extracellular micro-environmental factors, glucose, oxygen, TGFα, VEGF and fibronectin are included in this model. A set of reaction–diffusion equations describe the diffusion, penetration and uptake of glucose, oxygen and VEGF.
Glucose first penetrates blood vessels, and then diffuses in the extracellular microenvironment. After that, it is consumed by the tumor cells. This process is modeled by the following equation:
is the glucose concentration,
is the Laplace operator, DG
is the diffusivity of glucose. qG
, where pG
is the vessel permeability for glucose and r
is the blood vessels' average radius. In addition, Gblood
is the glucose concentration in blood and UG
is the cell’s glucose uptake rate. The time dependent characteristic function Xves
) is equal to 1, if a blood vessel is present at x
; otherwise it is equal to 0. Xtum
) is equal to 1 in the tumor region and is equal to 0 elsewhere. Xves
are updated at each simulation step according to the developing profile of the tumor and its micro-vascularity.
Oxygen also permeates the blood vessels' walls, diffuses in the surrounding and is consumed by tumor cells. This process is modeled by the following equation:
where C is the oxygen concentration, DC is the oxygen diffusivity, qC is the vessel permeability for oxygen, and UC is a cell’s uptake rate of oxygen.
TGFα, an analogue of EGF, is secreted by tumor cells and can be paracrine and juxtacrine [23
]. Equation (6) describes the diffusion and secretion of TGFα.
where T is the TGFα concentration, DT is its diffusivity, qT is vessel permeability to TGFα. ST is a cell’s net production rate of TGFα and δT is the natural decay rate of TGFα.
We applied homogeneous Neumann boundary conditions for all the above equations by assuming zero flux along the boundary of the considered domain. Additional file 5
: Table A5 and Additional file 6
: Equations A1–A5 list the parameters and initial conditions of the equations.
We solved these equations numerically with the finite difference method [24
Tissue scale: angiogenesis
Tumor induced angiogenesis is due to the secretion of VEGF by the tumor cells. VEGF diffuses into the surrounding corneal tissue and is also consumed by the endothelial cells [25
]. We model this process with the following equation:
where V is the VEGF concentration, DV is the diffusivity of VEGF, qV is the vessel permeability for VEGF and SV is a cell’s VEGF secretion rate. δV is the natural decay rate of VEGF.
Fibronectin is a component of the corneal tissue secreted by endothelial cells. In addition, tumor cells can consume fibronectin. This process is described by the following equation:
where F is the fibronectin concentration, β and γ are positive constants representing the production and uptake rates, respectively.
We assume that the motion of individual endothelial cell (EC) located at the tip of a capillary sprout governs the motion of the whole sprout. Chemotaxis in response to VEGF gradients and haptotaxis in response to fibronectin are the major factors that influence the motion of the endothelial cells at the capillary sprout tip [25
We defined the probability of migration of a tip endothelial cell (see Figure ) as
is the VEGF concentration, F
is fibronection concentration and lk
is the directional vector along the k
th direction. The term V
models chemotaxis in response to VEGF gradients [19
], where α is the chemotactic coefficient and kv
is a positive constant controlling the weight of VEGF concentration in chemotactic sensitivity. The term λF
models the haptotatic influence of fibronectin on the endothelial cells, where λ
is the haptotatic coefficient.
For a given tip of endothelial cells at location (i, j), the un-normalized migration probabilities can be calculated from equation (9) as follows:
The un-normalized probability P5, for a tip cell to remain stationary is the average of P1, P2, P3 and P4. After normalization, the above equations give the likelihood of the tip endothelial cell to move up, down, right, or left, or stay at its current position. The probability, P5, for a tip cell to remain stationary is the average of P1, P2, P3 and P4.
The algorithm related to angiogenesis is as follows:
1. Calculate the migration probabilities of ECs:
1.1 Solve the equations (7) and (8), and then calculate P1-P5 from (10);
1.2 Normalize the above numbers:
; define intervals
2. For every sprout tip cell, we check whether the age of vessel is greater than 18 hours and whether there are any free sites in its nearest neighborhood.
2.1 Sprout branching: If the above conditions are satisfied, two random numbers r1
between 0 and 1 are generated. If r1 I2
and r2 I3
, then we move two endothelial cells one below and one to the right of the spout tip endothelial cell.
2.2. Sprout migrating: If the above branching conditions are not satisfied, we generate another random number r
between 0 and 1. If r I3
, we move the tip endothelial cell to the right of spout tip endothelial cell.
3. Anastomosis: If two sprouts encounter each other, a new sprout continues to grow.
This study uses TKIs, particularly gefitinib, as inhibitors delivered by capillary vessels to treat brain cancer. The TKIs are modeled as a continuous concentration field. These molecules permeate the blood vessels and diffuse continuously in the microenvironment to produce an accumulative inhibitive effect on the tumor growth. We describe the evolution of the concentration of TKIs by the following equation:
where DTKi is the diffusivity of the TKIs, qTKi is the vessel permeability for TKIs, TKiblood is the blood TKIs concentration, UTKi is a cell’s uptake rate of TKIs, and δTki is the natural decay rate of TKIs.
The EGFR signaling pathway controls a cell’s phenotypic switch by binding TGFα molecules to the EGF receptors. Because the TKI molecules inhibit the autophosphorylation receptors, the downstream EGFR pathway responsible for a cell's phenotypic switch remains inactivated [26
]. The binding and unbinding processes of TKIs to EGFR with constant rates kb
, respectively, are described by the following equation:
Since the chemical reaction rate is much faster than the cells' phenotypic switch [26
], the concentration of EGFR: TKI complex in equation (12) can be obtained by applying Michaelis–Menten kinetics as follows:
where km is the Michaelis constant, km ≈ kb / ku, and [EGFR]0 is the initial concentration of the EGFR. We can then derive the effective amount of EGFR for the activation of downstream factors as follows:
Because of the TKI treatment, the effective amount of EGFR of some tumor cells will decrease. The decrease of the amount of effective EGFR results in a slow rate of change of PLCγ concentration. This in turn inhibits tumor progression by reducing the migration potential of these tumor cells (see equation 2 and detailed phenotype change of tumor cells in Figure ).
Finally, we summarize our computing algorithm at each step across multi-scales (Figure ) as follows. At micro-environmental scale, we solve the PDEs (equations 4–6) to obtain the spatial concentration distributions of glucose, oxygen and TGFα. At molecular scale, we use the calculated local TGFα concentration as the input for EGFR signaling pathway (equation 1) for each tumor cell. At cellular scale, tumor cells' migration potential (MP) (equation 2) is computed to determine their phenotypic switch (migration or proliferation); meanwhile, other phenotypic switches (quiescent or apoptosis) are associated with the current value of oxygen and glucose. At tissue scale, the spatial concentration distributions of VEGF and fibronectin (equations 7–8) will guide the tip endothelial cells' migration and sprout branching. In turn, the remodeled vasculature at tissue scale influences the spatial concentration distributions of glucose, oxygen and TGFα at micro-environmental scale. For TKI treatment, the TKI distribution is integrated into molecular scale by solving equation 11 along with aforementioned equations 4–6, and the initial value of EGFR is varied by equations 12–14 as well.