# Analytical solution to the radiotherapy fractionation problem including
dose bound constraints

Luis A. Fernández and Lucía Fernández  Dedicated to the memory of Juan Antonio
Fernández (1957-2018)  Dep. Mathematics, Statistics and Computation
University of Cantabria (SPAIN<EMAIL_ADDRESS>
###### Abstract.

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.

###### Key words and phrases:

Radiotherapy, fractionation, mixed and continuous optimization, linear
quadratic model.

###### 1991 Mathematics Subject Classification:

92C50, 90C20, 90C90

The work of the first author was supported by the Spanish “Ministerio de
Economía, Industria y Competitividad” under project MTM2017-83185-P

##  1\. 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

$S=\exp{\left(-\alpha_{T}d-\beta_{T}d^{2}\right)},$

where $\alpha_{T}$ and $\beta_{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\cdot 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^{N}=\exp{\left(-\sum_{i=1}^{N}\left(\alpha_{T}d_{i}+\beta_{T}d_{i}^{2}\right)\right)}.$

From here it is clear that the effect of radiation on the tumor is determined
by the quadratic function

(2)
$E_{T}(N,d)=\alpha_{T}\sum_{i=1}^{N}d_{i}+\beta_{T}\sum_{i=1}^{N}d_{i}^{2}.$

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 $\delta d$, with $\delta\in(0,1]$ being the so-called “sparing
factor”. The value of $\delta$ 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, $\delta$ will be about $0$; if not, it will be larger, until
reaching the value $\delta\approx 1$ at worst. Therefore, the effect of the
radiation on the OAR is determined by the following function

(3)
$E_{OAR}(N,d)=\alpha_{0}\delta\sum_{i=1}^{N}d_{i}+\beta_{0}\delta^{2}\sum_{i=1}^{N}d_{i}^{2},$

where $\alpha_{0}$ and $\beta_{0}$ are the parameters associated to the
healthy organs that we are trying to protect.

Typical values for $\alpha_{0},\beta_{0},\alpha_{T}$ and $\beta_{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}\leq d\leq
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.

##  2\. 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):

$(P_{1})\left\\{\begin{array}[]{ll}\begin{aligned} &\textrm{Maximize
}E_{T}(N,d),\\\ &\textrm{subject to }N\in\mathbb{N},d\in\mathbb{R}^{N}\textrm{
such that }\\\ &E_{OAR}(N,d)\leq\gamma_{OAR},\\\ &d_{min}\leq d_{i}\leq
d_{max},\ i=1,...,N,\end{aligned}\end{array}\right.$

where $E_{T}(N,d)$ is given by (2), $E_{OAR}(N,d)$ by (3) and
$d_{min},d_{max}$ and $\gamma_{OAR}$ are a priori known positive parameters,
that should be provided by the specialists. Roughly speaking, the restriction
$E_{OAR}(N,d)\leq\gamma_{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{(-\gamma_{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\in\mathbb{N}$, which corresponds to the
number of radiation doses, and $N$ continuous variables,
$d_{i}\in\mathbb{R},1\leq i\leq N$, which are the doses. In other words, this
problem has the peculiarity of having a variable number of unknowns.

Along this paper, it will be denoted

(4) $\varphi_{0}(r)=\alpha_{0}\delta r+\beta_{0}\delta^{2}r^{2},$ (5)
$\lambda_{0}=max\left\\{1,\dfrac{\gamma_{OAR}}{\varphi_{0}(d_{max})}\right\\},$

and

(6) $\rho_{0}=\dfrac{\gamma_{OAR}}{\varphi_{0}(d_{min})}.$

Also, we will denote by $\lfloor x\rfloor$ the greatest integer less than or
equal to $x$ and by $\lceil x\rceil$ the least integer greater than or equal
to $x$. Finally, the notation $d^{N}=(d_{0},\ldots,d_{0})$ means that
$d^{N}\in\mathbb{R}^{N}$ having all the $N$ components equal to $d_{0}$.

###  2.1. Existence of solution for $(P_{1})$

######  Theorem 1.

Let us assume $d_{min}>0$ and $\rho_{0}\geq 1$. Then, the problem $(P_{1})$
has (at least) one solution.

###### Proof.

Taking into account the restrictions for $(P_{1})$ and that $d_{min}>0$, we
have $N\leq\rho_{0}.$ Hence, the set of feasible values for $N$ is finite.

If $\rho_{0}=1$, the solution is $(N,d)=(1,d_{min})$, because this is the only
admissible pair for $(P_{1})$.

When $\rho_{0}\in(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 pair $(1,\min{\\{d_{max},\overline{d}_{0}\\}}),$ where
$\overline{d}_{0}=\dfrac{-\alpha_{0}+\sqrt{\alpha^{2}_{0}+4\beta_{0}\gamma_{OAR}}}{2\beta_{0}\delta}$.
Let us stress that $\varphi_{0}(\overline{d}_{0})=\gamma_{OAR}$.

If $\rho_{0}\geq 2$ we can reduce the problem $(P_{1})$ to a finite collection
of continuous optimization problems $(P_{1}^{N})$ with fixed $N$ given by:

$(P_{1}^{N})\left\\{\begin{array}[]{ll}\begin{aligned}
&\textrm{Maximize}\hskip
11.38092pt\tilde{E}^{N}_{T}(d)=\alpha_{T}\sum_{i=1}^{N}d_{i}+\beta_{T}\sum_{i=1}^{N}d_{i}^{2},\\\
&\textrm{subject to }d\in\mathbb{R}^{N}\textrm{ such that }\\\
&E_{OAR}(N,d)\leq\gamma_{OAR},\\\ &d_{min}\leq d_{i}\leq d_{max},\ \ \ \
i=1,...,N.\end{aligned}\end{array}\right.$

Firstly we will prove the existence of a solution for each problem
$(P_{1}^{N})$ (see Theorem 2), for $N$ running
$\left[1,\rho_{0}\right]\cap\mathbb{N}$ and denote it by $\overline{d}^{N}$.
Then, it is enough to take the pair
$\left(\overline{N},\overline{d}^{\overline{N}}\right)$ from the finite set

$\left\\{\left(N,\overline{d}^{N}\right):N\in\left[1,\rho_{0}\right]\cap\mathbb{N}\right\\},$

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_{1}^{N})$ is proved below:

######  Theorem 2.

Let us assume $d_{min}>0$, $\rho_{0}\geq 2$ and
$N\in\left[1,\rho_{0}\right]\cap\mathbb{N}$. Then the problem $(P_{1}^{N})$
has (at least) one solution.

###### Proof.

For small values of $N$, more precisely
$N\in\left[1,\lambda_{0}\right]\cap\mathbb{N}$, it is easy to verify that the
solution for $(P_{1}^{N})$ is the trivial one with maximum doses, that is,
$\overline{d}^{N}=(d_{max},...,d_{max}).$ For other values,
$N\in\left(\lambda_{0},\rho_{0}\right]\cap\mathbb{N}$, the existence of
solution for $(P_{1}^{N})$ 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_{1}^{N})$ 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\in\\{1,...,\overline{N}\\}$ such that
$\overline{d}_{i}\neq\overline{d}_{j}$ and interchange these coordinates to
generate a new solution.

  * c)

Under the condition $\rho_{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:

######  Example 1.

$(P_{10})\left\\{\begin{array}[]{ll}\textrm{Maximize }\displaystyle
E_{T}(N,d)=\sum_{i=1}^{N}d_{i}+\sum_{i=1}^{N}d_{i}^{2},\\\ \textrm{subject to
}N\in\mathbb{N},d_{i}\in\mathbb{R},\\\ \displaystyle
E_{OAR}(N,d)=\sum_{i=1}^{N}d_{i}+2\sum_{i=1}^{N}d_{i}^{2}\leq 10,\\\ 0\leq
d_{i}\leq 1,i=1,...,N.\end{array}\right.$

It is clear that for all feasible points we have

$E_{T}(N,d)\leq E_{OAR}(N,d)\leq 10.$

Let us stress that here $N$ can take any natural value, without restrictions.
Inspired by Theorem 13 below, let us consider the sequence given by

$d^{N}=(d_{0N},...,d_{0N}),\ \ \mbox{ with }\ \
d_{0N}=\dfrac{-1}{4}+\sqrt{\dfrac{1}{16}+\dfrac{5}{N}}.\ $

It is easy to check that it is feasible for $N\geq 4$,

$E_{OAR}(N,d^{N})=10,\ \ E_{T}(N,d^{N})=10-Nd_{0N}^{2}\longrightarrow 10,\ \
\mbox{ as }N\rightarrow+\infty.$

Hence, the problem $(P_{10})$ can not have solution
$(\overline{N},\overline{d})$: the supremum value $10$ can not be attained
since it should happen that

$10+\sum_{i=1}^{\overline{N}}\overline{d}_{i}^{2}=E_{T}(\overline{N},\overline{d})+\sum_{i=1}^{\overline{N}}\overline{d}_{i}^{2}=E_{OAR}(\overline{N},\overline{d})\leq
10,$

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$, $\rho_{0}\geq 2$ and
$N\in\left(\lambda_{0},\rho_{0}\right]\cap\mathbb{N}$. Then, the inequality
constraint of the problem $(P_{1}^{N})$ has to be active at
$\overline{d}^{N}$, with $\overline{d}^{N}$ being a solution for
$(P_{1}^{N})$.

###### Proof.

Arguing by contradiction, let us assume that the constraint is not active,
i.e.,

(7)
$\alpha_{0}\delta\sum_{i=1}^{N}\overline{d}_{i}+\beta_{0}\delta^{2}\sum_{i=1}^{N}\overline{d}_{i}^{2}<\gamma_{OAR}.$

Since $\lambda_{0}<N$, we know that there exists some index
$j\in\\{1,...,N\\}$ such that $\overline{d}_{j}<d_{max}$. Then, for
sufficiently small $\epsilon>0$, the point
$(\overline{d}_{1},...,\overline{d}_{j-1},\overline{d}_{j}+\epsilon,\overline{d}_{j+1},...,\overline{d}_{N})$
is feasible and satisfies

$\tilde{E}^{N}_{T}(\overline{d}^{N})<\tilde{E}^{N}_{T}((\overline{d}_{1},...,\overline{d}_{j-1},\overline{d}_{j}+\epsilon,\overline{d}_{j+1},...,\overline{d}_{N})),$

but this contradicts the fact that $\overline{d}^{N}$ is a solution for
$(P_{1}^{N})$. ∎

Hence, from now on, in this case we will consider the equality restriction

$\alpha_{0}\delta\sum_{i=1}^{N}d_{i}+\beta_{0}\delta^{2}\sum_{i=1}^{N}d_{i}^{2}=\gamma_{OAR}.$

Therefore, we deduce that

$\sum_{i=1}^{N}d_{i}^{2}=\dfrac{1}{\beta_{0}\delta^{2}}\left[\gamma_{OAR}-\alpha_{0}\delta\sum_{i=1}^{N}d_{i}\right],$

and the objective function can be written as

(8)
$\tilde{E}^{N}_{T}(d)=\left[\alpha_{T}-\dfrac{\beta_{T}\alpha_{0}}{\beta_{0}\delta}\right]\sum_{i=1}^{N}d_{i}+\dfrac{\beta_{T}\gamma_{OAR}}{\beta_{0}\delta^{2}}.$

Based on this identity, we can directly simplify the formulation of the
problem $(P_{1}^{N})$ as follows:

######  Proposition 1.

Let us assume $d_{min}>0$, $\rho_{0}\geq 2$,
$N\in\left(\lambda_{0},\rho_{0}\right]\cap\mathbb{N}$ and denote

(9)
$\displaystyle\omega_{\delta}=\frac{\alpha_{T}}{\beta_{T}}-\frac{\alpha_{0}}{\beta_{0}\delta}.$

  * i)

If $\omega_{\delta}>0$, then $(P_{1}^{N})$ is equivalent to

(10) $(P_{1}^{N,+})\left\\{\begin{array}[]{ll}\textrm{Maximize
}\displaystyle\sum_{i=1}^{N}d_{i},\\\ \textrm{subject to
}d\in\mathbb{K}_{1}^{N},\end{array}\right.$

where

(11)
$\mathbb{K}_{1}^{N}=\\{d\in\mathbb{R}^{N}:E_{OAR}(N,d)=\gamma_{OAR},d_{min}\leq
d_{i}\leq d_{max},1\leq i\leq N\\}.$

  * ii)

If $\omega_{\delta}<0$, then $(P_{1}^{N})$ is equivalent to

(12) $(P_{1}^{N,-})\left\\{\begin{array}[]{ll}\textrm{Minimize
}\displaystyle\sum_{i=1}^{N}d_{i},\\\ \textrm{subject to
}d\in\mathbb{K}_{1}^{N}.\end{array}\right.$

  * iii)

If $\omega_{\delta}=0,$ then every feasible point for $(P_{1}^{N})$ 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
$\alpha_{T}/\beta_{T}>\alpha_{0}/\beta_{0}$ and therefore, the case
$\omega_{\delta}>0$ is more frequent in clinical practice.

  * c)

As a consequence of Proposition 1, we can appreciate the great difference
between the cases $\omega_{\delta}>0$ and $\omega_{\delta}<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_{1}^{N})$

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:

$(P_{1}^{2,+})\left\\{\begin{array}[]{ll}\textrm{Maximize }\hskip
7.11317pt\displaystyle d_{1}+d_{2},\\\ \textrm{subject to
}(d_{1},d_{2})\in\mathbb{R}^{2},\\\ 2(d_{1}+d_{2})+d_{1}^{2}+d_{2}^{2}=12,\\\
1\leq d_{1},d_{2}\leq 3.\end{array}\right.\hskip
8.5359pt(P_{1}^{2,-})\left\\{\begin{array}[]{ll}\textrm{Minimize }\hskip
7.11317pt\displaystyle d_{1}+d_{2},\\\ \textrm{subject to
}(d_{1},d_{2})\in\mathbb{R}^{2},\\\ 2(d_{1}+d_{2})+d_{1}^{2}+d_{2}^{2}=12,\\\
1\leq d_{1},d_{2}\leq 3.\end{array}\right.$

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.

Figure 1.

Visually one can guess that the unique solution to $(P_{1}^{2,+})$ is located
on the diagonal (more precisely, it is given by
$(\overline{d}_{0},\overline{d}_{0})$ with $\overline{d}_{0}=\sqrt{7}-1$) and
there are two solutions of $(P_{1}^{2,-})$ lying on the boundary
(specifically, $(\overline{d}_{1},\overline{d}_{2})$ with
$\overline{d}_{1}=1$, $\overline{d}_{2}=\sqrt{10}-1$ and
$\overline{d}_{1}=\sqrt{10}-1$, $\overline{d}_{2}=1$).

####  2.2.1. Solving $(P_{1}^{N,+})$

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_{1}^{N,+})$ is a vector with equal coordinates:

######  Theorem 4.

Let us assume $d_{min}>0$, $\rho_{0}\geq 2$ and
$N\in\left(\lambda_{0},\rho_{0}\right]\cap\mathbb{N}.$ Then, the unique
solution to $(P_{1}^{N,+})$ has the form
$\overline{d}^{N}=(\overline{d}_{0},...,\overline{d}_{0})$ with

(13)
$\overline{d}_{0}=\dfrac{-\alpha_{0}N+\sqrt{(\alpha_{0}N)^{2}+4\beta_{0}N\gamma_{OAR}}}{2\beta_{0}\delta
N}.$

###### Proof.

By using the Cauchy-Schwarz inequality, we have

$\left(\sum_{i=1}^{N}d_{i}\right)^{2}\leq
N\left(\sum_{i=1}^{N}d_{i}^{2}\right).$

Therefore, for each feasible point it follows that

(14)
$\gamma_{OAR}=\alpha_{0}\delta\sum_{i=1}^{N}d_{i}+\beta_{0}\delta^{2}\sum_{i=1}^{N}d_{i}^{2}\geq\alpha_{0}\delta\sum_{i=1}^{N}d_{i}+\dfrac{\beta_{0}\delta^{2}}{N}\left(\displaystyle\sum_{i=1}^{N}d_{i}\right)^{2}.$

Defining $q(z)=\dfrac{\beta_{0}\delta^{2}}{N}z^{2}+\alpha_{0}\delta
z-\gamma_{OAR}$, previous inequality can be rewritten as

(15) $q\left(\sum_{i=1}^{N}d_{i}\right)\leq 0.$

Taking into account that the polynomial $q$ can be factorized in the form
$q(z)=\dfrac{\beta_{0}\delta^{2}}{N}(z-z_{1})(z-z_{2})$ with $z_{1}<0<z_{2}$,
we know that the relation (15) holds if and only if
$\sum_{i=1}^{N}d_{i}\in[0,z_{2}]$, because all the components $d_{i}$ have to
be positive.

Now it is clear that the maximum value is achieved when
$\sum_{i=1}^{N}\overline{d}_{i}=z_{2}$. Combining this fact with (14), we
deduce that

(16)
$\gamma_{OAR}=\alpha_{0}\delta\displaystyle\sum_{i=1}^{N}\overline{d}_{i}+\beta_{0}\delta^{2}\sum_{i=1}^{N}\overline{d}_{i}^{2}\geq\alpha_{0}\delta\sum_{i=1}^{N}\overline{d}_{i}+\dfrac{\beta_{0}\delta^{2}}{N}\left(\sum_{i=1}^{N}\overline{d}_{i}\right)^{2}=\gamma_{OAR}.$

Hence

(17)
$\left(\sum_{i=1}^{N}\overline{d}_{i}\right)^{2}=N\sum_{i=1}^{N}\overline{d}_{i}^{2}.$

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.,
$\overline{d}_{1}=...=\overline{d}_{N}.$ Therefore,
$\overline{d}^{N}=(\overline{d}_{0},...,\overline{d}_{0})$ with
$\overline{d}_{0}=\dfrac{z_{2}}{N}$ and (13) holds.

Let us emphasize that $\overline{d}_{0}$ satisfies
$d_{min}\leq\overline{d}_{0}\leq d_{max}$, thanks to the hypothesis
$N\in\left(\lambda_{0},\rho_{0}\right]$. ∎

####  2.2.2. Solving $(P_{1}^{N,-})$

Given $\overline{d}$ a solution of $(P_{1}^{N,-})$, 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
$\overline{\alpha}\in[0,+\infty),\overline{\lambda}\in\mathbb{R}$ and
$\\{\overline{\mu}_{i}\\}_{i=1}^{2N}\subset[0,+\infty)$ verifying

(18)
$\overline{\alpha}+|\overline{\lambda}|+\sum_{i=1}^{2N}\overline{\mu}_{i}>0,$
(19) $\overline{\alpha}\left(\begin{array}[]{c}1\\\ \vdots\\\
1\end{array}\right)+\overline{\lambda}\left(\begin{array}[]{c}\alpha_{0}\delta+2\beta_{0}\delta^{2}\overline{d}_{1}\\\
\vdots\\\
\alpha_{0}\delta+2\beta_{0}\delta^{2}\overline{d}_{N}\end{array}\right)+\left(\begin{array}[]{c}\overline{\mu}_{N+1}-\overline{\mu}_{1}\\\
\vdots\\\
\overline{\mu}_{2N}-\overline{\mu}_{N}\end{array}\right)=\left(\begin{array}[]{c}0\\\
\vdots\\\ 0\end{array}\right),$ (20)
$\overline{\mu}_{i}(d_{min}-\overline{d}_{i})=0,\ \
\overline{\mu}_{i+N}(\overline{d}_{i}-d_{max})=0,\ 1\leq i\leq N.$

Inspired by the 2D example, we will prove that $\overline{d}$ lies on the
boundary of $[d_{min},d_{max}]^{N}$. Let us argue by contradiction assuming
that $\overline{d}_{i}\in(d_{min},d_{max})$, for all $i\in\\{1,...,N\\}$.
Then, thanks to (20) we deduce that $\overline{\mu}_{i}=0,\forall
i\in\\{1,...,2N\\}.$ In this case, (19) reads:

(21) $\overline{\alpha}\left(\begin{array}[]{c}1\\\ \vdots\\\
1\end{array}\right)+\overline{\lambda}\left(\begin{array}[]{c}\alpha_{0}\delta+2\beta_{0}\delta^{2}\overline{d}_{1}\\\
\vdots\\\
\alpha_{0}\delta+2\beta_{0}\delta^{2}\overline{d}_{N}\end{array}\right)=\left(\begin{array}[]{c}0\\\
\vdots\\\ 0\end{array}\right).$

If $\overline{\lambda}=0$, the identity (21) implies that
$\overline{\alpha}=0,$ but this is not possible by (18). Therefore
$\overline{\lambda}\neq 0$ and from (21) we get

$\overline{d}_{1}=...=\overline{d}_{N}=\dfrac{-(\overline{\alpha}+\overline{\lambda}\alpha_{0}\delta)}{2\overline{\lambda}\beta_{0}\delta^{2}}.$

In other words, we arrive to the solution of problem $(P_{1}^{N,+})$,
contradicting our initial hypothesis about $\overline{d}$.

Consequently, there exists (at least) one index $j\in\\{1,...,N\\}$ such that
$\overline{d}_{j}\in\\{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_{1}^{N,-})$ by means for the following
auxiliary problem:

$(P_{1}^{N-1,-})\left\\{\begin{array}[]{ll}\begin{aligned}
&\textrm{Minimize}\hskip 14.22636pt\sum_{i=1}^{N-1}d_{i}+\overline{d}_{N},\\\
&\textrm{subject to }d\in\mathbb{R}^{N-1}\textrm{ such that }\\\
&E_{OAR}(N-1,d)=\gamma_{OAR}-\alpha_{0}\delta\overline{d}_{N}-\beta_{0}\delta^{2}\overline{d}_{N}^{2},\\\
&\hskip 2.84544ptd_{min}\leq d_{i}\leq d_{max},\hskip
11.38092pti=1,...,N-1.\end{aligned}\end{array}\right.$

######  Proposition 2.

Assume that $\overline{d}=(\overline{d}_{1},...,\overline{d}_{N})$ is a
solution of $(P_{1}^{N,-})$. Then $(\overline{d}_{1},...,\overline{d}_{N-1})$
is a solution of $(P_{1}^{N-1,-})$.

###### Proof.

Every feasible point $(d_{1},...,d_{N-1})$ for the problem $(P_{1}^{N-1,-})$
satisfies

$E_{OAR}(N-1,d)=\gamma_{OAR}-\alpha_{0}\delta\overline{d}_{N}-\beta_{0}\delta^{2}\overline{d}_{N}^{2}.$

This implies that $(d_{1},...,d_{N-1},\overline{d}_{N})$ is a feasible point
for $(P_{1}^{N,-}).$ Hence, using that $\overline{d}$ is a solution of
$(P_{1}^{N,-})$, we get

$\sum_{i=1}^{N-1}\overline{d}_{i}\leq\sum_{i=1}^{N-1}d_{i},$

which implies that $(\overline{d}_{1},...,\overline{d}_{N-1})$ is a solution
of $(P_{1}^{N-1,-})$. ∎

Arguing exactly in the same form as before with the problem $(P_{1}^{N-1,-})$,
we deduce that there must be an index $j\in\\{1,...,N-1\\}$ such that
$\overline{d}_{j}\in\\{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:

$(P_{1}^{1,-})\left\\{\begin{array}[]{ll}\begin{aligned}
&\textrm{Minimize}\hskip 14.22636ptd_{1}+\sum_{i=2}^{N}\overline{d}_{i},\\\
&\textrm{subject to }d_{1}\in\mathbb{R}\textrm{ such that }\\\
&\alpha_{0}\delta
d_{1}+\beta_{0}\delta^{2}d_{1}^{2}=\gamma_{OAR}-\alpha_{0}\delta\sum_{i=2}^{N}\overline{d}_{i}-\beta_{0}\delta^{2}\sum_{i=2}^{N}\overline{d}_{i}^{2},\\\
&\hskip 2.84544ptd_{min}\leq d_{1}\leq
d_{max}.\end{aligned}\end{array}\right.$

Clearly, it is enough to solve the quadratic equation to get the solution.

Summarizing previous results, given
$N\in\left(\lambda_{0},\rho_{0}\right]\cap\mathbb{N}$, the solution of
$(P_{1}^{N,-})$ has one of the following structures:

(22)
$\overline{d}^{N}=(\underbrace{d_{min},...,d_{min}}_{K},\underbrace{d_{max},...,d_{max}}_{N-K}),$

or

(23)
$\overline{d}^{N}=(\underbrace{d_{min},...,d_{min}}_{K},d^{*},\underbrace{d_{max},...,d_{max}}_{N-K-1}),$

with $d^{*}\in(d_{min},d_{max})$ being the unique positive root of the
quadratic equation

(24)
$\varphi_{0}(d^{*})=\gamma_{OAR}-K\varphi_{0}(d_{min})-(N-K-1)\varphi_{0}(d_{max}),$

with $\varphi_{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

(25)
$K=\dfrac{N\varphi_{0}(d_{max})-\gamma_{OAR}}{\varphi_{0}(d_{max})-\varphi_{0}(d_{min})}.$

Of course, this holds if and only if the right hand side is a natural number
or zero.

  * b)

In the case (23), since $\varphi_{0}$ is an strictly increasing function in
$[0,+\infty)$, we know that

$\varphi_{0}(d_{min})<\varphi_{0}(d^{*})<\varphi_{0}(d_{max}),$

and using (24) we get that

$K\in\left(\dfrac{N\varphi_{0}(d_{max})-\gamma_{OAR}}{\varphi_{0}(d_{max})-\varphi_{0}(d_{min})}-1,\dfrac{N\varphi_{0}(d_{max})-\gamma_{OAR}}{\varphi_{0}(d_{max})-\varphi_{0}(d_{min})}\right)\cap\mathbb{N},$

which means that

(26)
$K=\displaystyle\lfloor\dfrac{N\varphi_{0}(d_{max})-\gamma_{OAR}}{\varphi_{0}(d_{max})-\varphi_{0}(d_{min})}\rfloor.$

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:

######  Theorem 5.

Let us assume $d_{min}>0$, $\rho_{0}\geq 2$ and
$N\in\left(\lambda_{0},\rho_{0}\right]\cap\mathbb{N}.$ Then, a solution to
problem $(P_{1}^{N,-})$ is given by

  * a)

$\overline{d}^{N}=(\underbrace{d_{min},...,d_{min}}_{K},\underbrace{d_{max},...,d_{max}}_{N-K}),$
when $K$ defined by (25) belongs to $\mathbb{N}\cup\\{0\\};$ otherwise,

  * b)

$\overline{d}^{N}=(\underbrace{d_{min},...,d_{min}}_{K},d^{*},\underbrace{d_{max},...,d_{max}}_{N-K-1}),$
with $K$ defined by (26) and $d^{\star}$ satisfying (24).

######  Remark 3.

It is not difficult to show that a solution for $(P_{1}^{N,-})$ is also a
solution for the problem

(27) $\left\\{\begin{array}[]{ll}\textrm{Minimize
}\displaystyle\sum_{i=1}^{N}d_{i,}\\\ \textrm{subject to
}d\in\mathbb{R}^{N}\textrm{ such that }\\\ E_{OAR}(N,d)\geq\gamma_{OAR},\\\
d_{min}\leq d_{i}\leq d_{max}\textrm{, }1\leq i\leq N.\end{array}\right.$

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})$ will be the pair
$\left(\overline{N},\overline{d}^{\overline{N}}\right)$, where
$\overline{d}^{\overline{N}}$ denotes a solution of $(P_{1}^{\overline{N}})$
from the finite set

$\left\\{\left(N,\overline{d}^{N}\right):N\in\left[1,\rho_{0}\right]\cap\mathbb{N}\right\\},$

maximizing the value of $E_{T}(N,d)$. In fact, combining previous results, we
can avoid the calculation of most solutions for $(P_{1}^{N})$ 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
$\lfloor\lambda_{0}\rfloor=\lfloor\rho_{0}\rfloor$.

######  Theorem 6.

Let us assume $d_{min}>0,$ $\rho_{0}\geq 2$ and
$\lfloor\lambda_{0}\rfloor=\lfloor\rho_{0}\rfloor$. Then, the unique solution
to problem $(P_{1})$ is given by the pair
$(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}=\lfloor\lambda_{0}\rfloor$ and
$\overline{d}^{\overline{N}}=(d_{max},...,d_{max}).$

###### Proof.

In this case the set of feasible values for $N$ is
$\\{1,\ldots,N_{1}\\}\subset\mathbb{N}$ with
$N_{1}=\lfloor\rho_{0}\rfloor=\lfloor\lambda_{0}\rfloor.$ For those values of
$N$, the solution for $(P_{1}^{N})$ has the form
$\overline{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
$\lfloor\lambda_{0}\rfloor<\lfloor\rho_{0}\rfloor$. In the the trivial case
$\omega_{\delta}=0$, the function to be minimized and the one defining the
restriction are proportional. Therefore, we can be derive the following
result:

######  Proposition 3.

Let us assume $d_{min}>0,$ $\rho_{0}\geq 2,$
$\lfloor\lambda_{0}\rfloor<\lfloor\rho_{0}\rfloor$ and $\omega_{\delta}=0.$
Then any feasible pair $(\overline{N},d)$ with
$E_{OAR}(\overline{N},d)=\gamma_{OAR}$ is a solution to problem $(P_{1})$. In
particular, the pairs $(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}\in\\{\lceil\lambda_{0}\rceil,\ldots,\lfloor\rho_{0}\rfloor\\}$
and $\overline{d}^{\overline{N}}=(\overline{d}_{0},...,\overline{d}_{0}),$
where

(28)
$\overline{d}_{0}=\dfrac{-\alpha_{0}\overline{N}+\sqrt{(\alpha_{0}\overline{N})^{2}+4\beta_{0}\overline{N}\gamma_{OAR}}}{2\beta_{0}\delta\overline{N}},$

with $\overline{N}$ in the above set.

###### Proof.

Due to the hypothesis $\omega_{\delta}=0$, we deduce straightforwardly that
problem $(P_{1})$ is equivalent to

$(\tilde{P}_{1})\left\\{\begin{array}[]{ll}\begin{aligned} &\textrm{Maximize
}E_{OAR}(N,d),\\\ &\textrm{subject to
}N\in\mathbb{N},d\in\mathbb{R}^{N}\textrm{ such that }\\\
&E_{OAR}(N,d)\leq\gamma_{OAR},\\\ &d_{min}\leq d_{i}\leq d_{max},\
i=1,...,N,\end{aligned}\end{array}\right.$

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
$\overline{d}_{0}\in[d_{min},d_{max}]$ if and only if
$\overline{N}\in[\lambda_{0},\rho_{0}].$ ∎

######  Theorem 7.

Let us assume $d_{min}>0,$ $\rho_{0}\geq 2,$
$\lfloor\lambda_{0}\rfloor<\lfloor\rho_{0}\rfloor$ and $\omega_{\delta}>0.$
Then, the unique solution to problem $(P_{1})$ is given by the pair
$\left(\overline{N},\overline{d}^{\overline{N}}\right)$ with
$\overline{N}=\lfloor\rho_{0}\rfloor$ and
$\overline{d}^{\overline{N}}=(\overline{d}_{0},...,\overline{d}_{0})$, with
$\overline{d}_{0}$ given by (28).

###### Proof.

Here, the set of feasible values for $N$ is
$\\{1,\ldots,N_{2}\\}\subset\mathbb{N}$ with $N_{2}=\lfloor\rho_{0}\rfloor.$
Arguing as in previous theorem, among the small values (i. e.
$N\in\\{1,\ldots,N_{1}\\},$ with $N_{1}=\lfloor\lambda_{0}\rfloor,$ for
solving $(P_{1})$ we only retain $N=N_{1}$ and
$\overline{d}^{N_{1}}=(d_{max},...,d_{max})$. For the other values, i. e.
$N\in\\{N_{1}+1,\ldots,N_{2}\\},$ since $\omega_{\delta}>0$, the corresponding
solution for $(P_{1}^{N})$ is given by
$\overline{d}^{N}=(\overline{d}_{0},...,\overline{d}_{0})$ with
$\overline{d}_{0}$ defined in (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

$\psi(N)=N\overline{d}_{0}=\dfrac{-\alpha_{0}N+\sqrt{(\alpha_{0}N)^{2}+4\beta_{0}N\gamma_{OAR}}}{2\beta_{0}\delta}.$

Here, it follows easily that $\psi(N)$ is an strictly increasing function and
then, it will take its maximum value in $\\{N_{1}+1,\ldots,N_{2}\\}$ at
$N_{2}.$

Finally, we will derive that $(N_{2},\overline{d}^{N_{2}})$ is the unique
solution to problem $(P_{1})$ by showing that

(29) $E_{T}(N_{1},\overline{d}^{N_{1}})<E_{T}(N_{2},\overline{d}^{N_{2}}).$

To that end, let us consider the linear function

$H(x)=N_{2}(x\overline{d}_{0}+\overline{d}_{0}^{2})-N_{1}(xd_{max}+d_{max}^{2}),\
\ \ x\in[\dfrac{\alpha_{0}}{\beta_{0}\delta},+\infty),$

with

(30)
$\overline{d}_{0}=\dfrac{-\alpha_{0}N_{2}+\sqrt{(\alpha_{0}N_{2})^{2}+4\beta_{0}N_{2}\gamma_{OAR}}}{2\beta_{0}\delta
N_{2}}.$

Using that
$N_{1}\varphi_{0}(d_{max})\leq\gamma_{OAR}=N_{2}\varphi_{0}(\overline{d}_{0})$
by the admissibility, we get that
$H\left(\dfrac{\alpha_{0}}{\beta_{0}\delta}\right)\geq 0$. Also, taking into
account (30), it can also be checked that
$H^{\prime}(x)=N_{2}\overline{d}_{0}-N_{1}d_{max}>0,$ because $N_{1}<N_{2}$.
Then, from the assumption $\omega_{\delta}>0$ (see (9)), it follows that
$H\left(\dfrac{\alpha_{T}}{\beta_{T}}\right)>0$, which is equivalent to (29).
∎

In the case $\omega_{\delta}<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,$ $\rho_{0}\geq 2,$
$\lfloor\lambda_{0}\rfloor<\lfloor\rho_{0}\rfloor$ and $\omega_{\delta}<0$.
Then, a solution to problem $(P_{1})$ is given by one of the following pairs:

  * i)

$\left(\overline{N},\overline{d}^{\overline{N}}\right)$ with
$\overline{N}=\lfloor\lambda_{0}\rfloor$ and
$\overline{d}^{\overline{N}}=(d_{max},...,d_{max}),$

  * ii)

$\left(\overline{N},\overline{d}^{\overline{N}}\right)$ with
$\overline{N}=\lceil\lambda_{0}\rceil$ and
$\overline{d}^{\overline{N}}=(\underbrace{d_{min},...,d_{min}}_{K},\underbrace{d_{max},...,d_{max}}_{\overline{N}-K}),$
when $K$ defined by (25) with $N=\overline{N}$ belongs to
$\mathbb{N}\cup\\{0\\},$ or

  * iii)

$\left(\overline{N},\overline{d}^{\overline{N}}\right)$ with
$\overline{N}=\lceil\lambda_{0}\rceil$ and
$\overline{d}^{\overline{N}}=(\underbrace{d_{min},...,d_{min}}_{K},d^{*},\underbrace{d_{max},...,d_{max}}_{\overline{N}-K-1}),$
where $K$ is defined in (26) and $d^{\star}$ satisfies (24) with
$N=\overline{N}.$

######  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 $\omega_{\delta}>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 $\omega_{\delta}<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
$\dfrac{\alpha_{T}}{\beta_{T}}<\dfrac{\alpha_{0}}{\beta_{0}},$ because then
$\omega_{\delta}<0$ for all $\delta\in(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 $\dfrac{\alpha_{T}}{\beta_{T}}\approx 1.5\ Gy,$
while $\dfrac{\alpha_{0}}{\beta_{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 $\omega_{\delta}<0$) is very convenient
in the practice. Assuming that the other parameters are set, the condition
$\omega_{\delta}<0$ can always be achieved by taking $\delta$ 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 $\overline{N}=1$ is obtained when
$\omega_{\delta}<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.

Next, we illustrate the general process with a particular example:

######  Example 3.

Let us consider the following parameters taken from a typical clinical
situation: $\alpha_{T}=0.05\ Gy^{-1},\beta_{T}=0.005\ Gy^{-2},\alpha_{0}=0.04\
Gy^{-1},\beta_{0}=0.02\ Gy^{-2},$ see [10], together with $d_{min}=1\ Gy,$ and
$d_{max}=6\ Gy.$ Then, the problem reads

$(P_{11})\left\\{\begin{array}[]{ll}\textrm{Maximize
}E_{T}(N,d)=0.05\displaystyle\sum_{i=1}^{N}d_{i}+0.005\sum_{i=1}^{N}d_{i}^{2},\\\
\textrm{subject to }N\in\mathbb{N},d_{i}\in\mathbb{R},\\\ \displaystyle
E_{OAR}(N,d)=0.04\delta\sum_{i=1}^{N}d_{i}+0.02\delta^{2}\sum_{i=1}^{N}d_{i}^{2}\leq\gamma_{OAR},\\\
1\leq d_{i}\leq 6,i=1,...,N.\end{array}\right.$

  * i)

For $\delta=0.3$ and $\gamma_{OAR}=0.78$, we have $\lambda_{0}\approx 5.7$,
$\rho_{0}\approx 56.52$ and $\omega_{\delta}\approx 3.33>0$. Among the values,
$N\in\\{6,7,8,\dots,54,55,56\\}$, we have proved (see Theorem 7) that the
biggest one, $N=56$, and the solution of $(P_{11}^{56})$ (that here is the
hyperfractionated $\overline{d}^{56}=(1.008,...,1.008)$) provides the solution
of $(P_{11})$; in fact, $E_{T}(56,\overline{d}^{56})\approx 3.107$. We can
easily check that with the standard protocol $\tilde{d}_{S}^{25}=(2,...,2)$ we
get $E_{T}(25,\tilde{d}_{S}^{25})=3$, and therefore there is about $3.5\%$
gain in terms of effect on the tumor, while the efficiency regarding OAR is
the same
($E_{OAR}(56,\overline{d}^{56})=E_{OAR}(25,\tilde{d}_{S}^{25})=0.78$).

On the other hand, the hypofractionated radiotherapy given by
$\tilde{d}_{2}^{15}=(2.67,...,2.67)$ produces
$E_{T}(15,\tilde{d}_{2}^{15})\approx 2.54,$ although the damage on OAR is also
lower: $E_{OAR}(15,\tilde{d}_{2}^{15})\approx 0.67$. These last treatments are
mentioned in [12, pg. 16] in connection with breast cancer.

Of course, here we are only taking into account the mathematical point of
view. In clinical practice, other factors such as patient inconvenience and
additional cost may advise the use of fewer doses, if the difference in terms
of efficiency is considered small.

  * ii)

For $\delta=0.1$ and $\gamma_{OAR}=0.22$, we calculate $\lambda_{0}\approx
7.05$, $\rho_{0}\approx 52.38$ and $\omega_{\delta}=-10<0$. In this case, the
solution for $(P_{11})$ is given by $(\overline{N},\overline{d})$ with
$\overline{N}=8,\overline{d}=(1,d^{*},\underbrace{6,...,6}_{6})$ and
$d^{*}\approx 5.588\ Gy$ (see Theorem 8-iii)), having
$E_{T}(8,\overline{d})\approx 3.37$. Recall that this is a hypofractionated
type treatment. Just for comparison reasons, let us mention that the solution
of $(P_{11}^{7})$ is $\tilde{d}_{1}^{7}=(6,...,6)$ and the solution of
$(P_{11}^{9})$ is $\tilde{d}_{2}=(1,1,d_{2}^{*},\underbrace{6,...,6}_{6})$
with $d_{2}^{*}\approx 4.9\ Gy$ producing $E_{T}(7,\tilde{d}_{1}^{7})\approx
3.36$ and $E_{T}(9,\tilde{d}_{2})\approx 3.355$, that are smaller than
$E_{T}(8,\overline{d})$ as expected.

  * iii)

Let us emphasize that the difference between “few” and “many” doses is
relative to each particular problem and not an absolute classification. For
instance, in the problem $(P_{11})$ with $\delta=0.3$ and $\gamma_{OAR}=0.1$,
the solution is given by $(\overline{N},\overline{d})$ with
$\overline{N}=7,\overline{d}^{7}=(1.031,\ldots,1.031)$ which corresponds to
the hyperfractionated case (because $N\in\\{1,\ldots,7\\}$), although the
number of delivered doses is lower than in the previous hypofractionated
treatment, see $ii)$.

For $\omega_{\delta}<0$, in most practical situations the solution is the one
presented in Theorem 8-$iii)$, but the alternatives $i)$ and $ii)$ can also
appear as we show in the following example:

######  Example 4.

Let us take the following parameters: $\alpha_{T}=0.08\ Gy^{-1},$
$\beta_{T}=0.02\ Gy^{-2},$ $\alpha_{0}=0.01\ Gy^{-1},\beta_{0}=0.001\
Gy^{-2},\delta=1,$ $d_{min}=1\ Gy,$ and $d_{max}=6\ Gy.$ Here,
$\omega_{\delta}=-6<0$ and the problem under consideration is

$(P_{12})\left\\{\begin{array}[]{ll}\textrm{Maximize
}E_{T}(N,d)=0.08\displaystyle\sum_{i=1}^{N}d_{i}+0.02\sum_{i=1}^{N}d_{i}^{2},\\\
\textrm{subject to }N\in\mathbb{N},d_{i}\in\mathbb{R},\\\
0.01\displaystyle\sum_{i=1}^{N}d_{i}+0.001\sum_{i=1}^{N}d_{i}^{2}\leq\gamma_{OAR},\\\
1\leq d_{i}\leq 6,i=1,...,N.\end{array}\right.$

For $\gamma_{OAR}=0.961,$ we get that $\lambda_{0}\approx 10.01$ and
$\rho_{0}\approx 87.36$. For $N=\lfloor\lambda_{0}\rfloor=10$, the solution of
$(P_{12}^{10})$ is given by $\overline{d}^{10}=(6,...,6)$. While for
$N=\lceil\lambda_{0}\rceil=11$, the solution of $(P_{12}^{11})$ is given by
$\overline{d}^{11}=(1,d^{*},\underbrace{6,...,6}_{9})$, with $d^{*}\approx
5.53565\ Gy$. Since $E_{T}(10,\overline{d}^{10})=12>11.956\approx
E_{T}(11,\overline{d}^{11})$, the solution of $(P_{12})$ is given by
$\overline{N}=10$ and $\overline{d}^{10}$. This shows that sometimes the
option $i)$ is the valid one.

Finally, taking $\gamma_{OAR}=0.971,$ we have $\lambda_{0}\approx 10.11,$
$N=\lceil\lambda_{0}\rceil=11$ and
$K=\dfrac{N\varphi_{0}(d_{max})-\gamma_{OAR}}{\varphi_{0}(d_{max})-\varphi_{0}(d_{min})}=1.$
Then, the solution of $(P_{12}^{11})$ is given by
$\overline{d}^{11}=(1,\underbrace{6,...,6}_{10}).$ Now, the solution of
$(P_{12})$ is also given by $\overline{N}=11$ and $\overline{d}^{11},$ because
$E_{T}(11,\overline{d}^{11})=12.1$. So, in this case Theorem 8-$ii)$ holds.

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:

######  Corollary 1.

Let us assume $d_{min}>0,$ $\rho_{0}\geq 2$ and
$(\overline{N},\overline{d}^{\overline{N}})$ is a solution of $(P_{1})$. Then,
$\overline{N}$ is an increasing function of $\gamma_{OAR}$, decreasing with
respect to $\alpha_{0},\beta_{0}$ and $\delta$ and independent of $\alpha_{T}$
and $\beta_{T}$. Moreover,

  * a)

