A Hybrid Cellular Automata Model of Multicellular Tumour Spheroid Growth in Hypoxic Microenvironment

Abstract

A three-dimensional hybrid cellular automata (CA) model is developed to study the dynamic process of multicellular tumour spheroid (MTS) growth by introducing hypoxia as an important microenvironment factor which influences cell migration and cell phenotype expression. The model enables us to examine the effects of different hypoxic environments on the growth history of the MTS and to study the dynamic interactions between MTS growth and chemical environments. The results include the spatial distribution of different phenotypes of tumour cells and associated oxygen concentration distributions under hypoxic conditions. The discussion of the model system responses to the varied hypoxic conditions reveals that the improvement of the resistance of tumour cells to a hypoxic environment may be important in the tumour normalization therapy.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Yan Cai 1 and Jie Wu 2 and Shixiong Xu 3 and Zhiyong Li 1

Recommended by Martin Weiser

1, School of Biological Science and Medical Engineering, Southeast University, Nanjing 210096, China
2, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai 200240, China
3, Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China

Received 19 November 2012; Accepted 22 January 2013

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Multicellular tumour spheroids (MTS) are three-dimensional aggregates of malignant tumour cells, which can be grown in vitro under strictly controlled nutritional and mechanical conditions to mimic microtumour growth and metastases [ 1]. Sutherland et al. [ 2, 3] were the first to use MTS for the systematic study of tumour responses to therapy. Compared to monolayer cultures, a significant advantage of the MTS model is that they can more closely represent the in vivo microenvironment in solid tumours [ 4].

Hypoxia is one of the most important hallmarks of the abnormal metabolic microenvironment in solid tumours [ 5]. The imbalance of pro- and antiangiogenic factors results in the formation of layered capillary network in tumours and the abnormal structure and functions of tumour vessels [ 6]. Since the diffusion limit of oxygen in tissues is 100-200 μm [ 7], the regions far from the blood vessels become chronically hypoxic. In addition, due to the heterogeneous and disorganized vessel network in tumours, there is no clear relationship between the blood flow rate and oxygen tension (pO₂ ) [ 8], which may cause tumour tissue hypoxia despite the presence of blood flow.

Hypoxia can upregulate various angiogenic growth factors, including vascular endothelial growth factor (VEGF), thus triggering vascularization of tumours [ 9]. Furthermore, hypoxia also modulates the rate of progression through the cell cycle [ 10]. Whereas most normal cells undergo apoptosis when the hypoxic stress is too intense or persists for too long, tumour cells appear to have much more resistance to hypoxia [ 11]. One reason for this phenomenon is that the hypoxia inducible factor 1 α ( HIF 1 α ) upregulates the protein p27 to cause the arrest of the G1/S transition in the tumour cell cycle [ 12].

MTS have often been regarded as aggregates of cells with intimate cell-cell contacts; however, it demonstrated the existence of extracellular matrix (ECM) in the spheroids [ 13]. The ECM is a complex network composed of a variety of proteins and polysaccharides such as fibronectin, collagen, laminin, hyaluronate, heparan sulfate, and elastin, that are produced by the cells to form this matrix [ 14]. The role of ECM in the MTS growth has been investigated, including the cell-matrix adhesion through integrins [ 15] and ECM-dependent cell survival [ 16]. Although few studies have been conducted examining the cell-matrix interactions in MTS, the importance of these interactions in regulating cell growth, especially that in three-dimensional spheroids, should be addressed.

To investigate the physiological phenomena and pathological characteristics of MTS, mathematical models for MTS growth have been well developed in recent years [ 17- 20]. Growth models for MTS are usually classified into three groups: (a) continuum models, formulated through partial differential equations (PDEs) to describe the changes in distributions of certain substances [ 21]; (b) discrete models, often represented by cellular automata (CA) [ 22], agent-based [ 23], and Monte-Carlo inspired models [ 24]; (c) hybrid models, incorporating continuum models and discrete models [ 25]. CA models were originally introduced by von Neumann [ 26] as a possible idealization of biological systems and typically constructed by the setup of the transition rules (CA rules) generating the desired behaviour. In recent years, CA models for MTS growth have been studied to improve the understanding of the dynamic process of MTS growth and the changes of microenvironments in MTS. The benefit of using CA models is that they provide a description and allow a more realistic stochastic approach at both cellular and subcellular levels. However, one of the main disadvantages of CA models is that they require very large amounts of computational resources.

