Analytical solution to the radiotherapy fractionation problem including dose bound constraints

This paper deals with the classic radiotherapy dose fractionation problem for cancer tumors concerning the following goals: a) To maximize the effect of radiation on the tumor, restricting the effect produced to the organs at risk (healing approach). b) To minimize the effect of radiation on the organs at risk, while maintaining enough effect of radiation on the tumor (palliative approach). We will assume the linear-quadratic model to characterize the radiation effect and consider the stationary case (that is, without taking into account the timing of doses and the tumor growth between them). The main novelty with respect to previous works concerns the presence of minimum and maximum dose fractions, to achieve the minimum effect and to avoid undesirable side effects, respectively. We have characterized in which situations is more convenient the hypofractionated protocol (deliver few fractions with high dose per fraction) and in which ones the hyperfractionated regimen (deliver a large number of lower doses of radiation) is the optimal strategy. In all cases, analytical solutions to the problem are obtained in terms of the data. In addition, the calculations to implement these solutions are elementary and can be carried out using a pocket calculator.


Introduction
According to the World Health Organization [13], "radiotherapy is one of the major treatment options in cancer management. (...) Together with other modalities such as surgery and chemotherapy it plays an important role in the treatment of 40% of those patients who are cured of their cancer. Radiotherapy is also a highly effective treatment option for palliation and symptom control in cases of advanced or recurrent cancer. The process of radiotherapy is complex and involves understanding of the principles of medical physics, radiobiology, radiation safety, dosimetry, radiotherapy planning, simulation and interaction of radiation therapy with other treatment modalities".
Mathematical modelling has played an important role in understanding and optimizing radiation delivery for cancer treatment. Since its formulation more than 50 years ago, the linear-quadratic (LQ) model has become the preferred method for characterizing radiation effects. Usually, it is stated as follows: the survival probability S of a tumor cell after exposure to a single dose of radiation of d Gy is expressed as where α T and β T are two positive parameters describing the radiosensitivity of the cell, [9]. It is well known that these parameters depend on the type of radiation therapy chosen and also on the organ where the tumor is located [15]. More precisely, LQ model implies that if the initial size of the tumor is U , then it will be U · S after applying a d Gy dose. Let us recall that "Gray" (Gy) is the unit of ionizing radiation dose in the International System of Units.
LQ model has well documented predictive properties for fractionation/dose rate effects in the laboratory and "it is reasonably well validated, experimentally and theoretically, up to about 10 Gy per fraction and would be reasonable for use up to about 18 Gy per fraction", see [4]. Precisely, its range of validity is a key point of controversy; although there is a general consensus on the existence of this range, significant disagreements remain on the exact values of its limits. Let us illustrate this fact with other recent quotes: from [9], "in vitro (...) some authors suggesting significant discrepancies at doses of 5 Gy or above, while others report good agreement up to tens of Gy" and according to the French Society of Young Radiation Oncologists, "the dose / fraction must be between 1 and 6 Gy", see [11].
Given N doses, d 1 , ..., d N , eventually different, if we consider the stationary case (which means that neither the times of application of the doses nor the growth of the tumor produced between them are taken into account), the probability of accumulated survival is given by (1) S From here it is clear that the effect of radiation on the tumor is determined by the quadratic function On the other hand, radiation also affects healthy organs and tissues near the tumor (which we will denote by OAR, organs at risk, hereafter). In general, healthy organs and tissues receive less radiation than the tumor, which we will denote by δd, with δ ∈ (0, 1] being the so-called "sparing factor". The value of δ depends on factors such as the location and geometry of the tumor and also on the technology used to deliver the radiation, see [3]. It can be seen as a measure of the accuracy of the radiotherapy: if clinicians can keep the OAR almost unaffected by the radiation, δ will be about 0; if not, it will be larger, until reaching the value δ ≈ 1 at worst. Therefore, the effect of the radiation on the OAR is determined by the following function where α 0 and β 0 are the parameters associated to the healthy organs that we are trying to protect. Typical values for α 0 , β 0 , α T and β T can be found in the specialized literature such as [15]. These data come from conducting experiments and the corresponding adjustments (least squares regression) to achieve approximated values that best fit experimental data.
Let us now introduce the most common strategies for fractionating radiotherapy treatments: • Hypofractionation: Deliver higher doses of radiation on few occasions. This strategy results in a significant reduction of its duration. • Hyperfractionation: Deliver a large number of lower doses of radiation that are given more than once a day. In this paper we study the classic radiotherapy dose fractionation problem related to the following goals: a) To maximize the effect of radiation on the tumor, restricting the effect produced on the OAR (healing approach) in Section 2 and b) To minimize the effect of radiation on the OAR, maintaining enough effect of radiation on the tumor (palliative approach) in Section 3. The first novelty with respect to previous works in this framework concerns the presence of dose fraction bounds of the type 0 < d min ≤ d ≤ d max . On one hand, these restrictions are connected to the range of validity of the aforementioned LQ model and can be estimated for each particular tumor; on the other hand, they also take into account the minimum and maximum dose fraction that can be applied in practical situations in order to achieve a minimum effect and avoid undesirable side effects, respectively. It is well known that the dose per fraction value in most conventional treatments is around 2 Gy, see for instance [12]. Depending on the tumor type, the values of d min and d max can be tunned, but the reference values d min = 1 Gy and d max = 6 Gy could be a valid generic choice. In this sense one can not find in [12] a single treatment recomendation with a dose fraction less than 1 Gy and very few larger than 6 Gy.
The counterpart for imposing a positive minimum dose fraction is that the total number of radiations N should not be fixed a priori and this is the second important novelty of this work: N will also be considered another unknown of the problem and we will study the dependence of the solution with respect to N . Among others, this approach was followed by [8], but only for uniform dose treatments. Our approach here includes also the study for nonuniform protocols. A preliminary version of our results was presented by the second author as part of the academic project [7], except for the study of the dependence of the solution with respect to N which is new.
Summarizing, new analytical solutions in terms of the data are obtained for both problems, improving known results in the literature to the best of our knowledge, see for instance [10] and [14]. Moreover, the final calculations to implement these solutions are elementary and can be carried out using a pocket calculator.

