Abstract
Background
Reactiondiffusion based models have been widely used in the literature for modeling the growth of solid tumors. Many of the current models treat both diffusion/consumption of nutrients and cell proliferation. The majority of these models use classical transport/mass conservation equations for describing the distribution of molecular species in tumor spheroids, and the Fick's law for describing the flux of uncharged molecules (i.e oxygen, glucose). Commonly, the equations for the cell movement and proliferation are first order differential equations describing the rate of change of the velocity of the cells with respect to the spatial coordinates as a function of the nutrient's gradient. Several modifications of these equations have been developed in the last decade to explicitly indicate that the tumor includes cells, interstitial fluids and extracellular matrix: these variants provided a model of tumor as a multiphase material with these as the different phases. Most of the current reactiondiffusion tumor models are deterministic and do not model the diffusion as a local statedependent process in a nonhomogeneous medium at the micro and mesoscale of the intra and intercellular processes, respectively. Furthermore, a stochastic reactiondiffusion model in which diffusive transport of the molecular species of nutrients and chemotherapy drugs as well as the interactions of the tumor cells with these species is a novel approach. The application of this approach to he scase of nonsmall cell lung cancer treated with gemcitabine is also novel.
Methods
We present a stochastic reactiondiffusion model of nonsmall cell lung cancer growth in the specification formalism of the tool Redi, we recently developed for simulating reactiondiffusion systems. We also describe how a spatial gradient of nutrients and oncological drugs affects the tumor progression. Our model is based on a generalization of the Fick's first diffusion law that allows to model diffusive transport in nonhomogeneous media. The diffusion coefficient is explicitly expressed as a function depending on the local conditions of the medium, such as the concentration of molecular species, the viscosity of the medium and the temperature. We incorporated this generalized law in a reactionbased stochastic simulation framework implementing an efficient version of Gillespie algorithm for modeling the dynamics of the interactions between tumor cell, nutrients and gemcitabine in a spatial domain expressing a nutrient and drug concentration gradient.
Results
Using the mathematical framework of model we simulated the spatial growth of a 2D spheroidal tumor model in response to a treatment with gemcitabine and a dynamic gradient of oxygen and glucose. The parameters of the model have been taken from recet literature and also inferred from real tumor shrinkage curves measured in patients suffering from nonsmall cell lung cancer. The simulations qualitatively reproduce the time evolution of the morphologies of these tumors as well as the morphological patterns follow the growth curves observed in patients.
Conclusions
s This model is able to reproduce the observed increment/decrement of tumor size in response to the pharmacological treatment with gemcitabine. The formal specification of the model in Redi can be easily extended in an incremental way to include other relevant biophysical processes, such as local extracellular matrix remodelling, active cell migration and traction, and reshaping of host tissue vasculature, in order to be even more relevant to support the experimental investigation of cancer.
Background
As the name indicates, reactiondiffusion models consist of two components. For systems of molecules and atoms, the first component is a set of biochemical reactions which produce, transform or remove chemical species. The second component is a mathematical description of the diffusion process. At molecular level, diffusion is due to the motion of the molecules in a medium. If solutions of different concentrations are brought into contact with each other, the solute molecules tend to flow from regions of higher concentration to regions of lower concentration, and there is ultimately an equalization of concentration. The conceptual framework of a microscale reactiondiffusion system can be also adopted to describe the phenomenology of cellular proliferation in tumor growth. Indeed, reactiondiffusion equations based models have been widely used in the literature for modeling the tumor growth. A comprehensive review of reactiondiffusion models and spatial dynamics of tumor growth can be found in [14], while specific literature about different variant of reactiondiffusion models of tumor growth can be found in [2,513]. Recently, there have been also interesting approaches to the adaptation of general reactiondiffusion models to the specific patient [14,15].
In this application domain, the reactiondiffusion models describe the evolution of the tumors via proliferation of malignant cells and their infiltration into the surrounding healthy tissue (see [16] for a review). The building block of these models is the reactiondiffusion type partial differential equations expressing the rate of change of the tumor cell density as sum of two terms. These terms correspond to the two phenomena described by the model: the diffusion term models, via the first Fick's law, the migration of tumor cells within the tissue and the reaction term, that is polynomial in the tumor cell density, models the proliferation of tumor cells [15]. Different reactiondiffusion models proposed in the literature mostly differ by the construction of the diffusion tensor in the mathematical expression of the Fick's law and the form of the proliferation term [15].
In this study we present a multiscale reactiondiffusion model linking the proliferation of malignant cells to (i) the upshot of the interactions between the oncological drug and the tumor cell, (ii) the availability and the rate of uptake of nutrient by the tumor cells, (iii) and, finally the availability and the rate of consumption of oxygen. Moreover, unlike the majority of the existing models of tumor growth, our model is stochastic, i.e. the interactions between tumor cells and drugs, as well as the events of uptake and consumption of nutrients and oxygen are stochastic Markov events. All the reactiondiffusion events are parallel and concurrent. The probability of a given event to be executed is proportional to the number of substrate molecules (for biochemical reactions) or to the number of cells (for interactions like cell proliferation and tumor growth). Recent studies [17] support this approach and show how the competitive intracellular reactions for the uptake and consumption of nutrients and oxygen are crucial in determining the tumor morphology.
We developed a generalization of the Fick's laws to model diffusion of drugs, nutrients and oxygen in the tissue, whereas we use the standard Fick's law to model the tumor cell proliferation and invasion following the gradient of nutrients and oxygen. Namely, the number of tumor cells and their spatial proliferation depend on the diffusion of nutrients, oxygen and drugs through the space and on the results of the interaction of these cell with the anticancer drug.
Before proceeding with a detailed explanation of our model of reactiondiffusion system, we briefly introduce the motivations and the guidelines of our work.
The tumor size provides a measure that is useful for describing the time course of tumor response to the chemotherapic treatment. However, tumor growth changes can be observed only through repeat followingup visits and may require sophisticated and expensive hardware and software imaging techniques especially for monitoring the size of in deepseated tumors. Due to this reason, the measurements of tumor progression in time and space have yet to gain wide application as an end point for drug effects modelling in clinical trials. The measurements of tumor size are still principally used for tumor stage categorization, whereas in the earlyphase clinical trials the measurements of changes in hematologic variables have been used as pharmacodynamic targets [18]. Therefore, a pharmacodynamic model that describes the interactions of tumor cells with the oncological drug and with nutrients, as well as the drug effects may have practical potential as a midterm end point for decision making about drug administration schedule and treatment duration. In this article, we focus on the modelling and simulation of the growth of nonsmall cell lung cancer treated with gemcitabine. This kind of cancer and its pharmacological treatment have been recently studied and new experimental data concerning both the measurements of the progression of tumor size [18] and the pharmacokinetics and pharmacodynamics of the gemcitabine [19] are now available and enable us to use of the computational models presented in this study.
In order to build an accurate and predictive model of tumor growth, the physical and biochemical nonhomogeneous environments in which the tumor arises and progresses, a generalization of the current mathematical formalization of reactiondiffusion systems is needed both at the microscale of the intracellular phenomena and at the mesoscale of intercellular and tissue processes. A preponderance of reactiondiffusion models of intra and inter cellular kinetics is usually performed on the premise that diffusion is so fast that all concentrations are maintained homogeneous in space. However, recent experimental data on intra and intercellular diffusion constants, indicate that this supposition is not necessarily valid even for small prokaryotic cells [20,21]. If the system is composed of a sufficiently large number of molecules, the concentration, i. e. the number of molecules per unit volume, can be represented as a continuous and differentiable variable of space and time. In this limit a reaction diffusion system can be modeled using differential equations. In an unstructured solvent, ideally behaving solutes (i.e. solutes for which solutesolute interaction are negligible) obey the Fick's law of diffusion. However in biological systems even for purely diffusive transport phenomena classical Fickian diffusion is, at best, a first approximation [22,23]. The intra and the intercellular media are not homogeneous mixtures of chemical species, but highly structured environments partitioned into compartments in which the distribution of the biomolecules could be nonhomogeneous. The description of diffusion processes in this environment has to start from a model in which diffusion coefficient contains its dependency on the local concentrations of the solutes and solvent.
In order to tackle this problem, this paper presents a new model of diffusion coefficient for a nonhomogeneous nonwellstirred reactiondiffusion system. In this model the diffusion coefficient explicitly depends on the local concentration, frictional coefficient of the particles of the systems, and of the temperature of the reaction environment. In turn, the rates of diffusion of the biochemical species are expressed in terms of these concentrationdependent diffusion coefficients. In this study the purely diffusive transport phenomena of noncharged particles, and, in particular, the case in which diffusion is driven by a chemical potential gradient in the x direction only (the generalization to the threedimensional case easily ensues) are considered. The generalization of the Fick's law introduced in this work, consists of five main steps: 1. calculation of the local virtual force F per molecules as the spatial derivative of the chemical potential; 2. calculation of the particles mean drift velocity in terms of F and the local frictional coefficient f; 3. estimation of the flux J as the product of the mean drift velocity and the local concentration; 4. definition of diffusion coefficients as function of local activity and frictional coefficients and concentration; and 5. calculation of diffusion rates as the negative first spatial derivative of the flux J.
The diffusion events are modeled as reaction events and the spatial domain of the reaction chamber is divided into N_{s }subvolumes (or meshes) of size l, that from now on will be called meshes. The movement of a molecule A from box i to box j is represented by the reaction , where A_{i }denotes the molecule A in the box i and A_{j }denotes the molecule A in the box j. The reactiondiffusion system is thus modeled as a pure reaction system in which the diffusion events are first order reactions whose rate coefficients k are expressed in terms of statedependent diffusion coefficients. The time evolution of the system is computed by a Gillespielike algorithm [24] that selects at each simulation step in each mesh the fastest reaction, compares the velocities of the N_{s }selected reactions and finally executes the reaction that is by far the fastest.
The paper outlines as follows. First we describe our generalization of the Fick's law, then we briefly describe how it can be incorporated in a stochastic simulation framework. Finally, we present a model of tumor growth and the simulation results.
Methods
We summarize here the main passages of the generalization of the Fick's first law. We refer the reader to Lecca et al.[25] for a more comprehensive description of the mathematical structures and passages.
A generalization of Fick's law for modeling diffusion
In a chemical system the driving force for diffusion of each species is the gradient of chemical potential µ of this species. The chemical potential of any particular chemical species i is
where is the standard chemical potential of the species i (i.e. the Gibbs energy of 1 mol of species i at a pressure of 1 bar), R = 8.314 J · K^{}^{1 }· mol^{1 }is the ideal gas constant, and T the absolute temperature.
The quantity a_{i }is called chemical activity of component i, and it is given by
where γ_{i }is the activity coefficient, and c^{0 }a reference concentration, which, for example, could be set equal to the initial concentration. The activity coefficients express the deviation of a solution from ideal thermodynamic behavior and, in general, they may depend on the concentration of all the solutes in the system. For an ideal solution, the limit of γ_{i}, which is recovered experimentally at high dilutions is γ_{i }= 1. If the concentration of species i varies from point to point in space, then so does the chemical potential. For simplicity, here the case in which there is only a chemical potential gradient in the x direction is taken into account. Chemical potential is the free energy per mole of substance, free energy is the negative of the work, W, which a system can perform, and work is connected to force F acting on the molecules by dW = Fdx. Therefore an inhomogeneous chemical potential is related to a virtual force per molecule of
where N_{A }= 6.022 × 10^{23} mol^{1} is Avogadro's number, k_{B }= 1.381 × 10^{23} J · K^{1 }is the Boltzmann constant, and the sum is taken over all species in the system other than the solvent. This force is balanced by the drag force experienced by the solute (F_{drag}, _{i}) as it moves through the solvent. Drag forces are proportional to speed. If the speed of the solute is not too high in such a way that the solvent does not exhibit turbulence, the drag force can be written as follows
where f_{i }α c_{i }is the frictional coefficient, and v_{i }is the mean drift speed.
Moreover, if the solvent is not turbulent, the flux, defined as the number of moles of solute which pass through a small surface per unit time per unit area, can be approximated as in the following
i.e. the number of molecules per unit volume multiplied by the linear distance traveled per unit time.
Since the virtual force on the solute is balanced by the drag force (i. e. F_{drag,i }= F_{i}), the following expression for the mean drift velocity is obtained
so that Eq. (5) becomes
where
are the diffusion coefficients. The Eq. (7) states that, in general, the flux of one species depends on the gradients of all the others, and not only on its own gradient. However, here it is supposed that the chemical activity a_{i }depends only weakly on the concentrations of the other solutes, i.e. it is assumed that D_{ij }≈ 0 for i ≠ j and that Fick's laws still holds [26]. Let D_{i }denote D_{ii}. It is still generally the case that D_{i }depends on c_{i }in sufficiently concentrated solutions since γ_{i }(and thus a_{i}) has a non trivial dependence on c_{i }[26]. There is only one very special case, namely that of an ideal solution with γ_{i }= 1, in which the diffusion coefficient, D_{i }= k_{B}T/f_{i}, is constant. In order to find an analytic expression for the diffusion coefficient, D_{i}, in terms of the concentration, c_{i}, let us consider that the rate of change of concentration of the substance i due to diffusion is given by
Substituting Eq. (7) into Eq. (6), and then substituting the obtained expression for J_{i }into Eq. (8), give
so that
Let c_{i},_{k }denote the concentration of a substance i at coordinate x_{k}, and l = x_{k } x_{k}_{1 }the distance between adjacent mesh points. The derivative of c_{i }with respect to x calculated at is
By using Eq. (9) into Eq. (6) the diffusive flux of species i midway between the mesh points, is obtained:
where is the diffusion coefficient midway between the mesh points.
The rate of diffusion of substance i at mesh point k is
and thence
To determine completely the righthand side of Eq. (11) is now necessary to find an expression for the activity coefficient, γ_{i}, and the frictional coefficient, f_{i}, contained in the expression of the diffusion coefficient. In fact, by substituting Eq. (2) into Eq. (7) we obtain the diffusion coefficient in terms of activity coefficients γ_{i}:
where the frictional coefficient is assumed to be linearly dependent on the concentration of the solute like in sedimentation processes, i.e. in a mesh k, fi, k is
where k_{f }is an empirical constant, whose value can be derived from the ratio R = k_{f}/[η ]: Accordingly to the MarkHouwink equation [27], [η ] = kM^{α }is the intrinsic viscosity coefficient, α is related to the shape of the molecules of the solvent, and M is the molecular mass of the solute. If the molecules are spherical, the intrinsic viscosity is independent of the size of the molecules, so that α = 0. All globular proteins, regardless of their size, have essentially the same [η ]. If a protein is elongated, its molecules are more effective in increasing the viscosity and [η ] is larger. Values of 1.3 or higher are frequently obtained for molecules that exist in solution as extended chains. Longchain molecules that are coiled in solution give intermediate values of is α, frequently in the range from 0.6 to 0.75 [27]. For globular macromolecules, R has a value in the range of 1.4  1.7, with lower values for more asymmetric particles [28].
Although Eq. (13) is a simplified linear model of the frictional forces, it works quite well in many case studies and can be easily extended to treat more complex frictional effects (see [25,29]).
Let us focus now on calculation of the activity coefficients: a way to estimate the frictional coefficients will be presented in the next subsection. By using the subscript '1' to denote the solvent and '2' to denote the solute, it can be written that
where γ_{2} is the activity coefficient of the solute and c_{2} is the concentration of the solute. Differentiating with respect to c_{2} gives
The chemical potential of the solvent is related to the osmotic pressure (II) by
where V_{1} is the partial molar volume of the solvent and its standard chemical potential. Assuming V_{1} to be constant [30] and Differentiating μ_{1} with respect to c_{2} yield
Now, from the GibbsDuhem relation [31], the derivative of the chemical potential of the solute with respect to the solute concentration is
where M is molecular mass of the solute and is the partial molar volume of the solute divided by its molecular mass. The concentration dependence of osmotic pressure is usually written as
where B is the second virial coefficient, and thence the derivative of II with respect to the solute concentration is
Introducing Eq. (20) into Eq. (18) gives
From Eq. (15) and Eq. (21) it can be obtained that
so that
On the grounds that [32], solving the integral yields
The molecular mass M_{i,k }of the species i in the mesh k can be expressed as the ratio between the mass, m_{i,k}, of the species i in that mesh and the Avogadro number M_{i,k }= m_{i,k }/N_{A}. If p_{i }is the mass of a molecule of species i and c_{i,k }· l is the number of molecules of species i in the mesh k, then the molecular mass of the solute of species i in the mesh k is given by
Substituting the expression in Eq. (22) gives, for the activity coefficient of the solute of species i in the mesh k (γ_{i;k}), the following equation
Therefore, substituting Eq. (13) and Eq. (24) into Eq. (12), we obtain the following expression for a time and spacedependent diffusion coefficient
We finally estimated in the following way the second virial coefficient B. The statistical mechanics definition of the second virial coefficient is as follows
where u(r), which is given in Eq. (27), is the interaction free energy between two molecules, r is the intermolecular centercenter distance, k_{B }is the Boltzman constant, and T the temperature. In this work, it is assumed that u(r) is the LennardJones pair (12,6)potential (Eq. 27), that captures the attractive nature of the Van der Waals interactions and the very shortrange Born repulsion due to the overlap of the electron clouds:
By expanding the term exp into an infinite series, the Eq. (26) becomes
where T* ≡ 4/(k_{B}T ) and thus
An estimate of B is given by truncating the infinite series of functions to j = 4, since simulation results not shown here prove that taking into account the additional terms, obtained for j >4, does not significantly influence the simulation results [25].
Modelling the stochasticity
Both diffusion and reactions are modelled as reaction events whose dynamics is driven by the First Reaction Method of the Gillespie algorithm.
In particular, the diffusion events are modeled as firstorder reactions. namely, the movement of a molecule A from box i to box j is represented by the reaction , where A_{i }denotes the molecule A in the mesh i and A_{j }denotes the molecule A in the mesh j. In this way, the reactiondiffusion system is modeled as a pure reaction system.
The space domain of the system is divided into a number N_{s }of squared meshes of size l. The time evolution of the system is computed by the First reaction Method of Gillespie [24] that at each simulation step selects in each mesh the fastest reaction, compares the velocities of the N_{s }selected reactions and finally executes the reaction that is by far the fastest. The fastest reaction is defined as the reaction whose waiting time is the smallest.
The time at which each event is expected to occur is a random variable extracted by an exponential distribution [24]. Let R_{i }be the ith reaction channel expressed as
where l_{ij }is the stoichiometric coefficient of reactant S_{p}(_{i,j}), p(i, j) is the index that selects the species that participates in R_{i}, L_{i }is the number of reactants in R_{i}, and r_{i }is the rate constant. If the fundamental hypothesis of stochastic chemical kinetics [24] holds within a box, both diffusion and reaction events waiting times are distributed according to a negative exponential distribution, so that a typical time step has size
where R is the number of events. It is given by R = R_{diff }+ R_{react}, where R_{diff }is the number of possible diffusion events and R_{react }is the number of reaction events [33]. The diffusion and reaction propensities are given by the following expressions, respectively
where and are the number of chemical species that diffuse and the number of those the undergo to reactions, respectively. In general , since some species both diffuse and react. In Eq. (30), r_{i}^{(diff) }is the kinetic rate associated to the jumps between neighboring subvolumes, whereas in Eq. (31), r_{i}^{(react) }is the stochastic rate constants of the ith reaction.
From Eq. (11), the rate coefficient of the first order reaction representing a diffusion event is recognized to be as follows
Results
Here we describe the model of tumor growth implemented with the toll Redi. Redi is a software prototype Redi [25] that has been recently developed by Lecca et al. [25,34] to simulate the mathematical model of stochastic reactiondiffusion system that we have described in the previous section. We refer the reader to the references [25,3436] for technical details about the implementation and for a user manual of this software.
Model of tumor growth
The reactiondiffusion system modelling tumor growth involves four components: (i) the drug, gemcitabine, (ii) the tumor cell, (iii) oxygen, and (iv) glucose. The reaction events we modeled are the following:
R1. gemcitabine injection;
R2. gemcitabine diffusion;
R3. gemcitabine degradation (rate parameter k_{1});
R4. effective interaction of gemcitabine and death of tumor cell (rate parameter k_{2});
R5. ineffective interaction of gemcitabine: the tumor cell survives to the drug (rate parameter k_{3});
R6. tumor growth (rate parameter k_{4});
R7. glucose uptake (rate parameter k_{5});
R8. oxygen uptake (rate parameter k_{6});
R9. glucose diffusion;
R10. oxygen diffusion;
R11. tumor turnover (rate parameter k_{7}).
With regard to the dosage schedule of gemcitabine (event R1), we simulated the administration regime proposed by Tham et al. [18] and Soo et al. [37], i.e. gemcitabine was infused at a fixed dose rate of 1,000 mg/m^{2 }over 30 min on day 1 and 8 every three weeks. The diffusion coefficient of gemcitabine is automatically calculated by Redi as a function of space and time according to the formula (25). The efficacy of the gemcitabine (k_{2}), the rate constant for resistance appearance (k_{3}), and the tumor growth rate (k_{4}) have been inferred with KInfer (a maximum likelihood parameter estimator recently developed by Lecca et al. [38,39]) from the tumor size shrinkage curves observed in 56 patients treated with gemcitabine [18]. The 56 patients have been categorized by their age, sex and smoke history, and in our simulations we considered the average values of k_{2}, k_{3}, and k_{4} (see in Tables 1 and 2 the average values and the standard deviations of these three parameters).
Table 1. values of parameters and variables in the three models.
Table 2. categorization of patients and average values of gemcitabine efficacy.
Finally, the parameters of reactions R7 (k_{5}) and R8 (k_{6}) have been taken from [40,41] and [42,43], respectively.
According to reactions R9 and R10, tissues receive glucose and oxygen perfusing through the vessel wall and diffusing in the extracellular space.
Finally, the event R11 (tumor turnover) refers to the replacement of old tumor cells with newly generated ones from the existing ones. Tumor turnover is measured in units of sec · mm, and its value (k_{7}) for nonsmall lung cancer cell has been measured by Tham et al. [18].
We simulated the morphological changes of an irregular 2D spheroidal tumor having an initial diameter of 3 mm. The size of the computational space is 40 × 40 squared meshes each of which represents a squared portion of tissue having a side of 1 mm. If we assume that the cells have a diameter of 50 μm, a mesh of 1 mm^{2 }is approximately occupied by 2457 cells.
We assumed that the initial spatial distribution of gemcitabine exhibits a gradient pointing outside the tumor. Furthermore we assumed that the tumor as well as the surrounding healthy tissue are crossed by a vascular network of capillaries separated by a distance of 80 μm each from the other (see Figure 1). The glucose and oxygen are supplied by the capillaries and they diffuse through the tumor tissue with a rate of diffusion defined by Eq. (11). Their diffusion coefficients are calculated with the formula (25).
Figure 1. A simple model of the vascular network innervating the tumor. The distance between capillaries is 80 μm.
All the events R1R11 are modelled as reactionevents as in the following
R1. gemcitabine injection: zerothorder reaction ∅ → gemcitabine
R2. gemcitabine diffusion: first order reaction modelling the movement of gemcitabine molecules from mesh k to the mesh k': gemcitabine_{k } gemcitabine_{k', }where defined by Eq. (32);
R3. gemcitabine degradation (rate parameter k_{1});
R4. effective interaction of gemcitabine and death of tumor cell (rate parameter k_{2});
R5. ineffective interaction of gemcitabine: the tumor cell survives to the drug (rate parameter k_{3});
R6. tumor growth (rate parameter k_{4});
R7. glucose uptake (rate parameter k_{5}): zerothorder reaction glucose;
R8. oxygen uptake (rate parameter k_{6}): : zerothorder reaction ; oxygen;
R9. glucose diffusion: first order reaction modelling the movement of gemcitabine molecules from mesh k to the mesh k': glucose _{k } glucose_{k'}, where r^{(diff)} is defined by Eq. (32);
R10. oxygen diffusion: first order reaction modelling the movement of gemcitabine molecules from mesh k to the mesh k': oxygen_{k } oxygen_{k'}, where r^{(diff)} is defined by Eq. (32);
R11. tumor turnover (rate parameter k_{7}): a zerothorder reaction describes the generation of new tumor cells from the existing ones; a subsequent first order reaction specifies the replacement of old tumor cells with newly generated ones.
We developed three models, by changing the values of the glucose uptake, the dose of gemcitabine and the number of tumor cell per mesh. In vivo and in vitro experiments carried on in the last decade highlight the crucial role of these variables in governing the dynamics of tumor growth. Some reference experimental studies in this regards are reported in [4447]. Table 1 reports the initial values of the variables as well as the values of the parameters in the three models. The Model 1 is the reference model whose parameter have physiological values and the dose of gemcitabine is the one usually administered in vivo and in vitro clinical trials. In Model 2 we decreased the number of tumor cells per mesh and the rate of glucose uptake (100 cells per mesh instead of 2547 cells per mesh). In Model 3, we increased the dose of gemcitabine (10 μg instead of 1 μg). The orders of magnitude of changes of the parameter values in Model 2 and Model 3 are those that cause a significant change in the tumor growth progression.
The average of 100 simulations for each model (Model 1, 2 and 3) is showed in Figures 2, 3, and 4 respectively. Each subfigure is a screenshot of the state of the tumor size recorded each 10 weeks. Blue regions corresponds to areas occupied by more that 2000 tumor cells, yellow regions corresponds to areas of tissue with a number of tumor cell between 100 and 2000, and orange regions are those occupied by less that 100 tumor cells. The simulation of Model 1 in Figure 2 shows a progressive quasilinear growth of the tumor size at a rate of about 0.5 mm per week. The simulation of Model 2 shows that with a lower rate of glucose uptake the growth of the tumor is slowed down and the borders of the tumor ellipsoid are strongly irregular. Moreover, groups of cells originally belonging to the borders of the tumor proliferate in filaments in healthy parts of the tissue. The action of gemcitabine breaks these filaments but some tumor cells of the filaments still persist in isolated groups infiltrated into the healthy tissue.
Figure 2. Simulation of Model 1. The time unit is the week. The time separating a screenshot from the previous one is 10 weeks. The parameters of the model are listed in Table 1. The longitudinal initial size of the tumor spheroid is 3 mm. Screenshot number "0" is the state of the tumor after 10 weeks of treatment. In the spatial domain of tumor lesion each mesh hosts only tumor cells2. Blue regions are those occupied by more that 2000 tumor cells, yellow regions corresponds to areas of tissue with a number of tumor cells between 100 and 2000, and orange regions are those occupied by less that 100 tumor cells. The extension of the tumor increases linearly in time.
Figure 3. Simulation of Model 2. The time unit is the week. The time separating a screenshot from the previous one is 10 weeks. The parameters of the model are listed in Table 1. The initial diameter of the tumor ellipsoid is 3 mm. Screenshot number "0" is the state of the tumor after 10 weeks of treatment. In this model, in the spatial domain of tumor lesion, a mesh hosts both healthy and tumor cells. The number of tumor cells is 100 per mesh and the rate of glucose uptake is two order of magnitude smaller than in rate of glucose uptake in Model 1. As in Figure 1, blue regions are those occupied by more that 2000 tumor cells, yellow regions corresponds to areas of tissue with a number of tumor cells between 100 and 2000, and orange regions are those occupied by less that 100 tumor cells. The size of the tumor is approximately constant, but filaments of tumor cells propagate from the border of the tumor.
Figure 4. Simulation of Model 3. The time unit is the week. The parameters of the model are listed in Table 1. The initial diameter of the tumor ellipsoid is 3 mm. Screenshot number "0" is the state of the tumor after 10 weeks of treatment. In this model, in the spatial domain of tumor lesion, a mesh hosts both healthy and tumor cells. The number of tumor cells per mesh is 2457 as in Model 1 and the rate of glucose uptake is two order of magnitude smaller than the rate of glucose uptake in Model 1. As in Figure 1 and Figure 2, blue regions are those occupied by more that 2000 tumor cells, yellow regions corresponds to areas of tissue with a number of tumor cells between 100 and 2000, and orange regions are those occupied by less that 100 tumor cells. The size of the tumor is approximately constant, but filaments of tumor cells propagate from the border of the tumor, but are disrupted by the action of gemcitabine.
Figure 3 shows the simulation of Model 3, where we increased the dose of gemcitabine by a factor 10 to explore the behavior of the tumor mass for the extreme limit of an unrealistic dosage configuration. The rate of glucose uptake is the same as in Model 1. We observed a behavior similar to the one observed in the simulation of Model 1.
As expected, from these simulations we deduced that the effect of gemcitabine is stronger (i) at the early stage of the tumor (i.e. when the number of tumor cells is still low) and the rate of glucose uptake is also low (Model 2), or (ii) if the dose if greater them 1,000 mg.
At the best of our knowledge our study is the first attempt to model and simulate the tumor growth of nonsmall cell lung cancer in space and time. We validated our models by comparing the time behavior of the longitudinal size of the tumor ellipsoid with the theoretical and experimental results of Tham et al. [18]: we found a good agreement between the experimental data and the predictions of Model 2 and Model 3, such as the dosage of body gemcitabine necessary to slow down and arrest the growth of tumor (about 10 μg of body dose as in [18]), and the rate of tumor growth in the case of an insufficient amount of drug (between 0.5 an 1 mm per week as in the graph reported in [18]). Moreover these models confirm the correlation between glucose uptake and pharmacological treatment as reported by Duhaylongsod et al. in [44]. Namely, comparing the results of Model 2 and Model 3 with those of Model 1 we confirmed the necessity of a higher dose of gemcitabine or conversely of the reduction of the glucose uptake [18,48] for obtaining a significant increment of tumor shrinkage.
Conclusions
We have presented a computational framework for modeling and simulating the spatial dynamics of the diffusion of biological entities at micro and mesoscale in a nonhomogeneous medium. We use these mathematical and computational structure to model and simulate a nonsmall cell lung cancer treated with gemcitabine. The drug efficacy and the rate constant of resistance appearance have been estimated from real tumor growth curves recorded in 56 patients. The other parameters have been obtained from the literature reporting the in vitro experiments of the last decade. We explored the behavior of the model under different conditions concerning the rate of glucose uptake, the number of tumor cells and the dose of gemcitabine. The proposed models reproduce the expected tumor growth rate at the optimal body concentration of gemcitabine and confirm the correlation between glucose uptake and the response to the chemotherapy. At the best of our knowledge, this study is the first attempt to build a reactiondiffusion model of nonsmall cell lung cancer by integrating data from in vivo experiments and by inferring kinetic parameters from the tumor shrinkage curves of patients with the purpose to provide in silicogenerated dynamical images of the morphology of this kind of tumor.
Nonlinear models of cancer growth are needed to understand the phenomenon of realistic cancer growth. Simulations of such models conducted to determine the patterns of cancer growth and cancer response to drug and nutrient supply could support the design of the administration schedule and the duration of the therapy. Moreover, a computational model of a reactiondiffusion system taking into account the stochasticity of the interaction between drugs and tumor cells as well as the nonhomogeneity of the intra and intercellular medium may be a contribution toward this direction. Further extensions of this study are in progress and consider the opportunity to include immunological and angiogenic factors and interactions to make the current models more accurate, realistic and of greater medical interest.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
Each author contributed to this work in compliance with his/her expertise field. Paola Lecca developed the mathematical model of the stochastic dynamics of a nonhomogeneous reactiondiffusion systems. Paola Lecca also designed the in silico experiments and wrote the scripts for the simulations of the tumor growth with Redi. Daniele Morpurgo contributed to the literature referencing of the study, to the calibration of the models, and to the analysis and validation of the results.
Acknowledgements
The authors would like to thank Corrado Priami of the Microsoft Research  University of Trento Centre for Computational and Systems Biology of Trento (Italy), and University of Trento, for his valuable suggestions, R. A. Soo of the Department of HematologyOncology, National University Hospital of Singapore, and Lai San Tham of the LillyNUS Centre for Clinical Pharmacology for their indications.
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 14, 2012: Selected articles from Research from the Eleventh International Workshop on Network Tools and Applications in Biology (NETTAB 2011). The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/13/S14
References

