- Open Access
Mathematical modelling of spatio-temporal glioma evolution
Theoretical Biology and Medical Modelling volume 10, Article number: 47 (2013)
Gliomas are the most common types of brain cancer, well known for their aggressive proliferation and the invasive behavior leading to a high mortality rate. Several mathematical models have been developed for identifying the interactions between glioma cells and tissue microenvironment, which play an important role in the mechanism of the tumor formation and progression.
Building and expanding on existing approaches, this paper develops a continuous three-dimensional model of avascular glioma spatio-temporal evolution. The proposed spherical model incorporates the interactions between the populations of four different glioma cell phenotypes (proliferative, hypoxic, hypoglychemic and necrotic) and their tissue microenvironment, in order to investigate how they affect tumor growth and invasion in an isotropic and homogeneous medium. The model includes two key variables involved in the proliferation and invasion processes of cancer cells; i.e. the extracellular matrix and the matrix-degradative enzymes concentrations inside the tumor and its surroundings. Additionally, the proposed model focuses on innovative features, such as the separate and independent impact of two vital nutrients, namely oxygen and glucose, in tumor growth, leading to the formation of cell populations with different metabolic profiles. The model implementation takes under consideration the variations of particular factors, such as the local cell proliferation rate, the variable conversion rates of cells from one category to another and the nutrient-dependent thresholds of conversion. All model variables (cell densities, ingredients concentrations) are continuous and described by reaction-diffusion equations.
Several simulations were performed using combinations of growth and invasion rates, for different evolution times. The model results were evaluated by medical experts and validated on experimental glioma models available in the literature, revealing high agreement between simulated and experimental results.
Based on the experimental validation, as well as the evaluation by clinical experts, the proposed model may provide an essential tool for the patient-specific simulation of different tumor evolution scenarios and reliable prognosis of glioma spatio-temporal progression.
Introduction and background
Cancer is the second most fatal disease worldwide after heart disease . A cancer cell evolves from normal due to genetic mutations, which abnormally alter the cell proliferation rate. In particular, glioma is a rapidly evolving type of brain cancer, well known for its aggressive and diffusive behavior . This diffusive invasion has lead several research efforts to explore the glioma’s progression with the aid of mathematical diffusion equations [3–5], aiming to predict its spatial and temporal evolution. The high diffusion rate of glioma cells from the core tumor into the surrounding brain tissue often leads to treatment failure and tumor recurrence, even after the surgical resection. Gliomas vary from low- to high-grade, namely glioblastomas, which constitute the most malignant form of brain cancer, having an extremely poor prognosis.
Several diverse factors influence glioma progression. Tumor cell migration in brain can be stimulated by the surrounding extracellular matrix (ECM) that is a complex mixture of macromolecules, some of which like the collagens play a structural role and others, such as laminin, fibronectin and vitronectin are important for cell adhesion, spreading and motility. ECM is remodeled by glioma cells that are able to degrade the surrounding tissue, through the secretion of matrix-degradative enzymes (MDEs), such as the plasminogen activator and the family of matrix metalloproteinases. Essentially, all tumors go through an avascular stage of growth, where the nutrient (oxygen and glucose) supply is adequate. The nutrients are provided both by the intratumoral ECM and the surroundings of tumor via diffusion. In a future growth stage, the intratumoral ECM is completely degraded, while the surrounding extracellular material is not close enough to supply the entire tumor body with the necessary nutrients via diffusion. Thus the tumor area is divided into layers consisting of different cell populations, namely proliferating, hypoxic, hypoglycemic and necrotic regions of tumor cells. Fully developed tumors contain a necrotic core with absence of nutrients, surrounded by hypoxic and hypoglycemic rings with limited access to nutrients and an outward proliferating cell ring having full nutrient supply. Since the tumor dependents on nutrients diffusion, its avascular growth is limited. Further growth and proliferation leads to its vascular stage, where cancer cells enhance the existing vascular network through angiogenesis. Furthermore, the metastatic final tumor stage succeeds the vascular stage. Since cancer development comprises these distinct stages (avascular, vascular and metastatic), the majority of research efforts are concentrated on exploring each one of these stages separately. In particular, identifying the avascular tumor dynamics constitutes the first crucial step towards fully vascularized tumors investigation .
In parallel to identification of tumor characteristics, the prediction of tumor growth and diffusion can lead to useful insight into the disease dynamics, which may improve clinical outcomes. To this respect, several mathematical and computational models have appeared in the literature, which investigate the mechanisms that govern glioma’s progression and invasion, with the aim of predicting its future spatial and temporal evolution, with or without the effects of therapy . The models may constitute valuable tools for assisting the clinical practice towards the optimal individualized treatment, while facilitating medical research analysis. Current modelling techniques focus on either simulating individual cell behavior, or modelling tumor of clinically significant size, exploring volume expansion [7, 8]. Such models follow a discrete cell-based description that simulates cell growth, or a continuum framework that deals with the evolution of local cellular densities and is based on diffusion-reaction equations. Since both approaches reveal particular tumor aspects but also pose certain limitations, hybrid models have been developed to combine continuum and discrete variable descriptions, but they are still limited by model complexity and uncertainties in the mixing of models. Alternatively, the notion of multi-compartment continuum model appears to be effective in describing how subpopulations of various cell types proliferate and diffuse, while it remains computationally efficient .
Discrete models are able to incorporate cellular parameters and biological rules, as well as simulate the effects of chemotaxis, haptotaxis, cell to cell adhesion and address intracellular processes and intercellular communication issues. Cancer cells are usually represented by vectors of variables that determine their position and speed in space, their cell cycle phase, etc [10–13]. Moreover, these models are able to easily describe heterogeneous cell populations. However, they are computationally intensive and fail to be efficiently initialized with tumors of clinically significant size that already consist of over 1 billion cells. The discrete model of 1000 billion individual units is beyond the capability of most computers, even if simple interactions algorithms are exploited . Thus, such models are limited to the exploration of small tumors, below the level of clinical detection.
Continuum models represent alternative approaches in exploring tumor progression dynamics. Rather than considering each cell as a discrete individual, they assume tumors as spatial distributions of cell densities. As a result, computational requirements are substantially reduced. If the population of tumor cells is obtained as a continuum variable, the system can be described at macroscopic level as a biological tissue sample. Continuum models are also capable of accounting for the dynamic changes of chemical ingredients (such as glucose, oxygen, MDEs, chemotherapeutic drugs, etc) in the brain tissue. Additionally, they can take into account brain tissue heterogeneity and anisotropy. However, they lack a capability to simulate specific cellular factors, such as the genetic fingerprint of tumor cells and discontinuous changes, e.g. epithelial changes important for the invasion of tumor cells in adjacent tissue. Additionally, they are less efficient than discrete models in describing cellular-scale flux factors such as chemotaxis, haptotaxis and cell to cell adhesion.
Continuum approaches are most usually implemented by means of reaction-diffusion equations , where the cell proliferation factor is alternatively expressed by the Gompertz law, the second order polynomial equation, or the simple linear function [5, 16]. Recent developments in mathematical modelling of gliomas are extensively analyzed in . First, a model of untreated glioma is described followed by the presentation of models including chemotherapy, or surgical resection. The modelling assumptions vary from relatively simple considerations dealing with homogeneous brain tissue, to those involving heterogeneous brain tissue with different diffusion rates of glioma cells in grey and white matter on a geometrically complex brain domain. In the approach of , a three-dimensional diffusion model has been developed to address brain heterogeneity and anisotropy issues. The study introduces a mathematical framework incorporating alternative proliferation rate schemes and presents a solution of the diffusion equation through different numerical approximation of finite differences. Brain heterogeneity and anisotropy issues are dealt in the glioma growth model of , which is applied on real data. The proportions of white and grey matter, as well as the diffusion tensors are extracted by normal brain atlases, avoiding DTI processing. The research effort in  supports the hypothesis that the change in pH of tumor microenvironment may provide an important mechanism for cancer invasion. Another reaction-diffusion model presented in , deals with the diffusive behavior of glioma growth. This homogeneous model includes the Neumann boundary conditions imposed by the skull on gliomas and especially on glioblastoma multiforme (GBM). In  growth and diffusion issues of U87 glioblastoma tumor spheroids are investigated. The model includes diffusion, proliferation and cell to cell adhesion parameters, which differ for each cell line. U87 glioma cells are also explored in , where a model that includes the role of chemotaxis, haptotaxis and adhesion is presented. According to this model cell migration is strongly depend on adhesion, haptotactic and chemotactic parameters, apart from the random diffusion factor.
Expanding the potential of continuum mathematical models, some research efforts also engage the mechanical influence of tumor cells on the invaded tissue[5, 22–26]. Many continuum approaches exploit the exponential Gompertz law  in order to model the proliferation of cancer cells. According to this law, the growth rates of cell populations are high at early stages of growth, while they are slower at later stages. The Gompertz tumor growth models take into account the fact that as the tumor size increases, the growth rate is reduced because of the reduction of available space and the limited supply of nutrients. The total cell population follows the Gompertz growth model in , where a two-compartment model of cancer cells population dynamics is proposed. The model consists of a set of ordinary differential equations, which simulate proliferating and quiescent cell populations including the transition rates between them. Moreover, Gompertz law and Gompertzian parameters available in the literature are exploited in the models of  and .
Recently, hybrid models have been developed to overcome particular limitations of distinct approaches. These models exploit the continuum method in order to simulate the tumor microenvironment, while cell to cell interactions are implemented using the discrete approach [14, 31–37]. Some of these models are multiscale, incorporating the molecular and cellular level [6, 38–40]. However, hybrid models are still constrained by the discrete-model limitations in which they do not consider tissue-level associations and they still have a significant computational cost.
Multi-compartment continuum models form alternative continuum approaches that are expected to overcome some of the conventional modelling limitations and are able to include multiphase heterogeneous populations and heterogeneous extracellular matrix . According to this model philosophy, cells are grouped based on phenotype depending on their access to the necessary nutrients, or the level of differentiation . Such approaches are computationally efficient and effective in describing how subpopulations of various types of cells proliferate and diffuse. Invasion, proliferation, changes in phenotype and necrosis are readily expressed as additions or subtractions from the cell densities within each compartment. In  a continuum multi-compartment model is proposed investigating how tumor cells interact with their tissue microenvironment. The dispersal and interactive population changes of normoxic and hypoxic glioma cells, vascular endothelial cells, diffusible angiogenic factors and necrotic cells are implemented using diffusion equations. Furthermore, the study demonstrates how different proliferation and diffusion rates of glioma cells result to increasing degrees of mitoses, hypoxia-induced neoangiogenesis and necrosis, features that characterize the different glioma grades.
In this paper, a new continuum three-dimensional spherical and multi-compartmental mathematical model of avascular glioma growth in an isotropic and homogeneous medium is proposed. This model expands the multi-compartmental approach of  by integrating a new cell compartment that represents the hypoglycemic cell population. The initial model is further extended through the exploitation of the microenvironmental variables of the hybrid approach in , where apart from oxygen, we introduce a second vital nutrient, (i.e. glucose). The proposed model consists of heterogeneous tumor cell populations incorporating the interactions between four different glioma cell phenotypes into distinct cellular compartments, namely proliferative, hypoxic, hypoglycemic and necrotic, as well as their tissue microenvironment. Moreover, the model includes the effect of the host tissue, i.e. the extracellular matrix (ECM), the matrix-degradative enzymes (MDEs) and nutrients (oxygen and glucose) concentration on tumor cell proliferation and invasion, through different microenvironmental compartments. The proposed approach employs continuum variables (cell densities, concentrations of chemical ingredients) governed by diffusion principles. It incorporates the effects of concentration for both important nutrients in the tumor microenvironment, namely oxygen and glucose. Notice that the majority of existing approaches deal with a single nutrient, which most often stands for oxygen [14, 21, 24, 26, 32, 35, 37, 40], and causes the development of hypoxic cell regions. Furthermore, existing approaches employing both vital nutrients (e.g. [1, 6, 34, 39, 44]) consider only the combined effect of their availability, which accounts for the generation of only quiescent (non-proliferating) and necrotic cells. The novelty of our proposed model lies in the fact that apart from the combined effect, it takes into consideration the separate and independent impact of each nutrient in tumor growth, which leads to the formation of a new cell population (apart from hypoxic, or quiescent and necrotic), namely hypoglycemic cells, with a different metabolic profile. Moreover, this model includes innovative variation of particular factors, such as the local cell proliferation rate, the variable conversion rates of cells from one category to another and the nutrient-dependent thresholds of conversion. The model is validated through a comparison with specific experimental results of glioma models available in the literature and additional evaluation by clinical experts. The results demonstrate the model’s efficiency, providing an essential tool for the patient-specific simulation of different tumor evolution scenarios and reliable prognosis of glioma spatio-temporal progression.
The paper is organized as follows. In Section “Mathematical modelling”, the different cell and chemical compartments are described in detail as building blocks of the model, coupled with the appropriate differential equations. Section “Results and model validation” presents and discusses experimental results of tumor growth for different combinations of diffusion and proliferation rates and for various evolution times. Moreover, model validation is performed through comparisons with experimental results of glioma models available in the literature, as demonstrated in the same section. Finally, significant conclusions are drawn in Section “Discussion and conclusions”, where potential future work aiming to model improvement and completion is also discussed.
Current glioma modelling techniques, either apply discrete cell-based tumor growth simulations focusing in individual cell behavior, or follow continuum approaches dealing with the evolution of tumor cell densities. The discrete models are computationally feasible only if the initial number of tumor cells is relatively small, which is not possible for a tumor of clinically significant size that is detectable in current imaging techniques (MRIs, CTs, etc). Taking this into consideration, the rationale behind choosing to develop a continuum glioma evolution model in this paper relates to the need of dealing with virtual tumors of clinically visible size and density, consisting of large numbers of cancer cells. This model is initialized by a virtual spherical glioma tumor and simulates multiple cellular-microenvironmental interactions aiming to predict its three-dimensional spatio-temporal evolution.
where J represents the flux/motility of tumor cells, S stands for their net proliferation and T is the factor concerning the death of cancer cells due to microenvironmental conditions (e.g. absence of nutrients) and the beneficial effect of the applied treatment.
The motility of the tumor cells can be expressed as the result of four different fluxes J=J R a n d o m +J H a p t o t a x i s +J C h e m o t a x i s +J A d h e s i o n , namely the random cell diffusion, as well as the diffusion due to haptotaxis, chemotaxis and cell to cell adhesion [20, 21]. J R a n d o m =D R ∇C concerns the random cell diffusion, which is constant in a homogeneous medium and J H a p t o t a x i s =D H a p t ·C∇f stands for the haptotaxis that is the directed migratory response of cells to gradients of non-diffusible chemicals, such as the extracellular matrix. Moreover, J C h e m o t a x i s =D C h e m ·C∇G l is the tendency of cancer cells to move in the direction of the nutrients (glucose) gradient and J A d h e s i o n =C·K(C) corresponds to the movement due to cell to cell adhesion, according to which cells adhere to each other when they are close, while push apart when they are too compressed. In these equations D R , D H a p t , D C h e m , are the random diffusion, the haptotactic and the chemotactic coefficient respectively, C, f, G l the corresponding cancer cells, extracellular matrix and glucose concentrations and K represents the adhesive force of tumor cells.
The net proliferation rate of tumor cells S can be described by either a simple linear function S(C,t)=ρ·C, the second order polynomial equation , or the Gompertz law [5, 16], where C is the concentration of cancer cells, ρ corresponds to the proliferation rate constant and C m stands for the tissue cell maximum carrying capacity.
Cancer cells are capable of degrading the surrounding tissue, i.e. the extracellular matrix (ECM), which constitutes a crucial part of the invasive process. This degradation is achieved by means of the matrix-degradation enzymes (MDEs) that are produced by the cells and governed by the diffusion-reaction equation .
where m represents the concentration of MDEs, D m is the MDEs diffusion coefficient, C stands for the concentration of cancer cells, μ is the MDEs production rate and λ stands for the MDEs natural decay rate. The concentration change of ECM is simulated by the following equation .
where e represents the concentration of ECM and δ is the ECM degradation rate.
The tumor cells in order to survive, proliferate and diffuse need nutrients including oxygen and glucose. Many of them follow the anaerobic metabolism (glycolysis) instead of the normal aerobic metabolism that utilizes oxygen, even under of adequate oxygen presence [6, 45]. The cells that follow the glycolytic pathway do not use oxygen; however they need a much larger quantity of glucose, in order to produce energy, in comparison to those that have an aerobic metabolism . Hence, the tumor requires both oxygen and glucose, which are provided by the blood vessels and diffuse to the tumor area through the intratumoral and surrounding ECM. The nutrient production, diffusion and consumption are expressed by the diffusion-reaction equation .
where n represents the concentration of the vital nutrient, D n is the nutrient diffusion coefficient, β is the nutrient production rate through the ECM, γ stands for the nutrient consumption rate by the cancer cells and α is its natural decay rate.
Based on these equations, the proposed model expands the idea of glioma cell diffusion and proliferation taking into consideration different cell categories according to their metabolic profile, as well as the critical effect of two vital nutrients, which are consumed by the tumor cells.
Proposed model description: concepts and assumptions
The proposed continuum mathematical approach models the three-dimensional spherical spatio-temporal evolution of an avascular glioma in an isotropic and homogeneous medium. The model incorporates specific interactions between the heterogeneous populations of cancer cells and their microenvironment, which affect tumor growth and invasion. It takes into consideration the separate and independent effect of the concentration changes for both of the important vital nutrients in the tumor microenvironment, namely oxygen and glucose, leading to the formation of distinct cell populations with different metabolic profile and invasion properties. This innovative feature of the proposed cancer model can deal with different cell-profile changes for each nutrient, thus, allowing the formation of different areas in the cancer neighborhood. In other related approaches, when the concentration of one of the two nutrients is sufficiently reduced (no matter which) only one new cell population (e.g. quiescent cells) appears before leading to necrosis. The proposed model consists of continuous variables, described by diffusion equations. The invasion of tumor cells, which is governed by the diffusion principles is based on the following equation:
where C is the concentration of cancer cells at tumor point (r,θ,ϕ) (spherical coordinates) at time t, D represents the coefficient of random cell diffusion and f stands for the cell net proliferation rate. The proliferation rate in this model is based on the second order polynomial equation, since it takes in account the ratio of the local cell density and the maximum carrying capacity of the tissue, in order to limit the cell proliferation. Thus, when the tissue reaches its maximum carrying capacity, the tumor cells in this region, stop proliferating and remain quiescent, which is reasonable due to the lack of space. Furthermore, our model allows for variable local-cell proliferation rate due to the nutrients availability.
In order to decrease the complexity of the model, we express the invasion of tumor cells only due to random diffusion, i.e. haptotaxis, chemotaxis and cell to cell adhesion have been neglected. The latter mechanism, according to which cells exert adhesive forces to each other, is observed only on the cellular scale. Since our model deals with the macroscopic consideration of cell densities, it is reasonable to ignore these forces. Moreover chemotaxis, namely the tendency of cancer cells to move in the direction of the nutrients, has been neglected in the current model. This consideration does not significantly degrade the approximation of real tumor cell-motility, since the tendency of the cells to move in the direction of nutrients supply is implicitly taken into account in our model. More specifically, due to random diffusion, the tumor cells tend to move from the areas of large cell concentrations towards less crowded regions, which however hold higher nutrient concentrations because of their lower consumption rate in comparison to the regions of large cell populations; hence, tumor cells are indirectly forced move towards nutrients. Similarly, the effects of haptotaxis, which is the directed migratory response of cells towards non-diffusible chemicals, e.g. ECM, are implicitly considered; through random diffusion, cells tend to move towards regions of lower densities, which happen to hold higher concentrations of ECM.
The proposed continuum model of spherical shape is multi-compartmental. The evolution of tumor mass is modeled by means of coupled diffusion equations, which simulate the different cell densities. The compartmental approach is exploited to track glioma cell subpopulations based on viability and phenotype. Thus, this model includes a proliferative region, along with a hypoxic, hypoglycemic and a necrotic one. The cellular size of all phenotypes is assumed to be the same and the cell diffusion coefficients are considered constant, since the tumor is assumed to evolve inside an isotropic and homogeneous medium. Additionally, coupled continuum components are developed in order to account for the changes in the local environment that affect cells within each compartment. These components concern the diffusion and consumption of vital nutrients, concentration and degradation of the host tissue, i.e. ECM, as well as production and diffusion of MDEs. The diffusion coefficients of nutrients and MDEs are assumed constant within all tumor heterogeneous regions and outside the tumor as well. The incorporated nutrients that determine tumor cell viability and affect proliferation are the oxygen and glucose.
Oxygen and glucose are transported to the tumor region from blood vessels and diffuse through the ECM located within the tumor and its surroundings, which is destroyed by MDEs produced by the tumor cells as the tumor grows. Initially, the size and cell density of the tumor are small, thus all cells are sufficiently supplied with oxygen and glucose through the intratumoral and surrounding extracellular material [1, 6, 21, 26, 31, 36]. At this time all cancer cells are normoxic and proliferative, since the diffusion of nutrients is adequate for their survival and proliferation and most of them follow the aerobic metabolism. Furthermore, these cells compete for space with other neighboring cells; if the local tissue reaches its maximum carrying capacity, they stop proliferating and remain quiescent, while additionally they reduce their metabolism .
As time proceeds, the tumor radius increases and the concentrations of oxygen and glucose decrease in the central part of the tumor, closer to its core. At a certain time, each nutrient reaches a critical concentration value, below which it is not sufficient to meet the needs of all cells. Nutrients do not simultaneously fall below their critical concentrations, but are depleted independently depending on their consumption rates. Once the first nutrient, say oxygen, attains its critical concentration, the local tumor cells become hypoxic at a specific rate; a hypoxic region appears in the tumor’s center. Similarly, when the concentration of the other nutrient, i.e. glucose, falls below a critical level insufficient to meet energy needs of all cells, local cells turn to hypoglycemic and a hypoglycemic region appears in the central part of the tumor. The hypoxic and hypoglycemic cell populations neither proliferate nor die, while their metabolism changes. Specifically, the hypoxic cells perform increased glycolysis , which requires much less oxygen; the amount of oxygen consumed in hypoxic regions is not precisely known and could actually reach five times less the oxygen consumption in proliferative cells [32, 46]. However, since anaerobic metabolism is far less efficient in producing energy compared to the normal aerobic one, hypoxic cells eventually consume over 10 times more glucose for maintaining the same energy turn-over . As far as hypoglycemic cells are concerned, they reduce their metabolism, i.e. consumption of nutrients [6, 32, 35] and remain quiescent (they stop proliferating). Moreover, hypoxic cells are more migratory than normoxic proliferative cells , thus their diffusion coefficient is higher while they secrete larger amounts of degrading enzymes. When hypoxic cells migrate to an area with adequate oxygen, they convert back to normoxic at a specific constant rate; exactly the same applies to the case of the hypoglycemic cells, with respect to glucose . The cells that remain in a state of oxygen and glucose adequacy in a hypoxic or a hypoglycemic area continue proliferating; however, as time proceeds these cell populations decrease, since they increasingly convert to hypoxic or hypoglycemic respectively. Eventually, in the hypoxic and/or hypoglycemic areas the total proliferation rate is reduced to the appropriate levels. As the tumor evolves, both hypoxic and hypoglycemic populations grow in the central part of the tumor at the expense of proliferative cells, to finally form a hypoxic and/or a hypoglycemic core within the tumor.
When the vast majority of cells have turned to hypoxic or hypoglycemic in the central part of the tumor, one of the nutrients (oxygen or glucose) tends to vanish and tumor cells start dying, i.e. become necrotic. Thus a necrotic region appears close to the tumor core, where nutrients can no longer diffuse due to the long distance from the tumor surroundings. Additionally, necrotic cells appear in a tumor area where both nutrients have achieved their critical concentrations, regardless of hypoxic and/or hypoglycemic cell proportion. Actually, any of the nutrients can be the first to extinct at the tumor center depending on the initial concentrations, their diffusion coefficients and consumption rates. The necrotic cells do not have a metabolism, do not proliferate and do not move, since they are not alive. Moreover, necrosis can also result from the viable proliferative cells due to their contact with the necrotic ones, at a particular rate . In addition to necrosis, the natural cell death process is also present in the form of, apoptosis . However, apoptosis has been neglected in this model as to keep the complexity to a manageable level, since the apoptosis rate is low enough in comparison to the necrosis rate and does not significantly affect the cell proliferation and invasion processes.
Model equations for tumor-cell densities
Based on the above considerations of tumor progression, our model eventually considers tumor volume formed in four areas as in Figure 1, composed of the necrotic core, the hypoxic along with hypoglycemic rings around it, and the outer proliferating ring. Each of these areas consists mostly of cells belonging to the corresponding category. As long as the tumor radius increases, the hypoxic, hypoglycemic and necrotic radii increase as well. The tumor reaches a steady state when the hypoxic, hypoglycemic and necrotic regions occupy most of the entire tumor body.
In the model of Figure 1, four distinct glioma cellular compartments are developed. Each compartment constitutes a continuum equation, coupled with the remaining cellular compartments via common factors between them. The invasion of proliferative cancer cells follows the diffusion-reaction equation. The cell proliferation rate is expressed by the second order polynomial function, which takes into consideration the tissue maximum carrying capacity and diminishes the rate inversely according to the local cell density raise. As mentioned before, in a region of maximum cell density, the tumor cells apart from the fact that they do not proliferate, they reduce their metabolism (i.e. nutrient consumption) as well. Although these quiescent cells have temporarily stopped proliferating, they are included in the population of proliferative cells, since they still stay on a normoxic and normoglycemic state just as proliferative cells. The cell density equations are based on  and have been expanded to include the new cell compartment that represents the hypoglycemic population, as well as the effect of oxygen and glucose concentration in cell type (proliferative, hypoxic, hypoglycemic, necrotic) and proliferation rate. The equations governing the growth and invasion of proliferative cells are presented below.
This equation involves several terms describing the various cellular processes discussed above, where C, H, Q, N denote the concentrations of proliferative (normoxic/normoglycemic), hypoxic, hypoglycemic and necrotic cancer cells respectively at tumor point (r,θ,ϕ), at time t, D c represents the coefficient of random proliferative cell diffusion and T stands for the fraction of the local cell density to the maximum carrying capacity, given by:
Furthermore, b h , b q are the conversion rates of proliferative cells to hypoxic and hypoglycemic respectively, which unlike  are variable proportional to the oxygen or glucose relative concentrations. Additionally, a n stands for the constant of conversion rate of proliferative cells (as well as hypoxic and hypoglycemic) to necrotic due to the contact with the necrotic region and g h , g q , are the constant conversion rates of hypoxic and hypoglycemic cells to proliferative, respectively. Moreover, n h , g l q represent the thresholds for the proliferative cells to turn to hypoxic or hypoglycemic and in this model are related to the critical nutrients concentrations. The function f in Equation 7 stands for the local proliferation rate of tumor cells and corresponds to an innovative net proliferation rate variation, since in our model depends on the local oxygen and glucose availability. It can be modeled as:
where ρ corresponds to the proliferation rate constant and O c o n , G l c o n are coefficients that reduce the proliferation rate in relation to the local oxygen (n) and glucose (gl) concentration respectively and are given as in the following equations. The concentrations n and gl have been calculated including the current consumption by the cells at each time point, while n t m p and g l t m p are the corresponding concentrations before nutrient consumption. Additionally, n0 and g l0 represent the maximum concentrations of the nutrients which initially exist in the tumor surroundings. In particular,
The invasion of hypoxic cells is governed by diffusion, following the equation below:
where D h represents the coefficient of random hypoxic cell diffusion, a h is the variable (proportional to b h ) conversion rate of hypoxic cells to necrotic due to the lack of oxygen, while a g l h , is the conversion rate of hypoxic cells to necrotic, due to the absence of glucose, n h n , g l h n represent the thresholds for the hypoxic cells to turn to necrotic due to the absence of oxygen, or glucose respectively, modeled as:
Similarly the diffusion of hypoglycemic cells is expressed by the following equation:
where D q represents the coefficient of random hypoglycemic cell diffusion, a q is the conversion rate of hypoglycemic cells to necrotic due to the lack of oxygen, while a g l q is the variable (proportional to b q ) conversion rate of hypoglycemic cells to necrotic, due to the absence of glucose. Furthermore, n q n , g l q n represent the thresholds for the hypoglycemic cells to turn to necrotic due to the absence of oxygen or glucose respectively, described in the form:
Finally, necrotic cells do not diffuse and thus their concentration changes follow the rather simple equation:
The factors participating in Equation 14 have all been mentioned above.
Model equations for nutrient concentrations
The concentration-changes of oxygen and glucose are similarly governed by diffusion equations, which includes their production through the ECM, their natural decay and their consumption by the tumor cells. The oxygen concentration equation, as well as the new glucose concentration compartment is based on . Both of the nutrients equations have been expanded to involve the different consumption rates by each of the distinct cell types (proliferative, hypoxic and hypoglycemic). More specifically:
where n, gl represent the concentration of the oxygen, glucose at tumor point (r,θ,ϕ), at time t, D n , D g l are the respective diffusion coefficients, β n , β g l are the nutrients production rates through ECM, α n , α g l are the nutrients natural decay rates, γ c n , γ c g l stand for the nutrients consumption rates by the proliferative cancer cells, γ h n , γ h g l are the nutrients consumption rates by the hypoxic cells and respectively γ q n , γ q g l are the nutrients consumption rates by the hypoglycemic cells. In the tumor area where local cell-density has reached the tissue maximum carrying capacity, the normoxic/normoglycemic cells have been considered in this model to consume half of the nutrients’ amounts .
The concentration change of ECM because of its degradation by MDEs produced by the tumor cells, is demonstrated by the following equation (based on ):
where e is concentration of ECM at tumor point (r,θ,ϕ), at time t, m represents the concentration of MDEs and δ is the ECM degradation rate. MDEs are subjected to the principle of diffusion, produced by proliferative, hypoxic and hypoglycemic cells and are characterized by a natural decay rate. The MDEs equation is based on  and has also been expanded to include the distinct tumor cell categories:
where D m is the diffusion coefficient of MDEs, μ c , μ h , μ q are the production rates by the proliferative hypoxic and hypoglycemic cells respectively and λ represents the MDEs natural decay rate.
One of the major limitations of mathematical models concerns the parameter initialization using literature-based reference values. These values are often difficult to obtain, since most of them are patient specific and vary significantly among clinical and experimental cases. However, the existing wide range of values for each parameter allows models to incorporate different abnormal cases and examine many alternative outcomes.
A continuum glioma model can be initialized based either on the actual tumor geometry derived from medical images, or on a virtual spherical tumor , which approximates a real one, in terms of its imaging detectable size. In imaging modalities such as CT and MRI, the tumor is detectable above a certain size and density [3, 4]. According to the literature, the density threshold associated with the detectable boundary of a glioma in a medical image has been estimated to C2D=400c e l l s/m m2, which leads to , while the smallest detectable size corresponds to a radius of r=1 c m. According to medical experts tumors of radius smaller than 1 c m (e.g., r=5 m m) can be diagnosed. Yet, taking into consideration that oxygenated cells are limited to a distance of less than 0.1−0.2 m m from blood vessels, most of the initially detectable tumors (r≥5 m m) should already contain at least a hypoxic, even a necrotic region. Thus, in order to address tumorogenesis at early stages, our model is initialized with a virtual tumor of a smaller radius (i.e. r=1 m m), under the assumption that the initial tumor consists only of one cell category, namely proliferating cells without a significant obliquity. Accordingly, the initial tumor density has been considered to be slightly higher than the corresponding minimum detectable threshold, i.e. C0=104c e l l s/m m3. This constitutes the initial condition for the current model, which involves two boundaries; the inner one (r=r0) corresponding to the tumor border and the outer boundary (r=r1) taken far enough from the tumor and defining the area of the medium to be invaded by the tumor. The initial condition for tumor cell density is described as:
The invasion rate of glioma varies among patients depending on the grade of the disease. A wide range of invasion rates has been observed in clinical experience from serial MRIs and can be found in the literature in the form of several diffusion coefficient values. Similarly, for another patient-specific parameter namely the constant of tumor cell proliferation rate, a wide range of values has been reported. Different combinations of the diffusion and proliferation rates determine the different glioma grades. For instance, concerning a four grade malignant glioma (glioblastoma), a primary one is characterized by low diffusion and high proliferation, while a glioblastoma progressing from lower grade is determined by higher diffusion but lower proliferation rate . In this paper, various combinations of diffusion and proliferation rates have been considered, in order to explore different invasion and growth cases associated to various glioma grades. As far as the diffusion of hypoxic cells is concerned, since they are more migratory than normoxic cells , their invasion rate is higher; here their diffusion coefficient is taken one order of magnitude higher compared to the normoxic cells value. On the other hand, due to the lack of related data, the diffusion coefficient of hypoglycemic cells has been considered the same as the one of proliferative cells.
Initially, apart from outside the tumor, ECM is located within the tumor area as if it has not been completely degraded by MDEs. In our consideration, the initial intratumoral ECM concentration is considered one order of magnitude lower than the ECM maximum concentration at the tumor surroundings. As the tumor grows, intratumoral ECM is destroyed by MDEs, which are produced by the proliferative, hypoxic and hypoglycemic cells at specific rates. Since hypoxic cells are more invasive than normoxic cells, their MDEs secretion is considered larger and this is reflected in our model by accounting a twice higher MDEs production rate. Regarding hypoglycemic cells, they are associated to the same MDEs production rate as proliferative cells. The initial concentration of MDEs has been considered zero outside the tumor boundary and nonzero inside the tumor. More specifically, since the initial tumor density is lower than the tissue maximum carrying capacity, the initial intratumoral MDEs concentration is considered lower than the maximum reference value, equal to the ratio of the tumor density over the maximum carrying capacity.
The initial tumor is sufficiently provided with nutrients through the remaining intratumoral ECM, as well as from the surrounding ECM via diffusion. The maximum nutrients concentration exists in the area outside the tumor, which serves as the nutrient source. In the current model, nutrient concentrations are reduced towards the central part of the tumor  according to a normal distribution, depending on the distance from tumor periphery as demonstrated by the following equation:
where n0, g l0 represent the initial oxygen and glucose concentrations outside the tumor and D R is the distance of the current point from tumor periphery. The diffusion, as well as production, consumption and natural decay-rate parameters of the different chemicals (i.e. nutrients, ECM and MDEs) are taken from related research studies, such as [1, 3, 21, 32, 33, 35, 46–49]. Particularly, regarding the metabolic rates of the various cell phenotypes, hypoglycemic cells are assumed to consume nutrients at half the rate of proliferating cells. On the contrary, since the vast majority of hypoxic cells are considered to follow the glycolytic pathway, they are assumed to consume 10 times more glucose, but 1/5 less oxygen , compared to normoxic cells, the majority of which (approximately 60%) are considered to follow the aerobic metabolism in this model. Some parameters were hard to obtain and they were estimated according to related known values. The parameter values that were used in the equations of our model are summarized in Table A1 contained in an additional file (please see Additional file 1).
Boundary conditions for all tumor cell densities are determined by the first set of the following equation pairs. The same boundary conditions apply for MDEs concentration. Finally, nutrients and ECM are determined by the second equation pair.
Results and model validation
The compartmental equations of the proposed model along with its solutions have been numerically approximated by means of a finite differences method, namely Forward Euler. Several simulations have been performed using combinations of diffusion and proliferation rate constants, for the various evolution times of 1, 2, 3, 6, 9, 12, 15 and 18 months. More specifically, the results concern the combinations of [low diffusion, low proliferation], [high diffusion, low proliferation], [low diffusion, high proliferation], [high diffusion, high proliferation] and [medium diffusion, medium proliferation]. The model is initialized with an early-stage tumor of radius r0=1 m m, smaller than the minimum detectable size consisting only of proliferating cells, without any other assumptions on the initial radius of the various cell-population regions. Furthermore, the quiescent cells due to the lack of space have been included in the population of proliferative cells, since they still remain on normoxic and normoglycemic conditions. The simulation results of cell and nutrient densities with respect to the distance from the center of the tumor, for the combination of low diffusion (D=0.005 m m2/d a y) and high proliferation rate (ρ=0.04/d a y), corresponding to a primary high grade glioma, i.e. a glioblastoma from its first clinical detectability, are graphically illustrated in Figures 2, 3, 4 and 5 for the evolution times of 1, 6, 12 and 18 months. More results concerning the intermediate evolution times of 2, 3, 9 and 15 months are illustrated in an additional file (please see Additional file 2). This additional file also contains the graphical evolution of the low diffusion-high proliferation glioma evolution in the early stage of 1 day after its assumed detection, where all cancer cells remain proliferative since no hypoxia and/or hypoglycemia have appeared yet. Finally, all results concerning radius of regions for the various cell categories for the different times stages and the diffusion-proliferation pairs of interest are depicted in Table 1.
Figure 2 demonstrates the evolution of the primary glioblastoma 1 month after its assumed detection, when its radius has increased from r=1 m m to r=3.2 m m. It can be noticed that due to the high proliferation rate and the low invasion, cell density rapidly increases in the tumor center, which leads to the appearance of hypoxia (blue curve and region in Figures 2(a), 2(b)). The hypoxia is also verified in Figure 2(c) reflecting the oxygen inadequacy in the central part of the tumor. Concerning the other diffusion-proliferation combinations demonstrated in Table 1, the cells remain proliferative throughout the tumor, implying that there is still adequacy of oxygen and glucose. This proves that as D gets larger, the tumor is dominated by its diffuse extent that reduces cell density in its central part, thus it takes longer for the proliferative cells to fall into hypoxic and hypoglycemic conditions.
Simulation results of 6 months after the assumed tumor diagnosis (low diffusion, high proliferation) in Figure 3, demonstrate that the hypoxic cell density at the tumor center has increased. Additionally, a hypoglycemic cell population has appeared (magenta curve in Figure 3(a)) due to the glucose inadequacy, which is shown in Figure 3(d). Glucose concentration is further decreased since hypoxic cells consume much more glucose in order to survive. Apart from hypoxic and hypoglycemic cell populations, a necrotic core of approximately 3 mm, has been formed (red curve and region in Figures 3(a), 3(b)). This is followed by the hypoxic and hypoglycemic regions around it and the external proliferative region. Here the hypoglycemic zone is incorporated inside the hypoxic band. The hypoxic, hypoglycemic and necrotic populations have increased at the expense of proliferating cells (green curve in Figure 3(a)), which now constitute the minority for the tumor radius of r<5 m m. Moreover, at the same time hypoxic, hypoglycemic and necrotic regions have appeared for the high diffusion, high proliferation case (Table 1), where the tumor radius already exceeds the value of r=4 c m. This case could represent a progressive glioblastoma.
After 12 months (Figure 4), the hypoxic and hypoglycemic areas have grown while the separate zones have become more distinct, since the cells that compose each zone prevail on all other cells. It can be clearly observed that proliferating cells no longer exist in the central part of the tumor. Additionally, it is shown that the hypoglycemic region is still incorporated within the hypoxic one. Moreover, the glioma represented by medium diffusion, medium proliferation rates has progressed to grade IV, containing a necrotic central region, representing a secondary glioblastoma.
Finally, Figure 5 illustrates tumor growth after 18 months, where all the distinct regions have further increased in size and density from inside to outside. The formation of the extended necrotic core, surrounded by the hypoglycemic and hypoxic regions and the outward proliferating zone are clearly depicted.
Table 1 and Figure 6 demonstrate the relation between the diffusion coefficient and the tumor-growth rate. More specifically, it can be observed that in all cases the expansion rate is temporarily reduced after the first 3 months, while later it is increased proportionally to the diffusion-coefficient value concerning the medium and high diffusion cases, to finally reach a steady growth state. Regarding the low diffusion cases, the tumor growth rate is stabilized after its first reduction.
Model validation constitutes an essential part of mathematical modelling and it can be performed using either medical data stemming from clinical cases, or experimental model results available in the literature. The model in this paper is based on a virtual spherical tumor, so that its validation is more appropriate using experimental model results derived from similar existing research efforts. However, the same model can be initialized from a real tumor anatomy derived from MRI imaging. Nevertheless, the comparison with actual clinical data requires the use of nonhomogeneous media considerations, which is outside the scope of this paper. The results of the current study are compared to the corresponding outputs of the continuum model in  and the hybrid model in . In both cases, high correlation was revealed between simulated and experimental results.
First, the experimental results drawn from graphs presented in paper  are compared to simulations results of the proposed model. Grade association with size and time is depicted in Figure 7 for two different diffusion-proliferation combinations in comparison to respective results of . Specifically, Figure 7(a) demonstrates a primary glioblastoma (grade IV from its first detectability) and corresponds to low diffusion and high proliferation rate, while Figure 7(b) simulates a secondary glioblastoma, (progressing to grade IV from lower grade) and stands for medium diffusion and medium proliferation rate. For both cases, the initial tumor has a density equal to C0=104c e l l s/m m3 and a radius of r=1 m m (as considered in all the above-presented simulation results). This plot demonstrates the high agreement between the results of the proposed model and the experimental results of  concerning spatio-temporal evolution and grade association of glioma, even though the proposed model simulates an avascular tumor.
Table 2 depicts the time from the typical detectable size of r=1 c m to the typical fatal size of r=4 cm for all the five different combinations that were investigated in this paper. Additionally, the time intervals for the same D/ρ pairs are shown in this table for model in , which were approximated from the color-scale graph of predicted survival times for a spectrum of in silico patients. It can be observed that almost in all cases, the time that the tumor needs to reach the fatal radius of 4 c m according to the proposed model prediction is within the time-intervals of model  for the corresponding diffusion-proliferation tumor types.
In the model of  the simulated T2-detected tumor corresponds to a total cell density greater than the T2-threshold of detection, which according to this paper is approximated to 16% of the maximal tissue-carrying capacity. Following this rationale for comparison purposes, our model is also initialized with a spherical tumor of a density lower than the maximum carrying capacity, yet higher than the minimum detectable value of  (0.16·C m <C0<C m ), i.e. C0=105c e l l s/m m3. Figures 8(a), 8(c) and 8(e) illustrate the simulations of each grade of glioma, represented as a plot of cell density for the different cell populations (normoxic, hypoxic and necrotic), with respect to the distance from the center of the tumor. In Figure 8(e) the hypoglycemic cell population has been included in the normoxic cell population. Comparable curves derived from  corresponding to gliomas of similar grades and sizes are shown in Figures 8(b), 8(d) and 8(f). These figures demonstrate that a similar growth pattern is followed in both models, with respect to the different tumor regions in each glioma grade.
The deviations regarding cell-densities between the two models signify their differences and potential limitations. In particular, the proposed model, lacks a consideration of angiogenesis, but considers the impact of glucose gradual reduction on the cell proliferation rate. Moreover, the model of  is initialized with a tumor of a tenfold radius, namely the smallest T2 detectable size of r=1 c m, compared to the proposed model. This large initial size with high cell-density within its entire volume significantly differentiates the final cell densities, i.e. when the tumor has reached its fatal size.
In order to validate the proposed model for an initial tumor of a small radius r<1 m m, the results are compared to respective outcomes of model , which simulates avascular spherical tumor growth based on cell proliferation that depends on the concentrations of oxygen and glucose. Model results are demonstrated for a growing tumor far below its clinical detectable size. The initial concentration values of both nutrients outside the tumor are the same used in our model (n0=0.28·10−9m o l/m m3, g l0=16.5·10−9m o l/m m3). Since the initial tumor radius here is considered very small, r=0.05 m m, the concentrations of both nutrient inside the tumor are assumed to have their maximum values, equal to the reference concentration. The model in  does not include a diffusion factor for the tumor expansion, but its radius increases corresponding to cell proliferation. For comparison purposes, the diffusion coefficient value of our proposed model has been considered to be proportional to the mean area that each tumor cell occupies. The cell proliferation rate, simulates the tumor expansion due to the increasing space of proliferating cells. Additionally, the proliferation-rate constant in this case has been taken from  as ρ=0.12/d a y. These modified diffusion and proliferation rates and the tumor initialization for r=0.05 m m, lead to a similar tumor growth as the one in . The results are presented in Figure 9, where the tumor volume expansion for 28 days in comparison to the respective results of  are depicted. Concerning this growth, the radii of quiescent and necrotic regions in respect to the tumor radius are illustrated in Figure 10 and compared to the experimental results of . The quiescent radius corresponds to the hypoxic along with the hypoglycemic region of the proposed model. Figure 10 clearly demonstrates a very good agreement between the experimental data and the simulated growth curves, resulting from our model prediction.
Discussion and conclusions
Gliomas form an important class of brain cancer with high mortality rate. Mathematical models are increasingly used to understand and predict their behaviour. Identification of the avascular tumor dynamics constitutes the first crucial step towards the investigation of fully vascularized tumors. However, using current modelling techniques one must choose between simulating individual cell behavior and modelling tumors of clinically significant size. Both the continuum and discrete tumor models pose particular limitations. In order to overcome significant restrictions, a multi compartment continuum model appears to be effective in describing how subpopulations of various types of cells proliferate and diffuse, while it is computationally efficient. To this aim, a new continuum three-dimensional spherical mathematical model of avascular glioma growth in an isotropic and homogeneous medium is developed and proposed in this study. This model simulates malignant cell proliferation and invasion behavior, incorporating the interactions between the heterogeneous populations of four different glioma cell phenotypes, namely proliferative, hypoxic, hypoglycemic and necrotic, as well as their tissue microenvironment. Our model is capable of capturing intercellular interactions, spatial cell population heterogeneity and phenotype differentiation. It is validated through multiple simulations in order to determine sensitivity to changes in important model parameters, specifically, the proliferation and cell migration rates .
The different cell populations are implemented by means of distinct cellular compartments based on , which has been expanded towards the integration of tumor microenvironment with the inclusion of chemical compartments. Unlike the majority of related existing approaches that involve a single nutrient (oxygen) to guide tumor growth, our model incorporates the effects of the concentration changes of two nutrients (oxygen and glucose). The novelty of the proposed model lies in the fact that apart from the combined effect, each nutrient separately affects tumor cell proliferation and viability, while a new cell population is formed due to the lack of glucose, characterizing hypoglycemic cells by a different metabolic profile than hypoxic ones. Moreover, the effects of ECM, MDEs and the two vital nutrients, on cell survival, proliferation and invasion are simulated based on the equations of  (for oxygen, ECM and MDEs), which have been extended to include the different production and consumption rates corresponding to the distinct cell types. In essence, our model builds on relevant existing models by combining their capabilities and expanding their exposure and applicability to realistic glioma scenarios. All model variables, e.g. cell densities and concentrations of chemical ingredients are of continuous form described by diffusion-reaction equations. Several simulations have been performed using various diffusion and proliferation rate combinations corresponding to glioma reference values. The model results are presented for different evolution times in order to express the evolutionary characteristics of tumors. Furthermore, the model is validated through comparisons with specific experimental results of glioma models available in the literature, revealing high concordance of our model results with different aspects of glioma spatio-temporal evolution and tumor initialization phases.
Since the clinical cases of substantial significance are above the tumor size of 0.5−1 c m radius, our model is initialized with a much smaller tumor-radius r=1 m m containing only normoxic proliferative cells, as to consider the evolution of tumor composition. Our model predicts tumor radius and composition after its clinical diagnosis, while it also allows tumor initialization at a size far below the detectable in MRI images, so that it can approximate the actual tumor cell-densities since its infancy. As it was shown, the model enables identifying the size of the different cell regions, such as proliferative, hypoxic, hypoglycemic and necrotic, as well as their effect on the overall tumor growth. Therefore, the proposed model may prove useful for determining particular tumor growth parameters, which are hard to obtain from patient data but are crucial in the development of therapeutic strategies for cancer treatment. Moreover, the validation results demonstrate that our model can work effectively and provide the tumor-growth prediction initialized either from a very small or a larger tumor size. This enhances the applicability of the model irrespective of the stage of diagnosis.
In this study the effect of hypoxia as a hallmark of aggressive tumor behavior often met in glioblastomas that has been engaged in the model through nutrients availability/consumption has been associated with cell proliferation and invasion. According to medical experts this is of great significance in clinical practice, since it is associated with resistance to therapy, poorer survival and more malignant tumor phenotypes . Thus, one important utility of models is to simulate a quantitative link between tumor growth kinetics and the hypoxic burden of the tumor. To this aim, our mathematical model describes the spatial and temporal evolution of glioma in terms of concentration of malignant tumor cells. A major strength of this formalism is its potential prospective nature. Although it does not include the prediction of the tumor recurrence location, our method is potentially applicable to orient patient-specific definition of glioma margins.
Furthermore, since our model is able to detect cell densities of even below 1c e l l/m m3, it can effectively characterize the border line regions, i.e. the semi-cancerous regions, as well as the rest extra-tumoral suspicious area, which potentially consist of healthy tissue along with scattered cancer cells. In addition, the area that is occupied by an extra-tumoral edema is likewise considered as suspicious area by the radiologists, since it contains high concentrations of matrix-degradative enzymes (MDEs), which are known to degrade the extracellular material (ECM) and provide a favourable environment to cancer cells in order to migrate. Our model incorporates the variables of both the MDEs and ECM, thus it is able to detect those extra-tumoral suspicious areas containing high MDEs-concentrations and relatively low ECM-concentrations.
Despite its utility and effectiveness, the proposed model in its present form bears certain limitations. Tumor cells with various phenotypes initially develop asymmetric tumor morphologies, since they have different invasion rates and they migrate within an inhomogeneous medium; thus they are non-uniformly distributed. However, they eventually form a circular shape so that the spherical shape of the initial tumor used in our model constitutes a successive real tumor approximation. Furthermore, this model is homogeneous, i.e. it does not take into consideration brain heterogeneity regarding tumor cell invasion. As mentioned in , there are several indications that glioma diffuses faster in white than in grey matter. The cell diffusion coefficient considered in the current model is constant throughout the entire medium, while it should vary depending on the invaded tissue area. Along the same direction, model initialization is performed with a virtual tumor that lacks the specific spatial localization of a real tumor surrounded by particular brain anatomic structures (i.e. white and grey matter), where tumor cells would have different invasion rates. Thus, one limitation of the current model is the lack of support for medium heterogeneity, which however can be considered at a further stage of development. Moreover, a real tumor consists of at least three distinct regions of different texture properties (e.g. necrosis, proliferative, edema), where chemicals have different diffusion. The model presented does not consider variable diffusion properties for the nutrients due to intratumoral heterogeneity, since diffusion coefficient values are very difficult to obtain and additionally, such a differentiation would not offer a significant accession to model results.
Cell to ECM interaction (i.e. haptotaxis) and cell to cell adhesion play their specific roles in tumor growth and morphological dynamics. Cell to cell adhesion is important at the early stages of glioma development, while at later stages tumor growth is dominated by cell to ECM and nutrients interactions. Cell to cell adhesion could be incorporated in a macroscopic continuum approach, but has not been included in our model, since it concerns the cellular scale that is not explored in this work. Regarding haptotaxis, it has also been neglected as a distinct diffusion term, but it has been implicitly taken into account through random diffusion of tumor cells towards areas of lower cell concentration, where ECM concentration is higher. Similarly, chemotaxis has been implicitly incorporated through the factor of random diffusion. As another remark, the proposed model does not evaluate the impact of gene expression changes and mutations in cellular behavior, nor does it incorporate the effects of the epidermal growth factor receptor (EGFR) and the associated molecular pathways on tumor growth dynamics.
Overall, the major effort in this work is to demonstrate the results of the changes in the dominant cellular kinetics, i.e. diffusion and proliferation, considering the interactions between tumor cells and their microenvironment. Hence, the current effort forms a comprehensive modelling of the change in the tumor cellular proliferative and invasive phenotype due to the dynamic cell interactions with the microenvironment that contains the vital nutrients. As such, it may miss several key factors that determine cell phenotype, which nevertheless can be readily explored in future models.
Future efforts include the incorporation of several parameters aiming to overcome the limitations. A straight forward direction for model amendment is to address brain tissue heterogeneity and anisotropy issues, in order to take into account the different cell diffusion rates, depending on the invaded anatomic structure (i.e. white or grey matter). This can be easily performed by adopting the topography of imaging modalities (i.e. MRI) similar to other related studies . In the same direction we should consider the model initialization by a real tumor geometry extracted from corresponding medical images. Additionally, primary tumor formation in terms of proliferative, hypoxic and necrotic regions, derived from imaging texture, or biopsy data could also be integrated in this initialization. This primary intratumoral heterogeneity will allow the integration of variable nutrient diffusion properties. Subsequently, the implementation of additional flux factors, such as diffusion due to haptotaxis, chemotaxis and cell to cell adhesion should be addressed. The current model considers the effect of nutrient-concentration changes but it does not take into consideration the increase of waste concentrations, such as lactate and H+ ion, which act as cell proliferation inhibitors and threaten cell survival. The incorporation of those inhibitors is also considered essential in future developments. Furthermore, the effects of EGFR concentration and related pathways on tumor growth dynamics at a multicellular level will constitute a significant addition. Since the vascular morphology of gliomas has led to the hypothesis that the formation of new blood vessels is essential for tumor growth, future plans will target the investigation of angiogenesis within the proposed model.
According to the medical experts in our team, this model is successful in capturing avascular tumor growth observed clinically. It allows performing patient-specific simulation of different tumor evolution scenarios, towards reliable prognosis of glioma spatio-temporal progression. Since it explicitly incorporates the microenvironment interactions of proliferative, hypoxic and hypoglycemic glioma cells, along with the formation of necrosis, it allows evaluating clinically significant tumor sizes, while it is efficient in describing the dynamics of glioma tumors visualized with medical imaging. This mathematical model provides a means to simulate tumor development scenarios, which may lead to a better understanding of how altering fundamental parameters can influence brain tumor progression; thus it may constitute an important research tool for clinical assessment . In this regard, it can be exploited in cooperation with in-vivo and in-vitro models after establishing the initial diagnosis via a tumor biopsy or surgery. In-vitro invasion assays are important tools for investigating the tumor-matrix interactions and the effects of extracellular macromolecules on these interactions [51–53]. Since such studies are carried out through tightly-controlled experimental conditions, they allow the explicit determination of primary variables, e.g. nutrients initial concentrations [54, 55]. On the other hand, in vivo implanted xenographs (human brain tissue) in animal model variations can effectively investigate the reproducible cell migration and the tumor invasion into living non-neoplasmatic brain regions . Thus, both in vitro and in vivo studies of human tumor parts allow the examination of tumor behavior deriving important patient specific parameters, such as the invasion, the proliferation and the nutrients consumption rates . The proposed model can then be initialized, based on these parameters in order to attain a detailed and accurate short- and long-term patient specific prediction.
An important area of model exploration concerns the changes of tumor microenvironment in the context of a specific treatment. The treatment of tumors in the central nervous system represents a formidable challenge, further magnified by the fact that the brain is isolated by the blood brain barrier, rendering the delivery of high doses of chemotherapeutic agents, or gene vectors a very difficult task . This leads to unacceptable toxic systemic levels of drugs. The approaches proposing combinations of local tumor delivery, within the tumor mass and/or surrounding cavity, such as immune-stimulation and cytotoxic gene therapies, appear very promising in achieving the effective treatment of glioblastoma . Additionally, the introduction of new technologies, such as the microchip and convection enhanced drug delivery, will enable the local delivery of treatments, such as drugs and therapeutic gene vectors within the tumor mass, as well as the surrounding area, where the infiltrating tumor cells are localized [60, 61]. This provides the motivation for simulating and modelling the targeted cancer treatment scenarios. The multi-compartmental continuous nature of our model enables the integration of therapy-related issues in the form of either a simple factor, or a pharmakokinetic-pharmakodynamic model, which will investigate the beneficial impact of treatment in disease progression.
Kiran KL, Jayachandran D, Lakshminarayanan S: Mathematical modelling of avascular tumour growth based on diffusion of nutrients and its validation. Can J Chem Eng. 2009, 87: 732-740. 10.1002/cjce.20204.
Szeto MD, Chakraborty G, Hadley J: Quantitative metrics of net proliferation and invasion link biological aggressiveness assessed by MRI with hypoxia assessed by FMISO-PET in newly diagnosed glioblastomas. Cancer Res. 2009, 69 (10): 4502-4509. 10.1158/0008-5472.CAN-08-3884.
Swanson KR, Rockne RC, Claridge J, Chaplain MA, Alvord Jr EC, Anderson ARA: Quantifying the role of angiogenesis in malignant progression of gliomas: In Silico modeling integrates imaging and histology. Int Sys Tech: Math Onc Cancer Res. 2011, 71 (24): 7366-7375.
Swanson KR, Bridgea C, Murray JD, Alvord EC: Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. J Neurol Sc. 2003, 216: 1-10. 10.1016/j.jns.2003.06.001.
Clatz O, Sermesant M, Bondiau PY, Delingette H, Warfield SK, Malandain G, Ayache N: Realistic simulation of the 3D growth of brain tumors in MR images coupling diffusion with biomechanical deformation. IEEE Trans Med Imag. 2005, 24 (10): 1334-1346.
Jiang Y, Pjesivac-Grbovic J, Cantrell C, Freyer JP: A multiscale model for avascular tumor growth. Biophys J. 2005, 89: 3884-3894. 10.1529/biophysj.105.060640. [World Scientific Publishing Company]
Harpold HLP, Alvord EC, Swanson KR: The evolution of mathematical modeling of glioma proliferation and invasion. J Neuropathol Exp Neurol. 2007, 66: 1-9. 10.1097/nen.0b013e31802d9000.
Deisboeck TS, Stamatakos G: Multiscale Cancer Modeling. 2010, Boca Raton: Chapman and Hall/CRC Press
Roos T, Chapman SJ, Maini PK: Mathematical models of avascular tumor growth. SIAM. 2007, 49 (2): 179-208. 10.1137/S0036144504446291.
Wolfram S: Cellular Automata and Complexity: Collected Papers. 1994, USA: Addison-Wesley
Drasdo D: Coarse graining in simulated cell populations. ACS. 2005, 8 (2–3): 319-363.
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. 10.1006/jtbi.2000.2000.
Hatzikirou H, Deutsch A: Cellular automata as microscopic models of cell migration in heterogeneous environments. Curr Top Dev Biol. 2008, 81: 401-434.
Tanaka ML, Debinski W, Puri IK: Hybrid mathematical model of glioma progression. Cell Prolif. 2009, 42: 637-646. 10.1111/j.1365-2184.2009.00631.x.
Greenspan HP: Models for the growth of a solid tumor by diffusion. Stud Appl Math. 1972, 51: 317-340.
Roniotis A, Marias K, Sakkalis V, Tsibidis GD, Zervakis M: A complete mathematical study of a 3D model of heterogeneous and anisotropic glioma evolution. Proceedings of the 31st Annual International Conference of the IEEE Engineering in Medicine and Biology Society: 2–6 September 2009; Minneapolis; Minnesota; USA. 2009, New York: IEEE, 2807-2810.
Roniotis A, Manikis G, Sakkalis V, Zervakis M, Karatzanis I, Marias K: High grade glioma diffusive modeling using statistical tissue information and diffusion tensors extracted from atlases. IEEE Trans Inf Technol Biomed. 2012, 16 (2): 255-263.
Gatenby RA, Gawlinski ET: A reaction-diffusion model of cancer invasion. Cancer Res. 1996, 56: 5745-5753.
Giatili SG, Stamatakos GS: A detailed numerical treatment of the boundary conditions imposed by the skull on a diffusion–reaction model of glioma tumor growth. Clinical validation aspects. Appl Math Comp. 2012, 218: 8779-8799. 10.1016/j.amc.2012.02.036. [Elsevier]
Stein AM, Demuth T, Mobley D, Berens LM, Sander K: A mathematical model of glioblastoma tumor spheroid invasion in a three-dimensional In Vitro experiment. Biophys J. 2007, 92: 356-365. 10.1529/biophysj.106.093468.
Kim Y, Lawler S, Nowicki MO, Chiocca EA, Friedman A: A mathematical model for pattern formation of glioma cells outside the tumor spheroid core. J Theor Biol. 2009, 260: 359-371. 10.1016/j.jtbi.2009.06.025.
Nema S, Saxena VP: FEM Based study of concentration of proliferating cell in brain tumor. Appl Math. 2012, 3: 935-942. 10.4236/am.2012.38140.
May C, Kolokotroni E, Stamatakos G, Buechler SP: Coupling biomechanics to a cellular level model: an approach to patient-specific image driven multi-scale and multi-physics tumor simulation. Progr Biophys Mol Biol. 2011, 107: 193-199. 10.1016/j.pbiomolbio.2011.06.007.
Byrne H, Preziosi L: Modelling solid tumour growth using the theory of mixtures. Math Med Biol. 2003, 20 (4): 341-366. 10.1093/imammb/20.4.341.
Breward CJ, Byrne HM, Lewis CE: A multiphase model describing vascular tumour growth. Bull Math Biol. 2003, 65 (4): 609-640. 10.1016/S0092-8240(03)00027-2.
Bearer EL, Lowengrub JS, Frieboes HB, Chuang YL, Jin F, Wise SM, Ferrari M, Agus DB, Cristini V: Multiparameter computational modeling of tumor invasion. Cancer Res. 2009, 69 (10): 4493-4501. 10.1158/0008-5472.CAN-08-3834.
Preziosi L: Cancer Modelling and Simulation. 2003, Boca Raton: Chapman and Hall/CRC Press
Kozusko F, Bajzer Z: Combining Gompertzian growth and cell population dynamics. Math Biosci. 2003, 185 (2): 153-167. 10.1016/S0025-5564(03)00094-4.
Paltanea M, Tabirca S, Chen YJ, Tangney M: Cancer prediction modeling from volumetric data. Proceedings of the 11th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing: 26–29 September 2009; Timisoara; Romania. 2009, Washington: IEEE Computer Society, 162-167.
Simeoni M, Magni P, Cammia C, De Nicolao G, Croci V, Pesenti E, Germani M, Poggesi I, Rocchetti M: Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents. Cancer Res. 2004, 64: 1094-1101. 10.1158/0008-5472.CAN-03-2524.
Rejniak KA, Anderson AR: Hybrid models of tumor growth. Wiley Inter Rev. 2010, 3: 115-125.
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.
Anderson ARA, Chaplain MAJ, Newman EL, Steele RJC, Thompson AM: Mathematical modelling of tumour invasion and metastasis. J Theor Med. 2000, 2: 129-151. 10.1080/10273660008833042.
Gerlee P, Anderson ARA: A hybrid cellular automaton model of clonal evolution in cancer: The emergence of the glycolytic phenotype. J Theor Biol. 2008, 250: 705-722. 10.1016/j.jtbi.2007.10.038.
Jeon J, Quaranta V, Cummings PT: An off-lattice hybrid discrete-continuum model of tumor growth and invasion. Biophys J. 2010, 98: 37-47. 10.1016/j.bpj.2009.10.002.
Kim Y, Stolarska MA: A hybrid model for tumor spheroid growth In Vitro I: theoretical development and early results. Math Mod Meth Appl Sc. 2007, 17: 1773-1798. 10.1142/S0218202507002479. [World Scientific Publishing Company]
Perfahl H, Byrne HM, Chen T, Estrella V, Alarcon T, Lapin A, Gatenby RA, Gillies RJ, Lloyd MC, Maini PK, Reuss M, Owen MR: Multiscale modelling of vascular tumour growth in 3D: the roles of domain size and boundary conditions. PloS One. 2011, 6 (4): e14790-10.1371/journal.pone.0014790. [World Scientific Publishing Company]
Ribba B, Alcon T, Marron K, Maini PK, Agur Z: The use of hybrid cellular automaton models for improving cancer therapy. ACRI, LNCS. 2004, 3305: 444-453. [Springer]
Mallet DG, Pillis LGD: A cellular automata model of tumor-immune system interactions. J Theor Biol. 2006, 239 (3): 334-350. 10.1016/j.jtbi.2005.08.002. [World Scientific Publishing Company]
Gao X, Tangney M, Tabirca S: A multiscale model for hypoxia-induced avascular tumor growth. Proceedings of the International Conference on Bioscience, Biochemistry and Bioinformatics (IPCBEE): 26–28 February 2011; Singapore. 2011, Singapore: IACSIT Press, 53-58.
Cristini V, Lowengrub J: Multiscale Modeling of Cancer. 2010, Cambridge: Cambridge University Press
Martinez-Gonzalez A, Calvo GF, Romasanta LAP, Perez-Garcia VM: Hypoxic cell waves around necrotic cores in glioblastoma: a biomathematical model and its therapeutic implications. Bull Math Biol. 2012, 74: 2875-2896. 10.1007/s11538-012-9786-1. [World Scientific Publishing Company]
Li X, Cristini V, Nie Q, Lowengrub JS: Nonlinear three-dimensional simulation of solid tumor growth. Disc Cont Dyn Sys B. 2007, 7 (3): 581-604.
Casciari JJ, Sotirchos SV, Sutherland RM: Mathematical modelling of microenvironment and growth in EMT6/Ro multicellular tumour spheroids. Cell Prolif. 1992, 25: 1-22. 10.1111/j.1365-2184.1992.tb01433.x.
Smallbone K, Gavaghan DJ, Gatenby RA, Maini PK: The role of acidity in solid tumour growth and invasion. J Theor Biol. 2005, 235: 476-484. 10.1016/j.jtbi.2005.02.001.
Freyer JP, Tustanoff E, Franko AJ, Sutherland RM: In Situ oxygen consumption rates of cells in V-79 multicellular spheroids during growth. J Cell Physiol. 1984, 118: 53-61. 10.1002/jcp.1041180111.
Casciari JJ, Sotirchos SV, Sutherland RM: Glucose diffusivity in multicellular tumor spheroids. Cancer Res. 1988, 48: 3905-3909.
Freyer JP, Sutherland RM: A reduction in the In Situ rates of oxygen and glucose consumption of cells in EMT6/Ro spheroids during growth. J Cell Physiol. 1985, 124: 516-524. 10.1002/jcp.1041240323.
Mueller-Klieser WF, Sutherland RM: Oxygen consumption and oxygen diffusion properties of multicellular spheroids from two different cell lines. Adv Exp Med Biol. 1984, 180: 311-321. 10.1007/978-1-4684-4895-5_30.
Altman DA, Atkinson DS, Brat DJ: Best cases from the AFIP: glioblastoma multiforme. Radiograph Rev Publ Radiol Soc N A Inc. 2007, 27: 883-888.
Amberger VR, Avellana-Adalid V, Hensel T, Baron-van Evercooren A, Schwab ME: Oligodendrocyte-type 2 astrocyte progenitors use a metalloendoprotease to spread and migrate on CNS myelin. Euro J Neurosc. 1997, 9: 151-162. 10.1111/j.1460-9568.1997.tb01363.x.
Forsyth PA, Wong H, Laing TD, Rewcastle NB, Morris DG, Muzik H, Leco KJ, Johnston RN, Brasher PM, Sutherland G, Edwards DR: Gelatinase-A (MMP-2), gelatinase-B (MMP-9) and membrane type matrix metalloproteinase-1 (MT1-MMP) are involved in different aspects of the pathophysiology of malignant gliomas. Br J Cancer. 1999, 79: 1828-1835.
Amberger VR, Hensel T, Ogata N, Schwab ME: Spreading and migration of human glioma and rat C6 cells on central nervous system myelin in vitro is correlated with tumor malignancy and involves a metalloproteolytic activity. Cancer Res. 1998, 58: 149-158.
Friedlander DR, Zagzag D, Shiff B, Cohen H, Allen JC, Kelly PJ, Grumet M: Migration of brain tumor cells on extracellular matrix proteins in vitro correlates with tumor type and grade and involves alphaV and beta1 integrins. Cancer Res. 1996, 56: 1939-1947.
Friedl P, Noble PB, Walton PA, Laird DW, Chauvin PJ, Tabah RJ, Black M, Zanker KS: Migration of coordinated cell clusters in mesenchymal and epithelial cancer explants in vitro. Cancer Res. 1995, 55: 4557-4560.
Farin A, Suzuki SO, Weiker M, Goldman JE, Bruce JN, Canoll P: Transplanted glioma cells migrate and proliferate on host brain vasculature: a dynamic analysis. Glia. 2006, 53: 799-808. 10.1002/glia.20334.
Friedl P, Wolf K: Tumour-cell invasion and migration: diversity and escape mechanisms. Nat Rev Cancer. 2003, 3: 362-374. 10.1038/nrc1075.
Furnari FB, Fenton T, Bachoo RM, Mukasa A, Stommel JM, Stegh A, Hahn WC, Ligon KL, Louis DN, Brennan C, Chin L, DePinho RA, Cavenee WK: Malignant astrocytic glioma: genetics, biology, and paths to treatment. Genes Dev. 2007, 21: 2683-2710. 10.1101/gad.1596707.
Gajjar A, Chintagumpala M, Ashley D, Kellie S, Kun LE, Merchant TE, Woo S, Wheeler G, Ahern V, Krasin MJ, Fouladi M, Broniscer A, Krance R, Hale GA, Stewart CF, Dauser R, Sanford RA, Fuller C, Lau C, Boyett JM, Wallace D, Gilbertson RJ: Risk-adapted craniospinal radiotherapy followed by high-dose chemotherapy and stem-cell rescue in children with newly diagnosed medulloblastoma (St Jude Medulloblastoma-96): long-term results from a prospective, multicentre trial. Lanc Onc. 2006, 21: 813-820.
Fisher GH, Orsulic S, Holland E, Hively WP, Li Y, Lewis BC, Williams BO, Varmus HE: Development of a flexible and specific gene delivery system for production of murine tumor models. Oncogene. 1999, 18: 5253-5260. 10.1038/sj.onc.1203087.
Arora A, Scholar EM: Role of tyrosine kinase inhibitors in cancer therapy. J Pharmac Exp Therap. 2005, 315: 971-979. 10.1124/jpet.105.084145.
This work is supported by the YP.ER.THE.N. project, “Development of the research infrastructure of clinical computational tools and services for better diagnosis and prediction towards personalized treatment of oncological diseases” supported by the Interreg framework of Greek-Cyprus Co-operation, funded from European Union and National Sources.
The authors declare that they have no competing interests.
MP conceived the study, designed the model based on existing research efforts, developed its computational implementation, carried out the model simulations and drafted the manuscript. PK and XK participated in the design of the model providing guidelines from the clinical perspective, assisted in parameter estimation process, performed the evaluation of the model results, significantly contributed in the discussion section and suggested the future work extensions. MZ supervised the whole research, contributed to writing and improving the paper, suggested extensions and modifications and revised the manuscript critically. All authors read and approved the final manuscript.
Electronic supplementary material
Additional file 2: Figures of model results concerning intermediate tumor evolution times. Additional file 2 contains figures of model results (low diffusion-high proliferation glioma), concerning intermediate tumor evolution times, namely 2, 3, 9 and 15 months, as well as the simulation result of the same glioma in the early stage of 1 day after its assumed detection. (PDF 324 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Papadogiorgaki, M., Koliou, P., Kotsiakis, X. et al. Mathematical modelling of spatio-temporal glioma evolution. Theor Biol Med Model 10, 47 (2013). https://doi.org/10.1186/1742-4682-10-47
- Spatio-temporal evolution
- Cancer prediction
- Tumor cells
- Nutrients models
- Necrotic regions