Estimating Dose Painting Effects in Radiotherapy: A Mathematical Model

Tumor heterogeneity is widely considered to be a determinant factor in tumor progression and in particular in its recurrence after therapy. Unfortunately, current medical techniques are unable to deduce clinically relevant information about tumor heterogeneity by means of non-invasive methods. As a consequence, when radiotherapy is used as a treatment of choice, radiation dosimetries are prescribed under the assumption that the malignancy targeted is of a homogeneous nature. In this work we discuss the effects of different radiation dose distributions on heterogeneous tumors by means of an individual cell-based model. To that end, a case is considered where two tumor cell phenotypes are present, which we assume to strongly differ in their respective cell cycle duration and radiosensitivity properties. We show herein that, as a result of such differences, the spatial distribution of the corresponding phenotypes, whence the resulting tumor heterogeneity can be predicted as growth proceeds. In particular, we show that if we start from a situation where a majority of ordinary cancer cells (CCs) and a minority of cancer stem cells (CSCs) are randomly distributed, and we assume that the length of CSC cycle is significantly longer than that of CCs, then CSCs become concentrated at an inner region as tumor grows. As a consequence we obtain that if CSCs are assumed to be more resistant to radiation than CCs, heterogeneous dosimetries can be selected to enhance tumor control by boosting radiation in the region occupied by the more radioresistant tumor cell phenotype. It is also shown that, when compared with homogeneous dose distributions as those being currently delivered in clinical practice, such heterogeneous radiation dosimetries fare always better than their homogeneous counterparts. Finally, limitations to our assumptions and their resulting clinical implications will be discussed.


Introduction
Radiotherapy, the use of ionizing radiation to eliminate pathological tissues, is a treatment of choice for more than 50% of cancer patients diagnosed with solid tumors [1]. Technical and methodological advances have allowed radiation oncology to achieve local tumor control in a considerable number of patients. However, locoregional recurrence (LRR) remains a problem in many clinical settings. For example, a recent study in patients with Stage III lung cancer found a 5-year LRR rate of 31% [2]. In Glioblastoma Multiforme (GBM), the most common and aggressive malignant primary brain tumor, LRR approaches 90% [3]. In such critical cases, radiotherapy usually results in an initial shrinkage of malignancies, followed by a subsequent growth recovery that cannot be checked even by resorting to larger radiation doses.
The onset of radioresistance, and its resulting poor prognosis, is strongly correlated with the development of significant intratumoral heterogeneity. For that reason, there is growing interest in the clinical significance of tumor heterogeneity. In different works have been recently demonstrated extensive genetic variations in tumor cells due to intratumoral evolution [4], [5]. Moreover, tissue-level heterogeneity due to variations in vascular density and blood flow has been long since evident in clinical medical imaging. In recent years, accumulating evidence suggests that tumor heterogeneity is a key factor in the development of therapeutic resistance and therefore in radiation therapy outcomes [6], [7], [8]. As a consequence, increasing attention is being paid to ''dose painting'' (or ''dose sculpting''), a technique which consists in prescribing different radiation dosimetries to different regions within a given tumor, so that irradiation be boosted in more radioresistant (for instance, hypoxic, quiescent, etc.) regions [9], [10]. This strategy, which is in sharp contrast with the still prevailing homogeneous radiation dose delivery approach recommended by International Commission on Radiation Units (ICRU) reports (50, 1993; 62, 1999 and 83, 2010; see [11], [12] and [13] respectively), has been made possible by the availability of highprecision clinical particle accelerators, and looks particularly promising in those cases where current treatment techniques fail to provide sufficient tumor control.
However, in order for dose painting to show its full power, detailed information is needed about the internal structure of the tumor to be irradiated. Such information should ideally be provided by medical imaging techniques. These however are not yet able to distinguish different radiosensitivity regions except in a few cases, commonly related to hypoxia effects. On the other hand, even a modest miscalculation in the sizes of different radiosensitivity subvolumes has been suggested to produce serious consequences in clinical outcomes [14]. In view of current technical limitations, the question thus arises of providing tools to a) obtain as much information as possible about tumor heterogeneity before a radiation dosimetry plan is prescribed, b) simulate the effects of dose painting therapies which take into account whatever heterogeneity data are available, and c) compare such simulations with those corresponding to standard homogeneous radiation dose distributions currently delivered in clinical practice.
The work herein reported intends to yield some insight into these issues. More precisely, in the sequel a mathematical model for heterogeneous tumor growth is formulated, and the effects of various radiation dose distributions on it are investigated by means of computer simulations. Specifically, we consider a situation where two tumor cell phenotypes, cancer cell (CC) and cancer stem cell (CSC), are present at an early stage, when the tumor consists of about 10 5 cells in total. Concerning CCs and CSCs, we have assumed that i) CSCs represent only a small percentage of the total number of cells at that stage (say, about 15%), ii) CSCs have a significantly longer cell cycle duration than CCs and can replicate indefinitely, while CCs can perform only a limited number of cell divisions, and iii) CCs and CSCs show quite different resistance to radiation, CSCs being more radioresistant than CCs. These biological and radiobiological features have been reported in the literature, specifically for Glioblastoma Multiforme (GBM), where there is mounting evidence of CSCs presence in GBM tumors (cf. for instance [8], [15], [16], [17]). Growth of the heterogeneous tumor thus resulting is simulated by means of an agent-based model in which each cell is individually represented [18], [19]. Tumor growth is kept track of until a size of approximately 10 6 cells is attained, which roughly corresponds to a spheroid of about 1 cubic millimeter in size, a typical volume in multi-cellular spheroids (MCS) in vitro growth. At that stage, different (homogeneous and heterogeneous) radiation dose distributions are simulated using the Linear-Quadratic (LQ) model [20], [21], and their effects compared.
An interesting consequence of i) and ii) above is then shown to be that, as tumor grows, most of the CSCs concentrate themselves within the tumor core, irrespective of their initial distribution at an earlier stage. This fact, which will be described to be inversely correlated with cell migration rates when migration is not inhibited by cell-cell adhesion (which is the case, for instance, after cells undergo an epithelial-mesenchymal transition (EMT) [22]), is then used together with iii) to simulate the effects of different radiation dosimetries to achieve tumor control, or in the case where this cannot be obtained, to compare tumor heterogeneity (seen as an indicator of malignancy in terms of the proportion of CSCs) before and after treatment has been delivered. In this context it will be shown that for a given amount of radiation, heterogeneous dose distributions, where different radiation doses are delivered at different regions of the tumor according to the presence of more radioresistant cells there, invariably fare better than homogeneous ones when sufficient information about tumor spatial heterogeneity is available. In our case, such information will be shown to follow from assumptions i) and ii) above. It should be noticed that hypotheses i), ii) and iii) are amenable to experimental validation, at least in vitro.
Our work can be considered as a preliminary step towards analyzing preclinical models where larger tumors (of the order of cubic centimeters) should be dealt with, several tumor cell phenotypes would simultaneously be present (possibly as a consequence of mutations) as tumor expands, and vascular networks, immune response, and hypoxic and necrotic effects are also taken into account. While the case herein considered is still far from that situation, the simplicity of the setting selected allows us to stress the consequences derived from the minimal number of biological and radiobiological assumptions made on the tumor cell phenotypes involved. This last is particularly relevant in view of the scarcity of in vivo biological parameter measures available. Scaling results up to larger tumor sizes, as well as increasing phenotypic and anatomical complexity appear as feasible within the same approach, but only after key biological data retrievable by non-invasive probing had been identified, and their impact on tumor growth elucidated, an objective toward we intend to contribute with this work.
We conclude this introduction by observing that considerable attention is being currently paid to mathematical modeling as a tool towards designing patient-tailored and adaptive therapies; see for instance [22], [23], [24], [25], [26], [27], [28], [29], [30] and [31]. In particular, radiotherapy modeling and simulations have been addressed in [32], [33], [34], [35], [36], [37], [38], [39] and [40], as well as in [41], [42] and [43] where GBM cases are considered. Mathematical models and computer simulations on the impact of the presence of CSCs in tumor therapies have been discussed in [32], [44] and more recently in [43], where focus is made in a GBM case. It is worth to be stressed, however that in the cases previously mentioned, the total number of cells simulated (and thus the resulting structural complexity) remained way below that of the computer simulations arrived at in our current work.