Araujo RP, McElwain DLS: A history of the study of solid tumour growth: The contribution of mathematical modelling.
Bulletin of Mathematical Biology 2004, 66(5):10391091. PubMed Abstract  Publisher Full Text

Friedman A: A Hierarchy of Cancer Models and their Mathematical Challenges.
Discrete and Continuous Dynamical Systems  Series B, Mathematical Models in Cancer 2004, 4:147159.

Habiba S, MolinaParisb C, Deisboeckd TS: Compelx dynamics of tumors: modeling an emerging brain tumor with coupled reactiondiffusion equations.
Physica A: Statistical Mechanics and its Applications 2003, 327(34):501524. Publisher Full Text

Roose T, Chapman SJ, Maini PK: Mathematical Models of Avascular Tumor Growth.
SIAM REVIEW Society for Industrial and Applied Mathematics Research 2007, 49:179208.

Chaplain MAJ, Ganesh M, Graham IG: Spatiotemporal pattern formation on spherical surfaces: numerical simulation and application to solid tumour growth.
Journal of Mathematical Biology 2001, 42(5):387423. PubMed Abstract  Publisher Full Text

Clatz O, Sermesant M, ad H Delingette PYB, Warfield SK, Malandain G, Ayache N: Realistic simulation of the 3D growth of brain tumors in MR images coupling diffusion with biomechanical deformation.

