Optimization of Dose Fractionation for Radiotherapy of a Solid Tumor with Account of Oxygen Effect and Proliferative Heterogeneity

A spatially-distributed continuous mathematical model of solid tumor growth and treatment by fractionated radiotherapy is presented. The model explicitly accounts for three time and space-dependent factors that influence the efficiency of radiotherapy fractionation schemes—tumor cell repopulation, reoxygenation and redistribution of proliferative states. A special algorithm is developed, aimed at finding the fractionation schemes that provide increased tumor cure probability under the constraints of maximum normal tissue damage and maximum fractional dose. The optimization procedure is performed for varied radiosensitivity of tumor cells under the values of model parameters, corresponding to different degrees of tumor malignancy. The resulting optimized schemes consist of two stages. The first stages are aimed to increase the radiosensitivity of the tumor cells, remaining after their end, sparing the caused normal tissue damage. This allows to increase the doses during the second stages and thus take advantage of the obtained increased radiosensitivity. Such method leads to significant expansions in the curative ranges of the values of tumor radiosensitivity parameters. Overall, the results of this study represent the theoretical proof of concept that non-uniform radiotherapy fractionation schemes may be considerably more effective that uniform ones, due to the time and space-dependent effects.


Introduction
Approximately half of the patients, diagnosed with cancer, undergo radiotherapy (RT) [1]. The effect of irradiation on cancer and normal cells can be mathematically expressed via classical linear-quadratic model, which is known to fit experimental data well in a wide range of clinical parameters [2]. According to this model, the fraction of cells, which survive after a single radiation dose D, can be estimated as where α and β are radiosensitivity parameters of cells. Usually cancer cells have higher values of linear radiosensitivity parameter α that corresponding normal cells. However, normal tissues as a rule have greater α/β ratio, which restricts the use of high radiation doses [3]. One option to reduce normal tissue damage is to concentrate the radiation dose within the tumor mass. However, such option carries significant risks even for tumors with clear boundaries, due to the leakage radiation [4]. Other option to spare normal tissues is to fractionate the total dose, that is, to divide it into much smaller fractions, administered over a period of several weeks. The efficiency of a fractionation scheme depends on the effects that are widely referred to as the four "R"s of radiotherapy [5]. Two of the effects-reoxygenation and redistribution of cell cycle-indicate that the radiosensitivity of a cell depends on the surrounding concentration of oxygen and on the current stage of its cell cycle.
In particular, hypoxic and non-proliferating cells are more radioresistant [6]. The third effect is the repopulation of tumor cells that takes place between the irradiations. The fourth effect is the repair of sublethal damage, which can be neglected, unless the time interval between irradiations is as short as several hours [7]. The most typical RT fractionation schemes consist of fractions of 1.8 to 2.0 Gy, delivered once a day on weekdays within the period of several weeks [3]. However, different fractionation protocols were shown to lead to improvement in tumor cure and patient survival for most of the patients for some of the tumor types [8]. Importantly, in clinical practice optimization of RT fractionation is significantly complicated by a number of factors, including the great variability of tumors, belonging to the same type, and related ethical problems, associated with the fact that alternative protocols may worsen outcome for some of the patients. Of note, in the vast majority of the tested schemes the irradiation doses are distributed equally between the fractions. At that, the varied parameters of the schemes are the number of the fractions, the interval between them and the fractional dose, which are related through the constraint on total normal tissue damage.
Given the practical difficulties, mathematical modeling can be a powerful tool, providing insights into the problem of optimization of RT fractionation. Different approaches exist in this field that have their pros and cons. The use of non-spatially distributed phenomenological models of tumor growth, expressed in ordinary differential equations, often allows to obtain the globally optimal solutions for the considered problems via analytical methods [9][10][11]. However, these methods become complicated and even unsolvable under introduction of complex non-linear terms, aimed to account for timeand space-dependent effects [12] or for tumor-specific features [13]. In such cases, various heuristic approaches are used, ranging from direct comparison of the schemes [14] to more complex techniques like simulated annealing [13], which can not guarantee the global optimality of the solution. However, such methods can yield significant results, one of which has already been verified in a preclinical study [13].
Crucially, the use of ordinary differential equations can allow for only phenomenological consideration of the reoxygenation and redistribution of the cell cycle-the effects that result are the spatiotemporal variability of tumor cells' radiosensitivity [12,14]. Their explicit consideration is possible in spatially-distributed models, which can be divided in two types-a continuous model, expressed in partial differential equations (PDEs), and discrete models that usually treat every tumor cell as a separate agent but use PDEs for the consideration of dynamics of nutrients and other substances. However, to the best of our knowledge, the existing works on RT fractionation optimization that use continuous spatially-distributed models, only account for homogeneous and constant radiosensivity of tumor cells. That leads to the conclusion of optimality of standard radiation fractionation schemes [15,16]. Reoxygenation and redistribution of cell cycle can be straightforwardly incorporated into agent-based models. However, the complexity of such models and numerical costs of their simulations lead to practical impossibility of solving optimization tasks via them, at least under current level of computer technology. As a rule, in the corresponding works several fractionation schemes are compared directly, moreover, the considered numbers of tumor cells are several orders of magnitude less that the relevant numbers for the human tumors [17][18][19]. These factors significantly limit the usefulness of such models.
In this work, we present a spatially-distributed continuous mathematical model of solid tumor growth and treatment by fractionated RT that explicitly accounts for tumor cell repopulation, reoxygenation and redistribution of proliferative states. With the use of a specially-developed algorithm, we find the optimized fractionation schemes for varied radiosensitivity of tumor cells under the values of model parameters that correspond to different degrees of tumor malignancy. The resulting schemes lead to significant expansions in the curative ranges of the values of tumor radiosensitivity parameters.

