Skip to main content

Multi-scale agent-based brain cancer modeling and prediction of TKI treatment response: Incorporating EGFR signaling pathway and angiogenesis

Abstract

Background

The epidermal growth factor receptor (EGFR) signaling pathway and angiogenesis in brain cancer act as an engine for tumor initiation, expansion and response to therapy. Since the existing literature does not have any models that investigate the impact of both angiogenesis and molecular signaling pathways on treatment, we propose a novel multi-scale, agent-based computational model that includes both angiogenesis and EGFR modules to study the response of brain cancer under tyrosine kinase inhibitors (TKIs) treatment.

Results

The novel angiogenesis module integrated into the agent-based tumor model is based on a set of reaction–diffusion equations that describe the spatio-temporal evolution of the distributions of micro-environmental factors such as glucose, oxygen, TGFα, VEGF and fibronectin. These molecular species regulate tumor growth during angiogenesis. Each tumor cell is equipped with an EGFR signaling pathway linked to a cell-cycle pathway to determine its phenotype. EGFR TKIs are delivered through the blood vessels of tumor microvasculature and the response to treatment is studied.

Conclusions

Our simulations demonstrated that entire tumor growth profile is a collective behaviour of cells regulated by the EGFR signaling pathway and the cell cycle. We also found that angiogenesis has a dual effect under TKI treatment: on one hand, through neo-vasculature TKIs are delivered to decrease tumor invasion; on the other hand, the neo-vasculature can transport glucose and oxygen to tumor cells to maintain their metabolism, which results in an increase of cell survival rate in the late simulation stages.

Background

Brain cancer is a very complex and deadly disease. Traditional diagnoses and treatments of this disease are from in vitro experimental observations. Although biologists have developed many experimental data at the molecular, cellular, micro-environmental and tissue scales, only very few scientists have integrated these data into multi-scale models to study tumor response to treatment.

Cellular automata (CA) methods have been widely applied to model brain tumor growth [1, 2]. Although CA models are good at describing cell-cell and cell-microenvironment interactions, this type of discrete modelling approach falls short on investigating most fluid dynamic aspects of the tumor microenvironment. Alternatively, Continuum models employ systems of partial differential equations to simulate the solid tumor invasion by updating boundaries of different sub-domains of tumor based on the level-set method [3, 4]. It is, however, hard with this approach to describe cell-cell interactions, such as the competition among cells for nutrients. In general neither continuum nor discrete models can accurately simulate cancer spatio-temporal evolution with respect to the complexity of cancer. A hybrid discrete-continuum (HDC) model that couples a cellular automaton module with a continuum module was proposed by Anderson and co-workers [5, 6]. Although these models have not considered the effects of any gene-protein signaling pathways such as the EGFR pathway, computational cancer biologists have already studied them extensively. Recent studies showed that EGFR pathway plays an important role in the evolution of brain cancer [79]. Therefore, the existing multi-scale agent-based tumor models incorporated an EGFR signaling pathway [10, 11] at the molecular scale to enable individual cells to choose their phenotypic trait between proliferation and migration based on the pathway's state [10, 12, 13]. As indicated by Hanahan and Weinberg[14], angiogenesis is a significant transforming phase in tumor growth. During the angiogenesis phase, tumor cells secret vascular endothelial growth factors (VEGF) [15, 16] into the microenvironment to induce and sustain new capillary sprouts migrating from pre-existing vasculature towards the tumor. This in turn helps to maintain tumor cells’ metabolism by supplying them with glucose and oxygen and subsequently leads to metastasis. However, our previous agent-based models [1013, 17] did not explicitly take tumor-induced angiogenesis into consideration.

In this paper, we presented a novel multi-scale agent-based model to describe tumor growth with angiogenesis and study the response of brain cancer to EGFR tyrosine kinase inhibitors (TKIs) [18]. Several rates of changes of molecular species, such as PLCγ, CDh1 and cycCDK will determine a tumor cell’s phenotypic switch.

In order to integrate an angiogenesis module into the existing agent-based tumor growth models [10, 11, 13, 17], we developed a set of rules that underline the migration of endothelial cells and the branching of vessel sprouts. Compared to previous developed HDC rules [19], these rules which directly defined the probabilities of migration of endothelia cells are more suitable for implementation into an agent-based model. The tumor growth and angiogenesis are coupled through VEGF secreted by the tumor cells and through the glucose and oxygen permeated from the neo-vasculature. As the neo-vasculature develops, glucose and oxygen penetrate from the blood vessels and diffuse throughout the tumor microenvironment to promote further tumor growth, which in turn influences the concentration of VEGF.

The simulation results demonstrate that we can investigate the response of brain cancer to tyrosine kinase inhibitors (TKIs) and also can use the model to reveal the dual role of angiogenesis.

Implementation