Frieboes HB, Zheng X, Sun CH, Tromberg B, Gatenby R, Cristini V: An Integrated Computational/Experimental Model of Tumor Invasion.

Gatenby RA, Gawlinski ET: A ReactionDiffusion Model of Cancer Invasion.
Cancer Research 1996, 56:5745. PubMed Abstract  Publisher Full Text

Jiang Y, PjesivacGrbovic J, Cantrell C, Freyer JP: A Multiscale Model for Avascular Tumor Growth.
Biophysical Journal 2005, 89(6):38843894. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Khain E, Sander LM: Dynamics and pattern formation in invasive tumor growth.
Phys Rev Lett 2006, 96:188103. PubMed Abstract  Publisher Full Text

Konukoglu E, Sermesant M, Clatz O, Peyrat JM, Delingette H, Ayache N: A Recursive Anisotropic Fast Marching Approach to Reaction Diffusion Equation: Application to Tumor Growth Modeling.
Information Processing in Medical Imaging, Lecture Notes in Computer Science 2007, 20:687699.

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

Tracqui P: Biophysical models of tumour growth.
Rep Prog Phys 2009, 72:056701. Publisher Full Text

Hinow P, Gerlee P, McCawley LJ, Quaranta V, Ciobanu M, Wang S, Graham JM, Ayati BP, Claridge J, Swanson KR, Loveless M, Anderson ARA: A spatial model of tumorhost interaction: application of chemotherapy.