Equations for Tumor Growth
The mathematical model of tumor growth, considered herein, was based on our models, previously used for the investigation of various aspects of tumor growth and treatment [20][21][22][23]. There were five variables in this version of the model, which were the functions of space and time coordinates, r and t: the density of tumor cells n(r, t), the density of normal cells h(r, t), the fraction of necrotic tissue m(r, t), the concentration of glucose g(r, t) and the concentration of oxygen ω(r, t).
The main simplification of this version of the model was the absence of an explicit variable for the capillaries, from which the nutrients flow into the tissue. Instead it was assumed that their local density was proportional to the local fraction of the normal cells in the tissue.
The following set of equations governed the dynamics of the model variables under the absence of radiotherapy:

Dynamics of Cells and Necrotic Tissue
The term of tumor cells proliferation implied that the rate of this process was proportional to the rate of glucose consumption by tumor cells. This assumption was made on the basis that glucose is indispensable nutrient for biosynthesis [24]. Glucose is also a crucial energetic nutrient for tumor cells [25]; however, they are known to obtain energy under its depletion via multiple ways in order to increase their survival [26][27][28]. Therefore, the limiting nutrient in the model for the cell survival was oxygen, tumor cells being more or at least equally resistant to its depletion that the normal cells, that is, ≤ 1. The function of cell death rate M h (ω) was chosen to be smooth, tending to its maximum value under the full absence of oxygen, and equal to zero for the levels of oxygen which exceed the critical value ω * . Normal cells also died in the presence of tumor cells, which was introduced in the model to coarsely reflect two processes: the inability of normal cells to remain viable in acidic tumor microenvironment [29] and the degradation of capillary network inside the tumor [30].
Tumor cells were able to migrate throughout the tissue, which was governed by a diffusion-like term. Directed migration of tumor cells was neglected [31]. The convective terms described the bulk motion of the tissue elements, the velocity field I being determined by the dynamics of tumor cells. The drainage of necrotic tissue was neglected. Due to the assumption of the constancy of the total density of tumor cells, necrotic tissue, and normal cells, which was normalized to unity, the following expression for the gradient of I was obtained:

Dynamics of Nutrients
Model dynamics of both nutrients accounted for the same physiological processes-inflow of nutrients from the capillary network into the tissue, their consumption by tumor and normal cells and their diffusion within the tissue, the latter being much faster for oxygen. The form of the terms for the inflow of oxygen and glucose differed due to significant distinctions in the mechanisms of their blood and transvascular transport. The inflow of glucose is governed primarily by the process of passive diffusion through the pores in the walls of capillaries [32]. Therefore, the rate of the glucose inflow was taken to be proportional to the difference in glucose concentrations in blood and in tissue, and to the capillaries density, which, as it was mentioned above, was assumed to be in linear relationship with the density of normal cells. Glucose concentration in blood was considered to be constant and was normalized to unity.
Oxygen, as lipid-soluble substance with low molecular weight, passes directly through the capillary walls and flows into the tissue at much greater rate that glucose. Oxygen levels in arterial and venous blood differ more than twice even under normal conditions [33], which implies that its blood concentration should not be treated as constant. Moreover, the inflow of oxygen into the tissue is not proportional to the difference of its concentrations in capillary blood and tissue due to the complicated blood transport of oxygen, which molecules are carried in blood in two forms-bound to hemoglobin and unbound from it. Overall, the used term for the oxygen inflow assumed that the rate of this process is proportional to the difference between the fraction of oxygen-saturated hemoglobin under two values of unbound oxygen concentration-the one in arterial blood, which enters the capillaries, and the one in tissue. The function S(ω) represented oxygen-hemoglobin dissociation curve, the form of which is well-known in physiology [34]. For more detailed explanation of the assumptions, which underpinned the term of oxygen inflow, we refer the readers to our previous work [35].
The nutrient consumption was described via the terms of the classical Michaelis-Menten type. Tumor cells are known to consume nutrients much faster that normal cells, in order to support their proliferative activity, therefore, The rate of oxygen consumption by tumor cells fell down to the rate of oxygen consumption by normal cells under the decrease in tumor cells proliferation rate, caused by the glucose shortage [36].