When $\omega_{\delta}>0,$ $\overline{N}$ is also decreasing with respect to
$d_{min}$ and independent of $d_{max}$.

  * b)

When $\omega_{\delta}<0,$ $\overline{N}$ is also decreasing with respect to
$d_{max}$ and independent of $d_{min}$.

###### Proof.

It is just a consequence of the expressions
$\overline{N}=\lfloor\rho_{0}\rfloor$ with $\rho_{0}$ given by (6), when
$\omega_{\delta}>0,$ and $\overline{N}=\lfloor\lambda_{0}\rfloor$ or
$\lceil\lambda_{0}\rceil$ with $\lambda_{0}$ given by (5), when
$\omega_{\delta}<0.$ ∎

In Table 1, we summarize the resolution of problem $(P_{1})$ in algorithm
form, for the reader’s convenience.

ALGORITHM FOR SOLVING $(P_{1})$
---
DATA: $\alpha_{T},\beta_{T},\alpha_{0},\beta_{0},d_{min},d_{max},\delta$ and
$\gamma_{OAR}$
(all positive, $d_{min}<d_{max}$ and $\delta\leq 1$)
CALCULATE:
$\omega_{\delta}=\dfrac{\alpha_{T}}{\beta_{T}}-\dfrac{\alpha_{0}}{\beta_{0}\delta},$
$\lambda_{0}=max\left\\{1,\dfrac{\gamma_{OAR}}{\varphi_{0}(d_{max})}\right\\}$
and
$\rho_{0}=\dfrac{\gamma_{OAR}}{\varphi_{0}(d_{min})},$ with
$\varphi_{0}(r)=\alpha_{0}\delta r+\beta_{0}\delta^{2}r^{2}.$
IF $\rho_{0}<1,$ $(P_{1})$ has NO SOLUTION.
IF $\rho_{0}=1$, the pair
$(\overline{N},\overline{d}^{\overline{N}})=(1,d_{min})$ is the UNIQUE
SOLUTION of $(P_{1})$.
IF $\rho_{0}\in(1,2)$, the UNIQUE SOLUTION of $(P_{1})$ is the pair
$(\overline{N},\overline{d}^{\overline{N}})=(1,\min{\\{d_{max},\overline{d}_{0}\\}}),$
with
$\overline{d}_{0}=\dfrac{-\alpha_{0}+\sqrt{\alpha^{2}_{0}+4\beta_{0}\gamma_{OAR}}}{2\beta_{0}\delta}$.
IF $\rho_{0}\geq 2$ and $\lfloor\lambda_{0}\rfloor=\lfloor\rho_{0}\rfloor$,
the UNIQUE SOLUTION of $(P_{1})$
is the pair $(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}=\lfloor\lambda_{0}\rfloor$ and
$\overline{d}^{\overline{N}}=(d_{max},...,d_{max}).$
IF $\rho_{0}\geq 2,$ $\lfloor\lambda_{0}\rfloor<\lfloor\rho_{0}\rfloor$ and
$\omega_{\delta}>0,$ the UNIQUE SOLUTION of $(P_{1})$
is the pair $(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}=\lfloor\rho_{0}\rfloor$ and
$\overline{d}^{\overline{N}}=(\overline{d}_{0},...,\overline{d}_{0}),$
where
$\overline{d}_{0}=\dfrac{-\alpha_{0}\overline{N}+\sqrt{(\alpha_{0}\overline{N})^{2}+4\beta_{0}\overline{N}\gamma_{OAR}}}{2\beta_{0}\delta\overline{N}}$.
IF $\rho_{0}\geq 2,$ $\lfloor\lambda_{0}\rfloor<\lfloor\rho_{0}\rfloor$ and
$\omega_{\delta}<0,$ take $\overline{N}_{1}=\lceil\lambda_{0}\rceil$ and
CALCULATE:
$M=\dfrac{\overline{N}_{1}\varphi_{0}(d_{max})-\gamma_{OAR}}{\varphi_{0}(d_{max})-\varphi_{0}(d_{min})}$.
IF $M\in\mathbb{N}\cup\\{0\\}$, take $K=M$ and
$\overline{d}_{1}^{\overline{N}_{1}}=(\underbrace{d_{min},...,d_{min}}_{K},\underbrace{d_{max},...,d_{max}}_{\overline{N}_{1}-K})$.
IF $M\not\in\mathbb{N}\cup\\{0\\}$, take $K=\lfloor M\rfloor$ and
$\overline{d}_{1}^{\overline{N}_{1}}=(\underbrace{d_{min},...,d_{min}}_{K},d^{*},\underbrace{d_{max},...,d_{max}}_{\overline{N}_{1}-K-1}),$
with $d^{*}>0$ and
$\varphi_{0}(d^{*})=\gamma_{OAR}-K\varphi_{0}(d_{min})-(\overline{N}_{1}-K-1)\varphi_{0}(d_{max})$.
Also take $\overline{N}_{2}=\lfloor\lambda_{0}\rfloor$ and
$\overline{d}_{2}^{\overline{N}_{2}}=(d_{max},...,d_{max}).$
CALCULATE: $E_{T}(\overline{N}_{1},\overline{d}_{1}^{\overline{N}_{1}})$ and
$E_{T}(\overline{N}_{2},\overline{d}_{2}^{\overline{N}_{2}}).$
A SOLUTION of $(P_{1})$ is the pair
$(\overline{N},\overline{d}^{\overline{N}})$ that maximizes $E_{T}$ between
them.
IF $\rho_{0}\geq 2,$ $\lfloor\lambda_{0}\rfloor<\lfloor\rho_{0}\rfloor$ and
$\omega_{\delta}=0,$ ANY FEASIBLE PAIR $(\overline{N},d)$ such that
$E_{OAR}(\overline{N},d)=\gamma_{OAR}$ is a SOLUTION for $(P_{1})$.
In particular, the pairs $(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}\in\\{\lceil\lambda_{0}\rceil,\ldots,\lfloor\rho_{0}\rfloor\\}$
and $\overline{d}^{\overline{N}}=(\overline{d}_{0},...,\overline{d}_{0}),$
where
$\overline{d}_{0}=\dfrac{-\alpha_{0}\overline{N}+\sqrt{(\alpha_{0}\overline{N})^{2}+4\beta_{0}\overline{N}\gamma_{OAR}}}{2\beta_{0}\delta\overline{N}}.$
Table 1. Complete solution for problem $(P_{1})$ in algorithmic form.

