Algorithm of optimization of fractionated radiotherapy within its combination with antiangiogenic therapy by means of mathematical modeling

Antiangiogenic therapy is the anti-tumor treatment, that leads to the cessation of blood vessels formation and subsequent nutrient deprivation of the tumor. One of its frequently occurring effects is the transient improvement of tumor oxygenation, which lasts several days and leads to increased efficiency of successive exposure to radiation. It is followed by an escalation of hypoxia, and therefore the addition of antiangiogenic therapy to the whole course of fractionated radiotherapy has an ambiguous effect. Previously we have developed a mathematical model of such combined anti-tumor treatment. Herein we present an algorithm of optimization of radiotherapy fractionation, which aim is to find the most optimal distribution of irradiation doses in order to increase the efficiency of such combined treatment. We demonstrate an example of its work and discuss its further application.


Introduction
Radiotherapy (RT) is one of the classical options of cancer treatment. It acts by causing damage to cell DNA, and the presence of oxygen sufficiently enhances this process. In vitro experiments show, that radiosensitivity of cells in air and under hypoxia may vary up to three-fold [1]. This fact has led to the active search of methods for improvement of tumor oxygenation during RT. One of them is the combination of RT with antiangiogenic therapy (AAT), which is aimed at ceasing angiogenesis, i.e., the formation of new blood vessels, which supply tumor with nutrients. It was shown, that AAT can cause transient alleviation of intratumoral hypoxia, which lasts for several days at the beginning of the therapy [2]. This phenomenon was explained by the fact, that AAT causes normalization of structure of tumor microvessels and subsequent increase of their perfusion [3].
However, due to transient nature of this phenomenon, its manifestation does not guarantee that addition of AAT would enhance the overall efficiency of the whole course of fractionated RT. Moreover, the ultimate result -and primary goal -of AAT is the nutrient deprivation of a tumor, which means escalation of hypoxia in the long run [4]. Therefore, it is not surprising that preclinical and clinical studies demonstrate ambiguous results concerning the efficiency of such combined treatment [5,6]. In our recent work [7], we presented a physiologically justified mathematical model of tumor growth and treatment with RT and AAT and we demonstrated, that tumor radiosensitivity should be one of the factor which define the outcome of addition of AAT to RT: under its low values inclusion of AAT may prolong patient's survival; however, under high tumor radiosensitivity, AAT may compromise curative effect of RT. In that work we simulated the clinical scheme of RT fractionation, in which all doses of irradiation are equal. However, due to the influence of AAT on oxygen dynamics, it is possible, that increased efficiency of such combined treatment can be achieved by altering the distribution of irradiation doses over time. In order to investigate this possibility, we have developed an algorithm of optimization of radiotherapy fractionation with a limited total irradiation dose, which in practice is dictated by the negative effect of RT on normal tissues. Herein we present this algorithm, we demonstrate the result of its work within a model simulation under low tumor radiosensitivity and we discuss its further application.
There exist a number of works on optimization of antitumor tretament by means of mathematical modeling. However, mathematical modeling of combined radiotherapy and antiangiogenic therapy is nowadays by itself a little-studied field. The existing works on this topic consider only purely phenomenological dependencies, which significantly reduces their potential practical value [8,9]. To the best of our knowledge, there is only one work, in which an attempt to optimize such treatment was made [10]. It as well considers a phenomenological model wherein antiangiogenic therapy directly affects the parameters of tumor growth and radiosensitivity of tumor cells, and oxygen is not accounted for. The existing works on optimization of monoradiotherapy by the means of mathematical modeling mainly consider spatial redistribution of a single course of radiotherapy, that is dictated by clinical needs [9,11,12]. There exist several works on optimization of time scheduling of chemotherapy, which as well consider non-spatially distributed phenomenological models and rely on optimal control techniques [13,14]. Such approach would face tremendous difficulties when dealing with complex spatially-distributed physiologically-based models. In contrast, the approach presented herein can be integrated into the models of any reasonable complexity.