Tumor Cell Phenotypes Assumptions
For definiteness model parameter values corresponding to Glioblastoma Multiforme (GBM) cell lines have been used [45], [46], [47]. More precisely, we consider a tumor where two different phenotypes coexist at an early stage, when we assume a preponderant (approximately 85% of the total tumor volume) proportion of a tumor cell phenotype denoted as CC (cancer cell) coexisting with a second tumor cell phenotype CSC (cancer stem cell), randomly distributed, that roughly represents 15% of the total population at that stage. Both phenotypes CC and CSC are supposed to possess markedly different biological and radiobiological properties. In particular, we assume: P1.-The duration of cell cycle for CCs is significantly shorter than that of CSCs. In particular, CCs are assumed to divide every 26 hours. Then, for tumor cell phenotype CSC three cases are considered, corresponding respectively to a CSC cycle duration of 96 hours (four days), 72 hours (three days) and 48 hours (two days). Moreover, CCs are assumed to divide a maximum of 15 times, while CSCs are able to replicate indefinitely.
Concerning property P1, it is currently assumed that CSCs proliferate at a slower pace than ordinary cancer cells (see for instance [16], [48], [49], [50], [51], [52], [53], [54], [55] and [56]). Actually, as observed in the references previously quoted, slow-cycling is to be expected from CSCs since such cells belong to tumor phenotypes that are highly resistant to current therapies (radiotherapy, chemotherapy or combined) and these are targeted at killing cycling cells. On the other hand, recent in vivo experiments in a mouse model of Glioblastoma to identify and isolate CSCs through genetically engineered mice demonstrate the presence of a small pool of slow-cycling and highly tumorigenic cells that retain long-term self-renewal ability [55], [57]. We notice that cell cycle durations of 24 h-26 h for GBM have been reported [43], [45], [58], although considerably different cell cycle durations, which in particular include the values herein considered for CSCs, have been noticed as well [45], [59]. Concerning the assumption on the maximum number of CCs replications (see for instance [60]), we have selected the value 15 (cf. [32], [61], [62]), but our results continue to hold if this number is slightly changed. Actually, an arbitrary increase in CSC cycle duration is always compatible with our results, as long as CC cycle duration continues to be significantly faster.
In the course of tumor growth, each of the previous tumor cell phenotypes may transiently enter in a quiescent, non-proliferating stage, due to contact inhibition. Moreover, replication of CCs is always supposed to be symmetric. On the other hand, CSCs will be assumed to sustain either symmetric or asymmetric division, in which case one CSC and one CC will result from replication. Evidence for asymmetric division for CSCs, has been reported in [63], [64], [65]. Since reliable estimates about actual probabilities of asymmetric division p a do not seem to be available as yet, computer simulations will be performed for different choices of that model parameter, namely p a~0 :75, p a~0 :50 and p a~0 :25 (cf. for instance [32], [61]).
A second key assumption is: P2.-When irradiated, CSCs are significantly more resistant to radiation than CCs.
As a matter of fact, CSCs have been described as a comparatively small subpopulation that is highly radioresistant [15], [17], [61], [66]. Radioresistance and surviving cell fractions are estimated by means of the standard Linear-Quadratic (LQ) model [20], [21]. According to it, the surviving fraction of cells after a radiation dose D has been delivered, SF (D), is given by: where D is usually measured in Grays (Gy) (1 Gy is 1 Joule per Kilogram), (a,b) are the so-called radiosensitivity parameters, which depend of the cell phenotype considered, and j is a parameter introduced, as in [32], to distinguish the different radiosensitivities of the proliferating and quiescent states for CCs and CSCs. Actually, cells in a quiescent state (in the G0 cell cycle phase) are known to be more resistant to radiation than their nonquiescent counterparts [67]. It should be noticed that, when estimating the impact of radiation according to the LQ model, what matters is the particular combination of a and b that appears in (1), which provides the surviving cell fractions, rather than the separate values of a and b by themselves. For definiteness, we take in the sequel a~0:48 Gy {1 , b~0:02 Gy {2 and j p1~1 :00 for proliferating CCs. These radiosensitivity parameters have been reported in [47], where in vitro estimates on surviving cell fractions at 2.0 Gy, SF (2), can be found for different GBM cell lines (see also [68]); similar values for a and b have been recently proposed in [43]. In particular, SF (2)~0:36 for proliferating CCs in our case (to be compared to the value SF (2)~0:44 corresponding to a and b parameters proposed in [43]). For quiescent CCs we take j q1~0 :85, so that the corresponding value is SF (2)~0:42. For proliferating CSCs the value j p2~0 :30 (SF (2)~0:73) is taken, whereas for quiescent CSCs the value j q2~0 :20 (SF (2)~0:81) is selected. Also, such surviving cell fraction ranges at 2.0 Gy have been observed and reported in the literature (cf. [43], [68], [69]). We point out that the results to be described in this work continue to hold when the values selected for the radiosensitivity parameters a and b undergo considerable variations, which in particular include the ranges considered in the references quoted above. As a matter of fact, once assumptions P1 and P2 are made, our model is shown to be quite robust with respect to changes in its parameters.
A Three-dimensional (3D) Model of Stochastic Tumor Growth Different mathematical models of tumor growth and its radiation response have been reported in the literature. For instance, tumor growth models and radiation effects with continuous and discrete populations have been reviewed in [70], [71], [72] (see also [33], [34], [35], [36], [73], [74], [75], [76] and [77] for more details). On the other hand, the effects of different radiation dosimetries have been considered in [32], [33], [34], also in [43], [78] for fractionated radiotherapy and in [79] for a case of stereotactic radiosurgery. However, little seems to be known concerning mathematical modeling and computer simulations on the effects of heterogeneous radiation dose distributions on heterogeneous tumors, as in the case herein examined.
The model of tumor growth implemented in this work is as follows. Within the growing tumor, both tumor cell phenotypes, CC and CSC, will be subject to the same kinetic rules. More precisely, following [18], [80], [81], a three-dimensional (3D) cellular automata (CA) model for tumor growth is developed, where each cell is considered as an individual agent. In particular, each cell (whether CC or CSC) occupies a single node in a 3D unstructured lattice (a lattice with no rotational or translational symmetry [81]) thus avoiding symmetry artifacts. Cell division, migration, apoptosis (programmed death) and lysis (removal of debris) have been included and are represented by stochastic processes. Accordingly, each kinetic rule is characterized by a rate, and the governing equation to be solved is a multivariate master equation (see equation (2)). Figures 1 and 2 show a sketch of the processes included and a scheme describing the possible actions that a cell is able to perform in the mathematical model respectively (for further details, see Document S1 provided). Nutrient-limited growth is not accounted for in our model. This issue, as well as others, could be included at the expense of increasing complexity by adding degrees of freedom in the computer simulations, but they are not expected to play a significant role in a tumor cell colony of the size considered, which may be assumed to be fully oxygenated [80]. At any rate, for tumors of the size considered in this work our assumption is not unlikely. For instance, NIH3T3 cells form tumors of size larger than 1 cm 3 without apparent necrotic core even though microlesions may be observed [80].
As to the rules of the model of tumor growth, proliferation is only possible for cells located at a node having at least one free neighbor in the lattice. In the case that all neighbor sites are occupied, a cell enters in a quiescent state due to contact inhibition. However, quiescence is abandoned, and cells return to their normal state, as soon as one of the surrounding nodes becomes free. Proliferation, apoptosis, migration and lysis are modeled as stochastic processes occurring with certain rates. A Poisson process has been assumed for each individual stochastic process and a master equation for the change of the probability of the multi-cellular configuration at time t (denoted by the variable Z) in terms of the multi-cellular configuration in another state Z 0 at time t 0 is then used. It reads as follows: where p(Z,tjZ 00 ,t 00 ) denotes the conditional probability of finding the multi-cellular configuration Z at time t given the configuration was Z 00 at time t 00 , r Z?Z 0 being the transition rate from configuration Z to Z 0 . Notice that the master equation (2) is a balance equation. Indeed, the first term on the right of (2) is a gain term that summarizes all transitions that increase the probability of finding the corresponding multi-cellular system in configuration Z. On its turn, the second term in the right describes transitions that move the system away from Z, and thus represents a loss term. Equation (2) can be numerically solved if the initial condition p(Z,t~0)~d(Z 0 ) is given, where Z 0 denotes the initial multicellular (in our case, tumor cells) configuration. A configuration is determined by the spatial distribution of cells and the state of each cell (proliferating, quiescent, etc.). In our model both CCs and CSCs are able to migrate with the same rate. Migration is mimicked by a hopping process allowing any cell to move from one lattice site to a free neighbor lattice site. In case several free neighbor lattice sites exist, one of them is randomly chosen. In this work we have considered two different migration rates (k mig ), a comparatively low rate (0.025 h 21 ) in the range obtained from the cell diffusion constant [18] as outlined in the Document S1 provided, and a higher rate (1.75 h21 8) as estimated in vitro in [43] for a GBM cell line. These will be respectively referred to as low and high migration rates in the sequel. It should be stressed that we only consider here the case where the motion of a cell from one lattice site to another does not depend on the contact energy between neighboring cells, but only on the availability of space. In that case, the higher the migration rate, the stronger the cell dispersion is (see for instance Movies provided as supporting information). If however cell-cell adhesion would be considered, migrating cells would tend to fill holes and cavities [81], and migration will lead instead to tumor compactification. We assume that in our current setting this case is substantially included in very low migration cases. On the other hand, CCs and CSCs undergo programmed cell death (apoptosis) (see for instance [32]). Disposal of cellular debris resulting from apoptosis is carried out by a lysis process [82], for which a lysis rate k lys~0 :035 h {1 (about 30 h) has been assumed. This is about 10-fold less than phagocytosis (digestion of cellular debris by macrophages) observed in vivo in [83], but within the range reported for in vitro cultures 0.002 h 21 for Hybridoma VO 208 cell line [84] to 0.07 h 21 for Fibrobacter succinogenes [85].
The master equation (2) has been numerically solved by means of the so-called Gillespie algorithm [86], (also called kinetic Monte-Carlo algorithm or Bortz-Kalos-Lebowitz algorithm [87], see Document S1 for more details). Notice that one advantage of using a lattice model is the possibility of extending the same formalism at larger scales by permitting more than one cell to occupy a single lattice site [80]. In order to simulate the resulting biological effect when a radiation dose D is delivered, the surviving fraction is computed for each tumor cell phenotype according to the LQ model (1), taking into account the state of each cell, proliferating or quiescent. Surviving cells are randomly selected out of the total cell number involved. In the next section we will show a typical starting point (about 10 5 cells) in the computer simulations to be described below, as well as the resulting tumor once a size of about 10 6 cells has been reached under the kinetic rules just described, and before the radiotherapy treatment is applied. For convenience of the reader, we provide in Table 1 all parameter values used in our mathematical model to simulate the tumor growth and radiation response.

Results
As a reference case, we first consider the effect on a fully monoclonal tumor (containing only CCs) of a therapeutic irradiation protocol, consisting of 30 sessions of 2.0 Gy, each being homogeneously delivered on the tumor. According to standard radiotherapy scheduling, sessions are distributed into 6 weeks, each week including five sessions from Monday to Friday separated by 24 hours intervals, and with a 72 hours interval from Friday to next Monday in the following week (with weekend interruptions). The total radiation dose delivered with this treatment is thus 60 Gy. This is currently considered a standard radiotherapy treatment for most GBM tumors [88], [89], [90]. The corresponding process is illustrated in Figure 3, both for the cases of low and high migration, under our current assumption that migration is not inhibited by cell-cell adhesion. Figure 3 shows that tumor prior to treatment grows from week 1 to approximately week 5 for the low migration case, and just in one week in the case of high migration, until a size of about 10 6 cells is attained. Notice that the growth time of the tumor decreases with the migration rate due to the decreasing effect of contact inhibition inside the tumor. Then radiation therapy treatment starts, and accordingly tumor cell number diminishes during the first week (with small regrowth between each daily session and weekend interruptions, as represented by the knots in the straight line in the plot, see Figure 3). The pattern just described is reproduced until tumor eradication is achieved at the end of the radiotherapy treatment in these cases. A heterogeneous tumor containing the two tumor cell phenotypes is now considered. Starting from an initial configuration where 10 5 cells are present, out of which approximately 85% are CCs and 15% are CSCs, tumor growth is allowed until a size of about 10 6 cells is reached (see Figures 4 and 5). Then the impact of homogeneous and heterogeneous radiation dose distributions is modeled, and computer simulation results are compared in the cases where asymmetric division probabilities p a for CSCs are p a~0 :75, p a~0 :50 and p a~0 :25, the CSC cycle duration is taken to be 96 h, 72 h and 48 h, and the low and high migration rates are considered. The results obtained will show that the standard irradiation protocol described before fails now to achieve tumor control in any of the cases considered. To compare the dynamics of the tumor resulting after irradiation with respect to its pretreatment stage, computer simulations are stopped once the pretreatment population size of about 10 6 cells is again obtained.
As shown in Figure 4, for p a~0 :25, CSC cycle duration equal to 96 h and the low migration rate, the more radioresistant tumor cell phenotype CSC is confined within an inner, smaller region when irradiation is started. Such spatial CSCs distribution is neither a priori imposed, nor a consequence of the specific CSC cycle duration or the asymmetric division probability considered (see Table 2). It is due instead to the difference of the CSC and CC cycle durations. Indeed, a robust emerging feature is now observed. Namely, due to asymmetric division CSCs produce a certain fraction of CCs. Both CCs and CSCs then compete for resources including free space at the tumor border [18], [91]. For sufficiently small micro-motility, that competition is controlled by cell replication. As CCs proliferate faster than CSCs, they have a selective advantage in the competition for free space and will eventually outcompete the CSCs in the border region of the tumor, if (as it happens in our case) to achieve such dominance less replications are needed than the maximum number that CCs can perform. The precise timing depends on the relation of the cell cycle duration for CSCs vs. CCs, p a , and the fraction of CSCs in the initial population at 10 5 cells (notice that this fraction would itself be determined by p a and CSC cycle duration if the 10 5 cells would already have emerged by replication from a single initial CSC).
Therefore, for a low migration rate, CSCs will be contactinhibited by the fast proliferating CCs. As a consequence, CSCs will remain confined in an inner region in that case. Actually, on assuming a cell diameter of 20 mm, the diameter of the tumor in all cases is then of about 2680 mm (with a standard deviation of 56 mm over 20 simulations performed for each parameter set considered)   and the volume of this inner region where 100% of CSCs are located varies from the 15% to 25% of the total tumor volume, when asymmetric division probability and the CSC cycle duration are allowed to change (see Table 2, and Figures in the Document S1 provided). On the other hand, when the high migration rate is considered, CSCs are not fully concentrated in an inner region of the tumor (see Figure 5 for p a~0 :25 and CSC cycle duration equal to 48 h). However, we can define an inner region where at least 80% of CSCs are located (see Figure 6). When asymmetric division probability and CSC cycle duration are allowed to change in the parameter range considered, this inner region approximately represents between 21% to 40% of the volume where 90% of cells, both CCs and CSCs, are located (see Table 3, and Figures in the Document S1 provided). In this case, the diameter of the tumor for all cases is about 5294 mm (with a standard deviation of 778 mm over 20 simulations performed for each parameter set considered). In Tables 2 and 3 the number of CSCs just before treatment starts is shown, so that its dependence with migration rate, asymmetric division probability and CSC cycle duration can be observed. Actually, the number of CSCs existing before treatment starts is a key factor to estimate tumor resistance to radiation therapy, as we will recall below.
Bearing these facts in mind, it turns out that tumor control can be obtained in all cases when a radiation boost is applied at such internal regions. More precisely, in the case of low migration we observe that tumor control can be achieved for CSC cycle durations equal to 96 h and 72 h, when 2.5 Gy (for the case p a~0 :75), 2.9 Gy (for p a~0 :50) and 3.3 Gy (for p a~0 :25) are respectively delivered within the largest inner sphere containing 100% of CSCs, and 2.0 Gy is delivered in the rest of the tumor, according to the former standard fractionation protocol (5 days a week along 30 sessions at 24 hours intervals except for weekends). However, when CSC cycle duration is 48 h, tumor eradication is not possible with these heterogeneous therapies under the same conditions. In that case, to obtain tumor control, the dose delivered in the inner region has to be raised to 2.7 Gy, 3.4 Gy and 3.9 Gy respectively (see Table 4 and the Document S1 provided). Notice that these radiation doses increase as asymmetric division probability decreases.
Let us now consider the same heterogeneous therapies for the case of high migration. We now select an inner region where 80% of CSCs are located. Considering these heterogeneous therapies for each case as before, tumor control is now obtained only for p a~0 :75 with CSC cycle durations of 96 h and 72 h. This is due to the fact that i) the high migration rate permits less contact inhibition, which in turn allows for rapid re-growth, and ii) there are about 20% of CSCs which only receive a radiation dose of 2.0 Gy. Therefore, to obtain tumor control it is not only necessary to increase the radiation dose in the inner region, but also in the outer one (see Table 5 and the Document S1 provided). The radiation doses of the heterogeneous therapies required to obtain tumor control are provided for each case of migration rate considered in Tables 4 and 5. The respective temporal evolution of the number of each tumor cell phenotype is shown in Figure 7  It may appear at first glance that the successful results obtained for heterogeneous dosimetries could be a consequence of the overall radiation dose delivered over the tumor being larger than that administered according the standard irradiation protocol (2.0 Gy a day, 5 days a week at 24 hours intervals, with weekend interruptions and 60 Gy in total). However, heterogeneous dosimetry turns out to be crucial to achieve tumor control. In particular, tumor control fails to be attained when we deliver an averaged homogeneous dose (AD), corresponding to the same global radiation energy as in the heterogeneous dosimetry, carried out along a similar scheduling. The corresponding AD is given by: where V in , V out are the volume of the internal sphere and the remaining shell considered; D in , D out are the radiation doses delivered over the internal and external regions just described, and V tot is the total volume of the tumor. In the case of low migration the inner region is that where 100% of CSCs are located (see Table 2 for values of the diameter of this inner region for each case) and the volume of the outer region is computed with respect to the average diameter of the tumor at the beginning of the radiotherapy treatment (2680 mm). However, for the case of high migration the inner region is now selected as that where 80% of CSCs are located (see Table 3 for further details). Indeed, since some cells are now isolated far from the tumor bulk, instead of defining the tumor radius as the distance from its center of mass to the farthest cell, to compute the averaged dose the volume of the tumor is now considered as that of the region where 90% of total cells (CCs and CSCs) are located, where the diameter is about 3120 mm (with a standard deviation of 186 m over 20 simulations performed for each parameter set considered). The reason for this assumption is that it will yield a higher AD than that obtained when considering the maximum diameter of the tumor (5294 mm), which will extend to regions sparsely occupied by tumor cells.
Hence the averaged homogeneous therapies thus derived will deliver higher radiation doses than those that would be obtained if the outer shell were defined as that where all tumor cells are contained.
The averaged dose (AD) per session according to equation (3) is shown in Tables 4 and 5, for each of heterogeneous therapies described before. These AD vary from 2.10 Gy to 2.32 Gy for the low migration case and from 2.10 Gy to 2.52 Gy in the case of high migration for the asymmetric division probabilities and CSC cycle durations considered. Notice that the total radiation doses delivered by these averaged homogeneous therapies are higher than 60 Gy (the value corresponding to the standard irradiation protocol) for all cases. The total radiation doses corresponding to these new dosimetries range between 63.0 Gy to 69.6 Gy for the case of low migration and 63.0 Gy to 75.6 Gy for the high migration case. Some of these results are illustrated in Figure 7 (Figure 8 (B)) more cells remain isolated at the end of the treatment compared with the case of low migration (Figure 8 (A)). In that Figure, when tumor control is not achieved, computer simulations are performed until the surviving tumor reaches a size approximately equal to 10 6 cells, the number of cells it had before radiotherapy started.
Tables 4 and 5 reveal that tumor recurrence occurs in all cases for a homogeneous therapy delivering the corresponding average dose (AD). Besides, the number of CSCs in the tumor at the end of the radiotherapy treatment decreases with p a and CSC cycle duration (see Figure 9, and Tables in the Document S1 provided). In the case of low migration, for the heterogeneous therapies failing to achieve tumor control, the number of CSCs remaining alive at the end of the recurrence tumor stage is 107, 1785 and 4457 respectively, with the corresponding standard deviations being 8.53, 78.31 and 232.67 (see Figure 9 (A) to compare with the corresponding averaged homogeneous therapies). These values correspond to the cases p a~0 :75, p a~0 :50 and p a~0 :25 with a CSC cycle duration of 48 h. In Figure 9 (B), the number of CSCs at the end of the recurrence tumor stage is provided in the case of high migration for the heterogeneous therapies delivering 2.5 Gy (for the case p a~0 :75), 2.9 Gy (for p a~0 :50) and 3.3 Gy (for p a~0 :25) in the inner region, and 2.0 Gy in the rest of the tumor. Notice that, even when tumor control cannot be achieved with the heterogeneous therapies, the corresponding averaged homogeneous therapies always have more CSCs at the end of the recurrence tumor stage (see Figure 9 (C)). Moreover, in some cases that number of CSCs is larger than before the treatment started, resulting in more radioresistant tumors after treatment (see Document S1 for further details).
Thus, to achieve full eradication of a tumor consisting of two different tumor cell phenotypes, heterogeneous dosimetry is crucial. Actually, the choice of a minimal radiation dose sufficient to achieve tumor control depends on the value of p a , the CC and CSC cycle durations and on the internal spatial distribution of CSCs. In Tables 4 and 5, we describe the heterogeneous radiation therapies needed to achieve tumor control, and the corresponding averaged homogeneous therapies are also provided. Interestingly, the corresponding averaged homogeneous therapies in each case fail to obtain tumor control (see also Tables in the Document S1 provided). Moreover, homogeneous therapies needed to obtain tumor control are also provided in Tables 4 and 5. One readily sees that in all cases higher total radiation doses are needed for homogeneous than for heterogeneous therapies.
On the other hand, considering that for all choices of model parameters the AD is higher than 2.0 Gy, this implies that tumor recurrence will also occur for the standard irradiation protocol (2.0 Gy a day, 5 days a week at 24 hours intervals, with weekend interruptions and 6. Gy in total) for each case of p a , migration rate and CSC cycle duration considered. In terms of the number of remaining CSCs after treatment is completed, recurrence is certainly weaker when AD is delivered than for the standard fractionation protocol, as one could expect from the comparative increase in radiation delivery. Moreover, tumor control cannot be achieved for each case of p a , migration rate and CSC cycle duration considered even when the homogeneous therapy delivering the average radiation dose is rescheduled in 7 days a week along 30 sessions at 24 hours intervals, without weekend interruptions (see Figure 10 for some examples of averaged homogeneous therapies with this fractionation protocol). Since in many clinical scenarios radiation doses are mostly limited by damage inflicted at neighboring organs at risk and healthy tissues (see [92], [93] and [94]), it is important to estimate what amount of tumor control can be achieved when radiation dose distributions are kept as low as possible. In what follows, we shall consider a heterogeneous therapy for which the average radiation dose is approximately equal to 60 Gy, and a case of hyperfractionation (a type of scheduling consisting of comparatively many sessions, usually more than 1 per day, with low radiation doses [95]); cf. [96] for a specific study on GBM tumors. We shall see that in these cases a heterogeneous radiation dose distribution also yields better results than its averaged homogeneous equivalent, even when tumor control is not achieved.
Consider first the case of low migration, CSC cycle duration equal to 96 h and where the value of the total radiation dose is a bit less than 60 Gy for heterogeneous therapies consisting of 2.3 Gy for p a~0 :75, p a~0 :50 and p a~0 :25 within the largest inner sphere containing 100% of CSCs and 1.8 Gy in the rest of the tumor delivered 5 days a week along 30 sessions at 24 hours intervals with weekend interruptions. Computer simulations show that these radiation dosimetries fare better than their averaged homogeneous versions, where total AD lies between 56 Gy and 58 Gy for the values of p a considered (see Figure 11 (A), (B) and (C)). On the other hand, similar results can be obtained for a lower total radiation dose when the time lapse between sessions is also shortened. More precisely, consider the same cases but now for heterogeneous therapies consisting of 1.7 Gy within the largest inner sphere containing 100% of CSCs and 1.2 Gy in the rest of the tumor, delivered in two sessions per day, 5 days a week along 30 sessions 12 hours intervals, with weekend interruptions. While tumor control is not achieved, tumor radioresistance, measured in terms of the final proportion of CSCs, turns out to be lower for heterogeneous dosimetries than their averaged versions (see Figure 11 (D), (E) and (F)). Notice that in this case, the total averaged doses delivered by the heterogeneous dosimetries considered are much smaller than 60 Gy (between 39 Gy and 39 Gy for the values of p a considered). Figure 11 shows that there is tumor recurrence in all cases. However, it turns out that the number of CSCs remaining at the end of the recurrence tumor stage after radiation therapy is lower than that existing prior to therapy in all cases. Thus tumors surviving this therapy can be considered as less radioresistant than they were before radiation therapy started. An inspection of Figure 11 (C) and (F) quickly shows that in our current cases heterogeneous therapies yield better results than its averaged homogeneous counterparts previously discussed. For completeness, we provide estimates on the total number of CSCs after treatments are concluded and recurrence appears for each case of p a considered (see Figure 11).
We conclude this Section by remarking on the dependence of our model of tumor growth, and the results derived from its analysis, on data and parameters assumed.  To begin with, our results are not restricted to the figure selected (15%) for the proportion of CSCs within the tumor at the initial stage. In fact, they continue to hold as long as the more radioresistant tumor cell phenotype CSC represents a small percentage of the total tumor cell count. A particularly interesting limit case is that when tumor growth starts from a single CSC. Then for each value of p a (0.75, 0.50, 0.25), CSC cycle duration (96 h, 72 h and 48 h) and migration rate k mig (0.025 h 21 , 1.75 h 21 ) considered, the number of CSCs present when tumor has reached a size of about 10 6 cells (just before radiation treatment starts) is much smaller than that corresponding to the cases considered in this work. For instance, in the case of CSC cycle duration equal to 48 h, p a~0 :25, for the low and high migration cases the number of CSCs before the treatment starts is about 5956 and 14316 (with standard deviations of 129 and 530 over 20 simulations performed for each parameter set considered) respectively. The corresponding values for the case considered in this work are 21953 and 41629 respectively (see Tables 2 and 3). Moreover, the internal region where CSCs remain confined is smaller (1120 mm and 1840 mm with standard deviations of 39.7 mm and 71.4 mm over 20 simulations performed for each parameter set considered, respectively) than that reported in Tables 2 and 3. Hence any radiotherapy treatment that achieves tumor control in our case also does so for tumors staring from a single CSC under assumptions above.
On the other hand, we have made use of the assumption that the duration of cell cycle for CSCs is significantly longer than that of CCs, a hypothesis commonly assumed in the literature (cf. [16], [48], [49], [50], [51], [52], [53], [54], [55] and [56]). This fact notwithstanding, our model can be used to examine also the opposite situation, that is the case where CSC cycle duration is equal or smaller than that of ordinary CCs. As an example, we have considered the cases where CSC cycle lasts 26 hours (respectively 18 hours), which is equal to (respectively less than) the 26 hours cell cycle selected for CCs. As one can expect, the inner core where most CSCs remain concentrated is now larger than when slow-cycling CSCs is assumed. In particular, in the case of low migration and for a CSC cycle duration of 26 h, such internal regions (where 100% of CSCs are located) range from 20% to 83% of the total tumor volume for values of p a equal to 0.75, 0.50 and 0.25. Besides, when CSC cycle duration is taken to be 18 h that internal volume further expands, ranging now between 23% and 100% of the total tumor volume. Additional details, including the number of CSCs present when tumor size reaches about 10 6 cells and the case of high migration, are provided (see Document S1). A case which has not been addressed in this work is cancer cell plasticity, a hypothesis that has been advanced to better understand the onset of resistance after therapy; see for instance [97], [98] and [99]. According to this scenario, in addition to CSCs giving raise to CCs by asymmetric division, a (supposedly small) percentage of CCs may transform to a CSC phenotype, possibly as a reaction to radiation therapy. Although little quantitative information about this process seems to be available as yet, including such type of process in our model is possible. To support this statement, we have studied a particular example. Specifically, we have examined a model situation where a small percentage of CCs are transformed to CSCs along the radiation treatment. As expected, any increase in the number of CSCs results in increased malignancy, measured in terms of higher resistance to radiation therapy. However, our conclusion that heterogeneous, tumor-adapted radiation therapies fare better than their corresponding averaged homogeneous versions continues to hold. Details on this study can be found in the Document S1 provided.

Discussion
Tumor heterogeneity is being increasingly recognized as a key obstacle to achieve successful tumor control, either by means of radiotherapy, chemotherapy or through the use of combined therapies. Indeed, it is well known that tumors at an advanced stage contain different tumor subpopulations, which might have been generated as a consequence of sequential mutations of one initial clonogenic line, or could result from the presence of Cancer Stem Cells. Moreover, it is expected that such cell phenotypes may considerably differ in their biological and radiobiological properties, and in particular in their resistance to radiation (cf. for instance [15], [6], [60], [8], [100], [101]).
Accordingly, it has been proposed that the clinical prognosis of a given tumor would critically depend on the information that may be gathered about its internal heterogeneity, and more precisely, In all cases, treatment sessions were scheduled along 6 weeks separated by 24 hours intervals except for weekends, where a 72 hours interval is allowed. Data corresponding to 20 simulations (with different seeds of a random number generator) are presented. In the heterogeneous therapies, radiation doses are specified both for the outer (left) and inner (right) tumor regions, each case being indexed from (1) to (16). The averaged dose for any of the previous cases is labeled with the same number in the columns corresponding to homogeneous therapies. Within brackets the total dose of the radiation therapy treatment is also provided. See Tables   about the identification of regions within it with different sensitivities to a given therapy. In principle, once the spatial distribution of the various cell phenotypes coexisting in a tumor is known, and the resistance to therapy of each of these regions had been estimated, personalized strategies complementary to (or as a substitute to) surgery, could be designed to improve chances of clinical success. The latter can be understood either as total tumor eradication (the standard paradigm as of today) or as achieving instead a stable, chronically-controlled tumor burden where less aggressive lines keep at bay more resistant ones [102]. In either case, significant information towards a personalized treatment would be inferred from knowledge of the internal, non-homogeneous structure of a tumor and the resulting differences in therapy resistance corresponding to the regions thus identified. Unfortunately, to this day only partial information can be derived about tumor heterogeneity by means of non-invasive medical imaging modalities. Currently available information can be mainly used to distinguish various level of oxygen pressure within the tumor related to hypoxia processes [9], necrotic areas or highly proliferating regions detected by means of PET techniques [103]. While undoubtedly important, such information is not enough to design personalized therapies whose results could significantly improve those obtained by current procedures.
In this work, we have proposed a mathematical model of tumor growth to gain insight about two key issues: how heterogeneity unfolds in a growing tumor, and what type of radiation dosimetry is best suited to achieve control in heterogeneous tumors. Concerning the first issue, we noticed that substantial information about the evolution of spatial heterogeneity within a tumor can be retrieved from knowledge of a few key biological properties of the tumor cell phenotypes involved. In particular, we have shown in the first place that a difference in cell cycle duration between a majority of ordinary cancer cells (CCs) and a minority of comparatively slow-cycling cancer stem cells (CSCs) leads to a concentration of CSCs in regions that can be a priori estimated. In the cases just discussed, such regions consist in an internal core within an expanding tumor, but the result would apply to other geometries as well. In particular, it can be extended to larger tumors with corrugated shapes and boundaries.
We have already mentioned that our key assumption that CSCs have longer replication times than CCs is commonplace in the literature (see for instance [16], [48], [49], [50], [51], [52], [53], [54], [55] and [56]). As a matter of fact, such assumption is naturally associated to the consideration of CSCs as a subpopulation of tumor cells which is able to rescue tumor growth after therapies have been delivered. This is related to the fact that standard radiation therapies preferentially target dividing cells (which are more radiosensitive), and thus spare those that have slower cycles or remain quiescent. Notice that the cell cycle duration could in principle be estimated, at least in vitro, for all cell phenotypes known to appear in a given tumor. Importantly, the spatial heterogeneity pattern thus observed does not depend so much on the precise values of such biological parameters, but rather on the fact that they are significantly different for the tumor cell phenotypes involved. As a consequence, the result obtained is robust with respect to fluctuations in cell cycle duration due to systemic factors.
A second result obtained is that, once information about functional heterogeneity had been obtained, tumor-tailored radiation dosimetries can be designed to improve the treatment outcome. We have shown that heterogeneous radiation dosimetries do better than homogeneous ones when regions occupied by different radioresistant tumor subpopulations can be identified,  Table 4). (B) For heterogeneous therapies consisting of 2.5 Gy, 2.9 Gy and 3.3 Gy in the inner sphere and 2.0 Gy in the rest of the tumor (Top) and the corresponding averaged homogeneous therapies (Bottom) for the cases p a~0 :75, p a~0 :50 and p a~0 :25 (left, middle, right) with the high migration rate and CSC cycle durations equal to 48 h, 72 h and 93 h (see Table 5). In all cases (A, B, C), a standard scheduling (30 sessions along 6 weeks separated by 24 hours intervals except for weekends) was applied. Notice that the vertical coordinate is represented in a logarithmic scale. See Tables in the Document S1 for further details and  Movies S3 and this is in particular the case when more radioresistant phenotypes are assumed to replicate at a lower pace than less resistant ones. Interestingly, this result holds when the more radioresistant phenotype is allowed to sustain unlimited replication as in the case of CSCs, as opposed to the limited number of replications commonly assumed on CCs. The previous statement holds true, no matter the type of scheduling considered (with or without weekend interruptions) or the precise result pursued, being it total tumor eradication, controlled recurrence or palliative treatment. We believe that the comparative advantage of heterogeneous radiation dose distributions deserves some consideration, since to this day homogeneous dosimetries continue to be those being commonly implemented worldwide.
It is worth to stress that our model is quite robust with respect to changes in data and parameter values. In particular, our conclusions remain in force when CSC and CC cycle durations undergo considerable changes, as long as CSC cycle is significantly slower than that of CCs. Also, CCs and CSCs migration rates are allowed to undergo substantial changes (corresponding for instance to slow and fast migration processes) as far as both tumor cell phenotypes share a similar migration rate. Moreover, our results continue to hold when changes in the choice of the radiosensitivity parameters a and b in (1) are allowed, or when different fractions of the minority phenotype (CSC) are assumed. For instance, our results are not confined to the choice made for the assumed percentage (15%) of CSCs present at an early stage of tumor growth. They continue to hold if a different figure for that proportion is taken, as long as CSCs remain a small fraction of the total tumor population at that stage.
On the other hand, cancer cell plasticity has recently received considerable attention ( [97], [98] and [99]). We have studied a particular example where in addition to CCs being generated by CSCs with asymmetric division, a small percentage of CCs are transformed to CSCs as a consequence of radioresistance to therapy. Our conclusion that heterogeneous, tumor-adapted radiation therapies fare better than their corresponding averaged homogeneous versions continues to hold in this case.
We conclude by discussing on some of the limitations of this work, as well as on possible extensions thereof. To begin with, we are aware that more research is needed to understand the possible mechanisms that can be responsible for slow cycling of CSCs. Particularly relevant in this context would be to ascertain if slow cycling can, at least in some cases, be established as an intrinsic property of CSCs or if it could alternatively be induced by systemic feedback in the course of tumor growth. Interestingly, even if CSCs are assumed to cycle faster than CCs, our model still shows that heterogeneous dosimetries adapted to the resulting tumor heterogeneity continue to outperform standard homogenous therapies currently in use. A general conclusion that follows from our study is that detailed information about intratumoral heterogeneity is needed in order to implement efficient dose-painting techniques in clinical practice. In particular, in this work a clear dependence on tumor heterogeneity of the radiation doses needed to achieve tumor control has been obtained. In fact, the inner tumor regions where more radioresistant tumor cell phenotype remains confined are shown to strongly depend on CSC cycle duration and their probability of asymmetric division. In the particularly unfavorable assumption of fast CSCs cycling, this region may rank from 20% to 100% of the total tumor volume. In this latter situation, a worstcase scenario corresponding to a high and homogeneous radiation dose being prescribed and only limited by neighboring organs at risk tolerance, is recovered that corresponds to current clinical practice. Our results suggest that such situation could be considerably improved in many cases if and when sufficient information about key different biological and radiobiological properties of the tumor cell phenotypes present in a given tumor is available, be it either by estimating patient-specific parameters or by means of medical imaging techniques.
Finally, it looks feasible from a mathematical viewpoint to address within this framework situations where larger tumors are considered, a number of cell phenotypes coexist there due to mutations, and other effects (immune response, nutrient limitation, etc.) are accounted for. For example, the modeling framework selected in this work permits simulations to be scaled up to cubic centimeter sizes, though at the expense of lower spatial and functional resolution, and more computing resources. It can also be used to construct hybrid models, zooming in at the cell scale in regions of interest. In particular, we have chosen to represent each cell individually to exclude averaging effects when studying the relation between tumor heterogeneity and simulated radiation outcomes. We hope that this work could provide a starting point towards the study of the more general situations described above.

Supporting Information
Document S1 Details of computer simulations of the model of tumor growth and additional results. Movie S10 Time evolution of a heterogeneous tumor growth (where CCs and CSCs are present) for the low migration case with p a~0 :25 and CSC cycle duration equal to 96 h. Tumor growth is allowed unchecked from a size of about 10 5 cells until about 10 6 cells are present. Then, a homogeneous radiation therapy consisting of 1.30 Gy in the tumor is delivered. Treatment sessions are scheduled along 3 weeks, 5 days a week separated by 12 hours intervals except for weekends, where a 72 hours interval is allowed. Depicted in dark and light green (respectively, dark and light red) are proliferating and quiescent CCs (respectively, proliferating and quiescent CSCs). Dead cells are represented in black. See Figure 11 (E) in the article. (MP4)