We have established a 2D coupled mathematical model of tumour growth, angiogenesis, and blood perfusion in the previous work [ 27], in which the simulation of the early stage of tumour cell growth was simplified. At the avascular phase, the influence of hypoxic environment on the tumour growth is more significant than after angiogenesis process. Hypoxia is not only associated with the trigger of tumour angiogenesis, but also affects the tumour cell survival and proliferating. Moreover, some tumour cells can undergo a series of mutations in a certain hypoxic environment to become more aggressive and more resistant to chemotherapy and radiotherapy. Therefore, a mathematical model focuses on the hypoxia microenvironment, and associated MTS growth will be helpful in improving our understanding of the dynamic processes in early tumour growth.

The main aim of this study is to develop a 3D hybrid CA model to investigate the dynamic growth of MTS to the varied levels of hypoxic microenvironment. We define a continuum deterministic model (a system of conservation equations) to control the chemicals dynamics and a discrete CA model (individual-based processes) to control the individual cell migration and interactions. In addition, we focus on the influence of hypoxic microenvironment on the MTS growth by considering different cell phenotypes, which are caused by different hypoxia levels and may in turn induce different resistances of tumour cells to hypoxia.

2. Method

2.1. The Conservation Equations for Chemicals

The model is defined on a 3D region of 100 × 100 ×100 grids to cover a 2 mm × 2 mm ×2 mm space, for MTS cell proliferating, quiescent, and apoptosis. Each tumour cell occupies one individual element and competes for space and oxygen with others. The interactions of cells with the microenvironment are governed by the ECM concentration (denoted by _{c f} ) and the oxygen concentration (denoted by _{c o} ). The ECM can be degraded by the matrix-degrading enzymes (MDEs) (denoted by _{c m} ) which are produced by the tumour cells. The equations describing the interactions of tumour cells with ECM, MDE, and oxygen are [figure omitted; refer to PDF] The ECM concentration _{c f} is controlled by the degradation by MDE, while the MDE is governed by diffusion, production by the tumour cells, and decay of itself. Since we choose to focus on the MTS growth, the vascularization inside the tumour is not included in the current model. Therefore, we assume that oxygen is produced by the macromolecules in the ECM, transported by pure diffusion in interstitial space, and consumed by the tumour cells. For simplicity, oxygen production is assumed to be proportional to the ECM density [ 25]. The _{T i ,j ,k} term represents a tumour cell located at a node position (i ,j ,k ) . The value is either 1 if a cell is present or 0 if it is not. _{D m} and _{D o} are MDE and oxygen diffusion coefficients, respectively. δ , _{μ T} , λ , γ , and β are positive constants. The values of these parameters are shown in Table 1.

Table 1: Parameter values used in the simulation.

Parameter	Value
D m a	1 0 -9 c m 2 s -1
δ b	1.3 ×1 0 2 c m 3 M -1 s -1
μ T b	1.7 ×1 0 -18 M cell s -1 s -1
β a	5.8 ×1 0 -2 c m -3 s -1
D o a	1 0 -5 c m 2 s -1
λ a	1.7 ×1 0 -8 s -1
γ a	6.25 ×1 0 -17 M cell s -1 s -1

^a Pettet et al. [ 21].

^b Mansury and Deisboeck [ 23].

2.2. CA Model for an Individual Cell