Numerical Solving of Tumor Growth Model
The set of Equations (2) was solved numerically with assumption of the spherical symmetry of the tumor. The size of the computational region L was adjusted in order to be sufficiently small to spare computational time without imposing noticeable edge effects. The convective flow speed was set to zero at the left boundary, which represented the center of the tumor, where initially r = 0. For all the variables, the zero-derivative boundary conditions were used at both edges. The following initial conditions were used for tumor cells, normal cells and necrotic tissue: Equations for glucose and oxygen were considered in the quasi-stationary approximation, due to the fast dynamics of these variables, and were solved using the tridiagonal matrix algorithm. The equation for normal cells was not solved explicitly, the relation h = 1 − n − m being used for determining their density. For tumor cells and necrotic tissue the method of splitting into physical processes was used, that is, kinetic equations, migration equation and convective equations were solved successively during each time step. The implicit Crank-Nicholson scheme was used for the tumor cells migration equation. Convective equations were solved using the flux-corrected transport algorithm with the implicit anti-diffusion stage. The kinetic equations were solved by the simple explicit Euler method, which was justified by the relative smallness of the used time steps, adjusted for the solution of the transport equations. The flux-corrected transport algorithm was introduced in Reference [37], while other classical methods were described in many books (see, for example, Reference [38]). The choice of the time and space steps for different simulations is justified in Appendix A.
The computational code was implemented in C++ and can be found in the Supplementary Materials.

Equations for Radiotherapy
For the description of radiotherapy (RT), we relied on the classical linear-quadratic model, which was discussed in Section 1. We accounted for two effects that lead to spatiotemporal heterogeneity of the radiosensitivity of tumor cells. The first one is oxygen enhancement effect, which was introduced herein in the form presented in Reference [39], where the corresponding terms were deduced from the experimental data. The second effect is the decrease in radiosensitivity of quiescent cells, which was considered herein with the assumption that it fell along with the decrease in the cells' proliferation rate. We neglected the duration of each irradiation and assumed that the number of cells and the density of necrotic tissue changed in result of it instantaneously, which was realized in a code in a straightforward manner. We did not consider explicitly the death of normal cells due to RT, however, the total damage of the normal tissue was the crucial parameter for the optimization of RT fractionation, which is discussed in Section 2.3. Overall, the equations that expressed the densities of tumor cells and necrotic tissue after a single irradiation with the dose D through their values before it were as follows:

Optimization of Radiotherapy Fractionation
The task of finding the optimized fractionation of RT was formalized the following way. During all the simulations, the first irradiation was performed at the moment t = t 0 , when tumor radius, evaluated as the maximum space coordinate, at which n + m ≥ 0.1, reached 1 cm. We considered the RT schemes, which consisted of 42 doses of radiation, some of which could be zero, which were administered successively at 24 h interval. Therefore, each scheme D could be expressed as a vector of non-negative numbers, representing the values of doses, expressed in grays: As the standard reference scheme we used the following vector that corresponded to one of the typical courses in clinical practice, consisting of 30 doses of 2 Gy, delivered every weekday over six weeks [3]: However, we did not impose the condition of the obligatory presence of days-off in the tested fractionation schemes. All the considered schemes had to satisfy two constraints, related to the normal tissue damage: The first inequality was analogical to the condition that the biologically effective dose, delivered to the normal tissue, could not exceed its value for the standard fractionation scheme. The second inequality corresponded to the acute reactions and indicated that each dose could not exceed a certain threshold. The aim of the search was to find the scheme D, which would decrease the value of the following objective function F(D) as much as possible: wheren andr are the normalization parameters of the model for the number of tumor cells and length. Thus, the aim of the optimization procedure was to find the fractionation scheme, leading to the most efficient eradication of tumor cells, which should correspond to the increase in the tumor cure probability (TCP). The following formula was used for the estimation of TCP: which can be interpreted as the fraction of eradicated tumors among the identical tumors that have undergone the same treatment, under the assumption that the number of surviving cells throughout these tumors at the end of the treatment follows a Poisson distribution with the average of N.
For the search of the optimized RT fractionation schemes and the optimized values of the objective functions, Algorithm 1 was developed and implemented in the program code. Its repeating steps 2 and 3 represented an adaptation of the classical gradient descent method for the considered problem. Like the classical method, these steps could find only local optimum, and the aim of step 4 was to try to further optimize this result. The meaning of the actions, performed during step 4, is explained in Section 3.3. By themselves, steps 2-4 could yield different results depending on the initial scheme. By testing different initial schemes under various model parameters, we found out that these steps most often produced the best results with the use of the most optimal uniform fractionation scheme as the initial. The search for such scheme was performed during the step 1. The procedures of normalization of the schemes, with the aim of their compliance with the above-mentioned restrictions, were performed iteratively.

