Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

This article is part of the supplement: Selected articles from the 6th International Conference on Computational Systems Biology (ISB2012)

Open Access Research

A 3D multiscale model of cancer stem cell in tumor development

Fuhai Li1, Hua Tan12, Jaykrishna Singh1, Jian Yang1, Xiaofeng Xia1, Jiguang Bao2, Jinwen Ma13, Ming Zhan1* and Stephen TC Wong1

Author Affiliations

1 NCI Center for Modeling Cancer Development, Department of Systems Medicine and Bioengineering, The Methodist Hospital Research Institute, Weill Medical College of Cornell University, Houston, TX 77030, USA

2 School of Mathematical Sciences, Beijing Normal University, Laboratory of Mathematics and Complex Systems, Ministry of Education, Beijing 100875, China

3 Department of Information Science, School of Mathematical Sciences & LMAM, Peking University, Beijing 100871, China

For all author emails, please log on.

BMC Systems Biology 2013, 7(Suppl 2):S12  doi:10.1186/1752-0509-7-S2-S12


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/7/S2/S12


Published:17 December 2013

© 2013 Li et al.; licensee BioMed Central Ltd.

This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Recent reports indicate that a subgroup of tumor cells named cancer stem cells (CSCs) or tumor initiating cells (TICs) are responsible for tumor initiation, growth and drug resistance. This subgroup of tumor cells has self-renewal capacity and could differentiate into heterogeneous tumor cell populations through asymmetric proliferation. The idea of CSC provides informative insights into tumor initiation, metastasis and treatment. However, the underlying mechanisms of CSCs regulating tumor behaviors are unclear due to the complex cancer system. To study the functions of CSCs in the complex tumor system, a few mathematical modeling studies have been proposed. Whereas, the effect of microenvironment (mE) factors, the behaviors of CSCs, progenitor tumor cells (PCs) and differentiated tumor cells (TCs), and the impact of CSC fraction and signaling heterogeneity, are not adequately explored yet.

Methods

In this study, a novel 3D multi-scale mathematical modeling is proposed to investigate the behaviors of CSCsin tumor progressions. The model integrates CSCs, PCs, and TCs together with a few essential mE factors. With this model, we simulated and investigated the tumor development and drug response under different CSC content and heterogeneity.

Results

The simulation results shown that the fraction of CSCs plays a critical role in driving the tumor progression and drug resistance. It is also showed that the pure chemo-drug treatment was not a successful treatment, as it resulted in a significant increase of the CSC fraction. It further shown that the self-renew heterogeneity of the initial CSC population is a cause of the heterogeneity of the derived tumors in terms of the CSC fraction and response to drug treatments.

Conclusions

The proposed 3D multi-scale model provides a new tool for investigating the behaviors of CSC in CSC-initiated tumors, which enables scientists to investigate and generate testable hypotheses about CSCs in tumor development and drug response under different microenvironments and drug perturbations.

Background

The mechanisms of tumor initiation, progression, metastasis and drug resistance remain elusive due to the complex system of tumors. Recent studies have shown that a sub-population of tumor cells, named cancer stem cells (CSCs) or tumor initiating cells (TICs) tumor are responsible for tumor development and drug resistance [1-3]. The CSC concept is still controversial, as it is difficult to discover and validate cancer stem cells, particularly their unlimited self-renewal and differentiation capabilities [1]. However, CSCs have been being isolated from more and more cancers [2], since first discovered in the acute myeloid leukemia (AML) by using CD34++/CD38- biomarkers [3]. Recently, the breast CSCS were identified by using CD44+CD24-/low biomarkers in [4], and the colon CSCS were also reported [5]. CSCs are believed to have strong self-renewal capacity, could differentiate into heterogeneous tumor populations through asymmetric division [6,7], and are responsible for drug resistance and metastasis [8,9]. Reportedly, CSCs are heterogeneous with different self-renewal and tumor formation abilities, which might be caused by varying activation intracellular signaling (e.g., Wnt, Shh and Stat3) due to the diverse concentrations of external mE factors [10,11].

