 Research article
 Open Access
 Published:
Threedimensional tumor growth in timevarying chemical fields: a modeling framework and theoretical study
BMC Bioinformaticsvolume 20, Article number: 442 (2019)
Abstract
Background
Contemporary biological observations have revealed a large variety of mechanisms acting during the expansion of a tumor. However, there are still many qualitative and quantitative aspects of the phenomenon that remain largely unknown. In this context, mathematical and computational modeling appears as an invaluable tool providing the means for conducting in silico experiments, which are cheaper and less tedious than real laboratory experiments.
Results
This paper aims at developing an extensible and computationally efficient framework for in silico modeling of tumor growth in a 3dimensional, inhomogeneous and timevarying chemical environment. The resulting model consists of a set of mathematically derived and algorithmically defined operators, each one addressing the effects of a particular biological mechanism on the state of the system. These operators may be extended or readjusted, in case a different set of starting assumptions or a different simulation scenario needs to be considered.
Conclusion
In silico modeling provides an alternative means for testing hypotheses and simulating scenarios for which exact biological knowledge remains elusive. However, finer tuning of pertinent methods presupposes qualitative and quantitative enrichment of available biological evidence. Validation in a strict sense would further require comprehensive, casespecific simulations and detailed comparisons with biomedical observations.
Background
Introduction
Cancer is one of the main causes of mortality in the world. Statistics estimate that about one fifth of the population will suffer from cancer at some point in their lives [1]. Cancer is a category of diseases, which share several common features including sustained and uncontrolled cell proliferation, resistance to cell death, induction of angiogenesis, and activation of invasion and metastasis mechanisms [2].
The exact mechanisms that initiate cancer development remain largely unknown. However, it is widely accepted that cancer originates from cells which, due to various gene mutations, escape the body’s natural mechanisms of controlling the balance between cell proliferation and cell death [3]. These cells create a clump which grows faster than host cells. However, this small tumor grows with a decreasing rate; as the tumor grows, disorganization of the host vasculature and limited diffusion of nutrients to the center of the tumor lead to the formation of an internal necrotic core [4, 5]. Cells in the outer rim of the tumor proliferate, while cells in the interior die. For the tumor to grow large and become malignant, it needs to establish its own blood supply network, a process called angiogenesis. Angiogenesis results in a highly disorganized, tortuous and dilated vasculature, [4, 6,7,8] which however, provides the nutrients needed for further tumor growth. Evidently, during both the avascular and vascular phase of tumor development, the provision of nutrients to tumor cells through the blood supply network is highly inhomogeneous and time varying [4, 9,10,11]. In fact, many studies have shown that tumors contain hypoxic and hypoglycemic regions, particularly near the center, which affect local cell proliferation and death rates [4, 12, 13] and references therein.
Contemporary medical and biological literature on the subject shows that the vast majority of observations and results are conclusive only up to a point. In fact, there are still many qualitative and quantitative aspects of tumor progression that remain largely unknown. In this context, in silico modeling appears to be an invaluable tool for simulating scenarios and testing hypotheses pertaining to the aforementioned biological phenomena. To this end, this work describes an extensible and computationally efficient framework for in silico modeling of tumor growth in a 3dimensional, inhomogeneous and timevarying chemical environment.
Related literature
During the last few decades, mathematical and computational modeling of tumor growth has received a lot of attention from the scientific community. However, as noted in [14], there are no established first principle theories in celltissue modeling, and this seems to be the case even for today. Furthermore, there is still no generally agreed consensus on which modeling approach is the most suitable for modeling tumor growth. Scientists with different backgrounds have employed a variety of methods to attack the problem, in an effort to provide tools for conducting in silico experiments, which are significantly cheaper and less tedious than real laboratory experiments. Inspection of the literature shows that papers on the particular subject fall into two categories. Papers in the first category aim at an at least partial validation of a model against actual measurements. Such works include [15,16,17,18]. Papers in the second category aim at proposing new modeling methods or advancing already existing ones. Our work belongs in the second category.
In this section, we will review the main methods which have been most commonly employed in the pertinent literature.
Population models were probably among the first, simplest and yet effective of these approaches and utilize both deterministic and stochastic mathematics [19]. These models neglect tissue spatial structure and focus on the dynamics of the involved cell populations. They can address a variety of phenomena such as tumor clonal heterogeneity [20,21,22], tumorhost cell interactions [23,24,25] and response to therapy [26,27,28,29].
Despite the usefulness of population models, the spatial structure of tumors and the tissues they grow in appears to play an important role in tumor growth. To address this, several models have been proposed, with discrete entity, cellbased models forming a concrete class of such approaches. In these models, each tumor cell is treated as a discrete agent reacting to changes in its environment according to its own internal decision mechanism. Partial differential equations are most usually employed to model the background chemical environment. Most common approaches include latticebased [30,31,32,33,34,35,36,37,38,39,40,41,42], latticefree [14, 43,44,45] and Potts models [46,47,48,49,50]. Discrete agent models can address pertinent cellular, biochemical and biomechanical phenomena in considerable detail. However, they are computationally expensive and thus can simulate tumor sizes ranging in the order of at most 10^{6} cells, often considered in 2 dimensions.
Another popular approach is modeling the concentrations of both cells and chemical substances as continuous, spatially distributed quantities. Reactionadvectiondiffusion equations are most commonly employed to model multicellular tumor spheroids [51,52,53,54]. Mainly for (but not limited to) the case of gliomas, the reaction diffusion equation is invoked to model the infiltration of tumor cells in the surrounding healthy tissue [55,56,57,58,59,60,61].
In [62, 63] the authors simulated the temporal evolution of nonnecrotic, 2 and 3 dimensional tumors modeled as a continuum using a moving boundary. Using level set methods, this approach was extended in [64,65,66,67,68]. These papers demonstrate 2dimensional simulations which additionally considered angiogenesis, necrosis and features of the tumor microenvironment.
A different approach, employing diffuse interface, multiphase mixture models was taken in [69, 70]. In these works, tissue is modeled as a multiphase mixture of solid components (e.g. dead tumor cells, viable tumor cells, host tissue) and water. Their temporal evolution is derived by mass equations and thermodynamic constitutive laws. Further work on this approach includes [71,72,73] where the authors consider additional phenomena like angiogenesis and biomechanical effects.
Miscellaneous approaches include the spatially averaged cellular automata developed in [74,75,76] where space is discretized in voxels, each one containing a number of cells. The resulting cubic grid is treated as a cellular automaton, with specified rules governing the transition of cells through the cell cycle phases within each voxel, the expansion of the tumor and the effects of various treatment modalities. In [77, 78] a hybrid approach was taken; cell distributions were modeled as continuous quantities except for the proliferating regions at the tumor boundary, where cells were treated as discrete entities. For a recent collection of articles on multiscale cancer modelling we refer to [79].
Results
This paper develops a method for modeling 3dimensional tumor growth with an emphasis on macroscopic variables. More specifically, we focus on variables quantifying the spatial distribution of cells and molecules, the provision of nutrients by the vascular system, tumor cell proliferation and invasion, the effects of tumorinduced vascular remodeling and how they affect each other as the tumor grows. Inspired by the aforementioned literature, our work aims at providing an additional useful tool for in silico experimentation with the following properties:

a)
The resulting model is modular; that is, it consists of several discrete mathematical/algorithmic modules, each one addressing a particular biological phenomenon. This allows keeping track of the assumptions made for each module. Furthermore, it facilitates model readjustment in case a new set of hypotheses needs to be considered. As we will show in the following sections, this can be done by extending or even completely redesigning the modules pertaining to the new hypotheses.

b)
Although some of the models mentioned in the literature can address phenomena even at the single cell level, they are in general computationally intensive. For some of them, simulation times in the order of 10–24 h have been reported [39, 69]. This imposes restrictions on the shape and size of the simulated tissue; therefore, many models consider only 2dimensional tumors or spheroids with maximum size in the order of a few mm^{3}. However, realistic tumors can grow up to several cm^{3} in volume. Besides that, as demonstrated in [40], results obtained from simulated small tissue areas generally cannot be extrapolated to larger domains. This implies that a balance between consideration of microscopic details and the ability of simulating larger regions of tissue must be kept. The methods developed in this paper aim at a resulting model that can simulate large (in the order of cm^{3}) areas of tissue in 3 dimensions, with a spatial resolution in the order of 1–2 mm^{3}, i.e. the voxel size of contemporary imaging techniques.

c)
Tumor growth consists in a complex interaction of phenomena evolving in different time scales. From the macroscopic point of view we adopt, the shortest time scale concerns the diffusion of molecules (seconds) and the largest one the overall tumor expansion (months). It is definitely a challenge to choose an appropriate simulation time step, i.e. one that addresses all the involved mechanisms in a sufficient amount of temporal detail, while keeping the simulation computationally tractable. Therefore, some models focus solely on cell proliferation and neglect the local availability and diffusion of nutrients [50, 54,55,56,57,58,59,60,61, 74,75,76]. Another common approach is to choose a time step in the order of minutes or hours, and solve the resulting (quasi) steady state equations for the diffusion of nutrients [30,31,32, 41, 45, 52, 62,63,64, 66, 68, 70, 71, 73]. The methodology proposed in this work aims at time steps in the order of the shortest time scale, i.e. in the order of seconds.

d)
We also aim for computational efficiency, meant in a twofold sense.

First, in conjunction with (b) and (c) above. The simulations we present here consider tissue areas in the order of 4.2 × 4.2 × 4.2 cm^{3}, with a spatial resolution of 2 mm^{3} and a time step of 10 s for a time period of 3 months. The average simulation time on a standard desktop computer is about 10–12 min. The model is implemented in MATLAB; implementation in a precompiled language like C is expected to significantly decrease this time.