Algorithm 1: Optimization of dose fractionation for radiotherapy.
Data: Distribution of variables of tumor growth model, governed by Equations (2), at t = t 0 . Run simulation of radiotherapy (RT), governed by Equations (5), with fractional doses D st ; Remember the value of the objective function F st ; Step 1. Search for the optimal uniform fractionation scheme: j = 42; Step 2. Search for the "gradient": Normalize D j according to Equation (7), not altering the doses, equal to D max ; Simulate RT with the scheme D j ; Remember F j ; Step 3. Going down the "gradient": Normalize D n according to Equations (7) and (8); Simulate RT with the scheme D n ; Remember F n D ; if F n D < F n−1 D then if k n = 1 then n = n + 1; Step 4. Trying to improve the final part of the scheme: if D opt j max < D max then Result: D opt , F opt .

Parameters
The values of the model parameters, with which the simulations were performed, are listed in Table 1. For some parameters, three values are designated, which belonged to the parameter sets that were assumed to correspond to different levels of tumor malignancy-high, intermediate and low. The way the values of these parameters changed with the increase in tumor malignancy reflected the stronger manifestation of some of the hallmarks of cancer: self-sufficiency in growth signals and insensitivity to anti-growth signals, evading apoptosis, stimulating growth of new vessels and invasion into normal tissue [40]. The dimensionless model values of the parameters were the approximations of their normalized values, which were obtained with the use of the following normalization parameters-t = 1 h for time,r = 10 −2 cm for length,D = 1 Gy for radiation dose, g = 1 mg/mL for glucose concentration,ω = 1 mM for oxygen concentration,n = 3 × 10 8 cells/mL for maximum density of cells. The latter value was taken from the experimental work on the in vitro growth of multicellular tumor spheroids [36]. The values of the proliferation rate of tumor cells B and their nutrients consumption rates Q g n and Q ω n were also estimated according to the data of this work with the assumption that these values should be proportionally diminished during the growth of a relevant tumor in tissue. Furthermore, it was assumed that B, Q g n and Q ω n should proportionally increase with tumor malignancy. (α/β) h alpha-beta ratio for normal tissue 3 [3] D max maximum fractional dose 5 [11] δ S the amount of radiation dose added to each fraction 0.2 see the text during the search for the "gradient" The death rates parameters M and were assessed based on experimental data on cell behavior under extreme nutrient deprivation [41]. The death rate of tumor cells fell with the increase in tumor malignancy, reflecting the increased tolerance of malignant cells to nutrient deprivation. The coefficient of high malignant tumor cells' motility D n was an order of magnitude lower than the value which corresponds to high malignant glioma, one of the most invasive types of cancer [42]. The parameter of glucose inflow P g was estimated as the product of the experimental values of permeability of capillaries to glucose and normal capillary surface area density for human muscle [32]. The values of normal cells' rates of nutrients consumption for human muscle at rest were also used. The oxygen inflow parameter P ω was adjusted so that the initial oxygen concentration lied within its normal range for human muscle at rest [46]. The values of P g and P ω for low malignant tumor were obtained under the assumption of absence of tumor-induced angiogenesis, that is, the formation of new blood vessels. Their increase with tumor malignancy reflected the stimulation of angiogenesis by tumor, and followed different trends due to the following reasoning (see our previous work [35] for details). As was mentioned in Section 2.1.2, the inflow of glucose should be proportional to the density of capillaries and to their permeability, both of these parameters increasing due to the tumor-induced angiogenesis. On the contrary, oxygen inflow in tissue is not affected by the alterations in the number and sizes of capillaries' pores. Moreover, it should be at first approximation proportional to the blood flow rather than to the capillaries' density. Therefore, the inflow of oxygen should increase much slower than the inflow of glucose, in the result of angiogenesis.
Radiosensitivity parameters are known to vary dramatically between various tumor cell lines [50] and, moreover, can significantly differ even for tumors of the same type [51]. Therefore, the parameters of tumor cells' radiosensitivity were varied within a physiologically justified range in order to investigate the potential for optimization of RT fractionation under various response of tumor to the treatment. At that, the alpha-beta ratio for tumor cells was kept constant and equal to 10 [3]. The alpha-beta ratio for normal tissue (α/β) h was close to its experimental value for human muscle [3]. The value of the ratio of radiosensitivity of quiescent and proliferating tumor cells k for low malignant tumor was selected based on the fact that the radiosensitivity of quiescent and proliferating normal cells can differ five or more times [52]. Experimental data suggests that such difference should be leveled with the increase in tumor malignancy [53,54], therefore, the value of k increased up to unity with the increase in tumor malignancy. Other parameters of the optimization algorithm were adjusted manually in order to decrease the computational time as much as possible without noticeable distortion of the solution.  (6)) and the values of tumor cells' radiosensitivity α = 10β = 0.1. In the same way as it happened under other parameter values as well, after initial phase of exponential growth the living tumor cells concentrated at the tumor rim, closer to the source of nutrients, which in this model was considered to be proportional to the density of normal cells h = 1 − n − m. Due to the active consumption of glucose and oxygen by tumor cells, the proliferation rate of deeper located cells declined, and more deeper located cells began to die. Therefore, most of the volume of sufficiently large tumors was occupied by necrotic tissue, as Figure 1a demonstrates. Figure 1b corresponds to the day of the first irradiation. The radiosensitivity of tumor cells increased from the center of the tumor to its rim, where the most actively proliferating cells were situated and the concentration of oxygen was the highest throughout the tumor. During the first irradiation, approximately 79% of tumor cells survived in the nutrient-depleted regions, while only 12% of tumor cells survived in the outer layers of the tumor rim. The death of tumor cells due to this and following irradiations led to the gradual rise of the levels of nutrients that, in its turn, resulted in increase of radiosensitivity of tumor cells. Figure 1c demonstrates the 10 th day of RT, by which eight irradiations were performed, the number of tumor cells decreased by 20 times and the levels of glucose and oxygen in the tumor center reached correspondingly 32% and 74% of their values in the normal tissue. Therefore, the effective radiosensivity of tumor cells sufficiently increased and became almost constant throughout the tumor. During the next irradiation, ≈10.8% of tumor cells died in the outer layers of the tumor rim, and ≈11.5%-in the deeply located regions. Figure 1d shows the dynamics of the total number of tumor cells N(t) and the number of proliferating tumor cells N p (t), estimated by the following formulas:

Simulation of Tumor Growth and Radiotherapy
g(r, t) g(r, t) + g * r 2 dr.
In the considered simulation, the total number of tumor cells decreased during the RT course from ≈0.43 billion to ≈18 cells, at that the fraction of proliferating cells N p (t)/N(t) increased from ≈26% to ≈99.3%. Due to the nature of the model, tumor regrowth always happened after treatment, even if the total number of tumor cells that remained after RT, was less than one. Figure 1e shows the dynamics of the average oxygen pressure within the viable tumor rim, expressed in mmHg and estimated as where Ω(t) ≡ [nr 3 · 4π L 0 ω(r, t)n(r, t)r 2 dr]/N(t); (12) along with the dynamics of the oxygen pressure at the tumor center, that is, at the point r = 0. The former quantity cannot be straightforwardly measured in experiment, however, it can give a better estimation of the efficiency of the first irradiations. The latter quantity was up to five orders of magnitude smaller during the free tumor growth that pO 2 , but their values became almost equal during RT.  Figure 2 shows the numerically obtained values of tumor control probability (TCP), estimated by Equation (10), for standard and optimized fractionation schemes for different types of tumor. Some of the optimized schemes are also shown, which were obtained under the designated values of tumor radiosensitivity α = 10β. The solid lines interpolate the graphs of TCPs by the functions of the following form:

Optimization of Radiotherapy Fractionation
where γ and α cr are the fitting parameters. Let us introduce the following quantity as a high-level estimate of the effectiveness of RT fractionation optimization for different tumor types: where α st cr and α opt cr are the fitting parameters in Equation (13) for standard and optimized RT fractionation schemes of the considered tumor type. Roughly speaking, this quantity denotes the increase in the curative range of the values of tumor radiosensitivity parameters due to the RT fractionation optimization.     Table 1. The proliferation rates of these tumors' cells related as 4:2:1. This fact by itself stimulated increase in tumor cell number with the tumor malignancy. Other factors that contributed to the same effect were the decrease in tumor cells' death rate and the enhancement of nutrient inflow. However, at the moment, when these tumors reached the radii of 1 cm-at which the beginning of RT took place-their numbers of cells did not differ drastically and were equal to ≈0.51 billion, ≈0.43 billion and ≈0.43 billion correspondingly. This was due to the fact that tumor cells' rates of nutrients consumption also increased with the tumor malignancy, significantly reducing the pools of proliferating and alive cells. The presence of tumor cell motility in the cases of IM and HM tumors also led to a slight decrease in their numbers of tumor cells. However, it should be noted that in general case the variation of this parameter may have ambiguous effect on the total amount of tumor cells, which was discussed in our previous work [55].
The most prominent feature of all of the optimized RT fractionation schemes, obtained for these tumors for all values of α, is the fact that they could be clearly divided into two stages. The first stages were comprised of non-equal doses, which were noticeably less that the maximum fractional dose of D max = 5 Gy. In every case they lasted until the following two quantities became smaller than one by no more than several percent: the first quantity was the ratio of the oxygen level inside the tumor to its value for the normal tissue Ω(t)/ω(L, t) (see Equation (12)); the second quantity was the fraction of proliferating tumor cells N p (t)/N(t) (see Equation (11)). Thus, to the end of the first stages tumor cells radiosensitivity became close to its maximum level throughout all the pool of tumor cells. The aim of the second stages was to get advantage of the increased radiosensitivity of tumor cells. Therefore they represented a uniform sequence of doses, which were equal or close to the maximum fractional dose. The aim of the first stages was, therefore, to reach close to maximum sensitivity of tumor cells, reducing both the effective dose, delivered to normal tissues during the first stages (see Equation (7)), and their duration. The first aspect is crucial for the opportunity of increasing the number of more efficient irradiations during the second stage. The decrease of duration of the first stages, as well as of the whole courses, is crucial due to the process of tumor cells repopulation-that is, the shorter the treatment, the more effective it should be under the same amount of eradicated cells, since fewer acts of cell division should take place during its course.
For every considered value of α, the optimized RT schemes became longer with the decrease of tumor malignancy. This was due to the fact that quiescent tumor cells were more radioresistant in IM and especially LM tumors, because of the smaller value of the parameter k, which was equal to 1, 0.5 and 0.2 for HM, IM and LM tumors. Consequently, with the decrease of tumor malignancy it took longer time for radiation to eradicate enough tumor cells for the necessary increase in the level of glucose that would convert the remaining tumor cells into proliferative state, thus increasing their radiosensitivity. Therefore, the durations of the optimized RT courses were the shortest for HM tumors, which allowed to significantly decrease the influence of cell repopulation on the outcome of the treatments. Overall, the efficiency of RT fractionation optimization increased with tumor malignancy: the values of ∆α for LM, IM and HM tumors were ≈0.008, ≈0.014 and ≈0.028 correspondingly. Quite surprisingly, the three interpolated functions of optimized TCPs turned out to be very close to each other. For every value of α the optimized treatments for each of the three tumor types yielded very close TCP values, differing by no more than 7%, and the corresponding values of α cr for these tumors differed by less than 1%. On contrary, the efficiency of standard RT schemes significantly declined with the increase in tumor malignancy for every value of α. Of note, this happened despite the fact of greater radiosensitivity of HM tumor cells under glucose deficiency, and was due to their increased repopulation rate under normal level of glucose. It should be noted, however, that such qualitative outcome might change under different parameter values.
Since the presence of a radioresistant population within a malignant tumor may be of significant practical interest [56], we performed an analogical set of simulations for the fourth set of parameters, which corresponded to HM tumor with the only modification of decreased radiosensitivity of quiescent cells k = 0.2. The corresponding results are shown in Figure 2d. Expectedly, the standard fractionation schemes were the less effective for this case among the four considered parameter sets under any value of α. The length of the first stages of the schemes were significantly increased. Moreover, for the lowest values of α there were no second stages in the optimized schemes, since such weak therapies could not eradicate enough tumor cells for sufficient increase in the level of nutrients. The second stages became well-pronounced at α = 0.1 and their lengths increased with the increase of α, as it happened under another parameter sets as well. However, the optimization procedures led to the value of ∆α ≈ 0.028, close to its value for HM tumor with k = 1.