However the roles of CSCs in tumor development remain unknown because of tumor complexity in multiple levels, including signaling transduction, cell-cell communication, and cell-microenvironment interactions. The mathematical simulation models have been powerful tools for understanding the tumor systems [12]. In general, the existing mathematical models of the tumor development can be grouped into three major categories: discrete, continuous and hybrid. The discrete models, e.g., cellular automata [13] and Glazier and Graner model [14], simulate cell behaviors individually with a group of rules. The continuous models employ ordinary or partial differential equations to simulate the behaviors of tumor cell populations and dynamics of mE factors [15,16]. The hybrid models are the combination of the discrete (for modeling cells) and continuous (for modeling mE factors) models [12]. A few mathematical simulation studies have been developed to study CSCs functions in tumor development [17-20]. Whereas, the interactions between mE factors and CSCs, PCs and TCs, the impact of the heterogeneity and fraction of CSCs, have not been adequately considered in the mathematical modeling.

Herein, a 3D multi-scale mathematical modeling is proposed to study the functions of CSCs in tumor development. The model, extended our previous concept [21], enables us to study the roles of CSCs in tumor progression and chemo-drug resistance by simulating the tumor growth initiated by a set of heterogeneous CSC populations. The overview of the proposed multi-scale model is shown in Figure 1. In brief, mE factors are considered, including nutrients (e.g., oxygen and glucose), tumor angiogenesis factors (TAF), matrix degrading proteolytic enzyme (MDE), extracellular matrix (ECM), tissue pressure [22], and motility of cells are described in the molecular scale. The interactions between mE factors and cancer cells are shown in Figure 2. The behaviors of endothelial and cancer cells are described by the cellular automata in the cellular scale. The cancer cells are represented by a hierarchical organization of cellular subtypes [23,24], including CSC, a series of intermediate PCs and TCs, each endowed with different biological traits. At the tissue level, the tumor volume, shape and spatial distributions of tumor cells are investigated.

thumbnailFigure 1. The schematic of the 3D multiscale model and the implementation. The molecular scale describing the diffusion and reaction processes of mE factors with PDEs. The cellular level models the behaviours of cells with a 3D cellular automata. The tissue level visualizes the whole tumor morphology.

thumbnailFigure 2. The interactions between micro-environmental (mE) factors and tumor cells. The outer circle stands for mE factors being considered at the molecular level, the inner circle corresponds to the processes investigated at the cellular level. The interrelations are denoted by the solid arrows.

Methods

The PDE system characterizing reaction-diffusion process of mE factors

Five mE factors are considered in this model, including nutrients (n), tumor angiogenic factor (TAF) (c), matrix degrading proteolytic enzyme (MDE) (m), extracellular matrix (ECM) (f), and tissue pressure (p). A system of PDEs is used to delineate the diffusion and reactions mE factors. The χΩ is defined as follows.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M1">View MathML</a>

The diffusion and reaction of nutrients are modeled by the quasi-steady equation with non-zero Dirichlet boundary conditions [25,26]. Since they are much smaller comparing to cells, nutrient molecules diffuse quickly through ECM at each time point.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M2">View MathML</a>

(1)

where Dn is the diffusion rate of nutrient molecules, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M3">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M4">View MathML</a> denote the nutrient molecules transferring rates from pre-existing and neo-vasculature vessels. The <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M5">View MathML</a> is the nutrient molecule binding rate to fibronectin; <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M6">View MathML</a> is the uptake rate by cells, and it is different for specific types of cells. The χΣc is an indicator function that equals to 1 at the new generated vessels. The term (1-p) is used to indicate the difference of nutrient molecule transfer with different pressure, and (1-n) is to reflect nutrient molecules' saturation effect.