##  3\. 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:

$(P_{2})\left\\{\begin{array}[]{ll}\begin{aligned} &\textrm{Minimize
}E_{OAR}(N,d),\\\ &\textrm{subject to }N\in\mathbb{N},\\\
&d\in\mathbb{R}^{N}\textrm{ such that }\\\ &E_{T}(N,d)\geq\gamma_{T},\\\
&d_{min}\leq d_{i}\leq d_{max}\textrm{,
}i=1,...,N,\end{aligned}\end{array}\right.$

where $E_{OAR}(N,d)$ is given by (3), $E_{T}(N,d)$ is defined in (2) and
$\gamma_{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\in\mathbb{N}$, and the value of the $N$ doses, $d_{i}\in\mathbb{R},1\leq
i\leq N$.

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 $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

(31) $\varphi_{T}(r)=\alpha_{T}r+\beta_{T}r^{2},$ (32)
$\lambda_{T}=max\left\\{1,\dfrac{\gamma_{T}}{\varphi_{T}(d_{max})}\right\\},$

and

(33)
$\rho_{T}=max\left\\{1,\dfrac{\gamma_{T}}{\varphi_{T}(d_{min})}\right\\}.$

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\in[\lambda_{T},+\infty)\cap\mathbb{N}$. We will begin by showing that for
each fixed feasible value $N$, the associated problem $(P_{2}^{N})$ has a
solution, where

$(P_{2}^{N})\left\\{\begin{array}[]{ll}\begin{aligned}
&\textrm{Minimize}\hskip
11.38092pt\tilde{E}^{N}_{OAR}(d)=\alpha_{0}\delta\sum_{i=1}^{N}d_{i}+\beta_{0}\delta^{2}\sum_{i=1}^{N}d_{i}^{2},\\\
&\textrm{subject to }d\in\mathbb{R}^{N}\textrm{ such that }\\\
&E_{T}(N,d)\geq\gamma_{T},\\\ &d_{min}\leq d_{i}\leq d_{max}\textrm{,
}i=1,...,N.\end{aligned}\end{array}\right.$

For large values of $N$, specifically for $N\geq\rho_{T}$, the solution of
$(P_{2}^{N})$ is the trivial one with minimum values $d_{min}$. Among them
only the smallest value of $N$ could have practical interest, i.e.,
$\lceil\rho_{T}\rceil$. For the other values, when they exist, that is for
$N\in\left[\lambda_{T},\rho_{T}\right)\cap\mathbb{N}$, the existence of
solution for $(P_{2}^{N})$ 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_{2}^{N})$ that we will denote $\overline{d}^{N}$.
Again, it is enough to take the pair
$\left(\overline{N},\overline{d}^{\overline{N}}\right)$ from the finite set

$\left\\{\left(N,\overline{d}^{N}\right):N\in\left[\lambda_{T},\lceil\rho_{T}\rceil\right]\cap\mathbb{N}\right\\},$

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 $\overline{d}^{\overline{N}}$ are
equal, the solution will not be unique, because two different coordinates can
be permuted to generate a different solution.

  * b)

When $\omega_{\delta}>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})\left\\{\begin{array}[]{ll}\textrm{Minimize }\displaystyle
E_{OAR}(N,d)=\sum_{i=1}^{N}d_{i}+\sum_{i=1}^{N}d_{i}^{2},\\\ \textrm{subject
to }N\in\mathbb{N},d_{i}\in\mathbb{R},\\\ \displaystyle
E_{T}(N,d)=2\sum_{i=1}^{N}d_{i}+\sum_{i=1}^{N}d_{i}^{2}\geq 10,\\\ 0\leq
d_{i}\leq 1,i=1,...,N.\end{array}\right.$

Let us consider the sequence given by

$d^{N}=(d_{0N},...,d_{0N}),\ \ \mbox{ with }\ \
d_{0N}=-1+\sqrt{1+\dfrac{10}{N}}.\ $

It is easy to check that it is feasible for $N\geq 4$,

$E_{T}(N,d^{N})=10,\ \ E_{OAR}(N,d^{N})=10-Nd_{0N}\longrightarrow 5,\ \ \mbox{
as }N\rightarrow+\infty.$

Now, we can deduce that $(P_{20})$ can not have solution
$(\overline{N},\overline{d})$ because if
$E_{OAR}(\overline{N},\overline{d})=5$ together with the restriction
$E_{T}(\overline{N},\overline{d})\geq 10$ we get
$\sum_{i=1}^{\overline{N}}\overline{d}_{i}\geq 5,$ but this is incompatible
with $E_{OAR}(\overline{N},\overline{d})=5.$

  * c)