Second, for scalability reasons, the resulting model should be able to exploit multicore computation. We will show in part VII of the methods section that the collection of all model variables at each time instant (i.e. the state vector) essentially evolves in a dynamical systems fashion. Given the current state, the next state can be calculated by a sequential application of algorithmically defined operators. Each one of these operators is perfectly eligible for parallelization, thereby enabling implementations considering larger tissue areas with finer spatial and temporal resolutions.
Discussion
In this paper we present a novel methodology to approach the still open problem of modeling tumor growth. The presented modeling framework casts the problem in the realm of spatially distributed, stochastic dynamical systems by placing all pertinent spatial variables in a set of vectors, which collectively define the state vector of the overall system. At each time instant, a series of mathematically derived and algorithmically defined operators, each one corresponding to a particular biological mechanism, are applied to the state vector. Within the proposed framework, each one of these operators may be redesigned to consider different sets of starting assumptions, resulting generally in computationally efficient implementations. This facilitates the design of a large variety of hypothesis testing scenarios and corresponding in silico experiments, a process which, with the current limited qualitative and quantitative knowledge on the subject, seems inevitable. Several use cases are presented. Since we present simulations for a specific model, we do not attempt more detailed comparisons with biological data. Undoubtedly, there is still a lot of ground to be covered. Enrichment of both biological measurements and pertinent qualitative observations is necessary for further improvement of such methods. Future work should include much more casespecific and detailed comparisons between simulation results and available biological evidence.
Conclusions
We have developed an extensible and computationally efficient framework for modeling tumor growth in a threedimensional inhomogeneous and timevarying chemical environment, which constitutes an in silico alternative for testing different hypotheses and simulation scenarios. The model has been applied in the context of several use cases in order to visualize various aspects of tumour expansion and a multivariate analysis of the effects of model parameters on the number of live cancer cells of a growing tumor has been performed. Since many aspects of the pertinent biological mechanisms remain still largely unknown, finer tuning and validation of the simulation system in a strict sense presupposes qualitative and quantitative enrichment of the available biological evidence.
Methods
The rest of the paper is organized as follows. In section I we present the main ideas used to model the diffusion of particles. In section II, we discuss boundary conditions. Sections III and IV specialize the ideas of the previous sections in the cases of chemical and cellular diffusion. Section V discusses tumor cell metabolism and consumption of nutrients, and how they affect proliferation and necrosis. In section VI we model the macroscopic effects of tumorinduced vascular remodeling. Section VII presents the complete model architecture. In section VIII we present some use cases, including a multivariate study on the effects of various model parameters on the number of viable tumor cells after a period of free growth.
I. Modeling the diffusion of particles
The diffusion of particles is a natural phenomenon present in a vast variety of models regarding tumor growth. To model such phenomena, the diffusion partial differential equation is most commonly invoked:
where c(x, t) is the concentration of the species under consideration (cells or molecules) at time t and location x, and D the diffusion tensor of the species in the surrounding material. This equation has been widely used to model cell diffusion, particularly in the case of glioblastoma, as well as diffusion of molecules in tissue. In the case of isotropic diffusion D is a constant scalar, and equation (1) reduces to
where Δ_{x} is the Laplace operator in ℝ^{3}. In [80] we elaborated on the observation that (2) is the FokkerPlanck equation corresponding to the stochastic differential equation
where B_{t} denotes the standard Brownian motion in ℝ^{3}. Given the initial position x_{o} of a particle, the distribution of the random variable x_{t} (i.e. the solution of (3) at time t) provides c (x, t), i.e. the probability distribution over all possible locations of this particle at time t. The initial position of the particle may also be given in terms of a probability distribution c (x, 0). In this case, the probability distribution c (x, t) can be found by either of two equivalent ways: By solving (2) as a partial differential equation with initial value c (x, 0) to find the timely evolution of this distribution, or equivalently, by solving the stochastic differential equation (3) with initial distribution c(x, 0) to find the probability distribution of the random variable x_{t}.
To consider the collective movement of a population of particles, that is, molecules or cells located within a specified anatomic region, this notion of distribution can be utilized as follows: integration of c(x, t) over a region A of ℝ^{3} provides the fraction of the total particle population that is located in A at time t.
To model anisotropic diffusion, i.e. the preferential stochastic movement of particles along locally specific unit directions, the notion of the local diffusion ellipsoid (i.e. the diffusion tensor) is needed [59]. This notion corresponds to defining an ellipsoid in each point of the 3dimensional space under consideration. Mathematically, this is made explicit by defining, in each point, a 3 × 3 positive definite symmetric matrix:
This matrix can be decomposed in the following form:
where λ_{1}, λ_{2}, λ_{3} are the eigenvalues of D and u_{1,}u_{2}, u_{3} are the corresponding orthonormal eigenvectors. We note that D is positive definite, hence λ_{1}, λ_{2} and λ_{3} are positive numbers. The eigenvalues and eigenvectors of D define an ellipsoid whose principal axes lie on the directions of u_{1}, u_{2} and u_{3}. The principal axes have lengths \( 2\sqrt{\lambda_1} \), \( 2\sqrt{\lambda_2} \) and \( 2\sqrt{\lambda_3} \), respectively. The diffusion ellipsoid is depicted in Fig. 1.
The anisotropic diffusion of particles as dictated by the (local) diffusion ellipsoid can be modelled as follows; let p(x, x_{o}, dt) denote the probability of a particle starting at x_{o} to be at x after time dt. Then
Where
Note that the righthand side of (4) is essentially an anisotropic Gaussian in ℝ^{3}. The parameter α is a positive scalar, specific to the particles under consideration which in our case, will be tumor cells. This parameter rescales the eigenvalues, thereby rescaling conformally the axes of the diffusion ellipsoid. This reflects the fact that different kinds of particles may tend to move along the axes of same ellipsoid, but may do so with different velocities. Of note, since we are considering a local diffusion ellipsoid, in the most general case both eigenvalues and eigenvectors are functions of the position x. After some mathematical elaborations detailed in [80], equation (4) leads us to model anisotropic diffusion by
where b is a 3dimensional, normally distributed random vector with zero mean and covariance matrix the identity matrix times dt:
or equivalently, by the stochastic differential equation
For a given initial distribution c(x, 0) the solution of (5) provides the probability distribution of the random variable x_{t}, which can be interpreted exactly as described for the isotropic case. In the case of brain tumors, measurements concerning the diffusion tensor are obtained through the Diffusion Tensor Imaging (DTI) technique [56, 59].
In our approach, equations (3) and (5) will constitute the theoretical basis for modelling chemical and cellular diffusion.
Let us now assume that the diffusion of the particles (molecules or cells) under consideration takes place in a cubic lattice consisting of N × N × N geometrical cells (voxels). Each voxel is a cube of dimensions Δs × Δs × Δs. We fix a temporal discretization step equal to Δτ. Voxels in that cubic lattice can be classified into 4 categories, depending on the number of their neighboring voxels within the lattice. Voxels in the interior of the lattice have 26 neighbors. Voxels at the outer faces of the lattice have 17 neighbors. Voxels at the outer edges of the lattice have 11 neighbors and voxels at the outer vertices of the lattice have 7 neighbors. Furthermore, for each particular voxel, its neighboring voxels fall into 3 categories: the ones that share a common face, the ones that share a common edge and the ones that share a common vertex with the particular voxel.
For the remaining part of this section, let us adopt the assumption that, at each discrete time point, the distribution of the particles under consideration within each voxel is uniform. This is somewhat oversimplifying, and we will further elaborate on this assumption in the following sections, where we specifically consider chemical or cellular diffusion. For the moment, this assumption will render the presentation of the main ideas more straightforward.
Under the uniformity assumption, at any discrete time point t and for any pair of voxels A and B (not necessarily different) we can calculate the probability for a particle to lie within B at time t + Δτ, given that its position at time t is a uniformly distributed (u.d.) random variable supported in A. Numerical integration of equations (3) and (5) for a uniform initial distribution provides the means for a Monte Carlo calculation of this probability, which we will denote by Pr(A → B).
Let A denote a voxel not lying on the boundary of the lattice nor being adjacent to it, and B_{i}, i = 1, …, 26 its neighbors. A key step to our discretization process is to choose the voxels’ edge length Δs and the time step Δτ such that, from one time instant to the next one, the particles lying in each voxel diffuse at most into its neighbors. Mathematically, this means that Δs and Δτ should be chosen such that
If (6) holds then, for a spatially constant, isotropic diffusion like the one implied in (3), for each voxel A not lying on the boundary of the lattice nor being adjacent to it, these probabilities consist of essentially 4 numbers: one for the particles that are in A and will remain in A (i.e. Pr(A → A)) and three more, one for each of the common face, common edge and common vertice neighbors of A. Furthermore, for any two voxels A and B it holds that Pr(A → B) = Pr (B → A). For a spatially varying, anisotropic diffusion like the one implied in (5), these symmetries do not hold; the respective probabilities should be precalculated independently for any particular voxel.
Thus, under the uniformity assumption, for any voxel A not lying on the boundary of the lattice nor being adjacent to it, knowing the particle population of the voxel Q_{t}(A) and its neighbors Q_{t}(B_{i}), i = 1, …, 26 at a time instant t, allows us to calculate the population within A at the next time instant t + Δτ by
Apparently, this equation does not hold as such when the voxel under consideration lies on the boundary of the lattice. We will deal with these voxels (and their neighbors) in detail in the next section, where we discuss boundary conditions.
In what follows, it will come in handy to represent particle quantities within the voxels in vector form. Let us make this representation explicit by the following construction. Let the coordinates of each voxel in the lattice be given by a triad of integers, (i, j, k) where i, j, k = 1, …, N. We define a mapping L : ℕ^{3} → ℕ as follows: L(i, j, k) = i + (j − 1)N + (k − 1)N^{2}. This mapping is a bijection from the set of triads (i, j, k), i, j, k = 1, …, N to the integers from 1 to N^{3}. Let Q(i, j, k) denote the quantity of particles in the respective voxel. We define the vector q of N^{3} elements by q(L(i, j, k)) = Q(i, j, k). By this construction, the elements of the 3d matrix Q are explicitly mapped to the elements of the onedimensional vector q.
II. Dirichlet and Neumann boundary conditions
Knowing the particle population within each voxel of the lattice at a time instant, equation (7) allows us to calculate the population within each voxel at the next time instant for all the voxels, except the ones on the boundary of the lattice and their neighbors. Each of the boundary voxels has less than 26 neighbors. In fact, each of these voxels would have 26 neighbors in an infinite lattice, but not all of them are included in a bounded, N × N × N lattice. Unless we specify boundary conditions, the calculation in (7) cannot be carried out neither for those voxels and consequently, nor their neighbors. There are two types of boundary conditions that can be imposed in a classical diffusion problem, those of the Dirichlet type and those of the Neumann type. In this work, we consider two specific types of such boundary conditions, namely, the timeindependent Dirichlet and the homogeneous Neumann boundary conditions.
Timeindependent Dirichlet boundary conditions express the requirement that on the boundary of the region under consideration, the quantity of interest does not change with time. In the framework presented here, this is expressed mathematically by the following: If a voxel A lies on the boundary of the lattice, to calculate its population at the next time instant, instead of equation (7) simply apply Q_{t + Δτ}(A) = Q_{t}(A). Additionally, if a voxel is adjacent to the boundary, simply apply (7) by using the respective probabilities as they are calculated from the numerical integration of (3) or (5).
We note that in the case of tumor growth, some authors have also considered timedependent, periodic Dirichlet boundary conditions [40]. This is also feasible in the proposed framework, and can be implemented as follows. For any boundary voxel A, at any time instant t, apply Q_{t}(A) = g_{t}, where g_{t} is the desired periodic function.
Homogeneous Neumann boundary conditions express the requirement that the flux of the quantity of interest across the boundary of the region under consideration should be zero, i.e. the boundary is nonpermeable. In terms of calculus this is expressed by the requirement, at any time instant and at any point of the boundary, the projection of the gradient of the quantity on the outward normal of the boundary at that point to be zero. In the stochastics literature, a nonpermeable boundary within which a random motion takes place is often referred to as a reflecting boundary [81]. In the framework presented here, this is expressed mathematically as follows.
As previously mentioned, any voxel A lying on the boundary has 17, 11 or 7 neighbor voxels which we denote by B_{i}, where i is an integer from 1 up to 17, 11 or 7, depending on the position of the voxel. For each such voxel A, and its neighbors B_{i} we calculate the probabilities Pr(A → B_{i}) as described in the previous section, but only for the voxels that are contained in the lattice. Specifically, if A lies on the boundary, we calculate the probabilities Pr(A → A), Pr(A → B_{i}) where i is an integer from 1 up to 17, 11 or 7. We then normalize these probabilities to sum to one. By this calculation we acquire the probabilities we need, in order to apply equation (7) for any voxel in the lattice. Using these normalized probabilities when applying equation (7) for either boundary voxels, or their neighbors, ensures that every particle lying in a voxel on the boundary, will remain within the lattice, that is, within the region of interest, thereby capturing the notion of a reflecting boundary.
In our case, the region of interest is a cube. It is apparent that these methods of imposing Dirichlet or Neumann boundary conditions apply to more complex shapes, as long as space is properly discretized [82]. The case of brain tumors, where the skull naturally imposes a reflecting boundary to the diffusion of tumor cells is an example where this approach may be useful.
III. Diffusion of glucose and oxygen
In this section we will use the ideas presented previously to develop a model for the diffusion of chemical molecules (glucose or oxygen) in the region of interest, that is, the cubic lattice of dimensions N × N × N. For the moment, we will assume that the diffusion is isotropic and that the distribution of molecules within each voxel is uniform. Let q_{t} denote the N^{3} × 1 vector whose entries are the quantities of the molecule under consideration within each voxel at time t. Αssuming timeindependent Dirichlet boundary conditions, equation (7) implies that if we know q_{t}, we can calculate q_{t + Δτ} by performing a linear calculation. This means that there is a N^{3} × N^{3} square matrix T such that q_{t + Δτ} = Tq_{t}, .
We remind the reader that due to the symmetries holding for isotropic diffusion, for each voxel A not on the boundary, we can apply equation (7) by using only four numbers. We denote these numbers by Pr(A → A), Pr(F → A) (for common face neighbors), Pr(E → A) (for common edge neighbors), and Pr(V → A) (for common vertex neighbors). These numbers can be precalculated by numerically integrating (3) for a uniform initial distribution, where in each case, D is taken to be the diffusion coefficient of the respective molecule. We subsequently use these values and the mapping L defined in section I to construct the matrix T according to the following algorithm:
Note that, in view of the mapping L, each row i of T corresponds to a specific voxel A, i.e. to a specific position i of the vector q_{t + Δτ}. The nonzero entries of the particular row correspond to the probabilities Pr(A → A) and Pr(B_{i} → A), where B_{i} are the neighbors of A, as indicated by (7). Furthermore, each column j of T corresponds to a specific voxel B, i.e. to a specific position j of the vector q_{t}. The nonzero entries of the particular column correspond to the probabilities Pr(B → B) and Pr(B → A_{i}), where A_{i} are the neighbors of B. Probabilities of the type Pr(A → A) lie in the main diagonal of T.
In our model, we will use two such matrices, one for glucose and one for oxygen, denoted by T_{gl}, T_{o} respectively. Each of the respective N^{3} × 1 vectors will be denoted by gl_{t}, o_{t}.
The matrix constructed by Algorithm 1 has an interesting property. For sufficiently large N and due to the assumption implied in equation (6), it is a sparse matrix: in each row, at most 27 elements are nonzero. This provides a significant relief of the computational burden of the entire model.
Note that all arguments in this section, resulting in the simple model q_{t + Δτ} = Tq_{t} for molecular diffusion, rely on the uniformity assumption as it was stated in section I.
We mentioned that this assumption is somewhat oversimplifying. Indeed, in vivo measurements in tumor areas show that chemical gradients can be spatially nonuniform and time varying. Reportedly, oxygen profiles can vary locally up to 50% within an hourly time frame [83]. The highly irregular tumor vasculature may further complicate things and pertinent biological mechanisms remain largely unknown [6, 7]. Inherent stochasticity is also expected to play a role. Consequently, detailed quantification of the inflicted macroscopic effects such as time evolution of chemical fields seems currently infeasible. Therefore, we broaden our perspective as follows.
First, we note that our calculations showed that for both oxygen and glucose, it holds that Pr(A → A)> Pr(F → A)> Pr(E → A)> Pr(V → A) and that each of these numbers is one order of magnitude greater than the next one. Therefore, we will use the precalculated numbers Pr(A → A), Pr(F → A), Pr(E → A), and Pr(V → A) only as estimates for the (relative) orders of magnitude of these probabilities. To take the largely stochastic, collective effect of the aforementioned complex mechanisms and uncertainties into consideration, we will choose a time step Δτ in the order of seconds (i.e. the time scale of chemical diffusion) and introduce a certain degree of randomness to the matrices T_{gl} and T_{o}. Specifically, every some time steps, each column of these matrices corresponding to a voxel in the interior or in the vicinity of the tumor will be randomly perturbed, such that its nonzero entries retain their relative orders of magnitude and their sum remains one. We provide the respective implementation details in the appendix.
Introducing stochasticity implies that for each simulation scenario, several simulations will be required. Nonetheless, the results of these multiple simulations will enable a more comprehensive consideration of possible outcomes.
IV. Diffusion of tumor cells
To model the diffusion of cancer cells, we make the following assumptions:

i.
There are four types of cells within each voxel; live normal (host) cells, necrotic normal cells, live tumor cells and necrotic tumor cells.

ii.
For any voxel, there is an average cell population capacity (tumor +normal) which we denote by M and a maximum cell population capacity which we denote by M_{max}.

iii.
As the tumor grows, live normal cells become either dislocated by invading tumor cells or necrotic.

iv.
When the sum of living cancer, necrotic cancer and necrotic host cells in a voxel exceeds M_{max}, living cancer cells in excess of M_{max} invade neighboring voxels according to equation (5).

v.
Both tumor and normal necrotic cells remain in the voxel they became necrotic, that is, they do not “invade” neighboring voxels.
Our model will be based on equations (5) and (7). Since we assume that in each voxel, only tumor cells in excess of M_{max} diffuse to neighboring voxels, the probability Pr(A → A) is not needed. For each voxel A and its neighbors B_{i} we calculate the probabilities Pr(A → B_{i}) and normalize them to sum to 1. We then place these probabilities in the respective rows and columns of a matrix T_{c}, similar to the matrices T_{gl} and T_{o} of the previous section. Of note, the main diagonal of the matrix T_{c} consists of zeros, and not of probabilities Pr(A → A), as is the case with T_{gl} and T_{o}. Apparently, the matrix T_{c} is also sparse, lightening thereby the computational burden. Furthermore, by construction, the matrix T_{c} implies homogeneous Neumann boundary conditions for living tumor cells.
The tumor cell diffusion algorithm has a quite simple implementation. Let l_{t}, nc_{t}, nn_{t} denote the N^{3} × 1 vectors whose entries are respectively the live tumor, necrotic tumor and necrotic normal cells within each voxel at time t (we use this notation throughout the text, please see also the first paragraphs of the following section). Let w_{t} denote the sum of these vectors. Let \( \overline{M_{max}} \) denote the N^{3} × 1 vector whose entries are all M_{max}. We also adopt the following notation: for a vector x and a number μ, the vector (x ≥ μ) is the binary vector with elements set to 1 if the corresponding element of x is ≥μ and 0 otherwise. Finally, for any two vectors x and y let x. ∗ y denote their elementwise product. The resulting algorithm boils down to a simple vector algebraic representation:
Vector s_{1} contains, for each voxel, the total number of live (tumor and normal) cells that the particular voxel can hold additional to its necrotic cells. In view of assumptions (ii), (iv) and (v), s_{2} entries are numbers of live cancer cells which already exist in each voxel that can also remain in the respective voxel. In view of assumptions (iii) and (iv), the vector (l_{t} − s_{2}) contains the numbers of live tumor cells that lie within each respective voxel in excess of M_{max}, and therefore invade neighboring voxels by dislocating live normal cells. The populations of live tumor cells at the next time instant, i.e. after their diffusion to neighboring voxels is given by the sum of s_{2} and T_{c} ∗ (l_{t} − s_{2}).
Of note, the entire approach again relies on the uniformity assumption, but this time for tumor cells. Apparently, this is a simple, and definitely not the only way to model the invasion of tumor cells. For the simulations presented below, we will rely on it. More sophisticated approaches could introduce randomness in the matrix T_{c} exactly as described in section III; probabilities in T_{c} could be dynamically readjusted by taking into account nutrient quantities within each voxel, thereby introducing some at least rough, phenomenological notion of chemotaxis in the model. However, our aim for the moment is to keep the presentation of the main ideas as simple as possible. We will provide some suggestions for further elaboration and implementation of this module in the appendix.
V. Proliferation/necrosis according to cell metabolism and local availability of oxygen and glucose
In this section, we will study the proliferation of cells within a voxel A. For a time interval equal to the time step Δτ, we will neglect diffusion phenomena, and study the proliferation and necrosis of cells within A as they are dictated by the cells’ consumption needs and the local availability of oxygen and glucose. Although special effort has been made in order for the proposed approach to be founded on a consensus of biological evidence, it is definitely not the only one possible; somewhat different assumptions may lead to different approaches. In section VII, where we present the complete model architecture, it will become apparent that the approach adopted here can be completely readjusted, in order to consider different assumptions.
For any voxel A, we will need the following variables:
 l_{t}(A): number of living cancer cells within A at time t.
 nc_{t}(A): number of necrotic cancer cells within A at time t.
 nn_{t}(A): number of necrotic normal cells within A at time t.
The number of live normal cells within A at time t is given by s(M − l_{t}(A) − nc_{t}(A) − nn_{t}(A)) where s(∙) is the known function s(x) = x if x ≥ 0 and 0 if x < 0.
 o_{t}(A): Oxygen quantity (pmols) within A at time t. Initial value o_{0}(A) is the same for every voxel, denoted by \( \overline{o_0} \). \( \overline{o_0} \) is calculated such that the initial concentration of oxygen in each voxel equals the concentration of dissolved oxygen in the blood (see section VI).
 gl_{t}(A): Glucose quantity (pmols) within A at time t. . Initial value gl_{0}(A) is the same for every voxel, denoted by \( \overline{gl_0} \). \( \overline{gl_0} \) is calculated such that the initial concentration of glucose in each voxel equals the concentration of glucose in the blood (see section VI).
 o b_{t}(A): sec) by the local vascular network within/sec) by the local vascular network within A during the previous time interval t−Δτ→ t
 gl_b_{t}(A): Glucose supply rate (pmols/sec) by the local vascular network within A during the previous time interval−Δτ→t. The variables o_b_{t}(A) and gl_b_{t}(A) quantify macroscopically the role of the local vascular network in the provision of oxygen and glucose within each voxel. We will assume that their values remain constant during each time interval t → t + Δτ. In fact, these variables are too subject to a dynamic time evolution due to the effects of local vascular remodeling induced by the tumor. We will discuss this matter in detail in the next section. In this section, we will focus on the time interval t → t + Δτ, and aim at calculating the aforementioned cell populations and chemical quantities at the next time instant, i.e. t + Δτ. We note that the adopted notation implies that o_b_{t + Δτ}(A), gl_b_{t + Δτ}(A) denote the oxygen and glucose supply rates during the time interval under consideration, i.e. t → t + Δτ. In this section, these quantities are assumed to be (pre)calculated by the algorithm described in the next section.
We will further need the following parameters:
 M: average cell population capacity for each voxel.
 K_{o}: oxygen consumption rate (pmols/sec) of a normal cell.
 K_{gl}: glucose consumption rate (pmols/sec) of a normal cell. Typically, for a normal cell acquiring its energy mainly through combustion of glucose, K_{o} is 4 to 6 times larger than K_{gl} [36, 37, 84].
 K_{ATP}: ATP consumption rate (pmols/sec) of a normal cell.
 λ: The product λK_{ATP} defines the ATP consumption rate of an actively proliferating tumor cell. Since proliferating tumor cells consume much more resources than normal cells, λ should be >1 [39, 46, 85]. For quiescent tumor cells, λ is assumed to be 1.
 cc: Cell cycle duration (sec).
 a_{max}: maximum mitosis rate that tumor cells can achieve, when they are not limited by local oxygen and glucose levels. (mitoses/cell/ time step in secs).
Of note, concerning the initial values of o_b_{t}(A) and gl_b_{t}(A) i.e. o_b_{0}(A) and gl_b_{0}(A), we will make the following assumption. For every voxel at time t = 0 the provision of oxygen and glucose by the local vascular network and their consumption by normal cells should balance each other, such that the respective background concentrations remain constant, i.e. o_b_{0}(A) = MK_{o} and gl_b_{0}(A) = MK_{gl} [69, 71].
A normal cell’s ATP consumption is dictated by K_{o}, K_{gl} and the stoichiometry of clean combustion and glycolysis:
Thus, knowing K_{o}, K_{gl} allows us to explicitly determine K_{ATP}= \( \left(\raisebox{1ex}{$17$}\!\left/ \!\raisebox{1ex}{$3$}\right.\right) \) K_{o}+ 2K_{gl}.
The stoichiometry of the clean combustion of glucose requires that the glucose/oxygen uptake ratio is 1:6. It is well documented that for cancer cells, due to increased utilization of glycolysis, this is not the case [2]. Experimental measurements and estimations report that the ratio of glucose /oxygen consumption in tumors can vary up to 1:1 or even more [14, 46, 86,87,88]. Compared to clean combustion, glycolysis is 18 times less efficient in ATP production and cancer cells compensate this deficiency by upregulating glucose transporters, thereby increasing glucose import in the cytoplasm. Considerably increased glucose uptake and utilization has been reported in a variety of tumors by the use of positron emission tomography (PET) [2]. It has been also reported that local levels of oxygen and glucose have an effect on this ratio [86, 89]. Quantitative details of these phenomena are still unclear. The possibility that a single cell may employ both glycolysis and normal aerobic metabolism is not excluded. Qualitatively it seems evident that when oxygen falls below a certain threshold, cells tend to switch to a glycolytic phenotype. However, this observation does not tell the whole story, since cancer cells switch to glycolysis even when oxygen levels are abundant [90, 91].
To take account of this evidence from a modeling perspective, we take the following approach. We assume that tumor cells within a voxel may obtain the energy they need (i.e. λK_{ATP} pmols/sec) by acquiring a fraction β of it by glycolysis and the remaining fraction by combustion. This fraction will depend on the local availability of oxygen and glucose, plus, we will introduce a degree of randomness in it. Two additional, casespecific parameters are the minimum and maximum values of this fraction, which we respectively denote by β_{1} and β_{2}. Apparently, 0 ≤ β_{1} ≤ β ≤ β_{2} ≤ 1.
We have assumed that a quiescent tumor cell needs K_{ATP} pmols ATP/sec to stay alive. From this amount and according to the stoichiometry, βK_{ATP} pmols/sec should come from glycolysis of βK_{ATP}/2 pmols/sec glucose. The remaining (1 − β)K_{ATP} pmols/sec should come from combustion of (1 − β)K_{ATP}/36 pmols/sec glucose with (1 − β)K_{ATP}/6 pmols/sec oxygen. Thus, to stay alive, a quiescent tumor cell needs
On the other hand, an actively proliferating tumor cell needs λK_{ATP} pmols ATP/sec. From this amount, βλK_{ATP} pmols/sec should come from glycolysis of βλK_{ATP}/2 pmols/sec glucose. The remaining (1 − β)λK_{ATP} pmols/sec should come from combustion of (1 − β)λK_{ATP}/36 pmols/sec glucose with (1 − β)λK_{ATP}/6 pmols/sec oxygen. Thus, an actively proliferating tumor cell needs
With the above evidence and assumptions in mind, we devise the following algorithm, which will be executed for each voxel A at each time step. The involved calculations require a detailed analysis, consisting of several steps and subcases.
Step 1. In this step, we will calculate the following quantities:
O_{av}: available oxygen for tumor cells in A during Δτ.
Gl_{av}: available glucose for tumor cells in A during Δτ.
nn_{t + Δτ}(A): number of necrotic normal cells in A at the next time instant.
First, we calculate the amount of oxygen and glucose that will be available for tumor cells, by subtracting the consumption of normal cells:
Case 1.1: O_{1} ≥ 0 and Gl_{1} ≥ 0. In this case, oxygen and glucose suffice for all normal cells in A to stay alive, hence,
Case 1.2: O_{1} < 0 or Gl_{1} < 0. This means that either oxygen and/or glucose do not suffice for all normal cells in A to stay alive. We calculate
\( {N}_n=\mathit{\min}\left(\frac{o_t(A)+o\_{b}_{t+\varDelta \tau}(A)\varDelta \tau}{K_o\varDelta \tau},\frac{gl_t(A)+ gl\_{b}_{t+\varDelta \tau}(A)\varDelta \tau}{K_{gl}\varDelta \tau}\right) \), i.e. how many living normal cells in A will stay alive.
\( {\overline{N}}_n=s\left(M{l}_t(A){nc}_t(A){nn}_t(A)\right){N}_n \), i.e. how many normal cells in A will become necrotic.
In this case, the aforementioned quantities O_{av}, Gl_{av} and nn_{t + Δτ}(A) are
Step 2. In this step, by taking into account the results of Step 1, we will study the proliferation/necrosis of tumor cells. Eventually, we will calculate the quantities l_{t + Δτ}(A), nc_{t + Δτ}(A), o_{t + Δτ}(A), gl_{t + Δτ}(A). The quantity nn_{t + Δτ}(A) has been calculated in step 1.
Case 2.1: If l_{t}(A) = 0, i.e. there are no living tumor cells in the voxel, the calculation is simple:
Case 2.2. If l_{t}(A) > 0 the calculation is more elaborate. Let a be the mitosis rate of tumor cells per time step Δτ, i.e. the fraction of tumor cells within A that will divide during Δτ and cc the duration of their cell cycle. The duration of their cell cycle in time steps is cc/Δτ. Assuming a uniform distribution of proliferating cells at all time steps of the cell cycle, we estimate a total number al_{t}(A)(cc/Δτ) of actively proliferating tumor cells. Some stochasticity may be introduced in this estimate, but to keep the presentation simple we will not go into details.
In view of (8) and (9), we define
The quantities O_{c} and Gl_{c} are the amounts of oxygen and glucose that will be needed by tumor cells in A in order to proliferate with mitosis rate a. Those quantities should be limited respectively by O_{av} and Gl_{av}, calculated from Step 1. Mathematically, this is expressed by the inequalities O_{c} ≤ O_{av} and Gl_{c} ≤ Gl_{av}. These inequalities require a closer examination and can be equivalently written in the form
The inequalities (10) and (11) involve the up to now undetermined variables a and β. They reflect the fact that, for the tumor cells in A to proliferate with mitosis rate a, a number β ∈ [β_{1}, β_{2}] should exist, such that a ≥ 0, and (10), (11) are satisfied.
Investigation of (10): For each β ∈ [β_{1}, β_{2}] we define the function
For each given β ∈ [β_{1}, β_{2}], we observe the following:
10a) If a_{o}(β) ≥ 0: a_{o}(β) is the maximum mitosis rate that the tumor cells can achieve for the specific β, subject solely to the limitations imposed by the available oxygen in A.
10b) If a_{o}(β) < 0, it is implied that \( {O}_{av}{l}_t(A)\frac{1\beta }{6}{K}_{ATP}\varDelta \tau <0. \) This means that for the specific β, no positive mitosis rate can be achieved. In fact, available oxygen does not suffice for all tumor cells in A to stay alive.
Furthermore:
10c) β = 1 means that tumor cells can acquire the energy they need relying solely on glycolysis. Note that β → 1 implies a_{o}(β) → + ∞, reflecting the fact that in this case, the proliferation of tumor cells is not limited by the available oxygen.
10d) For each β ∈ [0, 1), a_{o}(β) is an increasing function of β.
10e) If a β ∈ [β_{1}, β_{2}] such that a_{o}(β) ≥ 0 exists, it should also satisfy
Since a_{o}(β) is increasing, \( \underset{\_}{\beta } \) is actually the lowest number for which a_{o}(β) ≥ 0. Hence,