Maximizing the effect of radiation on the tumor
The aim of this first problem is to determine the best strategy to maximize the effect of radiation on the tumor, while restricting the effect on the OAR (healing approach): where E T (N, d) is given by (2), E OAR (N, d) by (3) and d min , d max and γ OAR are a priori known positive parameters, that should be provided by the specialists. Roughly speaking, the restriction E OAR (N, d) ≤ γ OAR can be interpreted in the sense that the percentage of survival cells of the OAR should be greater than or equal to 100 exp (−γ OAR ). This is the classic fractionation problem that has been studied (with some variations) in several works, see for example the recent papers [2] and [14] (where more than one OAR is considered) and the references therein. The first novelty of our approach is that dose bound constraints are also included. Usually in the literature the lower bound 0 value is taken for d i and no upper bound is imposed; some exceptions are [5] and [6] where an upper bound is included, but not a positive lower bound. The danger of losing control of the tumor, due to the use of doses below a critical limit, has already been pointed out by [8]. In addition, our approach to the problem is more useful since the number of doses N is not initially set as in [2] and [14]. The case including repopulation was studied in [3], only assuming the non-negativity of d i .
From a mathematical point of view, this is a mixed optimization problem involving a discrete variable, N ∈ N, which corresponds to the number of radiation doses, and N continuous variables, d i ∈ R, 1 ≤ i ≤ N , which are the doses. In other words, this problem has the peculiarity of having a variable number of unknowns.
Also, we will denote by x the greatest integer less than or equal to x and by x the least integer greater than or equal to x. Finally, the notation d N = (d 0 , . . . , d 0 ) means that d N ∈ R N having all the N components equal to d 0 .
Proof. Taking into account the restrictions for (P 1 ) and that d min > 0, we have N ≤ ρ 0 . Hence, the set of feasible values for N is finite.
When ρ 0 ∈ (1, 2), the value N = 1 is still the only possible option. Consequently, we are faced with a maximizing problem of an increasing 1D function. Then, the solution will be given by the largest feasible value. In this case, it is quite easy to verify that the unique solution of (P 1 ) is the If ρ 0 ≥ 2 we can reduce the problem (P 1 ) to a finite collection of continuous optimization problems (P N 1 ) with fixed N given by: Firstly we will prove the existence of a solution for each problem (P N 1 ) (see Theorem 2), for N running [1, ρ 0 ] ∩ N and denote it by d N . Then, it is enough to take the pair N , d N from the finite set that maximizes the value of E T (N, d) as a solution to the problem (P 1 ).
The existence of a solution for each problem (P N 1 ) is proved below: Theorem 2. Let us assume d min > 0, ρ 0 ≥ 2 and N ∈ [1, ρ 0 ] ∩ N. Then the problem (P N 1 ) has (at least) one solution. Proof. For small values of N , more precisely N ∈ [1, λ 0 ] ∩ N, it is easy to verify that the solution for (P N 1 ) is the trivial one with maximum doses, that is, d N = (d max , ..., d max ). For other values, N ∈ (λ 0 , ρ 0 ] ∩ N, the existence of solution for (P N 1 ) follows from the classic Weierstrass Theorem, because we are maximizing a continuous objective function over a compact set.