Konukoglu E, Clatz O, Delingette H, Ayache N: Multiscale Cancer Modeling , CRC Press. Chapman and Hall/CRC Mathematical and Computational Biology 2010 chap. Personalization of ReactionDiffusion Tumor Growth Models in MR Images: Application to Brain Gliomas Characterization and Radiotherapy Planning;

Byrne HM: Modelling avascular tumour growth. Cancer Modelling and Simulation. Mathematical Biology and Medicine Series, 3, Preziosi ed., London: Chapman and Hall/CRC; 2003.

Ferreira SC, Martins ML, Vilela MJ: Reactiondiffusion model for the growth of avascular tumor.

Tham LS, Wang L, Soo RA, Lee SC, Lee HS, Goh BC, Holford NH: A pharmacodynamic model for the time course of tumor shrinkage by gemcitabine + carboplatin in nonsmall cell lung cancer patients.
Clinical Cancer Res 2008, 14(13):42134218. Publisher Full Text

Veltkamp SA, Pluim D, van Eijndhoven MA, Bolijn MJ, Ong FH, Govindarajan R, Unadkat JD, Beijnen JH, Schellens JHM: New insights into the pharmacology and cytotoxicity of gemcitabine and 2',2'difluorodeoxyuridine.
Mol Cancer Ther 2008, 7(8):24152425. PubMed Abstract  Publisher Full Text