If \( \underset{\_}{\beta }>{\beta}_2 \), we have that for each β ∈ [β_{1}, β_{2}] it holds that a_{o}(β) < 0. According to (10b), this means that the available oxygen in A does not allow proliferation and that not all tumor cells in A can stay alive. Since β is the percentage with which tumor cells rely on glycolysis, in this case, for any proliferation to happen, the available oxygen imposes greater reliance on glycolysis than the maximum, i.e. β_{2}.

If \( \underset{\_}{\beta}\le {\beta}_2 \), we have that for each \( \beta \in \left[\max \left(\underset{\_}{\beta },{\beta}_1\right),{\beta}_2\right] \) it holds that a_{o}(β) ≥ 0
Investigation of (11): For each β ∈ [β_{1}, β_{2}] we define the function
For each given β ∈ [β_{1}, β_{2}], we observe the following:
11a) If a_{gl}(β) ≥ 0: a_{o}(β) is the maximum mitosis rate that the tumor cells can achieve for the specific β, subject solely to the limitations imposed by the available glucose in A.
11b) If a_{gl}(β) < 0, it is implied that \( {Gl}_{av}{l}_t(A)\frac{17\beta +1}{36}{K}_{ATP}\varDelta \tau <0. \) This means that for the specific β, no positive mitosis rate can be achieved. In fact, available glucose does not suffice for all tumor cells in A to stay alive.
Furthermore,
11c) For each β ∈ [0, 1], a_{gl}(β) is a decreasing function of β.
11d) If a β ∈ [β_{1}, β_{2}] such that a_{gl}(β) ≥ 0 exists, it should also satisfy
Since a_{gl}(β) is decreasing, \( \overline{\beta} \) is actually the highest number for which a_{gl}(β) ≥ 0. Hence,

If \( \overline{\beta}<{\beta}_1 \), then for each β ∈ [β_{1}, β_{2}] it holds that a_{gl}(β) < 0. According to (11b), this means that the available glucose in A does not allow proliferation and that not all tumor cells in A can stay alive. Since β is the percentage with which tumor cells rely on glycolysis, in this case, for any proliferation to happen, the available glucose imposes lower reliance on glycolysis than the minimum, i.e. β_{1}.