Remark 1.
a) Let us point out that (P 1 ) is a nonconvex quadratically constrained quadratic optimization problem (even (P N 1 ) with fixed N ), because the objective is to maximize a convex function. Typically, this type of problems is computationally difficult to solve (see [14]), but here we will see that it can be done analytically. b) Unless all the components of the solution are the same, the uniqueness of solution fails: it is enough to take two indices i, j ∈ {1, ..., N } such that d i = d j and interchange these coordinates to generate a new solution. c) Under the condition ρ 0 < 1, it is apparent that the set of feasible points is empty and hence, the existence of solution for (P 1 ) fails. d) The hypothesis d min > 0 is also needed for proving the existence of solution for (P 1 ), as it can be shown through the following example: It is clear that for all feasible points we have Let us stress that here N can take any natural value, without restrictions. Inspired by Theorem 4 below, let us consider the sequence given by It is easy to check that it is feasible for N ≥ 4, Hence, the problem (P 10 ) can not have solution (N , d): the supremum value 10 can not be attained since it should happen that which is clearly impossible. e) The deficiency of this type of models to produce solutions (when d min = 0) that prescribe infinite doses with fractions tending to zero was pointed out by [8].
The following result provides a simpler version of the optimization problem for the bigger values of N : Theorem 3. Let us assume d min > 0, ρ 0 ≥ 2 and N ∈ (λ 0 , ρ 0 ] ∩ N. Then, the inequality constraint of the problem (P N 1 ) has to be active at d N , with d N being a solution for (P N 1 ). Proof. Arguing by contradiction, let us assume that the constraint is not active, i.e., Since λ 0 < N , we know that there exists some index j ∈ {1, ..., N } such that d j < d max . Then, for sufficiently small > 0, the point (d 1 , ..., d j−1 , d j + , d j+1 , ..., d N ) is feasible and satisfies but this contradicts the fact that d N is a solution for (P N 1 ). Hence, from now on, in this case we will consider the equality restriction Therefore, we deduce that and the objective function can be written as Based on this identity, we can directly simplify the formulation of the problem (P N 1 ) as follows: Proposition 1. Let us assume d min > 0, ρ 0 ≥ 2, N ∈ (λ 0 , ρ 0 ] ∩ N and denote iii) If ω δ = 0, then every feasible point for (P N 1 ) is a solution.

Remark 2.
a) The idea of this transformation can be found in [10] in the context of the problem (P 2 ) that we will study in the next section. b) Let us note that for the majority of tumors α T /β T > α 0 /β 0 and therefore, the case ω δ > 0 is more frequent in clinical practice. c) As a consequence of Proposition 1, we can appreciate the great difference between the cases ω δ > 0 and ω δ < 0: in the first one, to maximize the effect of radiation on the tumor we have to increase the total dose, while in the second the total dose remains minimum (see also the Subsection 3.3 below).