The TAF is selected by tumor cells, and during its diffusion, TAF will be either degrade naturally or captured by endothelial cells. This process is described as follows:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M7">View MathML</a>

(2)

where Dc is the diffusion rate of TAF, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M8">View MathML</a> is the boundary of necrotic and viable regions, the <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M9">View MathML</a> is the unit outer normal direction, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M10">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M11">View MathML</a> are the TAF secretion rates by dying and viable tumor cells, respectively; <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M12">View MathML</a> is the uptake rate by endothelial cells, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M13">View MathML</a> is the rate of natural degradation.

The ECM, such as fibronectin, represents a set of binding molecules that do not diffuse but increase the tumor cell adhesion. The concentration of ECM molecules is measured as follows:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M14">View MathML</a>

(3)

Where the <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M15">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M16">View MathML</a> denote fibronectin production rates by tumor and endothelial cells, respectively, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M17">View MathML</a> is the ECM degradation rate by MDE.

MDE is secreted by endothelial cells and viable tumor cells to degrade ECM. The MDE diffusion and reaction processes are defined as:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M18">View MathML</a>

(4)

In this equation, Dm is the diffusion coefficient, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M19">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M20">View MathML</a> denote MDE secretion rates by the endothelial cells and viable tumor cells, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M21">View MathML</a> is the MDE decay rate.

The cell velocity is resulted from the different cell proliferation at different regions, and is described by using the following Darcy-Stokes law [27].

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M22">View MathML</a>

(5)

And the velocity field follows the divergence equation:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M23">View MathML</a>

(6)

where, λa and λN are the volume loss rates caused by cellular apoptosis. The first term on the right of equation (6) is the source effect, while the second term can be considered as the sink effect. The diffusion of pressure is obtained by taking the divergence operation on both sides of equation (5), and combined with equation (6):

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M24">View MathML</a>

(7)

In the implementation, the finite element method is used to solve PDEs with diffusion item <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M25">View MathML</a>[28]. The ECM equation is solved by the 2nd order total variation Runge-Kutta method [29]. The time interval, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M26">View MathML</a>, is calculated to keep the stability of the PDEs:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M27">View MathML</a>

(8)

where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M28">View MathML</a> is the spatial interval, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M29">View MathML</a> is the function defined on the TAF, ECM, and cell velocity in [25], and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M30">View MathML</a> is the cell velocity at i-th spatial point [25].

3D Cellular Automata

For simplicity, cells are limited to interact with its six immediate orthogonal neighbor grids such that the cell can move or proliferate, based on probabilities calculated from the distribution of surrounding mE factors. Figure 3 shows a flowchart of cell behavior under interactions with mE factors and cell lineage in proliferation. Specifically, for a cell being at location k:

thumbnailFigure 3. The flowchart of control in the multiscale model of CSC-initiated tumor development. (A) Each cell evolves according to their life cycle under the conditions confined by the PDE system. Necrosis occurs when cell death is induced by hypoxia. (B) Hierarchical organization of different cell subtypes and their proliferative kinetics. Cancer stem cells (CSC) expand their own population by symmetric proliferation to two identical daughter cells still with CSC-like traits and expand the whole tumor through asymmetric differentiation to progenitor cells (PC) and terminally differentiated cells (TC) in this model. Similarly, PCs contribute to the constitution of the tumor in a similar way, but do not reversibly produce CSCs according to the hierarchical organization hypothesis. The TCs are assumed to either proliferate or be apoptotic at each time point without ability to divide into any other subtypes. The parameters above the arrows indicate the occurrence probability of the referred event at each time point.

1) Check available (empty) neighbor locations.

2) If there is no available neighbor location, go to step 8).

3) Calculate the motion probabilities to the m (m <= 6) available neighbor locations as: qi = ni/fi, i = 0, 1, ..., m, where q0 means the probability of staying at the same location.

4) Denote <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M31">View MathML</a>, and normalize them as: <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M32">View MathML</a>.