Elf J, Doncic A, Ehrenberg M: Mesoscopic reactiondiffusion in intracellular signaling.
Fluctuation and noise in biological, biophysical and biomedical systems. Procs. of SPIE 2003, 5110:114124. Publisher Full Text

Elowitz MB, M G Surette PEW, Stock JB, Leibler S: Protein mobility in the cytoplasm of Escherichia coli.
J Bacteriol 1999, 181:197203. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Agutter PS, Wheatley D: Random walks and cell size.
BioEssays 2000, 22(11):10181023. PubMed Abstract  Publisher Full Text

Agutter P, Malone P, Wheatley D: Intracellular transport mechanisms: a critique of diffusion theory.
J. Theor. Biol 1995, 176:261272. PubMed Abstract  Publisher Full Text

Gillespie D: Exact stochastic simulation of coupled chemical reactions.
Journal of Physical Chemistry 1977, 81(25):23402361. Publisher Full Text

Lecca P, Ihekwaba AEC, Dematté L, Priami C: Stochastic simulation of the spatiotemporal dynamics of reactiondiffusion systems: the case for the bicoid gradient.

Roussel CJ, Roussel MR: Reactiondiffusion models of development with statedependent chemical diffusion coefficients.
Progress in Biophysics & Molecular Biology 2004, 86:113160. PubMed Abstract  Publisher Full Text