2.2.
Solving (P N 1 ). Let us begin by showing a 2D-example of previous problems that will inspire the general results of this section.
Example 2. Let us consider the following optimization problems: In the Figure 1 the points on the blue surface are those that satisfy the equality constraint and the intersection of blue and orange surfaces gives the curve on which to maximize or minimize. Visually one can guess that the unique solution to (P 2,+ 1 ) is located on the diagonal (more precisely, it is given by and there are two solutions of (P 2,− 1 ) lying on the boundary (specifically, In fact, what happens in previous example can be extended to the N −dimensional case. More precisely, we will prove that the solution for (P N,+ 1 ) is a vector with equal coordinates: Theorem 4. Let us assume d min > 0, ρ 0 ≥ 2 and N ∈ (λ 0 , ρ 0 ] ∩ N. Then, the unique solution to (P N, Proof. By using the Cauchy-Schwarz inequality, we have Therefore, for each feasible point it follows that Defining q(z) = β 0 δ 2 N z 2 + α 0 δz − γ OAR , previous inequality can be rewritten as (15) q Taking into account that the polynomial q can be factorized in the form , because all the components d i have to be positive. Now it is clear that the maximum value is achieved when N i=1 d i = z 2 . Combining this fact with (14), we deduce that (16) In this case, Cauchy-Schwarz inequality becomes (in fact) an equality and this is true if and only if all the components are equal, i. e., . Given d a solution of (P N,− 1 ), since the objective function and the functions defining the restrictions are C 1 , we can apply the Lagrange Multipliers Rule [1] to deduce the existence of real numbers Inspired by the 2D example, we will prove that d lies on the boundary of [d min , d max ] N . Let us argue by contradiction assuming that d i ∈ (d min , d max ), for all i ∈ {1, ..., N }. Then, thanks to (20) we deduce that µ i = 0, ∀i ∈ {1, ..., 2N }. In this case, (19) reads: If λ = 0, the identity (21) implies that α = 0, but this is not possible by (18). Therefore λ = 0 and from (21) we get In other words, we arrive to the solution of problem (P N,+ 1 ), contradicting our initial hypothesis about d.
Consequently, there exists (at least) one index j ∈ {1, ..., N } such that d j ∈ {d min , d max }. Without loss of generality we can suppose that j = N . Let us see that in this case we can reduce the dimension of the optimization problem (P N,− 1 ) by means for the following auxiliary problem: . Hence, using that d is a solution of (P N,− 1 ), we get Arguing exactly in the same form as before with the problem (P N −1,− 1 ), we deduce that there must be an index j ∈ {1, ..., N − 1} such that d j ∈ {d min , d max } and we can reduce again the dimension of the problem, obtaining a new problem with N − 2 unknowns. Repeating this process several times we arrive to the final 1D problem: Clearly, it is enough to solve the quadratic equation to get the solution.
Summarizing previous results, given N ∈ (λ 0 , ρ 0 ] ∩ N, the solution of (P N,− 1 ) has one of the following structures: being the unique positive root of the quadratic equation with ϕ 0 defined in (4). We can characterize the unknown value K as follows: a) In the case (22), by using the equality restriction we derive that Of course, this holds if and only if the right hand side is a natural number or zero. b) In the case (23), since ϕ 0 is an strictly increasing function in [0, +∞), we know that and using (24) we get that which means that Taking into account the conditions (25) and (26), it is easy to conclude that the latter structure (23) is more frequently found in practice than (22). Previous argumentations lead us to the following result: ) is also a solution for the problem We will use this property in the proof of Theorem 8 (see Appendix 1).

2.3.
Analytical solution for (P 1 ). As we pointed out, a solution of (P 1 ) In fact, combining previous results, we can avoid the calculation of most solutions for (P N 1 ) by studying its dependence with respect to N . This is the goal of the next results. Let us start by studying the less frequent case: when λ 0 = ρ 0 . Proof. In this case the set of feasible values for N is {1, . . . , For those values of N , the solution for (P N 1 ) has the form d N = (d max , ..., d max ). Among them, it is clear that in order to solve (P 1 ) only the one with the largest number of components is of interest; this is attained at N 1 .
We will continue to analyze the most common case: when λ 0 < ρ 0 . In the the trivial case ω δ = 0, the function to be minimized and the one defining the restriction are proportional. Therefore, we can be derive the following result: with N in the above set.
Proof. Due to the hypothesis ω δ = 0, we deduce straightforwardly that problem (P 1 ) is equivalent to Obviously, the maximum value is reached when the restriction becomes an equality. This can be achieved in several ways, such as the treatments with equal doses described in the proposition statement. Let us emphasize that Then, the unique solution to problem (P 1 ) is given by the pair N , d (13). In order to study the dependence with respect to N for these values, thanks to Proposition 1, it is enough to consider the auxiliary function Here, it follows easily that ψ(N ) is an strictly increasing function and then, it will take its maximum value in Finally, we will derive that (N 2 , d N 2 ) is the unique solution to problem To that end, let us consider the linear function Using that N 1 ϕ 0 (d max ) ≤ γ OAR = N 2 ϕ 0 (d 0 ) by the admissibility, we get that H α 0 β 0 δ ≥ 0. Also, taking into account (30), it can also be the assumption ω δ > 0 (see (9)), it follows that H α T β T > 0, which is equivalent to (29).
In the case ω δ < 0, the situation is more complicated and it is detailed in the next result. The proof is a little bit technical and is postponed to the Appendix 1: Theorem 8. Let us assume d min > 0, ρ 0 ≥ 2, λ 0 < ρ 0 and ω δ < 0. Then, a solution to problem (P 1 ) is given by one of the following pairs: Remark 4. a) Let us strongly highlight that all expressions included in Theorems 6-8 and Proposition 3 can be explicitly calculated from the initial data of the problem. Moreover, the calculations to implement these solutions are elementary and can be carried out using a pocket calculator. b) On the other hand, when ω δ > 0, the optimal value of N is the largest one within its range of possibilities (i.e. it is a hyperfractionated type treatment) with equal doses, while in the case ω δ < 0 the optimal value is the smallest one (i.e. it is a hypofractionated type treatment). In this last case, let us stress that not all doses have to be equal or large; in fact, some of them may be minimum. As far as we know, this structure is not usually cited in the specialized literature. c) One interesting case appears when α T β T < α 0 β 0 , because then ω δ < 0 for all δ ∈ (0, 1] and the optimal regimen is always of hypofractionated type, independent of the technology used and the geometry of the tumor. In practice this condition holds in some special cases, such as the prostate tumor, where α T β T ≈ 1.5 Gy, while α 0 β 0 = 2 Gy, see [10] and [15]. d) After Remark 2-b) (see also the Subsection 3.3 below), it is clear that the hypofractionated case (associated with ω δ < 0) is very convenient in the practice. Assuming that the other parameters are set, the condition ω δ < 0 can always be achieved by taking δ close enough to 0.
This last fact is related to increasing the precision of the radiotherapy process (for instance, by using cutting-edge technology). e) A related problem to (P 1 ) is studied in [3] and [6], where the number of dose fractions N is also an unknown, jointly with d. The framework for that problems is more general, because a repopulation term is included in the objective function, but only the lower bound d min = 0 is assumed. Furthermore, the determination of the optimal value for N is carried out in [6] by means of numerical simulations, while in [3, Theorem 2] it is done explicitly and the value N = 1 is obtained when ω δ < 0. In this last case, the single dose could be too large in practice (remember that no upper bound is imposed in [3]) and then more fractions would have to be tried until an acceptable one is found.
Once we have obtained the analytical expressions for the solution of problem (P 1 ), we can deduce very easily its dependence with respect to the parameters defining the problem. This will help us to know how to adjust these parameters in order to achieve a desired solution. Let us show a result in the direction: is a solution of (P 1 ). Then, N is an increasing function of γ OAR , decreasing with respect to α 0 , β 0 and δ and independent of α T and β T . Moreover, a) When ω δ > 0, N is also decreasing with respect to d min and independent of d max . b) When ω δ < 0, N is also decreasing with respect to d max and independent of d min .
In Table 1, we summarize the resolution of problem (P 1 ) in algorithm form, for the reader's convenience.