6) Define <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M33">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M34">View MathML</a>, i = 1, 2, ..., m.

7) Pick up a number r randomly from [0, 1], then move the cell to the i-the neighbor location where r belongs to.

8) Update the cell age, and then check if the cell is mature to divide. If yes, add one new cell with right cell type to a neighbor location with the above rules. Check the mE conditions to determine if the cell should enter the quiescent or death status.

Simulation of chemotherapy

In simulation of chemotherapy, all cancer cells (CSCs, PCs, TCs) can be killed by chemo drugs with different doses. Since the diffusion and reaction processes of the drug molecules resemble that of nutrients, we use the same diffusion-reaction equation for drugs. We assume that TCs proliferation is reduced due to the effects of drugs, while CSC proliferation is accelerated as activated by the volume loss due to the drug effect, which is parallel to the normal stem cell functions [30]. We adjust the proliferation age of each cell subtype in such a way to stimulate the reaction of tumor cells to chemotherapy.

The Gompertz curve fitting

The tumor growth pattern is fitted by using the Gompertz curve, which shows slow change in the beginning and the end of tumor growth [31]. Though for tumors that undergo angiogenesis, a plateau in growth is not necessarily reached, the Gompertz cure is an appropriate choice in modeling the initial phase of the tumor growth with limited access to nutrients. Mathematically the tumor growth is represented as:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M35">View MathML</a>

(9)

where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M36">View MathML</a> is the volume of tumor at time <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M37">View MathML</a>, while positive constants k, b denote the axis displacement and growth speed. The least square is used to determine the optimal parameters as follows:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M38">View MathML</a>

(10)

Here, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M39">View MathML</a> is the tumor volume growth function defined in equation (9), <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M40">View MathML</a> is the measured volume of tumor at time point i, obtained either from biological experiments or by computer simulation, and N is the number of available measurements.

Measurements of tumor development

The tumor properties are evaluated by following measurements: proliferation potential (PP), time to reach potential (TtP), average aggressive index (AAI), and average fitting error (AFE). PP is defined as (9):

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M41">View MathML</a>

(11)

The PP value could not predict the final volume of a tumor because of many unforeseen contributing factors when tumors grow large. However, it can be used to compare the potential volume of tumors in a relative sense. TtP is estimated by solving an inverse problem of (9), that is, by searching the time when <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M39">View MathML</a> reaches PP for the first time:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M42">View MathML</a>

(12)

where 'inf' is the infimum operation. AAI is represented as:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M43">View MathML</a>

(13)

Here, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M44">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M45">View MathML</a> stand for the surface area and volume of the tumor at time point i, which are represented by the number of cells on the surface and the total number of cells composing the tumor respectively. The performance of the curve fitting is assessed by AFE, estimated by the least square method for (10):

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M46">View MathML</a>

(14)

where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M47">View MathML</a>.

Results

Simulation of tumor development under different CSCs contents

To investigate tumor development initiated from a set of tumor cells with different CSC contents (or fractions), herein, we simulated tumor growth with pure CSCs (CSCs only, initiated from ~20 CSCs) and unsorted (mixtures of CSCs and non-CSCs) tumor cells (200 tumor cells in which 4% are CSCs), respectively. Figure 4 shows the dynamics of tumor growth and concentration profiles of mE factors over the time under different CSC contents. As shown, the tumor initiated from the mixed tumor cells grows faster than that tumors initiated from pure CSCs. Whereas the tumors initiated from pure CSCs have more spiky morphology. This might be caused by biased migration of CSCs toward locations (without the limitation of other tumor cells) with a higher nutrient concentration and lower ECM concentration, which affect the tumor geometric morphology. We conducted growth curve fitting of the simulation results and found that the tumor initiated from pure CSCs has an elevated proliferation potential, though it takes a longer time to grow to its limit size (Figure 5). Tumors also exhibits a stronger aggressiveness compared to tumors initiated from the mixed tumor cells (Figure 5).