Figure 1shows a flowchart of the hybrid CA model under hypoxic microenvironment. Four phenotypes of tumour cells are defined in the model. Types P, Q, N, and M represent the proliferating cells, the quiescent cells, the necrotic cells, and the mutated cells, respectively. Initially, all tumour cells are proliferating cells. At each time step, a tumour cell will increase its age and be checked to see if there is enough space and oxygen ( _{c o} ...5; _{c thres _n} =0.5 ) for its proliferation; (a) if there is, the cell will be divided by mitosis to two daughter cells which are proliferating cells; one daughter cell will replace the parent cell, and the other one will move to the neighbouring element possessing the highest ECM concentration; (b) if there is no space for proliferation, the cell will become quiescent; for every quiescent cell, if the oxygen is inadequate, the cell will apoptose to become a necrotic cell with a predefined probability ( 0.1 % × _{T age} ). Tumour cells have been found to be able to survive in a very poorly oxygenated environment [ 9]. Therefore, we assume a phenotype of mutated cells, which has the most resistance to hypoxia. If the proliferating cell and quiescent cell have survived in a certain hypoxic environment ( _{c o} ...4; _{c thres _m} =0.8 ) for a long time ( t ...5; _{T m} =200 ), they are assumed to become mutated cells with a certain probability (10%). Since the process of ECM regeneration and the removal of dead tumour cells are not clear, the dead cell will still occupy the space but will cease any interaction with the environment. It needs to be mentioned here that there is no reverse change of cell from quiescent back to proliferating in this model.

Figure 1: Flowchart of the hybrid CA model showing the different phenotypes (proliferating, quiescent, necrosis, and mutation) under hypoxic microenvironment.

[figure omitted; refer to PDF]

Different coefficients of the MDE production and the consumption rate of oxygen are defined for different phenotypes of tumour cells, as shown in Table 2. The mutated tumour cells have the biggest production coefficients of MDE, which means they are more aggressive and more likely to invade into surrounding tissues. Based on the work of Anderson [ 25], each cell has its own internal adhesion value _{A i} . Only if the number of external neighbours of the cell is larger than _{A i} , the cell is allowed to migrate. The different values of _{A i} of different phenotypes are also shown in Table 2.

Table 2: Parameters of the four different phenotypes.

Phenotypes	MDE production	Oxygen consumption	A i
Proliferating cells (P)	μ T	γ	2
Quiescent cells (Q)	μ T /5	γ /2	4
Necrotic cells (N)	μ T /10	γ /4	6
Mutated cells (M)	4 /3 μ T	2 γ	0

2.3. Simulation

The simulation was carried out on a 3D domain of 2 mm ×2 mm ×2 mm , divided uniformly into 100 × 100 ×100 grids. A no-flux boundary condition was imposed on the square grid to restrict the tumour cells, MDE, ECM, and oxygen to the volume within the grid.

The initial distributions of oxygen and ECM were predefined as constant nondimensional values (they were 1 in this study) in the simulation space. Initially, 50 proliferating tumour cells were centered around grid element (50, 50, and 50) with age one. The simulation step is set to be 45 mins, according to the basic mitosis time for tumour cells [ 24].

The simulation process is as follows:

at each major iteration,

(1) MTS growth according to the CA rules (described in Section 2.2);

(2) solve ( 1) and update the chemicals (oxygen, MDE, and ECM) concentrations;

(3) update the distribution of tumour cells based on the changed microenvironment.

3. Results

In the results presented here, the simulation starts from the initial tumour and finishes at 800 time steps.

3.1. MTS Growth History

Figure 2shows the growth histories of MTS under hypoxic microenvironment and the cross-section at the plane z =50 . Different colours represent different cell phenotypes (proliferating cells: blue; quiescent cells: green; necrotic cells: red; and mutated cells: yellow). At the beginning, most of the tumour cells stay in proliferating and quiescent. However, the necrotic cells arise after t =200 and spread rapidly in the tumour interior, due to the limited oxygen supply and high oxygen consumption of proliferating cells. Once the necrotic core is formed at t =400 (see the cross-section pictures more clearly), the quiescent cells decrease continually, and the hypoxic environment in the MTS is satisfied for mutated cells. Eventually, some proliferating cells have survived and formed finger-like invasions to the surrounding matrix at the periphery area, while the mutated cells, which are the most aggressive, are scattered inside the MTS (also see Figure 5(a)).