Minimizing the effect of radiation on the organs at risk
In this section, we will consider a problem closely related to that of the previous section: the goal of this second issue will be to determine the best strategy to minimize the effect of radiation on the organs at risk, while maintaining a minimum effect of radiation on the tumor. It is clear that this approach can be interesting (at least) for palliative therapies. Mathematically we formulate it in the following way: (2) and γ T is a given positive parameter. Of course, this is also a mixed optimization problem with N + 1 unknowns: the number of radiation doses, N ∈ N, and the value of the N doses, This problem has recently been studied in the outstanding work [10], but with fixed N and only imposing the non-negativity constraint for the doses. Moreover, in [10] it is also remarked that "The real interest of the present approach would be the determination of the optimum solution for ALGORITHM FOR SOLVING (P 1 ) DATA: α T , β T , α 0 , β 0 , d min , d max , δ and γ OAR (all positive, d min < d max and δ ≤ 1) and , with ϕ 0 (r) = α 0 δr + β 0 δ 2 r 2 .
IF ρ 0 ≥ 2, λ 0 < ρ 0 and ω δ > 0, the UNIQUE SOLUTION of (P 1 ) IF ρ 0 ≥ 2, λ 0 < ρ 0 and ω δ < 0, take N 1 = λ 0 and IF M ∈ N ∪ {0}, take K = M and d Also take N 2 = λ 0 and d In particular, the pairs (N , d N in clinical practice". As an intermediate step, we have achieved here the expression of the optimal value for N in terms of the parameters of the problem in this particular setting, see Table 2 for a detailed description, depending on the case. As we will see, the study for problem (P 2 ) can be carried out following the same argumentation to that of (P 1 ) with minor differences.
3.1. Existence of solution for (P 2 ). In the sequel we will denote , .
Our first observation concerns the existence of solution for (P 2 ): Theorem 9. Let us assume d min > 0. Then, the problem (P 2 ) has (at least) one solution.
Proof. It is analogous to that of Theorem 1, although here there are infinite feasible values for N : combining the restrictions, those such that N ∈ [λ T , +∞) ∩ N. We will begin by showing that for each fixed feasible value N , the associated problem (P N 2 ) has a solution, where For large values of N , specifically for N ≥ ρ T , the solution of (P N 2 ) is the trivial one with minimum values d min . Among them only the smallest value of N could have practical interest, i.e., ρ T . For the other values, when they exist, that is for N ∈ [λ T , ρ T ) ∩ N, the existence of solution for (P N 2 ) is a consequence of Weierstrass Theorem, once more because the objective function is continuous and the feasible set is compact. Therefore, for each value of N in that interval, let us consider a global solution for the problem (P N 2 ) that we will denote d N . Again, it is enough to take the pair N , d N from the finite set that minimizes the value of E OAR (N, d) as a solution to the problem (P 2 ).

Remark 5.
a) Once more, except if all the coordinates of d N are equal, the solution will not be unique, because two different coordinates can be permuted to generate a different solution. b) When ω δ > 0 (see (9)), the hypothesis d min > 0 is necessary for proving the existence of solution for (P 2 ), as we can see through the following example: Example 5.
(P 20 ) Let us consider the sequence given by It is easy to check that it is feasible for N ≥ 4, When N ∈ [λ T , ρ T ) ∩ N, we know that d N = (d min , ..., d min ) is not a solution of (P N 2 ), because it is not even feasible. Hence, we can simplify the problem (P N 2 ) arguing in a similar way as in the proof of Theorem 3. Theorem 10. Let us assume d min > 0 and N ∈ [λ T , ρ T ) ∩ N. Then, the inequality constraint of the problem (P N 2 ) has to be active at d N , being d N a solution of (P N 2 ).
From now on, the restriction will be taken as one of equality, that is, Again, applying the same procedure as for (P 1 ) in Section 2 we have and the objective function will read Now, it is clear that we can simplify the formulation of the problem (P N 2 ), as follows: Proof. It is enough to take into account that Once we have seen that (P N,+ ii) a solution for (P N,− 2 ) has one of the following forms: ,