In contrast, when $\omega_{\delta}<0$, it is easy to show that problem
$(P_{2}),$ with only the lower bound constraints $d_{i}\geq 0$, has as
solution $(\overline{N},\overline{d})$ with $\overline{N}=1$ and
$\overline{d}=\dfrac{-\alpha_{T}+\sqrt{(\alpha_{T})^{2}+4\beta_{T}\gamma_{T}}}{2\beta_{T}},$
see [10].

  * d)

There are some particular cases in which the solution of $(P_{2})$ can be
determined from previous argumentations very easily. For instance, when
$\rho_{T}=1$, because then
$\left(\overline{N},\overline{d}^{\overline{N}}\right)=(1,d_{min})$ is the
only admissible pair. Also when $\rho_{T}>1$ and
$\lceil\lambda_{T}\rceil=\lfloor\rho_{T}\rfloor$, because only the large
values for $N$ are feasible (i.e., those verifying $N\geq\rho_{T}$) and
consequently $(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}=\lceil\rho_{T}\rceil$ and
$\overline{d}^{\overline{N}}=(d_{min},...,d_{min})$ is the solution of
$(P_{2})$.

When $N\in[\lambda_{T},\rho_{T})\cap\mathbb{N}$, we know that
$d^{N}=(d_{min},...,d_{min})$ is not a solution of $(P_{2}^{N})$, because it
is not even feasible. Hence, we can simplify the problem $(P_{2}^{N})$ arguing
in a similar way as in the proof of Theorem 3.