The model 2.1 Model of tumor growth and treatment
The model of tumor growth that we use is a continuous spatially-distributed model of reaction-diffusion-convection type, governed by partial differential equations. For full description of the model, including its equations and the values of parameters, we refer the reader to our work [7]. Herein we provide only the brief summary of the model.
The scheme of interactions between the model variables is presented in Fig. 1. The tumor grows in normal tissue, which consists of normal cells and includes preexisting normal microvasculature. Glucose and oxygen flow into tissue from microvasculature and are consumed by tumor and normal cells. Proliferating tumor cells push out surrounding tissues, which is one of the drivers of tumor growth. Microvasculature degrades within the tumor due to elevated solid pressure as well as chemical factors, which are not considered explicitly [15]. Under significant fall of glucose level, tumor cells pass into quiescent state, in which they, as well as normal cells, may die under long depletion of oxygen, turning into necrosis. Quiescent cells migrate throughout the tissue, which represents another driver or tumor growth. As these cells experience metabolic stress, they secrete vascular endothelial growth factor, or VEGF, which stimulates angiogenesis, i.e., formation of angiogenic capillaries [16]. They are described by a separate variable, which is introduced separately to account for "abnormality" of structure of tumor angiogenic capillaries, that results in their increased permeability for glucose as well as decreased blood flow and subsequent decreased oxygen inflow in tissue through a single capillary. Also, in presence of VEGF the structure of normal capillaries becomes abnormal and restores in its absence [17]. Antiangiogenic therapy (AAT) is modelled by instantaneous elimination of VEGF from the tissue till the end of a model simulation. Radiotherapy (RT) leads to direct instantaneous death of part of malignant cells -sensitivy of proliferating cells to it being an order of magnitude higher than that of quiescent cells. For RT description we rely on the classical linear-quadratic model, which is widely used in clinical applications [18]. The dependence of tumor cell radiosenstivity on oxygen concentration has been deduced from experimental data in [19]. The growth of tumor is simulated in one-dimensional region, with left boundary representing the center of the tumor, and right boundary being free and thus reflecting the assumption that the normal tissue does not hinder the tumor growth. Figure 2 Figure 2(c) refers to the same moment under combined RT and AAT, when the latter has already led to normalization of the absolute majority of capillaries, that has resulted in increase of oxygen inflow inside the tumor. In its turn, this has increased the efficiency of the first irradiation, which can be seen from the decrease in density of proliferating tumor cells.

Algorithm of optimization of radiotherapy fractionation
The task of finding the optimal protocol for fractionated RT under its combination with AAT was formalized as follows. The initial conditions are the distributions of the model variables at the time of the initiation of therapy, which takes place when the tumor radius has just reached 3 mm. Both RT and AAT start simultaneously, which is a standard clinical practice. The aim of the search is to find the distribution of fractional doses of irradiation, which would provide the maximum time of the tumor growth till the endpoint -radius of 7 mm -under the conditions, dictated by the clinical practice: 1) the total irradiation dosage is 60 Gy, 2) no more than one irradiation is carried out once every 24 hours, 3) the maximum number of fractions is 30 (or, equivalent to points 2 and 3, there are exactly 30 daily fractions, some of which may have zero dose).  For solution of this problem, the basic algorithm, summarized in Fig. 3, was developed and implemented within a single C++ code along with the code for the model of tumor growth and treatment. The utilized algorithm resembles the classical gradient descent method. For the initial scheme of RT, F actual = {F 1 , F 2 , ..., F 30 }, a scheme with 30 daily fractions with doses of 2 Gy is taken, which follows a standard clinical scheme with the exception that we neglect the fact that irradiations usually do not take place on days-off. The time, which it takes for the tumor to achieve the endpoint from the beginning of the therapy, T actual , is measured. Further the following steps are iteratively repeated. The first step is the search for "gradient", during which tumor growth times T i are measured for each of the 30 schemes of RT, which are obtained by adding the value of δ, which is the parameter of the algorithm, to one of the fractional irradiation doses, F i , i = 1, 2, ..., 30, and further normalizing the fractional doses so that their sum is 60 Gy. The second step is the descent along the "gradient", wherein a new scheme is introduced, in which each of the fractional doses is obtained by adding a value, proportional to T actual − T i to F i , i = 1, 2, ..., 30 (initial coefficient of proportionality is another parameter of the algorithm), and then all the fractional doses are once again normalized so that their sum is 60 Gy. Then the coefficient of proportionality is increased as long as the tumor growth time from the start of the treatment to the endpoint, T , increases. The scheme, which provides the longest T in this step, replaces F actual and these two steps repeat while the first scheme, obtained during the second step, does lead to additional delay in tumor growth time.