3.2.
Analytical solution for (P 2 ). As a consequence of previous results we arrive to the main theorems of this section that completely clarifies the situation concerning the problem (P 2 ). Recalling that ω δ = α T β T − α 0 β 0 δ (see (9)), we will see that ρ T and the sign of ω δ are the determinant factors in this analysis.
Again, the case ω δ = 0 is easily solved, because the function to be minimized and the one defining the restriction are proportional. Hence, the following result can be derived as Proposition 3: Proposition 5. Let us assume d min > 0, ρ T > 1, λ T < ρ T and ω δ = 0. Then any feasible pair (N , d) with E T (N , d) = γ T is a solution to problem (P 2 ). In particular, the pairs (N , d N ) with N ∈ { λ T , . . . , ρ T } and d N = Let us now continue by studying the more frequent case ω δ > 0. Here, we have to distinguish two different situations, depending on ρ T ∈ N or not: The proof is similar to that of Theorem 7 and postponed to Appendix 2. Finally, the case ω δ < 0 is studied in the next theorem and its proof is detailed in Appendix 3.
Theorem 13. Let us assume d min > 0, ρ T > 1, λ T < ρ T and ω δ < 0. Remark 6. a) Once more, let us emphasize that when ω δ > 0 the optimal value of N is the largest one within its range of possibilities (i.e. it is a hyperfractionated type treatment), while in the case ω δ < 0 the optimal value is the smallest one (i.e. it is a hypofractionated type treatment). This classification in terms of ω δ was described in [10], considering nonnegative doses. b) For the hypofractionated case, the single exposure is chosen in [10] as the preferred one. But this dose could be too large in practice and then two, three or more fractions would have to be tried until an acceptable one is found. This fact is remarked in [3] by saying that the case ω δ ≤ 0 "needs careful consideration since the validity of the model may be limited if N is small and the dose per fraction is large". Under our approach, this difficulty is overcome and we get the optimal number of dose fractions directly and its value, as in the other case. c) The uniqueness of solution fails when ω δ < 0 because (as it was said) any permutation of the coordinates of the indicated solution provides a new one.
In the following examples we will show that all the above possibilities mentioned in Theorems 12 and 13 can appear in practice: Example 6. Let us continue with the same parameters than in Example 3, which are: α 0 = 0.04 Gy −1 , β 0 = 0.02 Gy −2 , α T = 0.05 Gy −1 , β T = 0.005 Gy −2 , d min = 1 Gy and d max = 6 Gy. i) If δ = 1, then ω δ = 8 > 0. Hence, we are considering the problem When γ T = 4, we can easily calculate that λ T ≈ 8.33 and ρ T ≈ 72.72. Applying Theorem 12-b) the unique solution for (P 21 ) is one of these two pairs: For γ T = 4.35, we have that λ T ≈ 9.06 and ρ T ≈ 79.09. Therefore, applying Theorem 13-b), we know that a solution for (P 22 ) is given by (N , d)  iv) Again, it can be shown that, for some specific values of the parameter γ T , the alternative exposed in Theorem 13-a) is true. In particular, for γ T = 4.375, we have λ T ≈ 9.11 and ρ T ≈ 79.54 and a solution for (P 22 ) is given by (N , d) with N = λ T = 10 and d = (1, 6, ..., 6   9 ).
As we did for problem (P 1 ), here we can deduce quite easily how the solution depends with respect to the parameters defining the problem (P 2 ).
is a solution of (P 2 ). Then, N is an increasing function of γ T , decreasing with respect to α T and β T and independent of α 0 , β 0 and δ. Moreover, a) When ω δ > 0, N is also decreasing with respect to d min and independent of d max . b) When ω δ < 0, N is also decreasing with respect to d max and independent of d min .
For the reader's convenience, we have summarized the complete algorithm for the resolution of the problem (P 2 ) in Table 2. 3.3. Equivalent treatments. We will finish this section by mentioning another application of previous results. Let us begin by introducing biologically equivalent treatments. Two treatments of radiotherapy with doses d 1 , ..., d N andd 1 , ...,dÑ are said to be biologically equivalent for a certain tumor with characteristic parameters α T and β T when they have the same effect, that is, When all the doses are equal for both treatments, previous concept leads to biologically equivalent doses (BED) (see [3]), that can be calculated very easily (see [11] for instance) from the equality A clearly interesting question is to determine (among all the equivalent treatments) which one uses the lowest total dose. As in this case we are not paying attention to the effect of radiation on the OAR, this can be ALGORITHM FOR SOLVING (P 2 ) DATA: α T , β T , α 0 , β 0 , d min , d max , δ and γ T (all positive, d min < d max and δ ≤ 1) and , with ϕ T (r) = α T r + β T r 2 .
IF ρ T > 1 and λ T = ρ T , the UNIQUE SOLUTION of (P 2 ) is the pair (N , d N ) with N = ρ T and d N = (d min , ..., d min ).
with N = ρ T and d N = (d min , ..., d min ) is the UNIQUE SOLUTION of (P 2 ).
Also take N 2 = ρ T and d 2 ). A SOLUTION of (P 2 ) is the pair (N , d N ) that minimizes E OAR between them.
In particular, the pairs (N , d formulated in mathematical terms as the following optimization problem: After our study, we deduce that the solutions to this problem are treatments of type (37) or (39), see Theorem 11-ii). Let us mention that these are the same than in the case ω δ < 0, see Theorem 13, and consequently the hypofractionated protocols are also optimal in this sense.