######  Theorem 10.

Let us assume $d_{min}>0$ and
$N\in\left[\lambda_{T},\rho_{T}\right)\cap\mathbb{N}$. Then, the inequality
constraint of the problem $(P_{2}^{N})$ has to be active at
$\overline{d}^{N}$, being $\overline{d}^{N}$ a solution of $(P_{2}^{N})$.

From now on, the restriction will be taken as one of equality, that is,

$\alpha_{T}\sum_{i=1}^{N}d_{i}+\beta_{T}\sum_{i=1}^{N}d_{i}^{2}=\gamma_{T}.$

Again, applying the same procedure as for $(P_{1})$ in Section 2 we have

$\sum_{i=1}^{N}d_{i}^{2}=\dfrac{1}{\beta_{T}}\left[\gamma_{T}-\alpha_{T}\sum_{i=1}^{N}d_{i}\right],$

and the objective function will read

(34)
$\tilde{E}^{N}_{OAR}(d)=\left[\alpha_{0}-\dfrac{\beta_{0}\alpha_{T}\delta}{\beta_{T}}\right]\delta\sum_{i=1}^{N}d_{i}+\dfrac{\beta_{0}\delta^{2}\gamma_{T}}{\beta_{T}}.$

Now, it is clear that we can simplify the formulation of the problem
$(P_{2}^{N})$, as follows:

######  Proposition 4.

Let us assume $d_{min}>0$ and
$N\in\left[\lambda_{T},\rho_{T}\right)\cap\mathbb{N}$.

  * i)

If $\omega_{\delta}>0$, then $(P_{2}^{N})$ is equivalent to

$(P_{2}^{N,+})\left\\{\begin{array}[]{ll}\textrm{Maximize
}\displaystyle\sum_{i=1}^{N}d_{i},\\\ \textrm{subject to
}d\in\mathbb{K}_{2}^{N},\end{array}\right.$

where

(35)
$\mathbb{K}_{2}^{N}=\\{d\in\mathbb{R}^{N}:E_{T}(N,d)=\gamma_{T},d_{min}\leq
d_{i}\leq d_{max},1\leq i\leq N\\}.$

  * ii)

If $\omega_{\delta}<0$, then $(P_{2}^{N})$ is equivalent to

$(P_{2}^{N,-})\left\\{\begin{array}[]{ll}\textrm{Minimize
}\displaystyle\sum_{i=1}^{N}d_{i},\\\ \textrm{subject to
}d\in\mathbb{K}_{2}^{N}.\end{array}\right.$

  * iii)

If $\omega_{\delta}=0,$ then every feasible point for $(P_{2}^{N})$ is a
solution.

###### Proof.

It is enough to take into account that
$\alpha_{0}-\dfrac{\beta_{0}\alpha_{T}\delta}{\beta_{T}}=-\beta_{0}\delta\omega_{\delta},$
where $\omega_{\delta}$ is defined in (9). ∎

Once we have seen that $(P_{1}^{N,+})$ and $(P_{2}^{N,+})$ are essentially the
same problem (resp. $(P_{1}^{N,-})$ and $(P_{2}^{N,-})$), we can “translate”
the results obtained in section 2.2 to the current context as follows:

######  Theorem 11.

Let us assume $d_{min}>0$ and
$N\in\left[\lambda_{T},\rho_{T}\right)\cap\mathbb{N}$. Then,

  * i)

the unique solution to $(P_{2}^{N,+})$ is given by
$\overline{d}^{N}=(\overline{d}_{1},...,\overline{d}_{1})$ with

(36)
$\overline{d}_{1}=\dfrac{-\alpha_{T}N+\sqrt{(\alpha_{T}N)^{2}+4\beta_{T}N\gamma_{T}}}{2\beta_{T}N}.$

  * ii)

a solution for $(P_{2}^{N,-})$ has one of the following forms:

(37)
$\overline{d}^{N}=(\underbrace{d_{min},...,d_{min}}_{K},\underbrace{d_{max},...,d_{max}}_{N-K}),$

with

(38)
$K=\dfrac{N\varphi_{T}(d_{max})-\gamma_{T}}{\varphi_{T}(d_{max})-\varphi_{T}(d_{min})}\in\mathbb{N}\cup\\{0\\},\
\ \mbox{ or }$ (39)
$\overline{d}^{N}=(\underbrace{d_{min},...,d_{min}}_{K},d^{*},\underbrace{d_{max},...,d_{max}}_{N-K-1}),$

with

(40)
$K=\lfloor\dfrac{N\varphi_{T}(d_{max})-\gamma_{T}}{\varphi_{T}(d_{max})-\varphi_{T}(d_{min})}\rfloor,$

and $d^{*}\in(d_{min},d_{max})$ satisfying