Efficiency of the Optimization Algorithm
As was already mentioned, Algorithm 1 can provide at best the local optimal schemes. It should be noted that for some of the parameter values slightly more effective schemes that the ones, shown in Figure 2, were obtained by manual manipulations or by replacing the scheme, produced by step 1, with other initial schemes for the following steps. Of especial importance is the fact that these manipulations did not lead to any noticeable change of the TCP graphs. However, such cases, in which further optimizations can be provided, are worth being noticed. The simplest case is removing of initial zero fractions produced by the algorithm, like in the depicted optimized scheme for LM tumor with α = 0.07. Less trivial case took place when the search for the optimal uniform fractionation scheme, performed during step 1, yielded a bimodal distribution. In such case, setting another initial scheme for the following steps could optimize the result. For example, the optimized scheme for HM tumor with α = 0.19 is obtained with the initial uniform scheme with 9 irradiations of ≈4.47 Gy, produced by step 1. Its replacement by the uniform scheme with 13 irradiations of ≈3.53 Gy, which by itself results in 30% greater minimum number of tumor cells that the previous uniform scheme, allowed to produce an optimized scheme, resulting in halved minimum number of tumor cells, compared to the depicted optimized scheme. It did not lead to a noticeable change of TCP, since it was already close to 100%, however, it should be noted that such manipulation might turn out to be important under some other used parameter values.
Another possibility for the slight optimization of some of the schemes is the increase of several doses in the second stage of the scheme, if they are less than D max , by the expense of another doses. Such manipulation was not included in the algorithm for simplicity, since it was checked to provide only very slight improvements. However, the closeness of the doses of the second stages to the maximum fractional dose was by itself significant. For IM and LM tumors with sufficiently high values of α the first three steps of the algorithm by themselves yielded the schemes with doses of the second stage, equal to each other, but significantly less than D max . Further slight variations of the fractionation, performed during steps 2 and 3, were ineffective. This result corresponds well to the findings, described previously in other studies, which considered homogeneous and constant radiosensivity of tumor cells [15,16]. The aim of step 4 was to redistribute the fractions of the second stage, keeping them equal, which always resulted in close to maximum values of single doses. Two of the schemes, produced for IM and LM tumors with α = 0.22 without the use of step 4, are shown in Figure 2 via gray dots. The improvements, introduced by step 4, allowed to decrease the minimum number of tumor cells more than threefold in both cases.
Of note, under neglect of the constraint on the maximum fractional dose (see Equation (8)), Algorithm 1 generally produced the optimized schemes, consisting of longer first stages with smaller doses and second stages, consisting of a few strong irradiations, which most frequently were a couple of doses close to 10 Gy. Furthermore, easing the constraint on normal tissue damage (see Equation (7)) by increasing the value of (α/β) h resulted in more and more shorter optimized schemes with greater doses. These results are not surprising; moreover, Equation (1) immediately suggests that under absence of any time-dependent effects and any constraints increasing a single fraction should always be more efficient that its fractionation. Thus, in agreement with the clinical concepts, discussed in Section 1, the presented model indicates that dose fractionation is dictated by the constraints on normal tissue damage and by the alterations of tumor cells' radioresistance (which involves both reoxygenation and redistribution of proliferative states), while cell repopulation restricts the efficiency of the fractionation scheme.