The growth histories of MTS in hypoxic microenvironment and cross-sections at plane z =50 . Different colours represent different phenotypes of tumour cells (proliferating cells: blue; quiescent cells: green; necrotic cells: red; and mutated cells: yellow).

(a)

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 3shows the growth curves of MTS in the model. The number of total cells and the spheroid volume first grow exponentially at early days and then slow down after ~600 time steps (almost 18 days). The exponent growth and the saturation size are consistent with experimental observations by a series studies from Freyer et al. [ 29- 31].

The growth curves of MTS in the model. The number of total cells (a) and the volume of spheroid (b) in time.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

3.2. Oxygen Concentration Distribution

The integrated oxygen concentration ^{c o *} along the x -axis at the plane z =50 is shown in Figure 4, where x =0 is the centre of the simulation region and the MTS. The hypoxia inside the MTS becomes more and more extensive with the growth: (a) the hypoxic region is always at the MTS interior and spreads gradually to surrounding tissue; (b) the lowest value and the average value (shown as the green dotted line in Figure 3) of oxygen concentration decrease simultaneously with the growth time. However, the lowest oxygen concentration occurs in the centre area and may move towards the periphery of the MTS (especially at t =800 ), which is consistent with our previous work of coupled tumour growth model [ 27]. This is partly due to the fact that the quiescent and necrotic cells consume less oxygen than the proliferating ones. In addition, no assumptions about reverse changes from quiescent cells back to proliferating cells in this model also contribute this result.

Figure 4: Distribution of integrated oxygen concentration ^{c o *} along x -axis at z =50 plane. x =0 is the centre of the simulation region and the MTS. The green dotted line is the average value of ^{c o *} . All values are normalized to be 0-1.

[figure omitted; refer to PDF]

The growth history of MTS cells under three different hypoxia conditions. The baseline model with (a) _{c thres _n} =0.5 and _{c thres _m} =0.8 , (b) _{c thres _n} =0.7 , (c) _{c thres _m} =0.7 , and (d) _{c thres _m} =0.9 .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

At the end of the simulation process, the lowest value of oxygen concentration in the MTS is stable at around 0.75, which is between the two thresholds for necrotic cells ( _{c thres _n} =0.5 ) and mutated cells ( _{c thres _m} =0.8 ). This implies that the mutated cells will increase with the MTS growth according to the assumptions used in the current model.

3.3. Sensitivity to Hypoxic Microenvironment

To assess the sensitivities to the oxygen thresholds in the model, we changed the parameters _{c thres _n} , _{c thres _m} and performed three more simulations, as shown in Figure 5. Figure 5(a)shows the baseline model presented in the above sections. The numbers of proliferating cells and quiescent cells reach their peak values at t =300 and decrease gradually at the late period. As a consequence, a rapid increase of necrotic cells appears after t =300 due to the inadequate oxygen supply caused by the increased proliferating cells and quiescent cells propagated at an earlier stage. At the end of the simulation (after t =600 ), all the MTS cells of four different phenotypes almost stop growing and tend towards stability, which is one of the signals to stop the simulation.

The threshold for necrotic cells _{c thres _n} will influence the cell distribution by controlling the cell death in response to low oxygen concentration. A higher value of _{c thres _n} (=0.7) will reduce the resistance of tumour cells to the hypoxic environment and increase the necrotic cell number at the end of the simulation, as shown in Figure 5(b). At the same time, the proliferating cell number decreases due to the rapidly increasing number of necrotic cells which occupy most of the volume in the simulation region. This result suggests that the resistance of tumour cells to the hypoxic environment may play an important role in the development of MTS.