(41)
$\varphi_{T}(d^{*})=\gamma_{T}-K\varphi_{T}(d_{min})-(N-K-1)\varphi_{T}(d_{max}).$

###  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
$\omega_{\delta}=\frac{\alpha_{T}}{\beta_{T}}-\frac{\alpha_{0}}{\beta_{0}\delta}$
(see (9)), we will see that $\rho_{T}$ and the sign of $\omega_{\delta}$ are
the determinant factors in this analysis.

Again, the case $\omega_{\delta}=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,$ $\rho_{T}>1,$
$\lceil\lambda_{T}\rceil<\lfloor\rho_{T}\rfloor$ and $\omega_{\delta}=0.$ Then
any feasible pair $(\overline{N},d)$ with $E_{T}(\overline{N},d)=\gamma_{T}$
is a solution to problem $(P_{2})$. In particular, the pairs
$(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}\in\\{\lceil\lambda_{T}\rceil,\ldots,\lfloor\rho_{T}\rfloor\\}$
and $\overline{d}^{\overline{N}}=(\overline{d}_{1},...,\overline{d}_{1}),$
where

(42)
$\overline{d}_{1}=\dfrac{-\alpha_{T}\overline{N}+\sqrt{(\alpha_{T}\overline{N})^{2}+4\beta_{T}\overline{N}\gamma_{T}}}{2\beta_{T}\overline{N}}.$

Let us now continue by studying the more frequent case $\omega_{\delta}>0$.
Here, we have to distinguish two different situations, depending on
$\rho_{T}\in\mathbb{N}$ or not:

######  Theorem 12.

Let us assume that $d_{min}>0,$
$\lceil\lambda_{T}\rceil<\lfloor\rho_{T}\rfloor$ and $\omega_{\delta}>0$.

  * a)

If $\rho_{T}\in\mathbb{N}$, $\rho_{T}\geq 2,$ then the unique solution to
problem $(P_{2})$ is given by $\overline{N}=\rho_{T}$ and
$\overline{d}^{\overline{N}}=(d_{min},...,d_{min}),$

  * b)

If $\rho_{T}\not\in\mathbb{N},$ then the unique solution to problem $(P_{2})$
is given by $\left(\overline{N},\overline{d}^{\overline{N}}\right),$ where:

    * i)

$\overline{N}=\lceil\rho_{T}\rceil$ and
$\overline{d}^{\overline{N}}=(d_{min},...,d_{min}),$ or

    * ii)

$\overline{N}=\lfloor\rho_{T}\rfloor$ and
$\overline{d}^{\overline{N}}=(\overline{d}_{1},...,\overline{d}_{1})$, with
$\overline{d}_{1}$ given by (42).

The proof is similar to that of Theorem 7 and postponed to Appendix $2$.

Finally, the case $\omega_{\delta}<0$ is studied in the next theorem and its
proof is detailed in Appendix $3$.

######  Theorem 13.

Let us assume $d_{min}>0,$ $\rho_{T}>1,$
$\lceil\lambda_{T}\rceil<\lfloor\rho_{T}\rfloor$ and $\omega_{\delta}<0$.
Then, a solution to problem $(P_{2})$ is given by
$\left(\overline{N},\overline{d}^{\overline{N}}\right),$ with
$\overline{N}=\lceil\lambda_{T}\rceil$ and

  * a)

$\overline{d}^{\overline{N}}=(\underbrace{d_{min},...,d_{min}}_{K},\underbrace{d_{max},...,d_{max}}_{\overline{N}-K}),$
when $K$ defined by (38), with $N=\overline{N},$ belongs to
$\mathbb{N}\cup\\{0\\};$ otherwise,

  * b)

$\overline{d}^{\overline{N}}=(\underbrace{d_{min},...,d_{min}}_{K},d^{*},\underbrace{d_{max},...,d_{max}}_{\overline{N}-K-1}),$
with $K$ defined by (40) and $d^{\star}$ satisfies (41), both with
$N=\overline{N}.$

######  Remark 6.

  * a)

Once more, let us emphasize that when $\omega_{\delta}>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 $\omega_{\delta}<0$ the
optimal value is the smallest one (i.e. it is a hypofractionated type
treatment). This classification in terms of $\omega_{\delta}$ 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
$\omega_{\delta}\leq 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 $\omega_{\delta}<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:
$\alpha_{0}=0.04\ Gy^{-1}$, $\beta_{0}=0.02\ Gy^{-2},$ $\alpha_{T}=0.05\
Gy^{-1},$ $\beta_{T}=0.005\ Gy^{-2},$ $d_{min}=1\ Gy$ and $d_{max}=6\ Gy.$

  * i)

If $\delta=1$, then $\omega_{\delta}=8>0$. Hence, we are considering the
problem

$(P_{21})\left\\{\begin{array}[]{ll}\textrm{Minimize
}E_{OAR}(N,d)=0.04\displaystyle\sum_{i=1}^{N}d_{i}+0.02\sum_{i=1}^{N}d_{i}^{2},\\\
\textrm{subject to }N\in\mathbb{N},d_{i}\in\mathbb{R}\ \textrm{ such that }\\\
E_{T}(N,d)=0.05\displaystyle\sum_{i=1}^{N}d_{i}+0.005\sum_{i=1}^{N}d_{i}^{2}\geq\gamma_{T},\\\
1\leq d_{i}\leq 6,i=1,...,N.\end{array}\right.$

When $\gamma_{T}=4$, we can easily calculate that $\lambda_{T}\approx 8.33$
and $\rho_{T}\approx 72.72$. Applying Theorem 12-b) the unique solution for
$(P_{21})$ is one of these two pairs: $N_{1}=\lceil\rho_{T}\rceil=73$,
$\overline{d}^{N_{1}}=(1,...,1)$ with
$E_{OAR}(N_{1},\overline{d}^{N_{1}})=4.38$ or
$N_{2}=\lfloor\rho_{T}\rfloor=72$,
$\overline{d}^{N_{2}}=(\overline{d}_{1},...,\overline{d}_{1})$ with
$\overline{d}_{1}\approx 1.00926\ Gy$ and the objective function value
$E_{OAR}(N_{2},\overline{d}^{N_{2}})\approx 4.373.$ Clearly, we choose the
second pair.

  * ii)

Previous case is the most frequent in practice when $\omega_{\delta}>0$, but
for some specific values of $\gamma_{T}$ the alternative indicated in Theorem
12-b) occurs. For instance, taking $\gamma_{T}=4.014$ in $(P_{21})$, the
values become $\lambda_{T}\approx 8.36$ and $\rho_{T}\approx 72.98$. Now, the
two candidates for being the unique solution for $(P_{21})$ are:
$N_{1}=\lceil\rho_{T}\rceil=73$, $d^{N_{1}}=(1,...,1)$ with
$E_{OAR}(N_{1},\overline{d}^{N_{1}})=4.38$, as before, and
$N_{2}=\lfloor\rho_{T}\rfloor=72$,
$\overline{d}^{N_{2}}=(\overline{d}_{1},...,\overline{d}_{1})$ with
$\overline{d}_{1}\approx 1.0125\ Gy$ and $E_{OAR}(N_{2},d^{N_{2}})\approx
4.39.$ It is apparent that here the solution is the first one.

  * iii)

If $\delta=0.1$, then $\omega_{\delta}=-10<0$ and the problem under study is

$(P_{22})\left\\{\begin{array}[]{ll}\textrm{Minimize
}E_{OAR}(N,d)=0.004\displaystyle\sum_{i=1}^{N}d_{i}+0.0002\sum_{i=1}^{N}d_{i}^{2},\\\
\textrm{subject to }N\in\mathbb{N},d_{i}\in\mathbb{R}\ \textrm{ such that }\\\
E_{T}(N,d)=0.05\displaystyle\sum_{i=1}^{N}d_{i}+0.005\sum_{i=1}^{N}d_{i}^{2}\geq\gamma_{T},\\\
1\leq d_{i}\leq 6,i=1,...,N.\end{array}\right.$

For $\gamma_{T}=4.35$, we have that $\lambda_{T}\approx 9.06$ and
$\rho_{T}\approx 79.09$. Therefore, applying Theorem 13-b), we know that a
solution for $(P_{22})$ is given by $(\overline{N},\overline{d})$ with
$\overline{N}=\lceil\lambda_{T}\rceil=10$ and
$\overline{d}=(1,d^{*},\underbrace{6,...,6}_{8})$ with $d^{*}\approx 5.77\
Gy$. In this case, $E_{OAR}(10,\overline{d})\approx 0.2835$. As in Example 3,
for comparison reasons, we can calculate that the solution of solution of
$(P_{22}^{11})$ is $\tilde{d}_{1}=(1,1,d_{1}^{*},\underbrace{6,...,6}_{8})$
with $d_{1}^{*}\approx 5.247\ Gy$ producing $E_{OAR}(11,\tilde{d}_{1})\approx
0.2845$, that is slightly bigger than $E_{OAR}(10,\overline{d})$, as expected.

  * iv)

Again, it can be shown that, for some specific values of the parameter
$\gamma_{T}$, the alternative exposed in Theorem 13-a) is true. In particular,
for $\gamma_{T}=4.375$, we have $\lambda_{T}\approx 9.11$ and $\rho_{T}\approx
79.54$ and a solution for $(P_{22})$ is given by $(\overline{N},\overline{d})$
with $\overline{N}=\lceil\lambda_{T}\rceil=10$ and
$\overline{d}=(1,\underbrace{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})$.

######  Corollary 2.

Let us assume $d_{min}>0,$ $\rho_{0}\geq 2$ and
$(\overline{N},\overline{d}^{\overline{N}})$ is a solution of $(P_{2})$. Then,
$\overline{N}$ is an increasing function of $\gamma_{T}$, decreasing with
respect to $\alpha_{T}$ and $\beta_{T}$ and independent of
$\alpha_{0},\beta_{0}$ and $\delta$. Moreover,

  * a)

When $\omega_{\delta}>0,$ $\overline{N}$ is also decreasing with respect to
$d_{min}$ and independent of $d_{max}$.

  * b)

When $\omega_{\delta}<0,$ $\overline{N}$ is also decreasing with respect to
$d_{max}$ and independent of $d_{min}$.

###### Proof.

It follows from the expressions $\overline{N}=\lfloor\rho_{T}\rfloor$ or
$\lceil\rho_{T}\rceil$ with $\rho_{T}$ given by (33), when
$\omega_{\delta}>0,$ and $\overline{N}=\lceil\lambda_{T}\rceil$ with
$\lambda_{T}$ given by (32), when $\omega_{\delta}<0.$ ∎

For the reader’s convenience, we have summarized the complete algorithm for
the resolution of the problem $(P_{2})$ in Table 2.