thumbnailFigure 4. Simulation of tumor development under different CSC contents. (A) Time evolution of tumors initiated by pure CSCs and unsorted tumor cells. (B) The corresponding concentration profiles of micro-environmental (mE) factors.

thumbnailFigure 5. Quantitative comparison between tumor development initiated from pure CSCs and unsorted tumor cells. The proliferation potential, time needed to reach potential, average aggressive index, and average fitting error of tumors are compared between two kinds of tumors in five simulations.

Important parameters to tumor growth

In Table 1, all parameters of the proposed multi-scale model of tumor growth are listed. These parameters are determined mostly from either literature or experimental data. We performed sensitivity analysis to discover important parameters to tumor growth initiated from pure CSCs. The parameter values were perturbed in a range of 10%. Figure 6 shows the effects of parameters on the tumor growth. Some diffusion reaction related parameters, e.g., <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M3">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M10">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M19">View MathML</a>, and Dc are sensitive to the tumor growth, and some parameters proliferation abilities parameters, i.e., KCCP and KPP are also sensitive to tumor growth.

Table 1. Model parameters.

thumbnailFigure 6. Sensitivity analysis (A) The sensitivity analysis of the continuous parameters. (B) The sensitivity analysis of the proliferation related parameters.

Tumor response to drug treatment

The effect of chemo-drug treatment on tumor development and CSC fraction is also investigated. Figure 7 shows the dynamics of tumor development under chemo-drug treatment. As can be seen, the solid tumor shrinks during the drug treatment, whereas the CSC fraction increases. The simulation results also reveal that tumors will grow fast to half of its original volume after stop the drug treatment. This fast relapse may be because some CSCs reside in the interior of a tumor where drug molecules are not accessible. It might be due to that the fast proliferating non-CSCs dominate the outer rim of a tumor, and drug molecules will first kill the non-CSCs in order to reach the area where CSCs tend to gather. The quick relapse of tumor, once the treatment stops, as observed in our simulations, can be explained by the escape of CSCs from therapeutic interventions. The decreased percentage and lower proliferation rate of CSCs may be because these interior CSCs reside in a region where the access to nutrients is limited.

thumbnailFigure 7. Simulation of chemo-drug therapy response. (A) A 9-week treatment to a tumor initiated from pure CSCs. The tumor volume first reduces fast, then remains stable during drug treatment; whereas it grows back quickly after the treatment stops. (B) The evolution of the CSC fraction in the process. (C) The 3D morphology of the tumor at different stages of the treatment.

Simulation of CSC self-renewal heterogeneity in tumor development