Another threshold defined in this model for cell phenotype to local oxygen concentration is _{c thres _m} , which affects the mutated cell number. We assumed a phenotype of mutated cells to represent the most aggressive cells, which can survive in a very poorly oxygenated environment. The threshold for mutated cells _{c thres _m} has been changed to be 0.7 and 0.9, respectively, to test the system responses to this threshold. According to Figures 5(c)and 5(d), the changes of _{c thres _m} can affect the number of mutated cells. Although the influence on mutated cell number seems weak, the growth histories of other cells, especially the necrotic cells, have remarkable changes. This is due to the assumption that mutated cells have the strongest relationship with the hypoxic microenvironment. The varied mutated cell numbers may have a significant effect on the oxygen distribution in the simulation region and then influence the growth of the MTS. For example, a decreasing value of _{c thres _m} (=0.7) leads to a reduction of mutated cell number (Figure 5(c)), which will reduce the oxygen consumption in the MTS and improve the hypoxic environment. As a consequence, the growth speed of necrotic cell number decreases compared with the baseline model (Figure 5(a)). On the contrary, a higher value of _{c thres _m} ( =0.9 ) will induce the increase of the mutated cell number, as well as the growing speed of the necrotic cell number (Figure 5(d)).

4. Discussion

In this study, we established a hybrid CA model for 3D MTS growth by introducing hypoxia as an important microenvironment factor which influences cell migration and phenotype expression. It was assumed that (a) the distributions of the chemicals in the MTS (including ECM, MDE, and oxygen) are controlled by continuum equations and (b) the individual behaviour of single cell is determined by discrete CA rules. In addition, we introduced four different cell phenotypes (the proliferating cells, the quiescent cells, the necrotic cells, and the mutated cells) to reflect the different correlations between cells and environments.

The model has demonstrated the process of MTS growth and the spatial distribution of oxygen concentration in an MTS. The simulation results showed typical MTS growing features, such as a necrotic core in the centre, associated with the hypoxia region of the MTS, and more aggressive cells at the periphery area, which are consistent with the physiological knowledge (Figure 6).

Figure 6: An MTS section showing the proliferative rim and the necrotic core. This figure is reprinted from [ 28].

[figure omitted; refer to PDF]

To assess the model system responses to the varied hypoxic conditions, we examined the sensitivities to the oxygen thresholds in this model. The results showed that the resistance of tumour cells to hypoxic environment and the number of mutated cells influence the MTS growth significantly, which suggests that the improvement of the resistance of tumour cells to hypoxic environment may be important in the tumour normalization therapy.

Compared with our previous models of tumour development [ 27, 32, 33], the major improvements for the model presented here include the following: (a) 3D simulation instead of 2D to make the overall model more realistic; (b) more sophisticated assumptions for individual cell to let us go deep into the system responses of different environmental conditions; (c) the examinations of the sensitivities to the varied hypoxic environment show the importance of the resistance and adaption of tumour cells.

However, some physiological details were ignored or simplified during the model development. Firstly, we assumed phenotype M (mutated cells) to represent the most aggressive cells, having the most resistance to hypoxia, consuming the most oxygen, and being more likely invade to the surrounding tissues. This assumption was based on the fact that hypoxia can modulate the rate of progression of tumour cells through the cell cycle. However, the details of the regulation of cell cycle of certain cells were not included in the current model. The multiscale modelling considering the subcellular levels should be addressed in the future. Secondly, the possibility that quiescent cells revert to the proliferating stage was not included in the model. Many experimental investigations have revealed that quiescent cancer cells can reactivate their cell cycles, becoming proliferative again and thus restarting tumour growth, which is one common cause of chemotherapeutic treatment failure [ 34]. A model with CA rules for quiescent cells returning to the proliferating stage according to the microenvironment changing will be helpful in studying the impact of quiescent cells on MTS growth and chemotherapy effectiveness. Thirdly, the mechanical properties of ECM and the complex cell-matrix interactions were not investigated in the study. The inclusion of ECM in this modelling of MTS growth suggests that the model could be used to describe an in vivo tumour growth in the avascular phase. A deeper understanding of the interrelationships between the presence of ECM components and cell-cell contact may help to advance mathematical models to examine the migration and metastasis of in vivo cancer cells. These model potentialities will be the subjects of future studies.