Discussion
In this work, we presented a spatially-distributed continuous mathematical model of solid tumor growth and treatment by fractionated radiotherapy (RT). The model explicitly accounted for three factors that influence the efficiency of RT fractionation schemes-tumor cell repopulation, reoxygenation and redistribution of proliferative states. The main goal of this study was the search for optimized fractionation protocols that would increase the tumor cure probability under the constraints of maximum normal tissue damage and maximum fractional dose. For this goal, a special algorithm was developed. Its first step compared different uniformly fractionated RT schemes. By itself it showed that the length of an optimal treatment should grow with the decrease of the relative radiosensitivity of non-proliferating tumor cells, which occupied the main part of tumor at the beginning of RT. The next two steps of the algorithm represented an adaptation of the classical gradient descent method. They suggested dividing all the fractionation schemes in two stages. Fractionation during the first stages followed different trends, but their aim was always to spare the doses for the second stage, at that eradicating enough tumor cells for the levels of nutrients to increase close to their normal levels. Such approach brought the radiosensitivity of the remaining tumor cells close to the maximum level. If this goal was possible to achieve, the second stages began, which consisted of large equal doses. Of note, the qualitatively similar recommendation of a dose boost during the final part of the therapy was suggested in a previous theoretical study, which used another constitutive assumptions and considered treatment regimens, in which a patient is treated in several sessions, separated by weeks or months [57]. The aim of the fourth step of the algorithm was to optimize the number of fractions and the doses during the second stage, which previous steps could not do. Optimized RT fractionation schemes did not contain days-off, unlike standard clinical schemes. It should be noted that the current model neglected the change of radiosensitivity during the cell cycle for the proliferating cells, as well as the fact that cells do not die immediately due to irradiation [58]. The introduction of these aspects into the model may somehow alter the appearance of the optimized schemes, found by the used algorithm. However, they would hardly affect the main qualitative findings of this study.
We performed the optimization procedures using the objective function of minimum number of tumor cells during the treatment. Certainly, other objective functions can be incorporated within the introduced algorithm. One of them, which we implemented as well during the study, is the delay in tumor regrowth, which always happened during the simulations (see Section 3.1). The model simulations showed that only a rather moderate increase of it can be achieved for the considered parameter values under the restriction of maximum treatment duration of 6 weeks. However, the results suggested that sufficient tumor growth delays might be obtained for much longer treatments of slowly-proliferating tumors, which is in agreement with other theoretical studies [59,60]. Moreover, the presented model did not account for the drainage of necrotic tissue, which should be crucial for such problem. Furthermore, in this light, a very interesting augmentation of the model may be an introduction of concurrent antiangiogenic therapy, which not only influences the drainage of necrotic tissue [61], but also affects the intratumoral oxygen level in a complicated manner [62]. These factors should influence the outcome of combined radiotherapy and antiangiogenic therapy. Therefore, their consideration should provide insights into the ways of optimization of such treatment. This task lies within the scope of our future plans.
We tried to incorporate in the model the most basic features of malignant tumors, relevant for the considered task, and we varied the model parameters, assuming that some hallmarks of cancer should manifest themselves stronger with the increase in tumor malignancy. Certainly, this was a very general approach, and the results of this work are of purely qualitative nature. In our opinion, an important outcome of this study is the theoretical proof of concept that non-uniform RT fractionation schemes may be significantly more effective that uniform ones, due to the time and space-dependent effects. We hope that the presented algorithm would be useful for further, more specific, tasks. At that, one of the important aspects to be focused on is the consideration of a separate radioresistant population of cancer stem cells. Its determining role in optimization of RT treatment was already noticed in previous studies on mathematical modeling [18,63]. Funding: This work was supported by RFBR according to the research Project no. 19-01-00768. Numerical simulations were prepared with the support of the "RUDN University Program 5-100".

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Choice of Discretization
The choice of discretization for the numerical simulations was dictated by two goals. The first goal was sparing computing resources, the importance of which was increased due to the fact that a lot of simulations of different treatment courses should have been preformed for a single optimization task. Namely, the optimization tasks, discussed in Section 3.2, required an average of ≈882 treatment simulations, with their minimum and maximum numbers of 66 and 3560. The second goal was providing sufficient accuracy of the solution. Since the model described biological objects, for which significant variability is natural, there was no need to pursue the degree of accuracy that would be necessary for solving, for example, physical problems, and the main aim was to capture the qualitative behavior of the model properly.
The considered low malignant tumor had zero cell motility, therefore, the only transport term for the tumor cells was the convective term in its case. As it was mentioned in Section 2.1.3, the corresponding equation was solved by the flux-corrected transport algorithm [37]. This algorithm is of indeterminate order, and the crucial condition for its workability is |I(r, t) dt dr | < 1 2 ∀r ∀t, where I is the field of the convective flow speed, dt and dr are the time and space steps. The flux-corrected transport algorithm consists of two stages. The first stage solves the convective equation, maintaining the total number of tumor cells and the non-negativity of their density profile. However, it introduces erroneous diffusion, reduction of which is the aim of the second stage.
The simulations for the intermediate malignant and high malignant tumors included Crank-Nicholson method for the solution of tumor cell migration equation. Its main deficiency is the introduction of the spurious oscillations, which amplitude increases with the increase of D n dt/dr 2 , where D n is tumor cell motility. Obviously, the decrease of the time step under constant space step should increase the accuracy of the Crank-Nicholson method. However, such action would play an ambiguous role on the accuracy of the flux-corrected transport algorithm, since more frequent calculations within a single unit of time should amplify its total erroneous diffusion. Figure A1 shows the dependence of the low malignant, intermediate malignant and high malignant tumor growth speeds V, expressed in mm/week (which implies its multiplication by the factor of 16.8), on the time and space steps, equal to each other. The tumor growth speed was estimated via the method, described in our previous work [55], as the asymptotic value of the rate of change of the tumor radius, which was evaluated as the maximum space coordinate, at which n + m ≥ 0.1. The low malignant tumor growth speed changed non-monotonically with the refinement of discretization, reflecting the increase of the erroneous diffusion, introduced by the flux-corrected transport algorithm. This effect was not pronounced for the intermediate malignant tumor. This tumor had non-zero cell motility, therefore, the migration equation was solved in its case, and its accuracy fell under such refinement of discretization. Nevertheless, the tumor cell motility was sufficiently low in this case for this effect to remain unnoticed on this graph. However, this effect was strongly pronounced for the high malignant tumor with sufficiently high cell motility. Moreover, the utilized numerical approach turned out to be unstable for high malignant tumor under dr = dt = 0.01. Overall, based on these graphs and on the amount of computing resources, spent under different discretizations, time and space steps dr = dt = 0.1 were chosen to be used for all the simulations.