Our model encompasses four biological scales (Figure 1): 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 2). 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 3). 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 4). 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
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 cycle pathways (Figure 2). The cellular scale simulates the phenotypic switch of the tumor cells as determined by the molecular scale, the cell-cell and the cell-microenvironment interactions (Figure 3). The micro-environmental scale provides a connection between the micro-vasculature (angiogenesis) and the tumor cells. Tumor cells consume glucose and oxygen nutrients. The tumor growth factor alpha (TGFα) secreted by tumor cells triggers the EGFR signaling pathway, which in turn determines cell’s phenotype switch. Oxygen plays an important role in the cell cycle during the proliferation phase. The VEGF secreted by the tumor cells induces angiogenesis. The tissue scale focuses on angiogenesis (Figure 4). Blood vessel sprouts migrate and branch via tip endothelial cells' migration in response to VEGF chemotaxis and fibronectin haptotaxis.

Figure 2
figure 2

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.

Figure 3
figure 3

Illustrations of phenotype switch "decision" of tumor cell. Please see details in the main text (Cellular scale: Phenotype switch of tumor cell as "agent").

Figure 4
figure 4

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 4. 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, 2, 3, 4, and 5 of this paper.

Molecular scale: signaling pathway

The concentration of each component (Figure 2) in the EGFR signaling [11, 20] and cell-cycle pathways [21] are described by a system of coupled ordinary differential equations (ODE s) in the form

d X i / d t = v + v _
(1)

where v+ represents the production rate of X i and v _ the consumption rate. The parameters of ODE s are in Additional files 1, 2, 3, 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. 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. 2.

    Each active agent evaluates its migration potential (MP) by the following equation:

    MP PLC γ = d PLC γ / dt,
    (2)

    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. 3.

    If MP is less than σPLCγ, the agent starts to proliferate. If the concentration of CDh1 is less than a threshold thr 1 and the concentration of cycCDk is greater than the threshold thr 2 , 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. 4.

    Each agent chooses the "most attractive" location mentioned above according to the following probability:

    P j = ψ G j / F j + 1 ψ ε j ,
    (3)

    where Gj is the glucose concentration at location j, F j is the fibronectin concentration at j, ε 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:

G t = D G Δ G + X ves t , x q G G blood - G X tum t , x U G ,
(4)

where G is the glucose concentration, Δ 2 is the Laplace operator, D G is the diffusivity of glucose. qG = 2πrp G , where p G 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 X ves (t,x) is equal to 1, if a blood vessel is present at x; otherwise it is equal to 0. X tum (t,x) is equal to 1 in the tumor region and is equal to 0 elsewhere. X ves and X tum 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:

C t = D C Δ C + X ves t , x q C C blood C X tum t , x U C ,
(5)

where C is the oxygen concentration, D C is the oxygen diffusivity, q C is the vessel permeability for oxygen, and U C 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α.

T t = D T Δ T X ves t , x q T T + X tum t , x S T δ T T ,
(6)

where T is the TGFα concentration, D T is its diffusivity, q T is vessel permeability to TGFα. S T 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:

V t = D V Δ V X ves t , x q V V + X tum t , x S V δ V V ,
(7)

where V is the VEGF concentration, D V is the diffusivity of VEGF, q V is the vessel permeability for VEGF and S V 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:

F t = X ves t , x β X tum t , x γ F ,
(8)

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 4) as

P k α k v k v + V V + λ F . l k , k = 1 , 2 , 3 , 4 ,
(9)

where V is the VEGF concentration, F is fibronection concentration and l k is the directional vector along the k th direction. The term α k v / k v + V V models chemotaxis in response to VEGF gradients [19], where α is the chemotactic coefficient and k v 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:

P 1 = α k V / k V + V i , j ) · ( V i , j + 1 V i , j + λ ( F ( i , j + 1 ) F ( i , j ) ) , P 2 = α k V / k V + V i , j ) · ( V i , j 1 V i , j + λ ( F ( i , j 1 ) F ( i , j ) ) , P 3 = α k V / k V + V i , j ) · ( V i + 1 , j V i , j + λ ( F ( i + 1 , j ) F ( i , j ) ) , P 4 = α k V / k V + V i , j ) · ( V i 1 , j V i , j + λ ( F ( i 1 , j ) F ( i , j ) ) .
(10)

The un-normalized probability P 5 , 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, P 5 , 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. 1.

    Calculate the migration probabilities of ECs:

    1. 1.1

      Solve the equations (7) and (8), and then calculate P1-P5 from (10);

    2. 1.2

      Normalize the above numbers: P ˜ i = P i / j = 1 5 P j , . i = 1 , 2 , , 5 ; define intervals I 1 = ( 0 , P ˜ 1 ] , I i = ( j = 1 i - 1 P ˜ j , j = 1 i P ˜ j ] , i = 2 , , 5 .

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

    1. 2.1

      Sprout branching: If the above conditions are satisfied, two random numbers r1 and r2 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.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. 3.

    Anastomosis: If two sprouts encounter each other, a new sprout continues to grow.

TKI treatment

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:

T k i t = D Tki Δ T k i + X ves t , x q Tki T k i blood T k i X tum t , x U Tki δ Tki T k i ,
(11)

where D TKi is the diffusivity of the TKIs, q TKi is the vessel permeability for TKIs, TKiblood is the blood TKIs concentration, U TKi 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 k b and k u , respectively, are described by the following equation:

E G F R + T K I k u k b E G F R : T K I
(12)

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:

E G F R : T K I = E G F R 0 T K I k m + T K I
(13)

where k m is the Michaelis constant, k m k b / k u , 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:

E G F R eff = E G F R 0 E G F R : T K I
(14)

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 3).