Acknowledgments

This study was partially supported by the National Basic Research Program of China (973 Program) (no. 2013CB733800) and the National Natural Science Foundation of China (NSFC) (no. 11272091).

References

[1] G. Hamilton, "Multicellular spheroids as an in vitro tumor model," Cancer Letters , vol. 131, no. 1, pp. 29-34, 1998.

[2] W. R. Inch, J. A. McCredie, R. M. Sutherland, "Growth of nodular carcinomas in rodents compared with multi-cell spheroids in tissue culture," Growth, Development and Aging , vol. 34, no. 3, pp. 271-282, 1970.

[3] R. M. Sutherland, J. A. McCredie, W. R. Inch, "Growth of multicell spheroids in tissue culture as a model of nodular carcinomas," Journal of the National Cancer Institute , vol. 46, no. 1, pp. 113-120, 1971.

[4] J. M. Kelm, N. E. Timmins, C. J. Brown, M. Fussenegger, L. K. Nielsen, "Method for generation of homogeneous multicellular tumor spheroids applicable to a wide variety of cell types," Biotechnology and Bioengineering , vol. 83, no. 2, pp. 173-180, 2003.

[5] D. Fukumura, R. K. Jain, "Tumor microvasculature and microenvironment: targets for anti-angiogenesis and normalization," Microvascular Research , vol. 74, no. 2-3, pp. 72-84, 2007.

[6] D. Ribatti, B. Nico, E. Crivellato, A. Vacca, "The structure of the vascular network of tumors," Cancer Letters , vol. 248, no. 1, pp. 18-23, 2007.

[7] A. Krogh The Anantomy and Physiology of Capillaries , Yale University Press, New York, NY, USA, 1922.

[8] G. Helmlinger, F. Yuan, M. Dellian, R. K. Jain, "Interstitial pH and pO2 gradients in solid tumors in vivo: high-resolution measurements reveal a lack of correlation," Nature Medicine , vol. 3, no. 2, pp. 177-182, 1997.

[9] A. L. Harris, "Hypoxia--a key regulatory factor in tumour growth," Nature Reviews Cancer , vol. 2, no. 1, pp. 38-47, 2002.

[10] L. B. Gardner, Q. Li, M. S. Park, W. M. Flanagan, G. L. Semenza, C. V. Dang, "Hypoxia inhibits G1/S transition through regulation of p27 expression," Journal of Biological Chemistry , vol. 276, no. 11, pp. 7919-7926, 2001.

[11] T. Alarcón, H. M. Byrne, P. K. Maini, "A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells," Journal of Theoretical Biology , vol. 229, no. 3, pp. 395-411, 2004.

[12] J. J. Tyson, B. Novak, "Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions," Journal of Theoretical Biology , vol. 210, no. 2, pp. 249-263, 2001.

[13] T. Nederman, B. Norling, B. Glimelius, J. Carlsson, U. Brunk, "Demonstration of an extracellular matrix in multicellular tumor spheroids," Cancer Research , vol. 44, no. 7, pp. 3090-3097, 1984.

[14] M. T. Santini, G. Rainaldi, P. L. Indovina, "Apoptosis, cell adhesion and the extracellular matrix in the three-dimensional growth of multicellular tumor spheroids," Critical Reviews in Oncology/Hematology , vol. 36, no. 2-3, pp. 75-87, 2000.

[15] S. M. Frisch, E. Ruoslahti, "Integrins and anoikis," Current Opinion in Cell Biology , vol. 9, no. 5, pp. 701-706, 1997.

[16] F. Hirschhaeuser, H. Menne, C. Dittfeld, J. West, W. Mueller-Klieser, L. A. Kunz-Schughart, "Multicellular tumor spheroids: an underestimated tool is catching up again," Journal of Biotechnology , vol. 148, no. 1, pp. 3-15, 2010.