Laidler K, Meiser J, Sanctuary B: Physical Chemistry. Houghton Mifflin Company; 2003.

Harding S, Johnson P: The concentration dependence of macromolecular parameters.
Biochemical Journal 1985, 231:543547. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Solovyova A, Schuck P, Costenaro L, Ebel C: Non ideality of sedimantation velocity of halophilic malate dehydrogenase in complex solvent.
Biophysical Journal 2001, 81:18681880. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chang R: Physical Chemistry for the Biosciences. 3rd edition. Wiley; 2003.

Moran MJ, Shapiro HN: Fundamentals of Engineering Thermodynamics. 3rd edition. Wiley; 2003.

Tombs MP, Peacocke AR: The Osmotic Pressure of Biological Macromolecules. Monograph on Physical Biochemistry, Oxford University Press; 1975.

Bernstein D: Simulating mesoscopic reactiondiffusion systems using the Gillespie algorithm.

Lecca P, Dematté L: Stochastic simulation of reactiondiffusion systems.
Int J of Medical and Biological Engineering 2008, 1(4):211231.

Redi web page [http://www.cosbi.eu/index.php/research/prototypes/redi] webcite

Lecca P, nd M Lecca LD, Priami C: Stochastic modelling of diffusion systems: video image simulation of tubulin diffusion in cytoplasm: a case study.
II Eccomas thematic conference on computational vision and medical image processing 2009.

Soo RA, Wang LZ, Tham LS: A multicentre randomised phase II study of carboplatin in combination with gemcitabine at standard rate or fixed dose rate infusion in patients with advanced stage nonsmall cell lung cancer.
Ann Oncol 2006, 17:11281133. PubMed Abstract  Publisher Full Text

Lecca P, Palmisano A, Ihekwaba AEC, Priami C: Calibration of dynamic models of biological systems with KInfer.
European Jour of Biophysics 2010, 39(6):1019. Publisher Full Text

Lecca P, Kahramanoğullari O, Morpurgo D, Priami C, Soo RA: Modelling and estimating dynamics of tumor shrinkage with BlenX and KInfer.

Giovacchinia G, Picchio M, Schipani S, Landoni C, Gianolli L, Bettinardi V, Muzio ND, Gilardi MC, Fazio F, Messa C: Changes in glucose metabolism during and after radiotherapy in nonsmall cell lung cancer.
Tumori 2009, 95:177184. PubMed Abstract

Brown RS, Leung JY, Kison PV, Zasadny KR, Flint A, Wahl RL: Glucose Transporters and FDG Uptake in Untreated Primary Human NonSmall Cell Lung Cancer.

S KL, Shin YB, Pyo HB, Park SH, Ogura S, Okura I: Measurement of Oxygen Concentrations in Tumor Cells by the Phosphorescence Quenching Method.

Shen J, Khan N, Lewis LD, Armand R, Grinberg O, Demidenko E, Swartz H: Oxygen Consumption Rates and Oxygen Concentration in Molt4 Cells and Their mtDNA Depleted (ρ^{0}) Mutants.
Biophysical Journal 2003, 84:12911298. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Duhaylongsod FG, Lowe VJ, Patz EFJ, Vaughn AL, Coleman RE, Wolfe WG: Lung tumor growth correlates with glucose metabolism measured by fluoride18 fluorodeoxyglucose positron emission tomography.
Ann Thorac Surg 1995, 60(5):13481352. PubMed Abstract  Publisher Full Text

Pedersen MWB, Holm S, Lund EL, Hosfjgaard L, Kristjansen PEG: Coregulation of Glucose Uptake and Vascular Endothelial Growth Factor (VEGF) in Two SmallCell Lung Cancer (SCLC) Sublines In Vivo and In Vitro.
Neoplasia 2001, 3:8087. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Noguchi Y, Saito A, Miyagi Y, Yamanaka S, Marat D, Doi C, Yoshikawa T, Tsuburaya A, Ito T, Satoh S: Suppression of facilitative glucose transporter 1 mRNA can suppress tumor growth.
Cancer Letters 2000, 154(2):175182. PubMed Abstract  Publisher Full Text

Jones RG, Thompson CH: Tumor suppressors and cell metabolism: a recipe for cancer growth.
Genes & Dev 2009, 23:537548. PubMed Abstract  Publisher Full Text

Brú A, Albertos S, Subiza JL, GarciaAsenjo JL, Brú I: The Universal Dynamics of Tumor Growth.
Biophys J 2003, 85(5):29482961. PubMed Abstract  Publisher Full Text  PubMed Central Full Text