To investigate the heterogeneity of CSC self-renewal ability, we further extend our model as follows: 1) a tumor is initiated from a single CSC; 2) the initiating CSCs have different self-renewal abilities that are proportional to the Gaussian distribution with different means: <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M48">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M49">View MathML</a> is a Gaussian with mean μ and variance <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M51">View MathML</a>; and 3) the newly formed CSCs during tumor development will obtain a mean value randomly sampled from the <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M52">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M53">View MathML</a> is the mean of its parent CSC. In our simulation, we set <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M50">View MathML</a>, 0.5, and 0.75, respectively, to indicate the relative low, medium, and high self-renewal ability of initiating CSCs (with fixed <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/S2/S12/mathml/M54">View MathML</a>). The simulation starts from a single CSC to 30 days, then the same chemo-drug treatment is applied. The virtual chemo-drug treatment will be stopped when 85% of tumor cells are killed. Figure 8 shows the simulation results for the size and fraction of CSCs of tumors initiated from CSCs with different self-renewal abilities. As shown, before chemo-drug treatment, the tumor initiated from CSCs with high self-renewal ability has an average bigger tumor volume (about 1.5 fold), whereas the difference of the CSC fraction is about 2 fold. After chemo-drug treatment, the difference in the CSC fraction is significantly increased to about 3 fold, though the CSC fractions of all tumors are increased significantly comparing to those before chemo-drug treatment. The results indicated that though the pure chemo-drug treatment could reduce the size of tumor, it might also increase the aggressiveness of tumor due to the increased fraction of CSCs, especially for the tumor with the CSCs that have high self-renewal ability. Therefore, the combinations of chemo-drugs and anti-CSCs drugs are needed to achieve better treatment outcomes.

thumbnailFigure 8. Comparison of tumor development initiated from single CSCs with heterogeneous self-renewal ability. (A) Average size of tumors and (B) average CSC fractions of tumors initiated from a single CSC with relative low (0.25), medium (0.5) and high (0.75)self-renewal ability. (C) Average CSC fractions of tumors before and after chemo-drug treatment.

Discussion & conclusions

Here, a multi-scale and multi-factorial computational model is established in 3D space to study the behaviors and roles of CSCs in leading tumor development. The model is implemented at three hierarchical scales (molecular, cellular, tissue scales). The molecular subsystem characterizes the diffusion and reaction processes of mE factors by using PDEs. The cellular level subsystem simulates the proliferation and migration of all cancer cells and endothelial cells, considering the availability of mE factors, with a 3D cellular automaton. The tissue level subsystem evaluates the temporal and spatial variations of tumor morphology by, four indices. The model can be conveniently expanded to a particular application to generate testable hypotheses.

The simulation studies based on the multi-scale model could provide important insights into tumor development and treatment. For example, the simulation indicated that tumors in mice model initiated by the sorted CSC population had stronger aggressiveness and proliferation potential comparing to tumors from unsorted cancer cells. Also the simulation demonstrated that the neo-vasculature could grow into the interior of a tumor, suggesting a possibility of delivering drugs via neo-vasculature to target the CSCs in the interior of tumors, besides the anti-angiogenic therapy that elicits increased local invasion and distant metastasis of tumors [32,33]. In addition, the simulation indicated that pure chemo-drug treatment may increase the fraction of CSCs significantly, especially for the tumor with CSCs of high self-renewal ability, and consequently, the tumor residual will be more chemo-resistant and aggressive. Thus, a combination of chemo-drugs with CSC inhibition drugs would be more effective in cancer treatment without increasing tumor drug-resistance and aggressiveness.

Many parameters in our model were defined based on general understanding of tumor development through literature mining in this study, as in many other modeling studies. Despite limited experimental data used in defining the parameters, the model we proposed is still valid. It enables us to identify important CSCs behavior and interactions with interested mE factors, and to virtually test hypotheses that cannot be done in an animal model. The simulation studies based on the model can lead to new insight of CSCs in tumor development and shed light on the treatment. As more experimental data become available through our studies, the parameters can be better defined and calibrated, and the resulting model will be better predictable.

Several improvements of the proposed 3D tumor growth model will be conducted in the future work. Cancer is a complex and heterogeneous disease. CSCs from different types of cancer might have different functions and regulatory signaling pathways. With more data of cancer-specific signaling pathways, cell differentiation lineages, stromal cells, and detailed cell-cell interactions becoming available, the proposed model could be extended to study the CSCs in specific cancer types. Also the cell shape could be taken into account as it plays an important role in interactions between cells and mE factors, particularly when the cell density is high and causes shape deformation and cell-cell interaction through cell surface markers. On the other hand, relevant regulatory or signaling pathways could be integrated to refine the modeling. In addition, specific drug effects on different cell cycles and cell types could be considered.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

FL, MZ, SW conceived the study. FL, TH, MZ designed and coordinated the study and drafted the manuscript. FL, TH, JS, JY, XX, JB, JM, MZ participated in algorithm development and data analysis. All authors read and approved the final manuscript.

Acknowledgements

We would like to thank colleagues of the NCI-ICBP Center for Modeling Cancer Development (CMCD) at The Methodist Hospital Research Institute and Baylor College of Medicine for helpful discussions. Preliminary results of this study were published in the proceedings of IEEE ISB2012.

Declarations

The publication of this article has been funded by NIH U54CA149169.

This article has been published as part of BMC Systems Biology Volume 7 Supplement 2, 2013: Selected articles from The 6th International Conference of Computational Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/7/S2.

References

  1. Gupta PB, Chaffer CL, Weinberg RA: Cancer stem cells: mirage or reality?

    Nat Med 2009, 15(9):1010-1012. PubMed Abstract | Publisher Full Text OpenURL

  2. Nguyen LV, Vanner R, Dirks P, Eaves CJ: Cancer stem cells: an evolving concept.

    Nat Rev Cancer 2012, 12(2):133-143. PubMed Abstract | Publisher Full Text OpenURL

  3. Bonnet D, Dick JE: Human acute myeloid leukemia is organized as a hierarchy that originates from a primitive hematopoietic cell.

    Nat Med 1997, 3(7):730-737. PubMed Abstract | Publisher Full Text OpenURL

  4. Al-Hajj M, Wicha MS, Benito-Hernandez A, Morrison SJ, Clarke MF: Prospective identification of tumorigenic breast cancer cells.

    Proceedings of the National Academy of Sciences 2003, 100(7):3983-3988. Publisher Full Text OpenURL

  5. O'Brien CA, Pollett A, Gallinger S, Dick JE: A human colon cancer cell capable of initiating tumour growth in immunodeficient mice.

    Nature 2007, 445(7123):106-110. PubMed Abstract | Publisher Full Text OpenURL

  6. Pardal R, Clarke MF, Morrison SJ: Applying the principles of stem-cell biology to cancer.

    Nat Rev Cancer 2003, 3(12):895-902. PubMed Abstract | Publisher Full Text OpenURL

  7. Huntly BJ, Gilliland DG: Cancer biology: summing up cancer stem cells.

    Nature 2005, 435(7046):1169-1170. PubMed Abstract | Publisher Full Text OpenURL

  8. Rich JN: Cancer stem cells in radiation resistance.

    Cancer Res 2007, 67(19):8980-8984. PubMed Abstract | Publisher Full Text OpenURL

  9. Zhang M, Atkinson RL, Rosen JM: Selective targeting of radiation-resistant tumor-initiating cells.

    Proc Natl Acad Sci USA 2010, 107(8):3522-3527. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  10. Horst D, Chen J, Morikawa T, Ogino S, Kirchner T, Shivdasani RA: Differential WNT activity in colorectal cancer confers limited tumorigenic potential and is regulated by MAPK signaling.

    Cancer Res 2012, 72(6):1547-1556. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  11. Vermeulen L, De Sousa EMF, van der Heijden M, Cameron K, de Jong JH, Borovski T, Tuynman JB, Todaro M, Merz C, Rodermond H, et al.: Wnt activity defines colon cancer stem cells and is regulated by the microenvironment.

    Nat Cell Biol 2010, 12(5):468-476. PubMed Abstract | Publisher Full Text OpenURL

  12. Anderson AR: A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion.

    Math Med Biol 2005, 22(2):163-186. PubMed Abstract | Publisher Full Text OpenURL

  13. Kansal AR, Torquato S, Harsh GI, Chiocca EA, Deisboeck TS: Simulated brain tumor growth dynamics using a three-dimensional cellular automaton.

    J Theor Biol 2000, 203(4):367-382. PubMed Abstract | Publisher Full Text OpenURL

  14. Bauer AL, Jackson TL, Jiang Y: A cell-based model exhibiting branching and anastomosis during tumor-induced angiogenesis.

    Biophys J 2007, 92(9):3105-3121. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  15. Frieboes HB, Zheng X, Sun CH, Tromberg B, Gatenby R, Cristini V: An integrated computational/experimental model of tumor invasion.

    Cancer Res 2006, 66(3):1597-1604. PubMed Abstract | Publisher Full Text OpenURL

  16. Frieboes HB, Lowengrub JS, Wise S, Zheng X, Macklin P, Bearer EL, Cristini V: Computer simulation of glioma growth and morphology.

    Neuroimage 2007, 37(Suppl 1):S59-70. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  17. Zhu X, Zhou X, 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. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. Ganguly R, Puri IK: Mathematical model for the cancer stem cell hypothesis.

    Cell Prolif 2006, 39(1):3-14. PubMed Abstract | Publisher Full Text OpenURL

  19. Michor F: Mathematical models of cancer stem cells.

    J Clin Oncol 2008, 26(17):2854-2861. PubMed Abstract | Publisher Full Text OpenURL

  20. Sottoriva A, Verhoeff JJC, Borovski T, McWeeney SK, Naumov L, Medema JP, Sloot PMA, Vermeulen L: Cancer Stem Cell Tumor Model Reveals Invasive Morphology and Increased Phenotypical Heterogeneity.

    Cancer Research 2010, 70(1):46-56. PubMed Abstract | Publisher Full Text OpenURL

  21. Tan H, Li F, Singh J, Xia X, Cridebring D, Yang J, Zhan M, Wong S, Bao J, Ma J: A 3-dimentional multiscale model to simulate tumor progression in response to interactions between cancer stem cells and tumor microenvironmental factors.

    2012 IEEE 6th International Conference on Systems Biology (ISB) 18-20 Aug. 2012 2012 2012, 297-303. OpenURL

  22. Hori K, Suzuki M, Abe I, Saito S: Increased tumor tissue pressure in association with the growth of rat tumors.

    Jpn J Cancer Res 1986, 77(1):65-73. PubMed Abstract OpenURL

  23. Morrison SJ, Kimble J: Asymmetric and symmetric stem-cell divisions in development and cancer.

    Nature 2006, 441(7097):1068-1074. PubMed Abstract | Publisher Full Text OpenURL

  24. Dingli D, Traulsen A, Michor F: (A)symmetric stem cell replication and cancer.

    PLoS Comput Biol 2007, 3(3):e53. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  25. 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. PubMed Abstract | Publisher Full Text OpenURL

  26. Macklin P, McDougall S, Anderson AR, Chaplain MA, Cristini V, Lowengrub J: Multiscale modelling and nonlinear simulation of vascular tumour growth.

    J Math Biol 2009, 58(4-5):765-798. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  27. Truskey G, Yuan F, Katz D: Transport Phenomena in Biological Systems.

    Pearson Prentice Hall, Upper Saddle River, NJ 2004. OpenURL

  28. Gockenbach MS: Understanding and implementing the finite element method.

    Philadelphia: Society for Industrial and Applied Mathematics 2006. OpenURL

  29. Shu SGaC: Total variation diminishing runge-kutta schemes.

    Mathematics of computation 1998, 67:73-85. Publisher Full Text OpenURL

  30. Weissman IL: Stem cells: units of development, units of regeneration, and units in evolution.

    Cell 2000, 100(1):157-168. PubMed Abstract | Publisher Full Text OpenURL

  31. Norton L: A Gompertzian model of human breast cancer growth.

    Cancer Res 1988, 48(24 Pt 1):7067-7071. PubMed Abstract | Publisher Full Text OpenURL

  32. Paez-Ribes M, Allen E, Hudock J, Takeda T, Okuyama H, Vinals F, Inoue M, Bergers G, Hanahan D, Casanovas O: Antiangiogenic therapy elicits malignant progression of tumors to increased local invasion and distant metastasis.

    Cancer Cell 2009, 15(3):220-231. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  33. Ebos JM, Lee CR, Cruz-Munoz W, Bjarnason GA, Christensen JG, Kerbel RS: Accelerated metastasis after short-term treatment with a potent inhibitor of tumor angiogenesis.

    Cancer Cell 2009, 15(3):232-239. PubMed Abstract | Publisher Full Text OpenURL