If \( \overline{\beta}\ge {\beta}_1 \), we have that for each \( \beta \in \left[{\beta}_1,\mathit{\min}\left(\overline{\beta},{\beta}_2\right)\right] \) it holds that a_{o}(β) ≥ 0
We are now ready to complete Step 2:
Case 2.2.1. If \( \beta >{\beta}_2 \) or \( \overline{\beta}<{\beta}_1 \) or \( \min \left(\overline{\beta},{\beta}_2\right)<\max \left(\underset{\_}{\beta },{\beta}_1\right) \), the analysis above implies that for each β ∈ [β_{1}, β_{2}], the available resources in A (oxygen and/or glucose) do not suffice for all tumor cells in A to stay alive. We proceed as follows:
We pick a random \( \overset{\sim }{\beta } \) ∈[β_{1}, β_{2}].
If \( \overset{\sim }{\beta}\ne 1 \), the number of tumor cells that will remain alive is given by
If \( \overset{\sim }{\beta }=1 \), the number of tumor cells that will remain alive is given by
In any case, the number of tumor cells that will become necrotic is given by
Hence, for the next time instant we have
Case 2.2.2. If \( \beta \le {\beta}_2 \), \( \overline{\beta}\ge {\beta}_1 \) and \( \min \left(\overline{\beta},{\beta}_2\right)\ge \max \left(\underset{\_}{\beta },{\beta}_1\right) \), i.e. the complement of the condition in Case 2.2.1 holds, according to the preceding analysis we have that for each \( \beta \in \left[\max \left(\underset{\_}{\beta },{\beta}_1\right),\min \left(\overline{\beta},{\beta}_2\right)\right] \) there exists a nonnegative mitosis rate a such that the inequalities O_{c} ≤ O_{av} and Gl_{c} ≤ Gl_{av} are satisfied. Again, we pick a random \( \overset{\sim }{\beta } \) in \( \left[\max \left(\underset{\_}{\beta },{\beta}_1\right),\min \left(\overline{\beta},{\beta}_2\right)\right] \).
If \( \overset{\sim }{\beta}\ne 1 \), the corresponding mitosis rate is
If \( \overset{\sim }{\beta }=1 \), the corresponding mitosis rate is
Hence, for the next time instant we have
This concludes our analysis. Το summarize, tumor cells may be able to rely on glycolysis within certain limits, i.e. β_{1} and β_{2}, but these limits may become narrower by the available quantities of glucose and oxygen; in that case, the ability of tumor cells to proliferate with high mitosis rates is impaired. Mathematically, this is reflected by the decreased probability to attain high mitosis rates, or even the inability of tumor cells to stay alive.
We note that the calculation of O_{av} and Gl_{av} in Step 1 implies that normal cells are the first to fulfill their needs by the existing resources. This is certainly not accurate; again, a more realistic approach would be to introduce some randomness in the percentage of resources that would be available for normal versus tumor cells. Α random number roughly proportional to the ratio normal/tumor cells in the voxel may be a reasonable choice.
VI. The effects of tumorinduced vascular remodeling
In this section, we will propose a method to quantify the effects of tumorinduced vascular remodeling, in terms of how it affects the local provision of oxygen and glucose by the vascular system. The quantification we propose is based on the basic physiology of the vascular network plus additional biological observations regarding how it is affected by tumor growth.
Ideally, the vascular network works like a buffer of nutrients, in our case, oxygen and glucose. If, at a specific time instant, the concentration of a substance dissolved in the blood is higher than the respective concentration in the surrounding tissue, the substance diffuses through the vessel walls towards the tissue until the two concentrations are equal, and vice versa. The speed of this diffusion process as well as the capability of the local vascular network to quickly balance these concentrations is limited by the vessel density and total surface of the blood vessel walls in the region under consideration [92].
In our model, we will assume that the concentrations of dissolved oxygen and glucose in the blood are constant. In a voxel A of specific volume, these concentrations correspond to quantities of oxygen and glucose within A, which we denote by \( \overline{o_0} \) and \( \overline{gl_0} \). If at a specific time t, the quantity of say, oxygen in A i.e. o_{t}(A) is lower (higher) than \( \overline{o_0} \), this implies that the concentration of oxygen in A is lower (higher) than the concentration of dissolved oxygen in the blood. Hence, the provision of oxygen in A i.e. o_b_{t}(A) should increase (decrease) to level this imbalance. We model this increase (decrease) during each time step by a random fraction of the quantity \( \left(\left(\ \overline{o_0}{o}_t(A)\ \right)/\varDelta \tau \right) \). The respective quantity \( \left(\left(\ \overline{gl_0}{gl}_t(A)\ \right)/\varDelta \tau \right) \) is used for glucose. This results in the following equations:
where o_b_{t + Δτ}(A), gl_b_{t + Δτ}(A) are the oxygen and glucose supply rates (pmols/sec) by the local vascular network within A during the time interval t → t + Δτ. The quantities o_b_{t}(A) and gl_b_{t}(A) are the respective supply rates during the previous time interval, i.e. t − Δτ → t. The numbers r_{1} and r_{2} are random numbers uniformly distributed in the interval [0, 1].
We note that o_b_{t}(A) and gl_b_{t}(A) may take negative values. This simply reflects the fact that when the respective substance concentration in A is higher than the one in the blood and the consumption of the substance within A is low enough, diffusion may happen towards the vessels, decreasing thereby the substance quantity in the surrounding tissue.
It is clear, however, that o_b_{t}(A) and gl_b_{t}(A) cannot grow unboundedly neither towards positive nor towards negative values. As previously mentioned, they are limited by the local vessel density and total surface of the blood vessel walls in A. A detailed quantitative analysis addressing the pertinent mechanisms would render the model extremely complex. We therefore opt for a more macroscopic approach.
In [93] a maximum value of oxygen consumption for normal mammalian cells is given. From this, an upper bound for the absolute value of o_b_{t}(A) can be deduced. Assuming that in the tissue under consideration, normal cell metabolism does not change and utilizes oxygen and glucose in a steady ratio, we can deduce a similar bound for the absolute value of gl_b_{t}(A). We denote these two constant numbers by o_b_max and gl_b_max. In normal tissue, these numbers remain constant and are the same for each voxel. In our case, however, these bounds may be different for each voxel and are subject to a temporal evolution, inflicted by the tumorinduced vessel regression and angiogenesis. To address this in our model, we introduce the N^{3} × 1 vectors o_b_ max_{t} and gl_b_ max_{t}. Consistent with our previous notation, o_b_max_{t}(A) and gl_b_max_{t}(A) denote the aforementioned bounds for the voxel A during the time interval t − Δτ → t. For each voxel A, the initial values of o_b_max_{t}(A) and gl_b_max_{t}(A) at time t = 0 are respectively the constants o_b_max and gl_b_max. Essentially, o_b_max_{t}(A) and gl_b_max_{t}(A) quantify the capacity of the vascular network within A to provide/absorb molecules to/from the surrounding tissue, leveling thereby the concentration imbalances between blood and tissue. A disorganized and regressed vascular network in A implies lower values for o_b_max_{t}(A) and gl_b_max_{t}(A), as is usually the case in the interior of a tumor. A robust, dense vascular network in A implies higher values for o_b_max_{t}(A) and gl_b_max_{t}(A), as is the case for the outer proliferating rim of a tumor.
In what follows, we will use the temporal evolution of these quantities in each voxel to quantify the effects of tumorinduced vascular remodeling on the local provision/absorption of nutrients. We start from some basic biological background.
It is well documented, that as tumors grow, they coopt and affect the preexisting host vasculature by a number of ways, for which the collective term tumorinduced vascular remodeling is commonly used. Tumorinduced vascular remodeling consists in several mechanisms, including vessel occlusion, disintegration, and new vessel creation. The latter is most commonly referred as tumorinduced angiogenesis, and results in a tortuous, highly irregular vascular network [4, 6, 7].
In [94, 95] the authors observed that cancer cells initially coopt host existing vasculature and grow as well vascularized tumors for several days, up to 2 mm in diameter. No evidence of angiogenesis is observed during this period. Progressively, blood vessels near the tumor core start to regress and/or become occluded, while tumor periphery displays robust angiogenesis. Later work in [96, 97] demonstrated that this pattern repeats itself during later stages of tumor growth; once the tumor grows over a well vascularized region, local vasculature starts to regress. At the same time, tumor periphery displays high angiogenic activity, thereby further promoting tumor growth.
The exact biological mechanisms pertaining to these phenomena are not well understood. Generally, they are attributed to a variety of complex molecular and biomechanical interactions between existing vasculature and tumor cells. It is evident, however, that tumorinduced vascular remodeling affects the local supply of nutrients in tumors and this is where we are going to focus. Taking into consideration the aforementioned evidence, we will try to quantify the spatiotemporal evolution of local oxygen and glucose provision, i.e. the vectors o_b_{t} , gl_b_{t} , o_b_max_{t} and gl_b_max_{t}
Since tumorinduced vascular remodeling occurs either in the interior or in the close vicinity of a tumor, at each time step we consider only the voxels that have already been reached by the tumor, that is, voxels for which l_{t}(A) + nc_{t}(A) > 0. Let A denote such a voxel.
According to the aforementioned evidence, vessel regression should decrease o_b_max_{t}(A) and gl_b_max_{t}(A). As the occupation of A by tumor cells (live or necrotic) increases, the vessel regression rate should also increase, and hence, o_b_max_{t}(A) and gl_b_max_{t}(A) should decrease at a higher rate.
On the other hand, angiogenesis should increase o_b_max_{t}(A) and gl_b_max_{t}(A). A lower occupation of A by live tumor, necrotic tumor and necrotic normal cells, implies a higher angiogenesis rate and hence, a higher rate by which o_b_max_{t}(A) and gl_b_max_{t}(A) increase.
Furthermore, tumor angiogenesis is most commonly associated with nutrient deficit, i.e. when oxygen or glucose quantities fall below certain thresholds known respectively as hypoxia and hypoglycemia thresholds. Typical values for these thresholds are 0.30* o_{0} and 0.50* gl_{0} [5, 37, 98]. We introduce an additional N^{3} × 1 logical vector, denoted by sw_{t} with the following use: at the end of each time step, for each voxel A, the quantities o_{t}(A) and gl_{t}(A) are compared with their respective thresholds; if either of them is below its threshold and the voxel contains live tumor cells, sw_{t}(A) is set to 1, indicating that angiogenesis is on for this voxel. Otherwise, sw_{t}(A) is set to 0.
Let v_{r} and v_{e} denote the maximum rates by which the capacity of the vascular network in A to provide/absorb molecules to/from the surrounding tissue (as modeled by o_b_max_{t}(A) and gl_b_max_{t}(A)) decreases or increases, respectively. The orders of magnitude of the corresponding half and doubling times can be deduced from [94, 95] and are in the order of days.
According to the aforementioned biological evidence and assumptions, a general way to quantify the timely evolution of the macroscopic variables under consideration is
where f_{r}(l_{t}(A), nc_{t}(A)) and f_{e}(l_{t}(A), nc_{t}(A), nn_{t}(A)) are functions taking values in [0, 1]. Function f_{r}(l_{t}(A), nc_{t}(A)) is increasing in both of its arguments. Function f_{e}(l_{t}(A), nc_{t}(A), nn_{t}(A)) is decreasing in all three of its arguments.
These deterministic equations constitute only a rough approximation of the involved dynamics. For a more robust approach, we introduce randomness in them in the following way. Since v_{r} and v_{e} are the maximum rates of vessel regression/expansion, at each time step, and for each voxel A which has been reached by the tumor, we pick two random numbers r_{3} and r_{4} in the interval [0, 1] and introduce the more general, stochastic equations
It remains to choose the functions f_{r} and f_{e}. Since the monotonicity of these functions is determined, it remains to choose their shape, i.e. linear, convex or concave. For the simulations presented below, we have chosen functions of linear shape. Putting everything together, to calculate o_b_max_{t + Δτ}(A) and gl_b_max_{t + Δτ}(A) we will use the equations
We note that algorithmically, before each such calculation, a sanity check should be performed for the each of the quantities \( \frac{l_t(A)+{nc}_t(A)}{M} \) and \( \frac{M{l}_t(A){nc}_t(A){nn}_t(A)}{M} \) ensuring that their values stay respectively below 1 and above 0.
To summarize, to model the effects of tumorinduced vascular remodeling, we introduced two additional vectors, that is, two additional variables for each voxel A, namely o_b_max_{t}(A) and gl_b_max_{t}(A). These variables quantify the maximum values the local provision/absorption of oxygen and glucose may attain, i.e. the maximum absolute values of o_b_{t}(A) and gl_b_{t}(A) reflecting thereby the capacity of the local vascular network to provide/absorb molecules to/from the surrounding tissue. The spatiotemporal evolution of these variables reflects the effects of vessel regression and angiogenesis induced by the tumor. The resulting algorithm applied at each time step for each voxel A follows:
Note: For any positive number μ and xϵ ℝ we will make use of the function
Case 1: If l_{t}(A) = nc_{t}(A) = 0 that is, the tumor has not yet reached A. In this case, the maximum absolute values of oxygen and glucose provision/absorption during the time step t → t + Δτ simply equal the ones during the previous time step, t − Δτ → t:
We pick two random numbers r_{1}, r_{2}, uniformly distributed in [0, 1]. The actual provision/absorption for oxygen and glucose by the vascular system during time interval t → t + Δτ is calculated by
Case 2: If l_{t}(A) + nc_{t}(A) > 0 that is, the tumor has reached A. We first pick two random numbers r_{1}, r_{2}, uniformly distributed in [0, 1] and calculate
Again, we pick two random numbers r_{3}, r_{4} uniformly distributed in [0, 1]. The upper bounds for the absolute values of o_b_{t + Δτ}(A) and gl_b_{t + Δτ}(A) are given by
The provision/absorption for oxygen and glucose by the vascular system during time interval t → t + Δτ is calculated by
VII. The complete model architecture. Modularity and adjustability
The model we propose can be seen as a discrete time dynamical system. The state of the system consists of the nine N^{3} × 1 vectors l_{t}, nc_{t}, nn_{t}, o_{t}, gl_{t}, o_b_{t}, gl_b_{t}, o_b_max_{t} and gl_b_max_{t}
In sections IVI we have defined the following operators.

The operator defined in section VI, which we denote by F_{vr}. Applying this operator to the state vector consists in applying the algorithm described in the last section for each voxel. This operator calculates the supply rate of oxygen and glucose during the time interval t → t + Δτ from the respective values during the previous time interval, t − Δτ → t, as it is dictated by the effects of tumorinduced vascular remodeling.

The operator defined algorithmically in section V, which we denote by F_{pn}. Applying operator F_{pn} to the state vector, consists in checking all cases described in section V and performing the respective calculations for each voxel in the lattice. This operator calculates the proliferation/necrosis of cells in each voxel, as they are dictated by the voxels’ oxygen and glucose levels and supply rates.

The operators defined in section III, which we denote by F_{o}, F_{gl} . Applying each of these operators to the state vector consists in multiplying the matrices (T_{o} or T_{gl}) with the respective vector o_{t} or gl_{t}, thereby calculating how chemical fields change due to diffusion.