[17] P. P. Delsanto, C. Guiot, P. G. Degiorgis, C. A. Condat, Y. Mansury, T. S. Deisboeck, "Growth model for multicellular tumor spheroids," Applied Physics Letters , vol. 85, no. 18, pp. 4225-4227, 2004.

[18] G. Schaller, M. Meyer-Hermann, "Continuum versus discrete model: a comparison for multicellular tumour spheroids," Philosophical Transactions of the Royal Society A , vol. 364, no. 1843, pp. 1443-1464, 2006.

[19] C. A. Condat, S. A. Menchón, "Ontogenetic growth of multicellular tumor spheroids," Physica A , vol. 371, no. 1, pp. 76-79, 2006.

[20] B. Brutovsky, D. Horvath, V. Lisy, "Inverse geometric approach for the simulation of close-to-circular growth. The case of multicellular tumor spheroids," Physica A , vol. 387, no. 4, pp. 839-850, 2008.

[21] G. J. Pettet, C. P. Please, M. J. Tindall C.P., D. L. S. McElwain, "The migration of cells in multicell tumor spheroids," Bulletin of Mathematical Biology , vol. 63, no. 2, pp. 231-257, 2001.

[22] M. Aubert, M. Badoual, S. Féreol, C. Christov, B. Grammaticos, "A cellular automaton model for the migration of glioma cells," Physical Biology , vol. 3, no. 2, pp. 93-100, 2006.

[23] Y. Mansury, T. S. Deisboeck, "Simulating 'structure-function' patterns of malignant brain tumors," Physica A , vol. 331, no. 1-2, pp. 219-232, 2004.

[24] Y. Jiang, J. Pjesivac-Grbovic, C. Cantrell, J. P. Freyer, "A multiscale model for avascular tumor growth," Biophysical Journal , vol. 89, no. 6, pp. 3884-3894, 2005.

[25] A. R. A. Anderson, "A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion," Mathematical Medicine and Biology , vol. 22, no. 2, pp. 163-186, 2005.

[26] J. von Neumann Theory of Self-Reproducing Automata , University of Illinois Press, 1966.

[27] Y. Cai, S. Xu, J. Wu, Q. Long, "Coupled modelling of tumour angiogenesis, tumour growth and blood perfusion," Journal of Theoretical Biology , vol. 279, no. 1, pp. 90-101, 2011.

[28] T. Roose, S. J. Chapman, P. K. Maini, "Mathematical models of avascular tumor growth," SIAM Review , vol. 49, no. 2, pp. 179-208, 2007.

[29] J. P. Freyer, R. M. Sutherland, "Regulation of growth saturation and development of necrosis in EMT6/Ro multicellular spheroids by the glucose and oxygen supply," Cancer Research , vol. 46, no. 7, pp. 3504-3512, 1986.

[30] J. P. Freyer, R. M. Sutherland, "Proliferative and clonogenic heterogeneity of cells from EMT6/Ro multicellular spheroids induced by the glucose and oxygen supply," Cancer Research , vol. 46, no. 7, pp. 3513-3520, 1986.

[31] J. P. Freyer, "Role of necrosis in regulating the growth saturation of multicellular spheroids," Cancer Research , vol. 48, no. 9, pp. 2432-2439, 1988.

[32] Y. Cai, K. Gulnar, H. Zhang, J. Cao, S. Xu, Q. Long, "Numerical simulation of tumor-induced angiogenesis influenced by the extra-cellular matrix mechanical environment," Acta Mechanica Sinica/Lixue Xuebao , vol. 25, no. 6, pp. 889-895, 2009.

[33] Y. Cai, J. Wu, S. X. Xu, "Numerical simulation of inhibiting effects on solid tumour cells in anti-angiogenic therapy: application of coupled mathematical model of angiogenesis with tumour growth," Applied Mathematics and Mechanics , vol. 32, no. 10, pp. 1287-1296, 2011.

[34] S. A. Menchón, C. A. Condat, "Quiescent cells: a natural way to resist chemotherapy," Physica A , vol. 390, no. 20, pp. 3354-3361, 2011.

[]