Finally, we summarize our computing algorithm at each step across multi-scales (Figure 1) 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.

Results

We have implemented the above model into software "ABM-TKI" in the Matlab programming environment. “ABM-TKI” is a tool employing agent-based model (ABM) to simulate brain tumor growth. It includes an EGFR signaling pathway, a related cell-cycle, angiogenesis and TKIs treatment. We can employ this tool to predict the responses of brain cancer and reveal the dual roles of angiogenesis under TKI treatment.

Regarding software usage, the user can download and decompress the package from the project home page (https://sites.google.com/site/agentbasedtumormodeling/home). Then the user can run the program in Matlab (version R2007b or higher) with input as: angiog_tumor(time,isdrug).The input "time" is the period from the beginning of the simulation to the end. The "isdrug = 1" means TKI treatment and "isdrug = 0" means no TKI treatment. For example: "angiog_tumor(150,0)" will give the tumor growth profile without TKI treatment from 0 hour to 150 hours; "angiog_tumor(300,1)" will give the tumor growth profile with TKI treatment from 0 hour to 300 hours.

The output includes: (a) the vascular tumor growth pattern with or without TKI treatment; (b) the tumor growth visualization with the background of fibronectin; (c) the spatio-temporal evolution of the concentration of glucose, oxygen, TGFα and/or TKI; (d) various tumor cell numbers such as active cells, apoptotic cells, migratory cells, proliferative cells, quiescent cells and the number of endothelial cells; (e) the average change rate of PLCγ with or without TKI treatment.

Vascular tumor growth patterns

The vascularized tumor patterns at 60, 90, 120 and 150 hours are shown in Figure 5a, where different colors denote various tumor cell states: active (blue), quiescent (cyan), dead (black) and proliferative (pink). The endothelial cells are red. (See also Additional file 7: Figure A1 for more details). The active tumor cells are comprised of migratory cells, cells that have just completed their proliferation process, and cells that just switched back from quiescent state to the active state. Additional file 8: Figure A2 shows the growth of the tumor microvascularity and the fibronectin concentration at different time interval. Figure 5a showed that those tumor cells tend to migrate to locations near the vessels, where the glucose concentration level is the highest. Around t = 60 hours, the tumor was comprised of active cells and some quiescent cells. Around t = 90 hours some dead cells appeared in the tumor mass. From t = 120 hours to t = 150 hours, the number of both dead and active cells kept increasing. At t = 150 hours, we found that some quiescent cells near blood vessels switched back to active cells. At this time, the tumor developed a tree branching microvasculature, which is much denser near the tumor. Figure 5b shows the distributions of different micro-environmental factors (glucose, TGFα, oxygen, VEGF) at t = 150 hours. Given the space limitation, we show the results at a single time interval. The concentrations of these molecules at different times are shown in the Additional file 9: Figure A3, with different colors representing different levels of concentration. We found that the descending order of the diffusion rate for the four molecular species is as follows: oxygen, glucose, VEGF or TGFα.

Figure 5
figure 5

Simulation of vascular tumor growth without TKI treatment. (a) The evolution of vascular tumor patterns at 60, 90, 120 and 150 hours. (b) The distribution of glucose, oxygen, TGFα and VEGF in microenvironment at 150 hours. (c) Different types of cell numbers from 0 hour to 150 hours. (d) Proliferation rate of tumor cells derived from simulations (red line) and from in vitro experimental observations (blue line) for t=0-96 hours.

Figure 5c shows the numbers of different types of cells as a function of time. The number of active cells increased monotonically with time. The number of apoptotic cells increased abruptly at around t = 60 hours and kept increasing until the tumor microvasculature developed at t = 150 hours. The number of quiescent tumor cells kept increasing from t = 0 to 130 hours, and then began to decrease. The number of endothelial cells increased rapidly during the whole simulation time. The detailed evolutions of the numbers of various cells are shown in Additional file 10: Figure A4 separately.

In Figure 5d, we present the proliferation rate of tumor cells as a function of time from our simulation and from in vitro experimental results (at t = 96 hours) [27]. The in vitro experimental data are from human glioma tumor-initiating cells derived from 7 patients (GBM 1–7). The plotted experimental data are the mean and standard deviation values from the seven cell lines. We took the mean value of these data as the blue line in Figure 5d. The mean squared error of our prediction is 0.1421. The simulation and experimental data in Figure 5d are in very good agreement, which is an important validation of our model.

TKI treatment response

In our study we chose gefitinib as the EGFR TKI to treat brain cancer. In this section the term TKI refers to treating cell by gefitinib. Figure 6a shows the simulation profiles of the vascularized tumor growth with TKI treatment from t = 0 to 300 hours. The different colors represent cell types as defined in Figure 5a. The results in Figure 6a show that the tumor expansion slowed down obviously by the TKI treatment. Additional files 11- 12: Figure A5-6 show the detailed vascular tumor growth profile with TKI treatment as well as that in the presence of fibronectin at different time intervals.

Figure 6
figure 6

Simulation of vascular tumor growth with TKI treatment. (a) Vascular tumor growth patterns with TKI treatment at 60, 150, 240 and 300 hours; (b) The distribution of TKI in tumor vasculature at 300 hours; (c) Average change rate of PLCγ from hour 0 to 300 hours.

Figure 6b shows the distribution profile of TKIs concentration at t = 300 hours which is similar to the structure of tumor microvasculature. The Additional file 13: Figure A7 and Additional file 14: Figure A8 demonstrate the evolution of distributions of glucose, oxygen, TGFα and VEGF as well as TKIs during the treatment at different time intervals.

Figure 6c shows the average rate of change of PLCγ as a function of time. The PLCγ average rate of change increases from t = 0 to 60 hours and then it starts decreasing with an exception at around t = 75 hours where the data show a hump. Finally, the curve goes down after t = 125 hours. Since Additional file 10: Figure A4 shows the average rate of change of PLCγ always increases without TKI treatment, TKI treatment greatly affects the average rate of change of PLCγ. The numbers of various cells with TKI treatment are shown in Additional file 15: Figure A9.

Figure 7a shows the cell survival rate under TKI treatment from t = 0 to 300 hours which decreased until t = 150 hours and then increased again until t = 300 hours.

Figure 7
figure 7

Prediction and validation of cell survival percentage during TKI treatment. (a) Prediction of cell survival percentage during TKI treatment from 0 hour to 300 hours; (b) Cell survival percentage from simulations (purple line) and from in vitro experimental observations (black line) for t=0-96 hours.

Figure 7b demonstrates that the simulated cell survival rate has a trend similar with experimental results. The purple line in the figure represents the average from a hundred simulations. Human glioma tumor-initiating cells are derived from 7 patients (GBM 1–7) [27] for in vitro experiments observed for 96 hours, these experimental data are shown by the multiple lines in the figure. In the experiment human glioma tumor growth inhibitors of gefitinib are used as TKI at the concentration of 1 μM, while in our simulation we also chose gefitinib as TKI and its concentration near tumor region is also close to 1 μM from t = 0 to 100 hours. The relatively good agreement between the simulation prediction and the experimental results constitutes an important validation of our model.

Sensitivity analysis and model robustness

Parameter sensitivity analysis is to quantitatively discover sensitive parameters in the system. And robustness analysis is to examine whether the system is stable to modest fluctuations of these sensitive parameter values.

In this work, local parameter sensitivity analysis [2830] was employed to understand the relationship between the migration potential (defined above) of tumor cell and the variations in individual parameter values. The sensitivity coefficient (S) is calculated according to the following formula [31]:

S = Δ M P M P / Δ P P ,
(15)

where P is the parameter that is varied and MP is tumor cell's migration potential; ΔMP is the corresponding change in MP due to a small change in P denoted by ΔP. Each individual parameter is increased (or decreased) by 10% from its original value, and the resulting percentage change of the migration potential is computed. Figure 8 shows that the migration potential of the cells is most sensitive to the initial concentration of EGFR (X2) and TGFα (X1). The sensitivity analysis also demonstrates that kinetic rates, k1, k2, k3, k5, K4 and V4, are more critical than others in the system. Furthermore, it turns out that the developed intracellular pathway system is rather robust since all of the sensitivity coefficients are less than 1.8%.

Figure 8
figure 8

Sensitivity analysis of EGFR signaling pathway. The result shows that the migration potential of brain cancer cells is most sensitive to the initial concentration of EGFR (X2) and TGFα (X1). The result also shows us forward kinetic rates (k1, k2, k3, k5) are more important in the system and demonstrate that the developed intracellular pathway system is robust.

Based on above sensitivity analysis results, we selected 3 pathway components, TGFα (X1), EGFR (X2) and PLCγ (X6), and 6 kinetic rates including k1, k2, k3, k5, K4 and V4 as key parameters to investigate the robustness of the model to much larger ranges of parameter variations. For each of these parameters, we varied its value by 0.1-fold, 2.0 fold, 5.0 fold, 10.0 fold, 50.0-fold, and 100.0-fold of its original reference value respectively. Each time, only one of the above parameters was varied with a corresponding fold change, and all other parameters were kept fixed. We chose active cell number (ACN) at 100 hours in each simulation as the system outcome, and examined the system responses to the above parameter variations. We defined the robustness index (R) as follows:

R = A C N p A C N 0 A C N 0 ,
(16)

where ACN p is the active cell number with the varied parameter p, and ACN 0 is the average active cell number from the 100 simulations with unvaried reference parameters.

Figure 9 shows that for variations between 0.1-fold and 10.0-fold, the absolute values of the robustness index is less than 0.1894 with respect to most parameters (except TGFα and k1). Figure 10 shows the detailed analyses of the system robustness with respect to TGFα and k1 at variations between 5.0-fold to 10.0-fold. The results show that regarding TGFα and k1, at the variations between 0.1-fold and 7.0-fold, the absolute values of the robustness indexes are less than 0.1894. Furthermore, for some parameters, such as PLCγ, k2, k3, k5, K4, the robustness index is very small even with the variations up to 100.0-fold. Then we calculated the coefficient of variation (CV) of system outcomes which is defined as the ratio of the standard deviation to the mean of the above variational active cell numbers. At the parameter variations between 0.1-fold and 10.0-fold, the CV of system outcomes is 0.2007 which is much less than 1, indicating that the system is relatively stable to the parameter variations. These results reveal that our model is comparably robust with respect to relatively large ranges of parameter variations.

Figure 9
figure 9

Model robustness analysis. TGFα (X1), EGFR (X2), PLCγ (X6), and k1, k2, k3, k5, K4 and V4 were varied their values by 0.1-fold, 2.0 fold, 5.0 fold, 10.0 fold, 50.0-fold, and 100.0-fold of original reference values each time. The results show that the absolute values of the robustness indexes are less than 0.1894 with respect to most parameters (except TGFα and k1) for variations between 0.1-fold and 10.0-fold.

Figure 10
figure 10

Detailed robustness analyses for TGFα and k 1 . TGFα (X1), and k1 were varied by 5.0-fold, 6.0 fold, 7.0 fold, 8.0 fold, 9.0-fold, and 10.0-fold of original reference values each time. The results show the absolute values of the robustness indexes are less than 0.1894 with respect to TGFα and k1 for variations between 0.1-fold and 7.0-fold.

Discussion

We developed a multi-scale model by integrating a novel angiogenesis module into an agent-based tumor model based on a set of reaction–diffusion equations that describe the spatio-temporal evolution of the distributions of micro-environmental factors such as glucose, oxygen, TGFα, VEGF and fibronectin. These molecular species regulate tumor growth during angiogenesis. Each tumor cell is equipped with an EGFR signaling pathway linked to a cell-cycle to determine its phenotype.

Our simulations show several interesting findings. The first is that tumor cells tend to move towards blood vessels and gradually developed to a fan-shape as shown in Figure 5a. We interpret this result as follows. Since the nutrients (glucose and oxygen) concentrations are higher at locations near the blood vessels, they attract tumor cells.

The second interesting finding is that blood vessels tend to migrate to tumor and form a dense tree-branching vascular network. The reason is that a high VEGF gradient close to the tumor attracts endothelial cells, which in turn lead to branching of vessels in these regions.

The third interesting finding is that TKI treatment can inhibit tumor progression. The binding of TKI molecules to EGFR decreases the amount of effective EGFR, which results in low expression of PLCγ and low cell’s migration potential (Figure 6c). As a result, tumor invasion slows down.

The fourth interesting result is that the tumor cells' survival rate does not always decrease. This is due to the dual role of angiogenesis. Newly formed capillaries delivers a substantial amount of TKI molecules to tumor cells and blocks the EGFR signaling pathway, which lead to an inhibition of tumor growth. This in turn results in a decreased cell survival rate in the early stage of the tumor development. On the other hand, new capillaries transported a lot of glucose and oxygen to tumor cells which results in an increased cell survival rate at later stages (Figure 7a). The implications of the dual roles of angiogenesis reveal that clinical personnel should decrease cancer progression by using TKI treatment and inhibiting tumor-induced angiogenesis at the same time.

The sensitivity analysis reveals sensitive parameters in the EGFR signaling pathway. The robustness study confirms that our model is relatively robust and stable to fluctuations of these sensitive parameters.

Herein we used the two-dimensional in vitro experiments [27] to validate the effectiveness of the model. These experiments employed two-dimensional experimental protocol to isolate and plate human glioma tumor-initiating cells in Martrigel-coated culture flasks. Figure 5d and Figure 7b demonstrate that our in silico model does have strong predictive power and great potential for clinical work.

We are going to extend the model to three dimensions to simulate in vivo tumor growth for real clinical purposes. A three-dimensional lattice is indispensable for the simulation of in vivo tumor growth with angiogenesis, because tumor cells' activities, vasculature structure, and chemical cues' diffusion in three-dimensional heterogeneous tumor growth environment are very different from two-dimensional. Moreover, three-dimensional simulations require parallel computing techniques [32] to relieve the heavy computing request.

The potential of our model will further increase, by incorporating more realistic biological and physical features, such as blood flow and tumor growth-induced pressure [33] in the future.

Conclusions

This work presents a novel multi-scale agent-based brain tumor model encompassing an EGFR signaling pathway together with a related cell-cycle, an angiogenesis module and TKI treatment. It incorporates four relevant biological scales: the molecular scale, the cellular scale, the microenvironment scale and the tissue scale. At the molecular scale, a system of ordinary differential equations simulates the dynamics of the EGFR signaling pathway and the cell cycle to determine the cells' phenotypic switch at the cellular scale. We employed a set of partial differential equations to simulate the concentration changes of five extracellular chemical cues (glucose, oxygen, TGFα, VEGF and fibronectin) in the tumor micro-environmental scale. Angiogenesis was coupled into tumor growth through VEGF secreted by the tumor cells and through the glucose and oxygen permeated from the neo-vasculature at the tissue scale. Moreover, we integrated TKI treatment into EGFR signaling pathway to block the activation of EGFR.

Our simulations demonstrate that the entire tumor growth profile is a collective behaviour of its cells regulated by the EGFR signaling pathway and the cell cycle. We also discovered that angiogenesis has dual effects on TKI treatment: on one hand, neo-vasculature can deliver TKIs to decrease the tumor invasion, whereas on the other hand, it can transport a lot of nutrients ( glucose and oxygen) to tumor cells to maintain their metabolism, which results in an increase of cell survival rate at late simulation stage. There is a great similarity between the simulation results and existing in vitro experimental data. Further analyses show that our model has strong robustness regarding to the relatively large changes of the sensitive model parameters.

Availability and requirements

Project name: multi-scale agent-based brain tumor modeling project Project home page:http://www.methodisthealth.com/Softwarehttp://csysbio.org/Released%20Software.htmlhttps://sites.google.com/site/agentbasedtumormodeling/home Operating system(s): Platform independent Programming language: Matlab (R2007b) Other requirements: None License: GNU GPL, FreeBSD etc. Any restrictions to use by non-academics: license needed.

Abbreviations

EGFR:

Epidermal Growth Factor Receptor

VEGF:

Vascular Endothelial Growth Factor

TKIs:

Tyrosine Kinase Inhibitors

EC:

Endothelial Cell

MP:

Migration Potential.

References

  1. Patel AA, Gawlinski ET, Lemieux SK, Gatenby RA: A cellular automaton model of early tumor growth and invasion: the effects of native tissue vascularity and increased anaerobic tumor metabolism. J Theor Biol 2001, 213(3):315–331. 10.1006/jtbi.2001.2385

    Article  CAS  PubMed  Google Scholar 

  2. Kansal A, Torquato S, Harsh Iv G, Chiocca E, Deisboeck T: Cellular automaton of idealized brain tumor growth dynamics. Biosystems 2000, 55(1–3):119–127.

    Article  CAS  PubMed  Google Scholar 

  3. Zheng X, Wise SM, Cristini V: Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method. Bull Math Biol 2005, 67(2):211–259. 10.1016/j.bulm.2004.08.001

    Article  CAS  PubMed  Google Scholar 

  4. Macklin P, McDougall S, Anderson ARA, Chaplain MAJ, Cristini V, Lowengrub J: Multiscale modelling and nonlinear simulation of vascular tumour growth. J Math Biol 2009, 58(4):765–798. 10.1007/s00285-008-0216-9

    Article  PubMed Central  PubMed  Google Scholar 

  5. Anderson ARA: A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. Math Med Biol 2005, 22(2):163–186. 10.1093/imammb/dqi005

    Article  PubMed  Google Scholar 

  6. Rejniak KA, Anderson ARA: Hybrid models of tumor growth. Wiley Interdiscip Rev Syst Biol Med 2011, 3(1):115–125. 10.1002/wsbm.102

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  7. Chen P, Xie H, Sekar MC, Gupta K, Wells A: Epidermal growth factor receptor-mediated cell motility: phospholipase C activity is required, but mitogen-activated protein kinase activity is not sufficient for induced cell movement. J Cell Biol 1994, 127(3):847–857. 10.1083/jcb.127.3.847

    Article  CAS  PubMed  Google Scholar 

  8. Piccolo E, Innominato PF, Mariggio MA, Maffucci T, Iacobelli S, Falasca M: The mechanism involved in the regulation of phospholipase Cgamma1 activity in cell migration. Oncogene 2002, 21(42):6520–6529. 10.1038/sj.onc.1205821

    Article  CAS  PubMed  Google Scholar 

  9. Wang Z, Gluck S, Zhang L, Moran MF: Requirement for phospholipase C-gamma1 enzymatic activity in growth factor-induced mitogenesis. Mol Cell Biol 1998, 18(1):590–597.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  10. Zhang L, Athale CA, Deisboeck TS: Development of a three-dimensional multiscale agent-based tumor model: Simulating gene-protein interaction profiles, cell phenotypes and multicellular patterns in brain cancer. J Theor Biol 2007, 244(1):96–107. 10.1016/j.jtbi.2006.06.034

    Article  CAS  PubMed  Google Scholar 

  11. Zhang L, Chen LL, Deisboeck TS: Multi-scale, multi-resolution brain cancer modeling. Math Comput Simulat 2009, 79(7):2021–2035. 10.1016/j.matcom.2008.09.007

    Article  Google Scholar 

  12. Wang ZH, Zhang L, Sagotsky J, Deisboeck TS: Simulating non-small cell lung cancer with a multiscale agent-based model. Theor Biol Med Model 2007, 4: 50. 10.1186/1742-4682-4-50

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  13. Zhang L, Wang ZH, Sagotsky JA, Deisboeck TS: Multiscale agent-based cancer modeling. J Math Biol 2009, 58(4–5):545–559.

    Article  PubMed  Google Scholar 

  14. Hanahan D, Weinberg RA: The hallmarks of cancer. Cell 2000, 100(1):57–70. 10.1016/S0092-8674(00)81683-9

    Article  CAS  PubMed  Google Scholar 

  15. Graziano L, Preziosi L, Mollica F, Preziosi L, Rajagopal KR: Modeling of Biological Materials. In Mechanics in tumor growth. Edited by: Mollica F, Preziosi L, Rajagopal KR. Boston: Birkhäuser; 2007:263–321.

    Google Scholar 

  16. Risau W: Mechanisms of angiogenesis. Nature 1997, 386(6626):671–674.

    Article  CAS  PubMed  Google Scholar 

  17. Zhang L, Strouthos CG, Wang Z, Deisboeck TS: Simulating brain tumor heterogeneity with a multiscale agent-based model: Linking molecular signatures, phenotypes and expansion rate. Math Comput Model 2009, 49(1–2):307–319.

    Article  PubMed Central  PubMed  Google Scholar 

  18. Ellis P, Morzycki W, Melosky B, Butts C, Hirsh V, Krasnoshtein F, Murray N, Shepherd FA, Soulieres D, Tsao MS: The role of the epidermal growth factor receptor tyrosine kinase inhibitors as therapy for advanced, metastatic, and recurrent non-small-cell lung cancer: a Canadian national consensus statement. Curr Oncol 2009, 16(1):27.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  19. Anderson ARA, Chaplain MAJ: Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull Math Biol 1998, 60(5):857–899. 10.1006/bulm.1998.0042

    Article  CAS  PubMed  Google Scholar 

  20. Araujo RP, Petricoin EF, Liotta LA: A mathematical model of combination therapy using the EGFR signaling network. Biosystems 2005, 80(1):57–69. 10.1016/j.biosystems.2004.10.002

    Article  CAS  PubMed  Google Scholar 

  21. Alarcon T: A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells. J Theor Biol 2004, 229(3):395–411. 10.1016/j.jtbi.2004.04.016

    Article  CAS  PubMed  Google Scholar 

  22. Mansury YD: T: The impact of “search precision” in an agent-based tumor model. J Theor Biol 2003, 224: 25–337.

    Article  Google Scholar 

  23. Shvartsman SY, Wiley HS, Deen WM, Lauffenburger DA: Spatial range of autocrine signaling: modeling and computational analysis. Biophys J 2001, 81(4):1854–1867. 10.1016/S0006-3495(01)75837-7

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Morton KW, Mayers DF: Numerical solution of partial differential equations: an introduction. Cambridge: Univ Pr; 2005.

    Book  Google Scholar 

  25. Stokes CL, Lauffenburger DA: Analysis of the roles of microvessel endothelial cell random motility and chemotaxis in angiogenesis. J Theor Biol 1991, 152(3):377–403. 10.1016/S0022-5193(05)80201-2

    Article  CAS  PubMed  Google Scholar 

  26. Zhu XW, Zhou XB, Lewis MT, Xia L, Wong S: Cancer stem cell, niche and EGFR decide tumor development and treatment response: A bio-computational simulation study. J Theor Biol 2011, 269(1):138–149. 10.1016/j.jtbi.2010.10.016

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  27. Griffero F, Daga A, Marubbi D, Capra MC, Melotti A, Pattarozzi A, Gatti M, Bajetto A, Porcile C, Barbieri F: Different response of human glioma tumor-initiating cells to epidermal growth factor receptor kinase inhibitors. J Biol Chem 2009, 284(11):7138.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  28. Liu G, Swihart MT, Neelamegham S: Sensitivity, principal component and flux analysis applied to signal transduction: the case of epidermal growth factor mediated signaling. Bioinformatics 2005, 21(7):1194–1202. 10.1093/bioinformatics/bti118

    Article  CAS  PubMed  Google Scholar 

  29. Peng H, Wen J, Li H, Chang J, Zhou X: Drug inhibition profile prediction for NFkB pathway in multiple myeloma. PLoS One 2011, 6(3):e14750. 10.1371/journal.pone.0014750

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  30. Sun X, Su J, Bao J, Peng T, Zhang L, Zhang Y, Yang Y, Zhou X: Cytokine combination therapy prediction for bone remodeling in tissue engineering based on the intracellular signaling pathway. Biomaterials 2012, 33(33):8265–8276. 10.1016/j.biomaterials.2012.07.041

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  31. Rabitz H, Kramer M, Dacol D: Sensitivity analysis in chemical kinetics. Annu Rev Phys Chem 1983, 34(1):419–461. 10.1146/annurev.pc.34.100183.002223

    Article  CAS  Google Scholar 

  32. Petersen WP, Arbenz P: Introduction to parallel computing, vol. 9. USA: Oxford University Press; 2004.

    Google Scholar 

  33. Anderson ARA, Weaver AM, Cummings PT, Quaranta V: Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. Cell 2006, 127(5):905–915. 10.1016/j.cell.2006.09.042

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

This work was supported by Funding: NIH R01LM010185-03 (Zhou), NIH U01HL111560-01 (Zhou), NIH 1R01DE022676-01 (Zhou) and DoD TATRC (Zhou).

We would like to thank the members of Translational Biosystems Lab of Cornell Medical School for the valuable discussions.

Author information

Authors and Affiliations

Authors

Corresponding authors

Correspondence to Le Zhang or Xiaobo Zhou.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

XS participated to study conception, carried out the model programming, carried out the analysis of the model and drafted the manuscript. XZ participated to study conception and improved the manuscript. LZ participated to study conception, model analysis and improved the manuscript. HT participated to study conception and helped to initiate the model programming. JB helped to study conception. CS helped to improve the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material

12859_2012_5416_MOESM1_ESM.doc

Additional file 1:Table A1. Kinetic equations describing the reactions between the components of the simplified EGFR signaling pathway. (DOC 104 KB)

Additional file 2:Table A2. Coefficients of the simplified EGFR signaling pathway. (DOC 72 KB)

12859_2012_5416_MOESM3_ESM.doc

Additional file 3:Table A3. Kinetic equations describing the reactions between the components of the cell-cycle. (DOC 54 KB)

Additional file 4:Table A4. Parameter in cell-cycle pathway. (DOC 64 KB)

Additional file 5:Table A5. Parameters of microenvironmental PDEs in the model. (DOC 245 KB)

12859_2012_5416_MOESM6_ESM.doc

Additional file 6:Text A1. Equations describing the Initial distribution of glucose, oxygen, TGFα, VEGF and fibronectin. (DOC 178 KB)

Additional file 7:Figure A1. Tumor induced angiogenesis and vascular tumor growth without TKI treatment. (DOC 306 KB)

Additional file 8:Figure A2. Vascular tumor growths with concentration change of fibronectin. (DOC 517 KB)

12859_2012_5416_MOESM9_ESM.doc

Additional file 9:Figure A3. The concentration change of glucose, oxygen, TGFα and VEGF without TKI treatment. (DOC 722 KB)

Additional file 10:Figure A4. Various tumor cell numbers without TKI treatment. (DOC 82 KB)

Additional file 11:Figure A5. Vascular tumor growth with TKI treatment. (DOC 293 KB)

Additional file 12:Figure A6. Vascular tumor growth in the presence of fibronectin with TKIs treatment. (DOC 359 KB)

12859_2012_5416_MOESM13_ESM.doc

Additional file 13:Figure A7. The concentration change of glucose, oxygen, TGFα and VEGF with TKI treatment. (DOC 460 KB)

Additional file 14:Figure A8. The concentration change of TKIs shown at 60, 150, 240 and 300 hours. (DOC 398 KB)

Additional file 15:Figure A9. Various tumor cell numbers with TKI treatment. (DOC 218 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Sun, X., Zhang, L., Tan, H. et al. Multi-scale agent-based brain cancer modeling and prediction of TKI treatment response: Incorporating EGFR signaling pathway and angiogenesis. BMC Bioinformatics 13, 218 (2012). https://doi.org/10.1186/1471-2105-13-218

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1471-2105-13-218

Keywords