The operator defined algorithmically in section IV, which we denote by F_{c}. Applying the operator F_{c} to the state vector consists in executing Algorithm 2, calculating thereby how cancer cell populations within each voxel change due to cell diffusion.
Knowledge of the state vector at time t allows us to calculate the state vector at time t + Δτ, by applying the algorithm depicted in Fig. 2.
Phenomena pertaining to tumorinduced vascular remodeling, nutrient consumption, cell proliferation and cellular or molecular diffusion are modeled by separate operators (i.e. algorithmic modules), applied sequentially to the state vector. Within the proposed methodology, the algorithmic modules corresponding to operators F_{pn}, F_{vr} are completely readjustable. This facilitates the simulation of scenarios based on different hypotheses concerning the effects of tumorinduced vascular remodeling on nutrient supply rates, cell proliferation, necrosis and metabolism of chemical species. A variety of choices is also available for readjusting F_{c}; we provide some suggestions in the appendix. Introduction of additional diffusion operators like F_{o} and F_{gl} and extension of operators F_{pn} and F_{vr} enables the consideration of additional chemical species such as lactate, growth factors and chemotherapeutic agents. Gradual removal of necrotic cells from the tumor mass may (and should) also be considered. Introduction of additional cellular species is also feasible, by considering additional cellular diffusion operators and appropriate readjustment of F_{pn}.
Of note, there is a large disparity between the time scales of chemical and cellular diffusion, with the latter evolving much more slowly. The diffusion coefficients of oxygen and glucose are in the order of 10^{−5} cm^{2}/sec while the respective coefficient for tumor cells is in the order of 10^{−8} cm^{2}/sec. Furthermore, tumor cell diffusion in neighboring voxels is also affected by their proliferation. This allows applying F_{c} to the state vector every several time steps κ; for the simulations presented in the next section, we used a time step Δτ= 10 s and κ= 30 (i.e. 5 mins).
The methodology we described in the previous sections essentially casts the problem of modeling tumor growth in the realm of spatially distributed, stochastic dynamical systems. The state of the system evolves according to a law of the form x_{k + 1} = f(x_{k}) where the transition x_{k} → x_{k + 1} is stochastic and x has an additional spatial structure.
The model was implemented by making extensive use of MATLAB’s vectorized approach to coding. We note that apart from the operators F_{o}, F_{gl} and F_{c}, this vectorized implementation is feasible also for the operators F_{vr} and F_{pn}. However, to keep the code readable, in this work we have opted to implement F_{vr} and F_{pn} using loops. Although in the present work we did not exploit multicore computation, it is clear that each of the aforementioned algorithmic operators is eligible for parallel implementation. Furthermore, vectorized implementation of the resulting algorithmic modules opens the road for exploiting the capabilities of modern tools like Python’s Numba compiler or TensorFlow for computation on GPUs. This will facilitate simulations over larger tissue areas with finer spatial and temporal resolutions, plus, importantly, a comparative analysis of the numerical error induced by the discretization parameters. We note that such an analysis has not been performed yet, since it requires the consideration of more and smaller values for Δs and Δt. This, however, increases significantly the computational burden of each simulation, and requires a completely different implementation of the model in terms of programming. It is therefore left for future work.
VIII. Simulations and use cases
In this section, we use the model developed in the previous sections for a theoretical study of tumor growth, consisting of two parts. First, we use the model to visualize various aspects of tumor expansion. The resulting images are qualitatively compared with pertinent biological observations. We then perform a multivariate analysis regarding the effects of a subset of model parameters on the number of live cancer cells after a certain period of free growth.
The proposed model can be used for visualizing various phenomena encountered during the expansion of a tumor. Figures 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12, which we explain below, depict a series of such examples. These images are snapshots of a simulation with the following parameter values (see also sections about metabolism and tumorinduced vascular remodeling above). Maximum mitosis rate for cancer cells a_{max} was set such that their minimum doubling time is 5 days. Parameter lambda is set λ = 10.
Maximum vasculature regression rate (v_{r}) was set such that the corresponding minimum halftime time is 5 days. Maximum vasculature expansion rate (v_{e}) was set such that the corresponding minimum doubling time is 1 day. A spatially constant, synthetic diffusion tensor was used for the diffusion of live cancer cells, defined by the orthonormal vectors
The corresponding quantities \( \sqrt{\alpha {\lambda}_1},\sqrt{\alpha {\lambda}_2} \) and \( \sqrt{\alpha {\lambda}_3} \) appearing in equation (5) were calculated from the equation 2D = aλ derived in [80]. The value range for the cell diffusion coefficient D was chosen according to the estimations made in [99]. The quantity \( \sqrt{\alpha {\lambda}_1} \) was set much higher than \( \sqrt{\alpha {\lambda}_2} \) and \( \sqrt{\alpha {\lambda}_3} \), such that the diffusion of live cancer cells occurs primarily along the first vector. Table 1 summarizes parameter values used for the simulations of this section. An initial population of 5 ∙ 10^{5} live cancer cells was placed in the center of the simulated region.
Figs. 3 and 4 depict the tumor at the 70th and 90th day of the simulation, respectively. A general observation, holding for all performed simulations was that the chosen cell diffusion tensor affects the overall tumor shape in a noticeable way. This is depicted in Fig. 3, where the tumor shape is roughly similar to the ellipsoid of the diffusion tensor defined above. It can also be seen in Fig. 4, where the depicted tumor appears to be an approximately conformal expansion of the tumor in Fig. 3, although markedly distorted by the underlying stochasticity.
In Figs. 5 and 6, vertical sections of these tumors are shown, taken at the central (11th) voxel plane of the lattice. All colored voxels have been reached by the tumor, i.e. they contain a nonzero population of live plus necrotic cancer cells. We do not show voxels containing only host necrotic cells, something that was often observed in the near vicinity of the tumor periphery. We remind the reader that M is the average cell population capacity per voxel. The color code is as follows. Cyan voxels contain a total population of live tumor, necrotic tumor and necrotic host cells that does not exceed 50% of M. Where this quantity is above 50% of M, voxels are colored blue, gray or black, depending on the amount of necrotic cells they contain. Specifically, in blue voxels, the total necrotic (tumor+host) cell population is below 65% of M. In gray voxels, the total necrotic cell population is between 65 and 95% of M. In black voxels, necrotic cells are above 95% of M. Note that in both Figs. 5 and 6 necrosis is higher towards the tumor center. This is in agreement with the general observation that after some period of growth, due to vasculature disorganization and limited diffusion of nutrients near the tumor center, a necrotic core is formed; viable proliferating cells are located mainly at the outer rims of a tumor.
Figs. 7 and 8 depict a profile of the oxygen levels per voxel, corresponding to the vertical sections shown in Figs. 5 and 6, respectively. The darker the voxel shade, the lower the oxygen quantity in the specific voxel. The lightest shade indicates a voxel whose oxygen level is at least equal to normal. The darkest shade indicates a voxel whose oxygen level is below 15% of normal tissue level. The same color code holds for Figs. 9 and 10, which depict the respective glucose levels. Note that the tumor interior contains both hypoxic and hypoglycemic regions. From the simulations performed, hypoxic regions appeared to be much more spatially inhomogeneous and time varying than hypoglycemic ones, which generally tended to appear more congruent with necrotic areas of the tumor and overall tumor shape. Other than that, the shape of both hypoxic and hypoglycemic regions appeared to be largely random, and no pertinent spatial patterns were detected.
Similarly, Figs. 11 and 12 depict a similar profile of o_b_max_{t} for each voxel, corresponding respectively to the vertical sections shown in Figures 5 and 6. We remind the reader that, for a voxel A, o_b_max_{t}(A) quantifies the capacity of the vascular network within A to provide/absorb oxygen to/from the surrounding tissue, thereby reflecting the regressed or expanded vasculature in A (see also section VI). The color code for Figs. 11 and 12 is as follows: The darkest grayscale shade indicates that for the respective voxel A, the value of o_b_max_{t}(A) is below 50% of its value in normal tissue. The lightest grayscale shade indicates that for the respective voxel A, the value of o_b_max_{t}(A) is equal to its value in normal tissue. Magenta indicates that o_b_max_{t}(A) is between 100 and 125% of its value in normal tissue. Blue indicates that o_b _max_{t}(A) is between 125 and 150% of its value in normal tissue. We note that since o_b_max_{t} and gl_b_max_{t} evolve temporally in the same way (see section VI), using the same color code yields identical figures for gl_b_max_{t}. Figures 11 and 12 depict that towards the tumor center, o_b_max_{t} is lower, i.e. vasculature appears to be more regressed. Tumor periphery displays values larger than normal, reflecting the fact that at the outer rims of the tumor, angiogenesis takes place at a much faster rate than vessel regression.
As a more general use case, we performed a multivariate analysis on the effects of certain tumor growth related parameters on the number of viable tumor cells after a period of free growth. Specifically:

Maximum mitosis rate for cancer cells a_{max} was fixed such that the corresponding minimum doubling time is 5 days.

Parameter λ was varied in the set {2,4,6,8,10}. We remind the reader that λK_{ATP} is the ATP consumption rate of an actively proliferating tumor cell, where K_{ATP} is the ATP consumption rate of a normal host cell.

The parameters β_{1} and β_{2}, i.e. the minimum and maximum values of the energy fraction cancer cells acquire by glycolysis were varied in the set [β_{1} β_{2}] ={[0 0.1], [0.1 0.2], [0.2 0.3], [0.3 0.4], [0.4 0.5]}

Vasculature regression minimum halftime (in days) was varied in {1,2,3,4,5}.

Vasculature expansion minimum doubling time (in days) was varied in {1,2,3,4,5}. The case where no angiogenesis occurred throughout the entire simulation was also considered.
The initial tumor population was 5 ∙ 10^{5} live cancer cells. The evolution of these tumors was simulated for a time period of 90 days. Due to the stochasticity of the model, for each set of parameter values a total of 20 simulations was performed. The results of these simulations are shown in Figs. 13, 14, 15, 16, 17, and 18.
In each of the Figs. 13, 14, 15, 16, 17, and 18, the columns ‘lambda’,‘Vasculature regression minimum halftime’ and ‘Vasculature expansion minimum doubling time’ are selfexplanatory. Each row of the column ‘Live cancer cells after 90 days’ depicts the respective numbers of live tumor cells for each one of the 20 simulations performed for the parameters specified in the previous columns of the same row. In each row, these 20 numbers are drawn as horizontal lines, each one with length proportional to the resulting number of live cancer cells after 90 days. These 20 lines are drawn in sorted order with regard to their length, from longest (top) to shortest (bottom), and form the skewed bar observed in each row of the ‘Live cancer cells after 90 days’ column. We note that due to the smaller final populations observed absent angiogenesis, the aforementioned lines in Fig. 13 are drawn in a different scale than in Figs. 14, 15, 16, 17, and 18. In each row of the column ‘Probability of tumor survival’, we provide the fraction of simulations (out of the 20 simulations performed for the parameters of the specific row) that resulted to a nonzero population of live cancer cells. Each row of the column ‘Expected number of live cancer cells’ provides the expected number of live cancer cells after 90 days, conditioned on the survival of the tumor. This is essentially the mean of the resulting final populations, calculated by taking into account only the simulations that resulted in a nonzero final population of live cancer cells.
As explained above, each row of the ‘Live cancer cells after 90 days’ column depicts a skewed bar, formed by the final populations of live cancer cells for each of the 20 simulations performed for the specific row. Thus, the skewness of each such bar indicates the variance observed in the results of the respective simulations. For each such bar, a flatter right end indicates a lower variance; a more skewed right end indicates a higher variance. Visual inspection of these bars indicates that in most cases (i.e. rows), 20 simulations can provide a reasonable overview of the potential outcomes. However, there are cases where the skewness of the aforementioned bars is quite high, indicating a higher variance in the potential outcomes of the respective simulations. In such cases, like for example, the case in Fig. 13 where 0.1 ≤ β ≤ 0.2, lambda = 2 and Vasculature regression minimum halftime = 4 days or the same case but with lambda = 4, it is evident more simulations are needed. In general, our analysis showed that simulation parameters have an effect not only on the probability of survival and expected populations of cancer cells, but also, in several cases, on the variance of these populations. This was also observed for intermediate time points, i.e. cell populations calculated at time points within the overall time frame of 90 days.
Fig. 13 depicts the results of these simulations when tumors grow without angiogenesis. We see that the limits of the energy fraction β cancer cells can acquire from glycolysis play a crucial role on tumor growth and survival. For 0 ≤ β ≤ 0.1, i.e. when cells employ mainly combustion of glucose, tumors survive essentially only if they have minimal energy needs compared to host cells, and additionally, if vasculature regression evolves at a minimal rate. For 0.1 ≤ β ≤ 0.2, tumors survive in more cases than for 0 ≤ β ≤ 0.1 and in these cases they reach large end populations of viable cancer cells. For 0.2 ≤ β ≤ 0.3 tumors survive in even more cases, but don’t reach end populations as large as for 0.1 ≤ β ≤ 0.2. For 0.3 ≤ β ≤ 0.4 tumors survive almost like in the case where 0.2 ≤ β ≤ 0.3, although with lower probabilities and lower end populations. The same trend is observed when moving to the last case; for 0.4 ≤ β ≤ 0.5 tumors survive in slightly less cases than for 0.3 ≤ β ≤ 0.4, with lower probabilities and lower end populations. In each separate case, maximum vasculature regression rate (i.e. minimum halftime) and cancer cell energy requirements affect both probability of survival and final number of viable cells. In fact, in each case, the higher the maximum rate of vasculature regression, the lower are both the survival probability and the viable end population. Furthermore, in each case, for the same maximum rate of vasculature regression, higher energy requirements by tumor cells imply lower survival probabilities and viable end populations.
Figures 14, 15, 16, 17, and 18 consider additionally the effects of vasculature expansion. Visual inspection of these Figures shows that the limits of β affect the results in a way reminiscent of the one observed for the no angiogenesis case in Fig. 13. In Fig. 14, where 0 ≤ β ≤ 0.1, i.e. energy is acquired mainly through combustion, tumors grow only if vasculature regression evolves sufficiently slower than angiogenesis. In Fig. 15 (0.1 ≤ β ≤ 0.2), tumors survive and grow in many more cases, and in general they reach larger end populations of viable cells. In Fig. 16 (0.2 ≤ β ≤ 0.3), tumors survive in even more cases, however, end populations are lower than in Fig. 15. In terms of survival probabilities, tumors in Fig. 17 (0.3 ≤ β ≤ 0.4) tumors are slightly better than in Fig. 16, achieve, nevertheless, lower end populations. In Fig. 18 (0.4 ≤ β ≤ 0.5) tumors survive almost like in Fig. 17. Again, compared to Fig. 17, survival probabilities and end populations are lower.
A general observation is that, for the same limits of β and the same energy requirements λ, maximum rates of vasculature expansion and regression had monotonic effects on survival probabilities and viable end populations. Assuming other parameters equal, a higher maximum rate of vasculature expansion generally implies higher survival probability and viable end population. Respectively, a higher maximum rate of vasculature regression generally has the opposite effects. There are a few exceptions in these rules; they are marked by arrows at the left of each image and we will explain them below.
For all of the aforementioned exceptions, the expected number of live cancer cells was calculated for each day. These quantities were plotted in pairs for each case monotonicity did not hold; three such examples are shown in Figs. 19, 20, and 21. These Figures depict the pattern observed for all these cases. The expected behavior of tumor cells growing in more favorable conditions (blue curves) is to grow faster and larger than their respective counterparts (red curves) for most of the observed time period. However, they reach a maximum and start to regress sooner than tumor cells growing in less favorable conditions. Apparently though, on the 90th day, the sum of viable and necrotic tumor cells is larger for tumors growing in more favorable conditions.
The effects of λ (i.e. the parameter quantifying the energy requirements of cancer cells compared to host cells) are much more complicated. Using the abbreviations Vasculature Regression minimum Halftime (VRmHt) and Vasculature Expansion minimum Doubling time (VEmDt) defined for Figs. 19, 20 and 21, we observe the following. For 0 ≤ β ≤ 0.1, VRmHt = 5 days, VEmDt = 1 day, λ has an increasing effect on the viable end population. Apparently, in this case, high energy requirements induce hypoxia and hypoglycemia much more frequently, thereby triggering vessel expansion more often. Additionally, tumorinduced angiogenesis is much faster than tumorinduced vasculature regression; this results in higher final populations of viable cancer cells. The observation that, if angiogenesis occurs sufficiently faster than vascular regression, higher energy requirements act increasingly on the viable end population holds in general for 0 ≤ β ≤ 0.3, see Figs. 14, 15, and 16 . It is much less pronounced for 0.4 ≤ β ≤ 0.5 (Fig. 18), and does not hold for 0.3 ≤ β ≤ 0.4 (Fig. 17). As a counterexample, note that for 0.2 ≤ β ≤ 0.3, VRmHt = 5 days, VEmDt = 3 days, λ has a decreasing effect on the viable end population.
From the analysis performed so far, no general pattern providing a quantitative and concise interpretation of these observations was found. Evidently, the dependencies of the end population of viable cells on the parameters under consideration are quite complex, and display a rich structure of local maxima and minima. A more detailed study of these dependencies would require more simulations for each set of parameter values and additional consideration of standard deviations; we leave this for future work.
Availability of data and materials
Data sharing is not applicable to this article as no datasets were generated or analysed during the current study. All parameter values required to reproduce the presented simulations are included in this published article.
Abbreviations
 DTI:

Diffusion Tensor Imaging
 PET:

Positron Emission Tomography
 VEmDt:

Vasculature Expansion minimum Doubling time
 VRmHt:

Vasculature Regression minimum Halftime
References
 1.
Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P. Molecular Biology of the Cell. 5th ed. New York: Garland Science, Taylor & Francis Group; 2007.
 2.
Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144(5):646–74.
 3.
Michor F, Iwasa Y, Nowak MA. Dynamics of cancer progression. Nat Rev Cancer. 2004;4(3):197–205.
 4.
Vaupel P, Kallinowski F, Okunieff P. Blood flow, oxygen and nutrient supply, and metabolic microenvironment of human tumors: a review. Cancer Res. 1989;49(23):6449–65.
 5.
Brown JM, Wilson WR. Exploiting tumour hypoxia in cancer treatment. Nat Rev Cancer. 2004;4(6):437–47.
 6.
Carmeliet P, Jain RK. Angiogenesis in cancer and other diseases. Nature. 2000;407(6801):249–57.
 7.
Jain RK, di Tomaso E, Duda DG, Loeffler JS, Sorensen AG, Batchelor TT. Angiogenesis in brain tumours. Nat Rev Neurosci. 2007;8(8):610–22.
 8.
Anderson AR, Chaplain MA. Continuous and discrete mathematical models of tumorinduced angiogenesis. Bull Math Biol. 1998;60(5):857–99.
 9.
Hirst DG, Denekamp J. Tumour cell proliferation in relation to the vasculature. Cell Tissue Kinet. 1979;12(1):31–42.
 10.
Sutherland RM. Cell and environment interactions in tumor microregions: the multicell spheroid model. Science. 1988;240(4849):177–84.
 11.
Tannock I. F The relation between cell proliferation and the vascular system in a transplanted mouse mammary tumour. Br J Cancer. 1968;22(2):258–73.
 12.
Casciari JJ, Sotirchos SV, Sutherland RM. Glucose diffusivity in multicellular tumor spheroids. Cancer Res. 1988;48(14):3905–9.
 13.
Freyer JP, Sutherland RM. Proliferative and clonogenic heterogeneity of cells from EMT6/Ro multicellular spheroids induced by the glucose and oxygen supply. Cancer Res. 1986;46(7):3513–20.
 14.
Schaller G, MeyerHermann M. Multicellular tumor spheroid in an offlattice VoronoiDelaunay cell model. Phys Rev E Stat Nonlin Soft Matter Phys. 2005;71(5 Pt 1):051910.
 15.
Vavourakis V, Wijeratne PA, Shipley R, Loizidou M, Stylianopoulos T, Hawkes DJ. A Validated Multiscale InSilico Model for Mechanosensitive Tumour Angiogenesis and Growth. PLoS Comput Biol. 2017;13(1):e1005259. Published 2017 Jan 26. https://doi.org/10.1371/journal.pcbi.1005259.
 16.
López AG, Seoane JM, Sanjuán MA. A validated mathematical model of tumor growth including tumorhost interaction, cellmediated immune response and chemotherapy. Bull Math Biol. 2014;76(11):2884–906. https://doi.org/10.1007/s1153801400375 Epub 2014 Oct 28.
 17.
Collis J, Connor AJ, Paczkowski M, Kannan P, PittFrancis J, Byrne HM, Hubbard ME. Bayesian Calibration, Validation and Uncertainty Quantification for Predictive Modelling of Tumour Growth: A Tutorial. Bull Math Biol. 2017;79(4):939–74. https://doi.org/10.1007/s1153801702585 Epub 2017 Mar 13.
 18.
Portz T, Kuang Y, Nagy JD. A clinical data validated mathematical model of prostate cancer growth under intermittent androgen suppression therapy. AIP Advances. 2012;2:011002. https://doi.org/10.1063/1.3697848.
 19.
Wodarz D, Komarova NL. Dynamics of Cancer: Mathematical Foundations of Oncology. 1st ed. River Edge: World Scientific Publishing Co., Inc.; 2014.
 20.
Komarova NL, Burger JA, Wodarz D. Evolution of ibrutinib resistance in chronic lymphocytic leukemia (CLL). Proc Natl Acad Sci USA. 2014;111(38):13906–11.
 21.
RodriguezBrenes IA, Wodarz D. Preventing clonal evolutionary processes in cancer: Insights from mathematical models. Proc Natl Acad Sci USA. 2015;112(29):8843–50.
 22.
Beerenwinkel N, Schwarz RF, Gerstung M, Markowetz F. Cancer Evolution: Mathematical Models and Computational Inference. Syst Biol. 2015;64(1):e1–e25.
 23.
De Pillis LG, Radunskaya AE, Wiseman CL. A validated mathematical model of cellmediated immune response to tumor growth. Cancer Res. 2005;65(17):7950–8.
 24.
Leon K, Garcia K, Carneiro J, Lage A. How regulatory CD25(+)CD4(+) T cells impinge on tumor immunobiology? On the existence of two alternative dynamical classes of tumors. J Theor Biol. 2007;247(1):122–37.
 25.
RobertsonTessi M, ElKareh A, Goriely A. A mathematical model of tumorimmune interactions. J Theor Biol. 2012;294:56–73.
 26.
De Pillis LG, Radunskaya A. A Mathematical Tumor Model with Immune Resistance and Drug Therapy: An Optimal Control Approach. J Theor Med. 2001;3(2):79–100.
 27.
Castorina P, Carcò D, Guiot C, Deisboeck TS. Tumor growth instability and its implications for chemotherapy. Cancer Res. 2009;69(21):8507–15.
 28.
Stura I, Venturino E, Guiot C. A twoclones tumor model: Spontaneous growth and response to treatment. Math Biosci. 2016;271:19–28.
 29.
Forouzannia F, Enderling H, Kohandel M. Mathematical Modeling of the Effects of Tumor Heterogeneity on the Efficiency of Radiation Treatment Schedule. Bull Math Biol. 2018;80(2):283–93.
 30.
Patel AA, Gawlinski ET, Lemieux SK, Gatenby RA. A cellular automaton model of early tumor growth and invasion. J Theor Biol. 2001;213(3):315–31.
 31.
Alarcón T, Byrne HM, Maini PK. A cellular automaton model for tumour growth in inhomogeneous environment. J Theor Biol. 2003;225(2):257–74.
 32.
Alarcón T, Byrne HM, Maini PK. A Multiple Scale Model for Tumor Growth. Multiscale Model Simul. 2005;2:440–75.
 33.
Anderson AR. A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. Math Med Biol. 2005;22(2):163–86.
 34.
Anderson AR, Weaver AM, Cummings PT, Quaranta V. Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. Cell. 2006;127(5):905–15.
 35.
Wang Z, Zhang L, Sagotsky J, Deisboeck TS. Simulating nonsmall cell lung cancer with a multiscale agentbased model. Theor Biol Med Model. 2007;4:50.
 36.
Gerlee P, Anderson AR. An evolutionary hybrid cellular automaton model of solid tumour growth. J Theor Biol. 2007;246(4):583–603.
 37.
Gerlee P, Anderson AR. A hybrid cellular automaton model of clonal evolution in cancer: the emergence of the glycolytic phenotype. J Theor Biol. 2008;250(4):705–22.
 38.
Wang Z, Birch CM, Deisboeck TS. Crossscale sensitivity analysis of a nonsmall cell lung cancer model: linking molecular signaling properties to cellular behavior. Biosystems. 2008;92(3):249–58.
 39.
VitalLopez FG, Armaou A, Hutnik M, Maranas CD. Modeling the effect of chemotaxis on glioblastoma tumor progression. AIChE J. 2011;57:778–92.
 40.
Perfahl H, Byrne HM, Chen T, et al. Multiscale Modelling of Vascular Tumour Growth in 3D: The Roles of Domain Size and Boundary Conditions. Secomb TW, ed. PLoS One. 2011;6(4):e14790.
 41.
Li F, Tan H, Singh J, Yang J, Xia X, Bao J, Ma J, Zhan M, Wong TCS. A 3D multiscale model of cancer stem cell in tumor development. BMC Syst Biol. 2013;7(Suppl 2):S12.
 42.
Haridas P, Browning AP, McGovern JA, Sean McElwain DL, Simpson MJ. Threedimensional experiments and individual based simulations show that cell proliferation drives melanoma nest formation in human skin tissue. BMC Syst Biol. 2018;12:34.
 43.
Kansal AR, Torquato S, Harsh GR IV, Chiocca EA, Deisboeck TS. Cellular automaton of idealized brain tumor growth dynamics. Biosyst. 2000;55(1–3):119–27.
 44.
Macklin P, Edgerton ME, Thompson AM, Cristini V. Patientcalibrated agentbased modelling of ductal carcinoma in situ (DCIS): from microscopic measurements to macroscopic predictions of clinical progression. J Theor Biol. 2012;301:122–40.
 45.
Kempf H, Bleicher M, MeyerHermann M. SpatioTemporal Dynamics of Hypoxia during Radiotherapy. PLoS One. 2015;10(8):e0133357.
 46.
Jiang Y, PjesivacGrbovic J, Cantrell C, Freyer JP. A Multiscale Model for Avascular Tumor Growth. Biophys J. 2005;89(6):3884–94.
 47.
Rubenstein BM, Kaufman LJ. The role of extracellular matrix in glioma invasion: a cellular Potts model approach. Biophys J. 2008;95(12):5661–80.
 48.
Shirinifard A, Gens JS, Zaitlen BL, Popławski NJ, Swat M, Glazier JA. 3D multicell simulation of tumor growth and angiogenesis. PLoS One. 2009;4(10):e7190.
 49.
Szabó A, Merks RM. Cellular potts modeling of tumor growth, tumor invasion, and tumor evolution. Front Oncol. 2013;3:87.
 50.
Jeanquartier F, JeanQuartier C, Cemernek D, Holzinger A. In silico modeling for tumor growth visualization. BMC Syst Biol. 2016;10(1):59.
 51.
Ward JP, King JR. Mathematical modelling of avasculartumour growth. IMA J Math Appl Med Biol. 1997;14(1):39–69.
 52.
Venkatasubramanian R, Henson MA, Forbes NS. Incorporating energy metabolism into a growth model of multicellular tumor spheroids. J Theor Biol. 2006;242(2):440–53.
 53.
Schaller G, MeyerHermann M. Continuum versus discrete model: a comparison for multicellular tumour spheroids. Philos Trans A Math Phys Eng Sci. 2006;364(1843):1443–64.
 54.
Stein AM, Demuth T, Mobley D, Berens M, Sander LM. A Mathematical Model of Glioblastoma Tumor Spheroid Invasion in a ThreeDimensional In Vitro Experiment. Biophys J. 2007;92(1):356–65.
 55.
Swanson KR, Bridge C, Murray JD, Alvord EC Jr. Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. J Neurol Sci. 2003;216(1):1–10.
 56.
Jbabdi S, Mandonnet E, Duffau H, Capelle L, Swanson KR, PélégriniIssac M, Guillevin R, Benali H. Simulation of anisotropic growth of lowgrade gliomas using diffusion tensor imaging. Magn Reson Med. 2005;54(3):616–24.
 57.
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 Imaging. 2005;24(10):1334–46.
 58.
Rockne R, Rockhill JK, Mrugala M, Spence AM, Kalet I, Hendrickson K, Lai A, Cloughesy T, Alvord EC Jr, Swanson KR. Predicting efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach. Phys Med Biol. 2010;55(12):3271–85.
 59.
Painter KJ, Hillen T. Mathematical modelling of glioma growth: the use of Diffusion Tensor Imaging (DTI) data to predict the anisotropic pathways of cancer invasion. J Theor Biol. 2013;323:25–39.
 60.
Patel V, Hathout L. Imagedriven modeling of the proliferation and necrosis of glioblastoma multiforme. Theor Biol Med Model. 2017;14(1):10.
 61.
Swan A, Hillen T, Bowman JC, Murtha AD. A PatientSpecific Anisotropic Diffusion Model for Brain Tumour Spread. Bull Math Biol. 2018;80(5):1259–91.
 62.
Cristini V, Lowengrub J, Nie Q. Nonlinear simulation of tumor growth. J Math Biol. 2003;46(3):191–224.
 63.
Li X, Cristini V, Nie Q, Lowengrub JS. Nonlinear threedimensional simulation of solid tumor growth. Discrete Control Dyn Syst. 2007;7(3):581–604.
 64.