ALGORITHM FOR SOLVING $(P_{2})$
---
DATA: $\alpha_{T},\beta_{T},\alpha_{0},\beta_{0},d_{min},d_{max},\delta$ and
$\gamma_{T}$
(all positive, $d_{min}<d_{max}$ and $\delta\leq 1$)
CALCULATE:
$\omega_{\delta}=\dfrac{\alpha_{T}}{\beta_{T}}-\dfrac{\alpha_{0}}{\beta_{0}\delta},$
$\lambda_{T}=max\left\\{1,\dfrac{\gamma_{T}}{\varphi_{T}(d_{max})}\right\\}$
and
$\rho_{T}=max\left\\{1,\dfrac{\gamma_{T}}{\varphi_{T}(d_{min})}\right\\},$
with $\varphi_{T}(r)=\alpha_{T}r+\beta_{T}r^{2}.$
IF $\rho_{T}=1$, the pair
$(\overline{N},\overline{d}^{\overline{N}})=(1,d_{min})$ is the UNIQUE
SOLUTION of $(P_{2})$.
IF $\rho_{T}>1$ and $\lceil\lambda_{T}\rceil=\lfloor\rho_{T}\rfloor$, the
UNIQUE SOLUTION of $(P_{2})$
is the pair $(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}=\lceil\rho_{T}\rceil$ and
$\overline{d}^{\overline{N}}=(d_{min},...,d_{min}).$
IF $\rho_{T}\in\mathbb{N}$, $\rho_{T}\geq 2$,
$\lceil\lambda_{T}\rceil<\lfloor\rho_{T}\rfloor$ and $\omega_{\delta}>0,$ the
pair $(\overline{N},\overline{d}^{\overline{N}})$
with $\overline{N}=\rho_{T}$ and
$\overline{d}^{\overline{N}}=(d_{min},...,d_{min})$ is the UNIQUE SOLUTION of
$(P_{2})$.
IF $\rho_{T}\not\in\mathbb{N},$
$\lceil\lambda_{T}\rceil<\lfloor\rho_{T}\rfloor$ and $\omega_{\delta}>0,$
take $(\overline{N}_{1},\overline{d}_{1}^{\overline{N}_{1}}),$ with
$\overline{N}_{1}=\lfloor\rho_{T}\rfloor$ and
$\overline{d}_{1}^{\overline{N}_{1}}=(\overline{d}_{1},...,\overline{d}_{1})$
where
$\overline{d}_{1}=\dfrac{-\alpha_{T}\overline{N}_{1}+\sqrt{(\alpha_{T}\overline{N}_{1})^{2}+4\beta_{T}\overline{N}_{1}\gamma_{T}}}{2\beta_{T}\overline{N}_{1}}$.
Also take $\overline{N}_{2}=\lceil\rho_{T}\rceil$ and
$\overline{d}_{2}^{\overline{N}_{2}}=(d_{min},...,d_{min})$.
CALCULATE: $E_{OAR}(\overline{N}_{1},\overline{d}_{1}^{\overline{N}_{1}})$ and
$E_{OAR}(\overline{N}_{2},\overline{d}_{2}^{\overline{N}_{2}}).$
A SOLUTION of $(P_{2})$ is the pair
$(\overline{N},\overline{d}^{\overline{N}})$ that minimizes $E_{OAR}$ between
them.
IF $\rho_{T}>1,$ $\lceil\lambda_{T}\rceil<\lfloor\rho_{T}\rfloor$ and
$\omega_{\delta}<0,$
take $\overline{N}=\lceil\lambda_{T}\rceil$ and CALCULATE
$M=\dfrac{\overline{N}\varphi_{T}(d_{max})-\gamma_{T}}{\varphi_{T}(d_{max})-\varphi_{T}(d_{min})}$.
IF $M\in\mathbb{N}\cup\\{0\\}$, take $K=M$ and
$\overline{d}^{\overline{N}}=(\underbrace{d_{min},...,d_{min}}_{K},\underbrace{d_{max},...,d_{max}}_{\overline{N}-K})$.
IF $M\not\in\mathbb{N}\cup\\{0\\}$, take $K=\lfloor M\rfloor$ and
$\overline{d}^{\overline{N}}=(\underbrace{d_{min},...,d_{min}}_{K},d^{*},\underbrace{d_{max},...,d_{max}}_{\overline{N}-K-1}),$
with $d^{*}>0$ and
$\varphi_{T}(d^{*})=\gamma_{T}-K\varphi_{T}(d_{min})-(\overline{N}-K-1)\varphi_{T}(d_{max})$.
A SOLUTION of $(P_{2})$ is the pair
$(\overline{N},\overline{d}^{\overline{N}}).$
IF $\rho_{T}>1$, $\lceil\lambda_{T}\rceil<\lfloor\rho_{T}\rfloor$ and
$\omega_{\delta}=0,$ ANY FEASIBLE PAIR $(\overline{N},d)$ such that
$E_{T}(\overline{N},d)=\gamma_{T}$ is a SOLUTION for $(P_{2})$.
In particular, the pairs $(\overline{N},\overline{d}^{\overline{N}})$ with
$\overline{N}\in\\{\lceil\lambda_{T}\rceil,\ldots,\lfloor\rho_{T}\rfloor\\}$
and $\overline{d}^{\overline{N}}=(\overline{d}_{1},...,\overline{d}_{1}),$
where
$\overline{d}_{1}=\dfrac{-\alpha_{T}\overline{N}+\sqrt{(\alpha_{T}\overline{N})^{2}+4\beta_{T}\overline{N}\gamma_{T}}}{2\beta_{T}\overline{N}}.$
Table 2. Complete solution for problem $(P_{2})$ in algorithmic form.

###  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}$ and
$\tilde{d}_{1},...,\tilde{d}_{\tilde{N}}$ are said to be biologically
equivalent for a certain tumor with characteristic parameters $\alpha_{T}$ and
$\beta_{T}$ when they have the same effect, that is,

$\displaystyle\alpha_{T}\sum_{i=1}^{N}d_{i}+\beta_{T}\sum_{i=1}^{N}d_{i}^{2}=\displaystyle\alpha_{T}\sum_{i=1}^{\tilde{N}}\tilde{d}_{i}+\beta_{T}\sum_{i=1}^{\tilde{N}}\tilde{d}_{i}^{2}.$

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

$N(\alpha_{T}d+\beta_{T}d^{2})=\tilde{N}(\alpha_{T}\tilde{d}+\beta_{T}\tilde{d}^{2}).$

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 formulated
in mathematical terms as the following optimization problem:

$(P_{3})\left\\{\begin{array}[]{ll}\begin{aligned} &\textrm{Minimize
}\sum_{i=1}^{N}d_{i},\\\ &\textrm{subject to
}N\in\mathbb{N},d\in\mathbb{R}^{N}\textrm{ such that }\\\
&\alpha_{T}\sum_{i=1}^{N}d_{i}+\beta_{T}\sum_{i=1}^{N}d_{i}^{2}=\gamma_{T},\\\
&d_{min}\leq d_{i}\leq d_{max}\textrm{,
}i=1,...,N.\end{aligned}\end{array}\right.$

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 $\omega_{\delta}<0$, see Theorem 13, and consequently
the hypofractionated protocols are also optimal in this sense.

##  4\. 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

$\omega_{\delta}=\dfrac{\alpha_{T}}{\beta_{T}}-\dfrac{\alpha_{0}}{\beta_{0}\delta}.$

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 $\omega_{\delta}>0$, we have shown that the optimal number of doses $N$
are $\lfloor\rho_{0}\rfloor$ for $(P_{1})$ and $\lfloor\rho_{T}\rfloor$ or
$\lceil\rho_{T}\rceil$ for $(P_{2})$, the upper values of their ranges of
interest (i.e. hyperfractionated type treatments) with equal doses; while in
the case $\omega_{\delta}<0,$ the optimal values of $N$ are
$\lfloor\lambda_{0}\rfloor$ or $\lceil\lambda_{0}\rceil$ for $(P_{1})$ and
$\lfloor\lambda_{T}\rfloor$ 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.

## Acknowledgement

The authors would like to express their gratitude to Prof. Cecilia Pola
(University of Cantabria) for fruitful discussions and helpful comments.

## References

  * [1] Bertsekas DP (2003) Nonlinear Programming. Athena Scientific, Belmont, Massachusetts.
  * [2] Bertuzzi A, Bruni C, Papa F, Sinisgalli C (2013) Optimal solution for a cancer radiotherapy problem, J Math Biol 66, 311-349.
  * [3] Bortfeld T, Ramakrishnan J, Tsitsiklis JN, Unkelbach J (2015) Optimization of radiation therapy fractionation schedules in the presence of tumor repopulation, INFORMS J Comput 27 (4), 788-803.
  * [4] Brenner DJ (2008) The linear-quadratic model is an appropriate methodology for determining iso-effective doses at large doses per fraction. Semin Radiat Oncol 18 (4), 234-239.
  * [5] Bruni C, Conte F, Papa F, Sinisgalli C (2015) Optimal weekly scheduling in fractionated radiotherapy: effect of an upper bound on the dose fraction size, J Math Biol 71, 361-398.
  * [6] Bruni C, Conte F, Papa F, Sinisgalli C (2019) Optimal number and sizes of the doses in fractionated radiotherapy according to the LQ model, Math Med Biol 36, 1-53.
  * [7] Fernández L (2020) Problemas de optimización asociados a los tratamientos de radioterapia. Final Undergraduate Project, University of Cantabria.
  * [8] Jones B, Tan LT, Dale RG (1995) Derivation of the optimum dose per fraction from the linear quadratic model, Br J Radiol 68, 894-902.
  * [9] McMahon SJ (2019) The linear quadratic model: usage, interpretation and challenges, Phys Med Biol 64, 1-24.
  * [10] Mizuta M, Takao S, Date H, Kishimoto N, Sutherland KL, Rikiya O, Shirato H (2012) A mathematical study to select fractionation regimen based on physical dose distribution and the linear-quadratic model. Int J Radiat Oncol Biol Phys 84, 829-833.
  * [11] Radiation therapy dose calculator of the French Society of Young Radiation Oncologists (2021) http://www.sfjro.fr/ilq/en/
  * [12] Radiotherapy dose fractionation (2019) The Royal College of Radiologist, London.
  * [13] Radiotherapy risk profile (2008). Technical Manual. World Health Organization, Geneva.
  * [14] Saberian F, Ghate A, Kim M (2015) A two-variable linear program solves the standard linear-quadratic formulation of the fractionation problem in cancer radiotherapy. Oper Res Lett 43, 254-258.
  * [15] Van Leeuwen CM, Oei AL, Crezee J, Bel A, Franken NAP, Stalpers LJA, Kok HP (2018) The alfa and beta of tumours: a review of parameters of the linear-quadratic model, derived from clinical radiotherapy studies. Radiat Oncol 13, 1-11.

## Appendix $1$: Proof of Theorem 8

The expression given in $i)$ is derived exactly as in Theorem 6 for the values
$N\leq\lambda_{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. To that end, it is enough to show that if
$(\underbrace{d_{min},...,d_{min}}_{K},d^{*},\underbrace{d_{max},...,d_{max}}_{N-K-1})$
is a solution for $(P_{1}^{N})$ and
$(\underbrace{d_{min},...,d_{min}}_{\tilde{K}},\tilde{d}^{*},\underbrace{d_{max},...,d_{max}}_{N-\tilde{K}})$
is a solution for $(P_{1}^{N+1})$, with $N>\lambda_{0}$, then the following
relation holds

(43)
$Kd_{min}+d^{*}+(N-K-1)d_{max}\leq\tilde{K}d_{min}+\tilde{d}^{*}+(N-\tilde{K})d_{max}.$

Together with (8) and the assumption $\omega_{\delta}<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 $\overline{N}=\lceil\lambda_{0}\rceil$, the lowest value of $N$ in
the set $(\lambda_{0},\rho_{0}]\cap\mathbb{N}$.

Comparing their expressions in the form (26) with $N+1$ and $N$, resp., we
conclude that $\tilde{K}\geq K+1$. Hence, if we denote
$K_{0}=\tilde{K}-K\in\mathbb{N}$, the inequality (43) can be written as

(44) $d^{*}-\tilde{d}^{*}\leq K_{0}d_{min}+(1-K_{0})d_{max}.$

Let us recall that $d^{\star}$ satisfies (24) and $\tilde{d}^{*}$ verifies

(45)
$\varphi_{0}(\tilde{d}^{*})=\gamma_{OAR}-\tilde{K}\varphi_{0}(d_{min})-(N-\tilde{K})\varphi_{0}(d_{max}).$

We will show that (44) holds dividing the argumentation in three cases:

Case 1.- Suppose that $K_{0}d_{min}+(1-K_{0})d_{max}\leq 0$. We choose the
point

$(d_{1},...,d_{N})=(\underbrace{\dfrac{K_{0}}{K_{0}-1}d_{min},...,\dfrac{K_{0}}{K_{0}-1}d_{min}}_{K_{0}-1},\underbrace{d_{min},...,d_{min}}_{\tilde{K}-K_{0}},\tilde{d}^{*},\underbrace{d_{max},...,d_{max}}_{N-\tilde{K}}),$

that under the assumption satisfies the bounds restrictions and

$E_{OAR}(N,d)=\alpha_{0}\delta\left(\tilde{K}d_{min}+\tilde{d}^{*}+(N-\tilde{K})d_{max}\right)+$
$+\beta_{0}\delta^{2}\left(\dfrac{(K_{0})^{2}}{K_{0}-1}d_{min}^{2}+(\tilde{K}-K_{0})d_{min}^{2}+(\tilde{d}^{*})^{2}+(N-\tilde{K})d_{max}^{2}\right)\geq$
$\geq\tilde{K}\varphi_{0}(d_{min})+\varphi_{0}(\tilde{d}^{*})+(N-\tilde{K})\varphi_{0}(d_{max})=\gamma_{OAR}.$

This means that it is feasible for the problem $(\ref{PEQ-})$. Taking into
account Remark 3, we get (44).

Case 2.- Suppose now that $K_{0}d_{min}+(1-K_{0})d_{max}>0$ and furthermore

$K_{0}\varphi_{0}(d_{min})+(1-K_{0})\varphi_{0}(d_{max})\leq 0$. We can argue
similarly choosing

(46)
$(d_{1},...,d_{N})=(\underbrace{d_{min},...,d_{min}}_{K},\tilde{d}^{*},\underbrace{d_{max},...,d_{max}}_{N-K-1}).$

Due to (45) and the hypothesis we have

$E_{OAR}(N,d)=K\varphi_{0}(d_{min})+\varphi_{0}(\tilde{d}^{*})+(N-K-1)\varphi_{0}(d_{max})=$
$=\gamma_{OAR}-K_{0}\varphi_{0}(d_{min})+(K_{0}-1)\varphi_{0}(d_{max})\geq\gamma_{OAR}.$

Again, we have a feasible point for the problem $(\ref{PEQ-})$ and therefore
we deduce $d^{*}\leq\tilde{d}^{*}$ and hence (44), because

$d^{*}-\tilde{d}^{*}\leq 0<K_{0}d_{min}+(1-K_{0})d_{max}.$

Case 3.- Finally, suppose that $K_{0}d_{min}+(1-K_{0})d_{max}>0$ and moreover

$K_{0}\varphi_{0}(d_{min})+(1-K_{0})\varphi_{0}(d_{max})>0$. Here, we
introduce the auxiliary function defined for $s\in[0,1]$ by

$G(s)=\sqrt{\alpha_{0}^{2}+4\beta_{0}\left(\gamma_{OAR}-(\tilde{K}-sK_{0})\varphi_{0}(d_{min})-(N-\tilde{K}+s(K_{0}-1))\varphi_{0}(d_{max})\right)}.$

Solving the quadratic equations (24) and (45), it is easy to derive that

$d^{*}=\dfrac{-\alpha_{0}+G(1)}{2\beta_{0}\delta},\ \
\tilde{d}^{*}=\dfrac{-\alpha_{0}+G(0)}{2\beta_{0}\delta}.$

Using the Mean Value Theorem, we deduce that there exists $\theta\in(0,1)$
such that

$d^{*}-\tilde{d}^{*}=\dfrac{G(1)-G(0)}{2\beta_{0}\delta}=\dfrac{G^{\prime}(\theta)}{2\beta_{0}\delta}=\dfrac{K_{0}\varphi_{0}(d_{min})+(1-K_{0})\varphi_{0}(d_{max})}{G(\theta)\delta}.$

Therefore, the inequality (44) is equivalent to

(47)
$\dfrac{K_{0}\varphi_{0}(d_{min})+(1-K_{0})\varphi_{0}(d_{max})}{K_{0}d_{min}+(1-K_{0})d_{max}}\leq
G(\theta)\delta.$

Under the present hypotheses, the function $G$ is strictly increasing and,
since $\theta$ is an unknown value in $(0,1)$, we will verify that (47) is
valid if it holds for $\theta=0$. On the other hand, the value
$K_{0}\in\mathbb{N}$ is also unknown, but we can verify that the function

$F(m)=\dfrac{m\varphi_{0}(d_{min})+(1-m)\varphi_{0}(d_{max})}{md_{min}+(1-m)d_{max}}=$
$=\alpha_{0}\delta+\beta_{0}\delta^{2}\left(\dfrac{md_{min}^{2}+(1-m)d_{max}^{2}}{md_{min}+(1-m)d_{max}}\right),$

is strictly decreasing, because

$F^{\prime}(m)=\beta_{0}\delta^{2}\dfrac{d_{min}d_{max}(d_{min}-d_{max})}{(md_{min}+(1-m)d_{max})^{2}}<0.$

Hence, the inequality (47) will be true if $F(1)\leq G(0)\delta.$ We conclude
by noting that

$\tilde{d}^{*}\in[d_{min},d_{max}]\Longleftrightarrow\varphi_{0}(\tilde{d}^{*})\in[\varphi_{0}(d_{min}),\varphi_{0}(d_{max})]\Longleftrightarrow$
$\stackrel{{\scriptstyle(\ref{E620})}}{{\Longleftrightarrow}}\gamma_{OAR}-\tilde{K}\varphi_{0}(d_{min})-(N-\tilde{K})\varphi_{0}(d_{max})\in[\varphi_{0}(d_{min}),\varphi_{0}(d_{max})].$

Then,

$G(0)=\sqrt{\alpha_{0}^{2}+4\beta_{0}\left(\gamma_{OAR}-\tilde{K}\varphi_{0}(d_{min})-(N-\tilde{K})\varphi_{0}(d_{max})\right)}\geq$
$\geq\sqrt{\alpha_{0}^{2}+4\beta_{0}\varphi_{0}(d_{min}})\geq\alpha_{0}+\beta_{0}\delta
d_{min}=\dfrac{F(1)}{\delta},$

as asserted.

## Appendix $2$: Proof of Theorem 12

It follows the same lines to that of Theorem 7.

Case a).- Assume $\rho_{T}\in\mathbb{N}$, $\rho_{T}\geq 2.$