Results and discussion
The results of the work of the algorithm under a basic set of parameters of tumor model (see [7]) and parameters of optimization model δ = = 0.05 is illustrated in Fig. 4. The first figure represents the initial RT scheme. The others refer to the schemes, obtained at the end of the first, second, fourth, tenth and twentieth iterations, the latter scheme being the final result.
During the simulation, the utilized algorithm proposes the following modifications of the RT scheme. At the first iteration it suggests to enhance the first three doses at the expense of the middle ones and thus to take advantage of the transient hypoxia alleviation window. The several final doses also increase, and the next iteration suggests their further sharp growth, while the alterations in the first doses are on the contrary smoothed out. These two iterations provide approximately 7% increase in the tumor growth time T . During the consequent two iterations, the very first fractional dose isolates by achieving the value of 4 Gy, while the values of its neighbors fall down close to their initial value of 2 Gy. The twelfth, thirteenth and three of the last doses also rise up above their neighbors, and all these modifications add up about 3.5% of increase in T . Each of the following sixteen iterations provides only slight alteration in the RT schedule, however, in total they add up more than 4% to T , leading to overall almost 15% gain in tumor growth time. In the RT schedule, eventually suggested by the algorithm, the first several doses line up in a smooth curve, which begins with the first irradiation dosage of about 3 Gy, with the following doses slightly falling close to the initial value of 2 Gy during the first five days. In the middle of the course, starting from the fifteenth day, four days of absence of irradiation are situated, which are preceded by one strong and one mild irradiation of about 3.6 and 1.5 Gy correspondingly, and are followed by four days of intermittent schedule, in which each of the two irradiations of more than 3 Gy is paired by a following day of rest. The end of the course is marked by three fractions of increasing dose, the latter one achieving almost 6 Gy. However, it should be noted, that such dosage for fractionated RT can be excessive from a biological point of view.
These results of an implementation of the developed algorithm serve as its approbation rather than a ready proposal of optimization of investigated type of combined anti-tumor treatment, due to the following reasons. Firstly, the major disadvantage of similar algorithms is the fact that they are able to indicate only locally optimal solution. In attempt to find the globally optimal solution, we are planning to introduce the technique of simulated annealing within the scope of the algorithm, and to start its work with various initial RT schedules and various values of its parameters, δ and . However, it should be noted, that it will nevertheless be practically impossible to prove that the most optimal of the proposed schemes will indeed be the globally optimal. Secondly, we assume that the considered combined treatment can be further optimized by one promising modification of an algorithm, that is the possibility of the shift of the beginning of AAT, by which the escalation of hypoxia, that inhibits the RT efficiency in the end of its course, can be delayed. Thirdly and most importantly, investigation of the simulation of the tumor growth under only a single set of parameter values is of little use by itself. Although some qualitative recommendations can be derived from the presented optimized scheme -i.e., enhancement of the first and the last fractional doses and intermittent schedule in the middle of the course -it is unclear, which of them, and what new recommendations, will show up in the simulations under other values of model parameters. The solutions to the posed issues can not be quickly obtained due to the bulkiness of calculations, which is caused by the complexity of the tumor growth model, -e.g., the above discussed simulation, being performed in parallel on eight cores, took about a day on an average home PC.
Certainly, we understand that the results of mathematical modeling can be regarded as only preliminary ones. However, if the effectiveness of some qualitative recommendations for optimizing the RT fractionation will be demonstrated within a wide range of both tumorspecific and patient-specific model parameters, then, in our opinion, it will serve as a good reason for the necessity of their experimental verification.