Conclusions
In this work, we have derived the analytical expressions for the optimal total number of radiations N and their specific doses d for problems (P 1 ) and (P 2 ). They are presented in Tables 1 and 2 in algorithmic form (see also Theorems 6,7,8,12 and 13) and there exists a clear parallelism between the structure of the solutions for both problems. We have proved that they essentially depend on the sign of the quantity For fixed N , this fact is well known in the literature and it has been reported several times in different frameworks (see for instance [10], [2] and [6]). Moreover, this is consistent with some clinical findings as noted in [10].
When ω δ > 0, we have shown that the optimal number of doses N are ρ 0 for (P 1 ) and ρ T or ρ T for (P 2 ), the upper values of their ranges of interest (i.e. hyperfractionated type treatments) with equal doses; while in the case ω δ < 0, the optimal values of N are λ 0 or λ 0 for (P 1 ) and λ T for (P 2 ), the lower values of those ranges (i.e. hypofractionated type treatments). In this last case, let us stress that not all doses have to be maximum; in fact, some of them may be minimum and at most one of them can take an intermediate value. For non-uniform protocols, the lack of uniqueness for the solution can be used to our benefit, because the doses can be administered in any order depending on various external factors such as the condition of the patient. The study concerning the derivation of the optimal number of doses N had already been performed for example in [8] in the hyperfractionated case, but (as far as we know) it is completely new for the hypofractionated case.
Let us emphasize again that the calculations to apply all these results are elementary and can be carried out using a pocket calculator from the initial data. Of course, the algorithms described in Tables 1 and 2 can be implemented quite straightforwardly in any platform using any programming language to make them more accessible.
We hope that these theoretical results may provide useful insights to address more complete models (including repopulation terms and multiple OAR) and that, ultimately, will lead to some improvement (however small) in clinical practice, due to the impact it would have on the large number of patients who could benefit.
Appendix 1: Proof of Theorem 8 The expression given in i) is derived exactly as in Theorem 6 for the values N ≤ λ 0 . Taking into account that d * can be very close to d max or d min , item ii) can be seen as a kind of special case of iii). So, we will focus on proving iii) that it is the most complicated case.
Together with (8) and the assumption ω δ < 0, this implies iii), because (43) means that the values of the objective function E T at the solutions are decreasing with N and therefore, the maximum value will be attained at N = λ 0 , the lowest value of N in the set (λ 0 , ρ 0 ] ∩ N. Comparing their expressions in the form (26) with N + 1 and N , resp., we conclude thatK ≥ K + 1. Hence, if we denote K 0 =K − K ∈ N, the inequality (43) can be written as Let us recall that d satisfies (24) andd * verifies We will show that (44) holds dividing the argumentation in three cases: that under the assumption satisfies the bounds restrictions and ).
Due to (45) and the hypothesis we have Again, we have a feasible point for the problem (27) and therefore we deduce d * ≤d * and hence (44), because Case 3.-Finally, suppose that K 0 d min + (1 − K 0 )d max > 0 and moreover Here, we introduce the auxiliary function defined for s ∈ [0, 1] by Solving the quadratic equations (24) and (45), it is easy to derive that Using the Mean Value Theorem, we deduce that there exists θ ∈ (0, 1) such that Therefore, the inequality (44) is equivalent to Under the present hypotheses, the function G is strictly increasing and, since θ is an unknown value in (0, 1), we will verify that (47) is valid if it holds for θ = 0. On the other hand, the value K 0 ∈ N is also unknown, but we can verify that the function Then, as asserted.
Appendix 2: Proof of Theorem 12 It follows the same lines to that of Theorem 7. Case a).-Assume ρ T ∈ N, ρ T ≥ 2.
As usual, we divide the interval for feasible values of N in two parts: [λ T , ρ T ) ∩ N and [ρ T , +∞) ∩ N.
In order to study the dependence with respect to N in the interval [λ T , ρ T ), thanks to Proposition 4 (with ω δ > 0) and (36), it is enough to consider the auxiliary function Once more, it follows easily that ψ 1 (N ) is an strictly increasing function.
Since we are assuming ρ T ∈ N and ρ T ≥ 2, then ψ 1 will take its maximum value in the set [λ T , ρ T ) ∩ N at N 1 = ρ T − 1. Therefore, the candidate for solution to problem (P 2 ) is given by the pair (N 1 , d N 1 ) with d N 1 = (d 1 , ..., d 1 ), where d 1 is given by On the other hand, in the interval [ρ T , +∞), we know that the other candidate for solution to problem (P 2 ) is given by the pair (N 2 , d N 2 ) with N 2 = ρ T and d N 2 = (d min , ..., d min ).
To derive that (N 2 , d N 2 ) is the unique solution to problem (P 2 ), we will show that Following the same idea to that of the proof of Theorem 7, we introduce the auxiliary function Taking into account (48) and that N 2 ϕ T (d min ) = γ T (by the definition of ρ T ), it can be checked that H 1 (x) = N 1 d 1 − N 2 d min < 0, since N 1 < N 2 .
In the interval [ρ T , +∞), the other candidate is N 2 = ρ T with d N 2 = (d min , ..., d min ). When ρ T ∈ N, any of them can provide the unique solution to problem (P 2 ) (see for instance Example 6).
Arguing as in Appendix 1, the candidate when N runs [λ T , ρ T ) ∩ N is N 1 = λ T with d N 1 given by Theorem 13-a) or b) and N = N 1 , thanks to Theorem 11−ii). We will conclude by showing that Let us argue with the expression b) for d N 1 , because (as we have pointed out before) the value d * can be very close to d min or d max and hence item a) can be seen as a special case of b). Therefore, the inequality (50) is equivalent to (51) Kϕ 0 (d min ) + ϕ 0 (d ) + (N 1 − K − 1)ϕ 0 (d max ) ≤ N 2 ϕ 0 (d min ).
If H 2 is an increasing function, since ω δ < 0, we will have which gives (51). So, taking into account that let us finish the proof by showing that H 2 (x) ≥ 0. If N 2 d min > N 1 d max , this is true straightforwardly, because we know that N 1 d max > Kd min + d * + (N 1 − K − 1)d max .
When N 2 d min ≤ N 1 d max , we can argue as in the proof of Theorem 8, taking the pointd = (d 1 , ..., d N 1 ) = ( N 2 N 1 d min , ..., N 2 N 1 d min ), that satisfies the bounds restrictions and This means that it is feasible for the problem (P N 1 ,− 2 ). Taking into account that (N 1 , d N 1 ) is a solution for that problem, see Proposition 4-ii) and Remark 3, we get H 2 (x) ≥ 0.