As usual, we divide the interval for feasible values of $N$ in two parts:
$[\lambda_{T},\rho_{T})\cap\mathbb{N}$ and $[\rho_{T},+\infty)\cap\mathbb{N}.$

In order to study the dependence with respect to $N$ in the interval
$[\lambda_{T},\rho_{T})$, thanks to Proposition 4 (with $\omega_{\delta}>0$)
and (36), it is enough to consider the auxiliary function

$\psi_{1}(N)=N\overline{d}_{1}=\dfrac{-\alpha_{T}N+\sqrt{(\alpha_{T}N)^{2}+4\beta_{T}N\gamma_{T}}}{2\beta_{T}}.$

Once more, it follows easily that $\psi_{1}(N)$ is an strictly increasing
function. Since we are assuming $\rho_{T}\in\mathbb{N}$ and $\rho_{T}\geq 2$,
then $\psi_{1}$ will take its maximum value in the set
$[\lambda_{T},\rho_{T})\cap\mathbb{N}$ at $N_{1}=\rho_{T}-1$. Therefore, the
candidate for solution to problem $(P_{2})$ is given by the pair
$(N_{1},\overline{d}^{N_{1}})$ with
$\overline{d}^{N_{1}}=(\overline{d}_{1},...,\overline{d}_{1}),$ where
$\overline{d}_{1}$ is given by

(48)
$\overline{d}_{1}=\dfrac{-\alpha_{T}N_{1}+\sqrt{(\alpha_{T}N_{1})^{2}+4\beta_{T}N_{1}\gamma_{T}}}{2\beta_{T}N_{1}}.$

On the other hand, in the interval $[\rho_{T},+\infty)$, we know that the
other candidate for solution to problem $(P_{2})$ is given by the pair
$(N_{2},\overline{d}^{N_{2}})$ with $N_{2}=\rho_{T}$ and
$\overline{d}^{N_{2}}=(d_{min},...,d_{min}).$

To derive that $(N_{2},\overline{d}^{N_{2}})$ is the unique solution to
problem $(P_{2})$, we will show that

(49)
$E_{OAR}(N_{2},\overline{d}^{N_{2}})<E_{OAR}(N_{1},\overline{d}^{N_{1}}).$

Following the same idea to that of the proof of Theorem 7, we introduce the
auxiliary function

$H_{1}(x)=N_{1}(x\overline{d}_{1}+\overline{d}_{1}^{2})-N_{2}(xd_{min}+d_{min}^{2}),\
\ \ x\in[\dfrac{\alpha_{0}}{\beta_{0}\delta},+\infty).$

Taking into account (48) and that $N_{2}\varphi_{T}(d_{min})=\gamma_{T}$ (by
the definition of $\rho_{T}$), it can be checked that
$H_{1}^{\prime}(x)=N_{1}\overline{d}_{1}-N_{2}d_{min}<0,$ since $N_{1}<N_{2}.$

Using that also $\gamma_{T}=N_{1}\varphi_{T}(\overline{d}_{1})$, we get that
$H_{1}\left(\dfrac{\alpha_{T}}{\beta_{T}}\right)=0$ and from the assumption
$\omega_{\delta}>0$ (see (9)), it follows that
$H_{1}\left(\dfrac{\alpha_{0}}{{\beta_{0}\delta}}\right)>0$, which is
equivalent to (49).

Case b).- Assume $\rho_{T}\not\in\mathbb{N}$. Here, the optimal value of $N$
in the interval $[\lambda_{T},\rho_{T})$ is $N_{1}=\lfloor\rho_{T}\rfloor$ and
$\overline{d}^{N_{1}}=(\overline{d}_{1},...,\overline{d}_{1})$ with
$\overline{d}_{1}$ given by (48). In the interval $[\rho_{T},+\infty)$, the
other candidate is $N_{2}=\lceil\rho_{T}\rceil$ with
$\overline{d}^{N_{2}}=(d_{min},...,d_{min}).$

When $\rho_{T}\not\in\mathbb{N},$ any of them can provide the unique solution
to problem $(P_{2})$ (see for instance Example 6).

## Appendix $3$: Proof of Theorem 13

When $\omega_{\delta}<0$, it is still true that $N_{2}=\lceil\rho_{T}\rceil$
and $\overline{d}^{N_{2}}=(d_{min},...,d_{min}).$ Arguing as in Appendix $1$,
the candidate when $N$ runs $[\lambda_{T},\rho_{T})\cap\mathbb{N}$ is
$N_{1}=\lceil\lambda_{T}\rceil$ with $\overline{d}^{N_{1}}$ given by Theorem
13-$a)$ or $b)$ and $\overline{N}=N_{1}$, thanks to Theorem 11$-ii)$. We will
conclude by showing that

(50) $E_{OAR}(N_{1},\overline{d}^{N_{1}})\leq
E_{OAR}(N_{2},\overline{d}^{N_{2}}).$

Let us argue with the expression $b)$ for $\overline{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\varphi_{0}(d_{min})+\varphi_{0}(d^{\star})+(N_{1}-K-1)\varphi_{0}(d_{max})\leq
N_{2}\varphi_{0}(d_{min}).$

For proving (51), we consider again a linear function such as

$H_{2}(x)=(N_{2}-K)(xd_{min}+d_{min}^{2})-(xd^{*}+(d^{*})^{2})-(N_{1}-K-1)(xd_{max}+d_{max}^{2}).$

By construction, we know that

$K\varphi_{T}(d_{min})+\varphi_{T}(d^{\star})+(N_{1}-K-1)\varphi_{T}(d_{max})=\gamma_{T}\leq
N_{2}\varphi_{T}(d_{min}).$

This is equivalent to say that
$H_{2}\left(\dfrac{\alpha_{T}}{\beta_{T}}\right)\geq 0$.

If $H_{2}$ is an increasing function, since $\omega_{\delta}<0$, we will have

$H_{2}\left(\dfrac{\alpha_{0}}{\beta_{0}\delta}\right)\geq
H_{2}\left(\dfrac{\alpha_{T}}{\beta_{T}}\right)\geq 0,$

which gives (51).

So, taking into account that

$H^{\prime}_{2}(x)=(N_{2}-K)d_{min}-d^{*}-(N_{1}-K-1)d_{max},$

let us finish the proof by showing that $H^{\prime}_{2}(x)\geq 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}\leq N_{1}d_{max},$ we can argue as in the proof of Theorem
8, taking the point

$\tilde{d}=(d_{1},...,d_{N_{1}})=(\dfrac{N_{2}}{N_{1}}d_{min},...,\dfrac{N_{2}}{N_{1}}d_{min}),$

that satisfies the bounds restrictions and

$E_{T}(N_{1},\tilde{d})=\alpha_{T}N_{2}d_{min}+\beta_{T}\dfrac{(N_{2})^{2}}{N_{1}}d_{min}^{2}\geq
N_{2}\varphi_{T}(d_{min})\geq\gamma_{T}.$

This means that it is feasible for the problem $(P_{2}^{N_{1},-})$. Taking
into account that $(N_{1},\overline{d}^{N_{1}})$ is a solution for that
problem, see Proposition 4-$ii)$ and Remark 3, we get $H^{\prime}_{2}(x)\geq
0$.