Zheng X, Wise SM, Cristini V. Nonlinear simulation of tumor necrosis, neovascularization and tissue invasion via an adaptive finiteelement/levelset method. Bull Math Biol. 2005;67(2):211–59.
 65.
Hogea CS, Murray BT, Sethian JA. Simulating complex tumor dynamics from avascular to vascular growth using a general levelset method. J Math Biol. 2006;53(1):86–134.
 66.
Macklin P, Lowengrub J. Nonlinear simulation of the effect of microenvironment on tumor growth. J Theor Biol. 2007;245(4):677–704.
 67.
Macklin P, Lowengrub JS. A New Ghost Cell/Level Set Method for Moving Boundary Problems: Application to Tumor Growth. J Sci Commun. 2008;35(2–3):266–99.
 68.
Macklin P, McDougall S, Anderson AR, Chaplain MA, Cristini V, Lowengrub J. Multiscale modelling and nonlinear simulation of vascular tumour growth. J Math Biol. 2009;58(4–5):765–98.
 69.
Wise SM, Lowengrub JS, Frieboes HB, Cristini V. Threedimensional multispecies nonlinear tumor growthI Model and numerical method. J Theor Biol. 2008;253(3):524–43.
 70.
Cristini V, Li X, Lowengrub JS, Wise SM. Nonlinear simulations of solid tumor growth using a mixture model: invasion and branching. J Math Biol. 2009;58(4–5):723–63.
 71.
Frieboes HB, Jin F, Chuang YL, Wise SM, Lowengrub JS, Cristini V. Threedimensional multispecies nonlinear tumor growthII: Tumor invasion and angiogenesis. J Theor Biol. 2010;264(4):1254–78.
 72.
Sciumè G, Shelton S, Gray W, Miller C, Hussain F, Ferrari M, Decuzzi P, Schrefler B. A multiphase model for threedimensional tumor growth. New J Phys. 2013;15:015005.
 73.
Chen Y, Wise SM, Shenoy VB, Lowengrub JS. A stable scheme for a nonlinear, multiphase tumor growth model with an elastic membrane. Int J Numer Method Biomed Eng. 2014 Jul;30(7):726–54.
 74.
Dionysiou DD, Stamatakos GS, Uzunoglu NK, Nikita KS, Marioli A. A fourdimensional simulation model of tumour response to radiotherapy in vivo: parametric validation considering radiosensitivity, genetic profile and fractionation. J Theor Biol. 2004;230(1):1–20.
 75.
Kolokotroni EA, Dionysiou DD, Uzunoglu NK, Stamatakos GS. Studying the growth kinetics of untreated clinical tumors by using an advanced discrete simulation model. Math Comput Model. 2011;54:1989–2006.
 76.
Kolokotroni E, Dionysiou D, Veith C, Kim YJ, Sabczynski J, Franz A, Grgic A, Palm J, Bohle RM, Stamatakos G. In Silico Oncology: Quantification of the In Vivo Antitumor Efficacy of CisplatinBased Doublet Therapy in NonSmall Cell Lung Cancer (NSCLC) through a Multiscale Mechanistic Model. PLoS Comput Biol. 2016;12(9):e1005093.
 77.
Kim Y, Stolarska MA, Othmer HG. A Hybrid Model for Tumor Spheroid Growth in vitro I: Theoretical Development and Early Results. Math Models Methods Appl Sci. 2007;17:1773–98.
 78.
Kim Y, Othmer H. Hybrid models of cell and tissue dynamics in tumor growth. Math Biosci Eng. 2015;12(6):1141–56.
 79.
Deisboeck TS, Stamatakos GS. Multiscale Cancer Modeling. Boca Raton: Taylor & Francis; 2011.
 80.
Antonopoulos M, Stamatakos G. In Silico NeuroOncology: Brownian MotionBased Mathematical Treatment as a Potential Platform for Modeling the Infiltration of Glioma Cells into Normal Brain Tissue. Cancer Inform. 2015;14(Suppl 4):33.
 81.
Grigoriu M. Stochastic Calculus: Applications in Science and Engineering (Chap. 7). Birkhäuser: Switzerland; 2003.
 82.
Stamatakos G, Giatili S. A Numerical Handling of the Boundary Conditions Imposed by the Skull on an Inhomogeneous Diffusion Reaction Model of Glioblastoma Invasion Into the Brain: Clinical Validation Aspects. Cancer Informat. 2017;16(16):1–16.
 83.
Helmlinger G, Yuan F, Dellian M, Jain RK. Interstitial pH and pO2 gradients in solid tumors in vivo: highresolution measurements reveal a lack of correlation. Nat Med. 1997;3(2):177–82.
 84.
Freyer JP, Sutherland RM. Regulation of growth saturation and development of necrosis in EMT6/Ro multicellular spheroids by the glucose and oxygen supply. Cancer Res. 1986;46(7):3504–12.
 85.
Freyer JP, Tustanoff E, Franko AJ, Sutherland RM. In situ oxygen consumption rates of cells in V79 multicellular spheroids during growth. J Cell Physiol. 1984;118(1):53–61.
 86.
Casciari JJ, Sotirchos SV, Sutherland RM. Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular pH. J Cell Physiol. 1992;151(2):386–94.
 87.
KunzSchughart LA, Doetsch J, MuellerKlieser W, Groebe K. Proliferative activity and tumorigenic conversion: Impact on cellular metabolism in 3d culture. Am J Physio Cell Physiol. 2000;278:765–80.
 88.
Wehrle JP, Ng CE, McGovern KA, Aiken NR, Shungu DC, Chance EM, Glickson JD. Metabolism of alternative substrates and the bioenergetic status of EMT6 tumor cell spheroids. NMR in Biomed. 2000;13:349–460.
 89.
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(3):516–24.
 90.
Gatenby RA, Gillies RJ. Why do cancers have high aerobic glycolysis? Nat Rev Cancer. 2004;4(11):891–9.
 91.
Vander Heiden MG, Cantley LC, Thompson CB. Understanding the Warburg effect: the metabolic requirements of cell proliferation. Science. 2009;324(5930):1029–33.
 92.
Despopoulos A, Silbernagl S. Color Atlas of Physiology (6th edition). Stuttgart: Thieme; 2003.
 93.
Wagner BA, Venkataraman S, Buettner GR. The rate of oxygen utilization by cells. Free Radic Biol Med. 2011;51(3):700–12.
 94.
Holash J, Maisonpierre PC, Compton D, Boland P, Alexander CR, Zagzag D, Yancopoulos GD, Wiegand SJ. Vessel cooption, regression, and growth in tumors mediated by angiopoietins and VEGF. Science. 1999;284(5422):1994–8.
 95.
Holash J, Wiegand SJ, Yancopoulos GD. New model of tumor angiogenesis: dynamic balance between vessel regression and growth mediated by angiopoietins and VEGF. Oncogene. 1999;18(38):5356–62.
 96.
Lee DS, Rieger H, Bartha K. Flow correlated percolation during vascular remodeling in growing tumors. Phys Rev Lett. 2006;96(5):058104.
 97.
Bartha K, Rieger H. Vascular network remodeling via vessel cooption, regression and growth in tumors. J Theor Biol. 2006;241(4):903–18.
 98.
Ganong W. Review of Medical Physiology. 19th ed. New York: Appleton & Lange; 1999.
 99.
Murray JD. Mathematical Biology II: Spatial Models and Biomedical Applications (3rd edition). New York: SpringerVerlag; 2011.
Acknowledgements
Fruitful discussions with Christos Kyroudis, Eleftherios Ouzounoglou, Norbert Graf and Stefaan Van Gool are dully acknowledged.
Funding
This work was funded in part by the CHIC project (grant agreement No 600841). The funding body did not play any role in the design of the study, the collection, analysis, and interpretation of data, and in writing the manuscript.
Author information
Affiliations
Contributions
MA conceived the ideas, developed the respective mathematical arguments and algorithms and wrote the first draft of the manuscript. DD, GS, NU provided feedback, reviewed and corrected the manuscript. MA, DD, GS, NU have read and approved the final manuscript.
Corresponding author
Correspondence to Markos Antonopoulos.
Ethics declarations
Ethics approval and consent to participate
Not applicable
Consent for publication
Not applicable
Competing interests
The authors declare that they have no competing interests of any kind.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix
Appendix
IX.i The random perturbation of matrices T_{o} and T_{gl} described in part III of the methods section is implemented as follows. There are two parameters involved

κ, that is, the number of time steps every which the perturbation is applied to the matrices T_{o} and T_{gl} calculated by algorithm 2.

μ, which is a number ranging from 0 to 9.
Every κ time steps, the respective columns of matrices T_{o} and T_{gl} are perturbed in the following way. Each nonzero probability of the column to be perturbed is multiplied by a uniformly distributed random number ranging from 1 to 1 + μ. There are 27 nonzero probabilities in each column, hence, we use 27 such random numbers. We then normalize these 27 products to sum to 1. Note that every time we apply this process on a column, the resulting ratio of any two of the column’s 27 probabilities could be anywhere between 1/(1 + μ) smaller and 1 + μ larger than its original value. Thus, their relative order of magnitude is preserved.
Sample simulations were performed with time step Δτ = 10 s, κ = 30, 60 (i.e. perturbations applied every 5 or 10 min) and μ = 1, 4, 9. The resulting cell populations did not vary significantly. In fact, they displayed an almost exact agreement on the numbers corresponding to their greatest order of magnitude and slight variations on the number corresponding to their second greatest order of magnitude. For the simulations described in section VIII, the values κ = 30 and μ = 1 were used.
IX.ii In this last part, we provide some suggestions on how the cellular diffusion module (section IV) may be readjusted, in order to consider slightly different sets of starting assumptions. We do this by discussing four specific examples, and show that each of them can be implemented using the sparse matrix, vector algebraic framework we propose. The resulting algorithmic modules were not used for a comprehensive analysis like the one presented in the main text. In this paper, they serve mainly as ideas on how to further elaborate on the specific module.

a)
For the simulations in the main text M_{max} was treated as a constant, and was set 2% larger than M for all voxels. A more general assumption would be to assume that for each voxel A the corresponding number M_{max}(A) is a random number ranging close to M, that also changes with time; The corresponding implementation is straightforward: each time before applying Algorithm 2, generate this vector of random numbers; Then simply apply Algorithm 2, by first substituting \( \overline{M_{max}} \) with this vector.

b)
Additionally to the assumptions (i)(v) in section IV consider the following:
(vi) Live cancer cells do not invade voxels whose total necrotic (cancer+normal) cell population is above M.
This is an additional scenario that can be quite easily implemented in the sparse matrix framework we propose; each time before applying Algorithm 2 find all (new) voxels whose total necrotic population is above M. For each one of these voxels, do the following:

Find the row of the matrix T_{c} corresponding to the voxel.

Find the row’s nonzero elements;

Set all row elements to zero;

Normalize the columns of T_{c} corresponding to the nonzero elements of the row;
Subsequently, apply Algorithm 2 by using the resulting form of the matrix T_{c}

c)
Alternatively, instead of the previous additional assumption, consider the following:
(vi) Live cancer cells invading voxels whose total necrotic (cancer+normal) cell population is above M turn immediately to necrotic.
The implementation of the corresponding algorithm is a slight variation of Algorithm 2. We remind that for a vector x and a number μ, the vector (x ≥ μ) is the binary vector with elements set to 1 if the corresponding element of x is ≥μ and 0 otherwise. For vectors x, y, x. ∗ y denotes their elementwise product. The resulting algorithm is:
Note that in the last two lines, invading cells are added to the live cancer cells of the invaded voxel if its total necrotic population is below M. If its total necrotic population is above M, invading cells are added to the necrotic cancer cells of the voxel.

d)
Last, we informally discuss some ideas on how to introduce some at least crude, phenomenological notion of chemotaxis in this module. Chemotaxis refers to the phenomenon where cells tend to move towards higher concentrations of specific molecules, say, for example, oxygen. We remind that each column of the matrix T_{c} corresponds to a voxel A. Each nonzero element of the column corresponds to a specific neighbor of A. In fact, each nonzero element is the fraction of live cancer cells in A that are in excess of M_{max} that will invade the respective neighbor. How can we quantify the tendency of cells to move towards higher concentrations of oxygen in this framework? The simplest, although admittedly crude way is the following. Each time before executing Algorithm 2, for each voxel that has been reached by the tumor or is adjacent to it, multiply elementwise the respective column with the current vector of oxygen quantities o_{t}. Then normalize its resulting nonzero elements (which remain in the same place) to sum to one. This process rescales the probability of invading each neighbor according to its respective oxygen content; We stress that this is merely a rough idea, which however, we wanted to share. The quantitative details and overall soundness of this idea remain to be furtherly analyzed and improved.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Tumor growth
 Cellular diffusion
 Chemical diffusion
 In silico modeling
 Computational modeling