Analytical Solution to the Radiotherapy Fractionation Problem Including Dose Bound Constraints

According to the World Health Organization, [11], “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, [8]. 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, [14]. 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 [8], “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 http://www.sfjro.fr/ilq/en/.

Hence, given N doses, \(d_{1},...,d_{N}\), possibly different, the probability of accumulated survival can be expressed by

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

It should be noted that throughout this work, we will not consider tumor repopulation between doses.

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

where \(\alpha _{0}\) and \(\beta _{0}\) are the parameters associated with 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 [14]. 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 two common strategies for fractionating radiotherapy treatments:

• Hypofractionation: Deliver higher doses of radiation in fewer sessions. This strategy results in a significant reduction in 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:

1. (a)

   To maximize the effect of radiation on the tumor, restricting the effect produced on one OAR (healing approach) in Sect. 2 and

2. (b)

   To minimize the effect of radiation on an OAR, maintaining enough effect of radiation on the tumor (palliative approach) in Sect. 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} \le d \le 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 [10]. Depending on the tumor type, the values of \(d_{min}\) and \(d_{max}\) can be tuned, but the reference values \(d_{min} = 1 \ Gy\) and \(d_{max} = 6 \ Gy\) could be a valid generic choice (one exception is Example 4, see as follows). In this sense one cannot find in [10] a single treatment recommendation 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 [7], but only for uniform dose treatments. Our approach here includes also the study for nonuniform protocols. The origin of this work was an academic project, which has been extended here with the detailed study of the dependence of the solution with respect to N.

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 [9] and [12].

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 (1), \(E_{OAR}(N,d)\) by (2) 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) \le \gamma _{OAR}\) can be interpreted in the sense that the percentage of surviving cells of the OAR should be greater than or equal to \(\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 [12] (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 [7]. In addition, our approach to the problem is more useful since the number of doses N is not initially set as in [2] and [12]. The case including repopulation was studied in [3], only assuming the nonnegativity of \(d_i\).

From a mathematical point of view, this is a mixed-integer 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 \le i \le 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

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 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} \ge 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 \le \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 no other pair is feasible. When \(\rho _0 \in (1,2)\), the value \(N =1\) is still the unique 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}=\psi _0(1)\). Let us stress that \(\varphi _0(\overline{d}_{0}) = \gamma _{OAR}\).

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

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})\). \(\square\)

The existence of a solution for each problem \((P_{1}^{N})\) is proved below:

Theorem 2

Let us assume \(d_{min}>0\), \(\rho _{0} \ge 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. \(\square\)

Remark 1

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, but here we will see that it can be done analytically.

2. (b)

   Unless all the components of the solution are equal, the uniqueness of solution fails: it is enough to take two indices \(i,j \in \lbrace 1,...,\overline{N} \rbrace\) such that \(\overline{d}_{i} \ne \overline{d}_{j}\) and interchange these coordinates to generate a new solution.

3. (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.

4. (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

It is clear that for all feasible points we have \(E_{T}(N,d) \le E_{OAR}(N,d) \le 10.\) Let us stress that here N can take any natural value, without restrictions. Inspired by Theorem 4, let us consider the sequence given by

It is easy to check that it is feasible for \(N \ge 4\),

Hence, the problem \((P_{10})\) has no solution \((\overline{N},\overline{d})\): the supremum value 10 cannot be attained since it should happen that

which is clearly impossible.

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} \ge 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 any solution \(\overline{d}^{N}\) of \((P_{1}^{N})\).

Proof

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

Since \(\lambda _{0} < N\), we know that there exists some index \(j \in \lbrace 1,...,N \rbrace\) 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

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

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

Therefore, the objective function can be written as

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} \ge 2\), \(N \in \left( \lambda _{0}, \rho _{0} \right] \cap \mathbb {N}\) and denote

1. (i)

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

   $$\begin{aligned} (P_{1}^{N,+}) \text { Maximize } \displaystyle \sum _{i=1}^{N} d_{i}, \text { subject to } d \in \mathbb {K}_{1}^{N}, \end{aligned}$$

   where

   $$\begin{aligned} \mathbb {K}_{1}^{N} = \lbrace d \in \mathbb {R}^{N} : E_{OAR}(N,d) = \gamma _{OAR}, d_{min} \le d_{i} \le d_{max}, 1 \le i \le N \rbrace . \end{aligned}$$
2. (ii)

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

   $$\begin{aligned} (P_{1}^{N,-}) \text { Minimize } \displaystyle \sum _{i=1}^{N} d_{i}, \text { subject to } d \in \mathbb {K}_{1}^{N}. \end{aligned}$$
3. (iii)

   If \(\omega _{\delta }=0,\) then every feasible point for \((P_{1}^{N})\) is a solution.

Remark 2

1. (a)

   The idea of this transformation can be found in [9] in the context of the problem \((P_{2})\) that we will study in the next section.

2. (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.

3. (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 cumulative dose, while in the second the cumulative dose remains minimum.

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:

In Fig. 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.

2D Example

Full size image

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 general \(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} \ge 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 \(\overline{d}_{0}=\psi _0(N).\)

Proof

By using the Cauchy-Schwarz inequality, we have

Therefore, for each feasible point it follows that

Defining \(q(z)=\beta _{0} \delta ^{2}z^{2}/N+\alpha _{0} \delta z - \gamma _{OAR}\), previous inequality can be rewritten as

Taking into account that the polynomial q can be factorized in the form \(q(z)= \beta _{0} \delta ^{2}(z-z_{1})(z-z_{2})/N\) with \(z_1< 0 < z_2\), we know that relation (8) 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 (7), we deduce that

Hence

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}=z_{2}/N\) and this is the desired expression, see (4).

Let us emphasize that \(\overline{d}_{0}\) satisfies \(d_{min} \le \overline{d}_{0} \le d_{max}\), thanks to the hypothesis \(N \in \left( \lambda _{0}, \rho _{0} \right]\). \(\square\)

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, see [1], to deduce the existence of real numbers \(\overline{\alpha } \in [0,+\infty ), \overline{\lambda } \in \mathbb {R}\) and \(\lbrace \overline{\mu }_{i} \rbrace _{i=1}^{2N} \subset [0,+\infty )\) verifying

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 \lbrace 1,...,N \rbrace\). Then, thanks to (11) we deduce that \(\overline{\mu }_{i}=0, \forall i \in \lbrace 1,...,2N \rbrace .\) In this case, (10) reads:

If \(\overline{\lambda } = 0\), identity (12) implies that \(\overline{\alpha } = 0,\) but this is not possible by (9). Therefore \(\overline{\lambda } \ne 0\) and from (12) we get \(\overline{d}_{1}=...=\overline{d}_{N}\). In other words, we arrive to the solution of problem \((P_{1}^{N,+})\), contradicting our initial hypothesis.

Consequently, there exists (at least) one index \(j \in \lbrace 1,...,N \rbrace\) such that \(\overline{d}_{j} \in \lbrace d_{min}, d_{max} \rbrace\). 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:

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

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

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

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 \lbrace 1,...,N-1 \rbrace\) such that \(\overline{d}_{j} \in \lbrace d_{min},d_{max} \rbrace\) 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 \in \left( \lambda _{0}, \rho _{0} \right] \cap \mathbb {N}\), the solution of \((P_{1}^{N,-})\) has one of the following structures:

or

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

with \(\varphi _0\) defined in (3).

We can characterize the unknown value K as follows:

1. (a)

   In case (13), by using the equality restriction we derive that

   $$\begin{aligned} K = \dfrac{N \varphi _{0}(d_{max})-\gamma _{OAR}}{\varphi _{0}(d_{max})-\varphi _{0}(d_{min})}. \end{aligned}$$

   (16)

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

2. (b)

   In case (14), since \(\varphi _{0}\) is a strictly increasing function in \([0,+\infty )\), we know that

   $$\begin{aligned} \varphi _{0}(d_{min})< \varphi _{0}(d^{*}) < \varphi _{0}(d_{max}), \end{aligned}$$

   and using (15) we get that

   $$\begin{aligned} 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}, \end{aligned}$$

   which means that

   $$\begin{aligned} K = \displaystyle \lfloor \dfrac{N\varphi _{0}(d_{max})-\gamma _{OAR}}{\varphi _{0}(d_{max})-\varphi _{0}(d_{min})} \rfloor . \end{aligned}$$

   (17)

Taking into account conditions (16) and (17), it is easy to conclude that the latter structure (14) is more frequently found in practice than (13). Previous argumentations lead us to the following result:

Theorem 5

Let us assume \(d_{min}>0\), \(\rho _{0} \ge 2\) and \(N \in \left( \lambda _{0}, \rho _{0} \right] \cap \mathbb {N}.\) Then, a solution to problem \((P_1^{N,-})\) is given by

1. (a)

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

2. (b)

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

Remark 3

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

We will use this property in the proof of Theorem 8.

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} \ge 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\). \(\square\)

We will continue to analyze the most common case: when \(\lfloor \lambda _0 \rfloor < \lfloor \rho _0 \rfloor\). In the trivial case \(\omega _{\delta } = 0\), the function to be minimized and the one defining the restriction are proportional. Therefore, we can derive the following result:

Proposition 3

Let us assume \(d_{min}>0,\) \(\rho _{0} \ge 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 \(\overline{d}_{0}=\psi _0(\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

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]\). \(\square\)

Theorem 7

Let us assume \(d_{min}>0,\) \(\rho _{0} \ge 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})\), where \(\overline{d}_{0}=\psi _0(\overline{N}).\)

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} = \psi _0(N)\). 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 \(\phi _0(N)= N \psi _0(N).\) Here, it follows easily that \(\phi _0\) is a 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

To that end, let us consider the linear function

with \(\overline{d}_0= \psi _0(N_2).\)

Using that \(N_1 \varphi _0(d_{max}) \le \gamma _{OAR} = N_2\varphi _0(\overline{d}_0)\) by the admissibility, we get that \(H_1(\alpha _0/(\beta _0 \delta )) \ge 0\). Also, it can be checked that \(H_1'(x) = N_2\overline{d}_0-N_1 d_{max} > 0,\) because \(N_1 < N_2\). Then, from the assumption \(\omega _{\delta } >0\) (see (6)), it follows that \(H_1(\alpha _T/\beta _T) > 0\), which is equivalent to (19). \(\square\)

In the case \(\omega _{\delta } < 0,\) the situation is more complicated and it is detailed in the next result:

Theorem 8

Let us assume \(d_{min}>0,\) \(\rho _{0} \ge 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:

1. (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}),\)

2. (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 (16) with \(N = \overline{N}\) belongs to \(\mathbb {N} \cup \{0\},\) or

3. (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 by (17) and \(d^\star\) satisfies (15) with \(N = \overline{N}.\)

Proof

The expression given in i) is derived exactly as in Theorem 6 for the values \(N \le \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

Together with (5) and the assumption \(\omega _{\delta } < 0\), this implies iii), because (20) 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 form (17) with \(N+1\) and N, resp., we conclude that \(\tilde{K} \ge K+1\). Hence, if we denote \(K_0 = \tilde{K} - K \in \mathbb {N}\), inequality (20) can be written as

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

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

\(\underline{\text{Case 1.- Suppose that} \ K_0 d_{min}+ (1-K_0)d_{max} \le 0}\). We choose the point

that under the assumption satisfies the bounds restrictions and

This means that it is feasible for problem (18). Taking into account Remark 3, we get (21).

\(\underline{\text{Case 2.- Suppose now that} \ K_0 d_{min}+ (1-K_0)d_{max} > 0}\) and furthermore

\(\underline{K_0 \varphi _0(d_{min})+ (1-K_0)\varphi _0(d_{max}) \leq 0}\). We can argue similarly choosing

Due to (22) and the hypothesis we have

Again, we have a feasible point for problem (18) and therefore we deduce \(d^{*} \le \tilde{d}^{*}\) and hence (21), because \(d^{*} - \tilde{d}^{*} \le 0 < K_0 d_{min}+ (1-K_0)d_{max}.\)

\(\underline{\text{Case 3.- Finally, suppose that} \ K_0 d_{min}+ (1-K_0)d_{max} > 0}\) together with

\(\underline{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

Solving the quadratic Eqs. (15) and (22), it is easy to derive that

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

Therefore, inequality (21) is equivalent to

Under the present hypotheses, the function G is strictly increasing and, since \(\theta\) is an unknown value in (0, 1), we will verify that (23) 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

is strictly decreasing, because

Hence, inequality (23) will be true if \(F(1) \le G(0)\delta\). We conclude by noting that

Then,

as asserted. \(\square\)

Remark 4

1. (a)

   Note that all expressions included in Theorems 6–8 and Proposition 3 can be explicitly calculated from the problem data.

2. (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.

3. (c)

   One interesting case appears when \(\alpha _T/\beta _T < \alpha _0/\beta _0,\) because then \(\omega _{\delta } < 0\) for all \(\delta \in (0,1]\) and the optimal regimen is always of hypofractionated type, regardless 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 \(\alpha _T/\beta _T \approx 1.5 \ Gy,\) while \(\alpha _0/\beta _0 = 2 \ Gy\), see [9] and [14].

4. (d)

   After Remark 2-b), 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).

5. (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 problem includes a repopulation term 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, it is clear that 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 [9], together with \(d_{min}=1 \ Gy,\) and \(d_{max}=6 \ Gy\). Then, the problem reads

1. (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 [10] (see 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.

2. (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.

3. (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.

We have summarized the resolution of problem \((P_1)\) in algorithmic form in Table 2.

Example 4

Radiotherapy treatment of spinal tumors is complicated by the proximity to the major nerve tracts and the risk of radiation myelopathy (RM). For this reason, doses cannot be easily increased. Stereotactic Body Radiation Therapy (SBRT) is a modern technique allowing to deliver much higher radiation doses in fewer sessions and with great precision. Due to its growing interest, there are recent studies on the determination of maximum doses and the number of fractions for spinal SBRT. This case fits into our framework because the spinal cord is a single OAR.

For illustrative purposes only, we present a numerical example by choosing the following values for the parameters: \(\alpha _T/\beta _T = 3.8 \ \ Gy^{-1}\) (corresponding to a meningioma, see [14]), \(\alpha _0/\beta _0 = 2 \ \ Gy^{-1}\) and \(\delta = 0.16\). In the first three columns of the following table we have summarized the treatments with the same BED (biologically effective dose) as those ones cited in [13] for spinal SBRT and our \(\alpha _T/\beta _T\); in the last two columns we have calculated the effects on the tumor (these values coincide with BED\(\times \alpha _T\) and are all very close, as expected by construction) and on the spinal cord (associated with a risk of up to \(5\%\) of RM):

For these values, \(\omega _{\delta } = -8.7 < 0,\) which corresponds to the hypofractionated case. For a better comparison with the results in Table 1, we have chosen \(d_{min} = 6.55\), \(d_{max} = 16.62\) and \(\gamma _{OAR} = 49.91 \alpha _0 \delta\) (associated with the most harmful treatment). As mentioned in the Introduction, there is controversy over the validity of the LQ model for large doses as \(d_{max}\). This issue has also been noted in [13], but their authors decided to use it due to the lack of clinical validity of the alternative models. Here we will follow the same approach.

Using our algorithm, we arrive to the solution for \((P_{1})\) given by \(\overline{N}=2\) and non-equal doses \(\overline{d}=(7.13,16.62)\). The corresponding effects are \(E_T = 109.81\alpha _T\) and \(E_{OAR} = \gamma _{OAR}.\) We can appreciate that only two fractions (with a cumulative dose of 23.75) are needed to achieve an effect on the tumor approximately \(23\%\) higher than the other treatments in Table 1, reaching the maximum allowed effect on OAR.

Remark 5

Although this work deals with the case of only one OAR, when multiple OARs are involved, there is still some hope to take advantage of our study. If one OAR is clearly more important than the others, we can prioritize it, apply our algorithm with the corresponding constraint and then check if the damage to the others OARs is acceptable. Clearly, this approach may work sometimes, but not always.

When all the OARs have (more or less) the same relevance, assuming that the constraints are written in the form

none of them being deductible from the rest, and denoting

we can construct two different global “fictitious” OARs that combine (24) into one as follows:

1. 1.

   “Adding” all the constraints to obtain

   $$\begin{aligned} K^N_0 = \{(d_1,\ldots ,d_N) \in [d_{min},d_{max}]^N: \alpha _0 \sum _{i=1}^N d_i + \beta _0 \sum _{i=1}^N d^2_i \le \gamma _{OAR}\}, \end{aligned}$$

   with \(\alpha _0 = \sum _{j=1}^m \alpha _j \delta _j,\) \(\beta _0 = \sum _{j=1}^m \beta _j \delta ^2_j\) and \(\gamma _{OAR} = \sum _{j=1}^m \gamma _{OAR_j}.\) This set is the smallest one defined by only one of this type of relations and containing \(K^N\). We can assure that any vector in \(K^N_0\) will satisfy (at least) one of the original constraints (24) and (hopefully) more than one.

2. 2.

   On the other hand, the largest set of this type contained in \(K^N\) is defined by

   $$\begin{aligned} \tilde{K}^N_0 = \{ (d_1,\ldots ,d_N) \in [d_{min},d_{max}]^N: \tilde{\alpha }_0 \sum _{i=1}^N d_i + \tilde{\beta }_0 \sum _{i=1}^N d^2_i \le \tilde{\gamma }_{OAR}\}, \end{aligned}$$

   with \(\tilde{\alpha }_0 = \min \{\frac{\gamma _{OAR_j}}{\beta _j \delta ^2_j}:j=1,\ldots ,m\}\), \(\tilde{\beta }_0 = \min \{\frac{\gamma _{OAR_j}}{\alpha _j \delta _j}:j=1,\ldots ,m\}\) and \(\tilde{\gamma }_{OAR} = \tilde{\alpha }_0 \tilde{\beta }_0\).

Applying our theory to \(K^N_0\) and \(\tilde{K}^N_0\) may provide some insights (upper and lower bounds) about what happens in \(K^N\).

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 one organ at risk (OAR), 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:

where \(E_{OAR}(N,d)\) is given by (2), \(E_{T}(N,d)\) is defined in (1) and \(\gamma _{T}\) is a given positive parameter. Of course, this is also a mixed-integer optimization problem in which the number of radiation doses, \(N \in \mathbb {N}\), is an unknown, as well as the value of the N doses, \(d_{i} \in \mathbb {R}, 1 \le i \le N\).

This problem has recently been studied in the outstanding work [9], but with fixed N and only imposing the nonnegativity constraint for the doses. Moreover, in [9] 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.

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

For large values of N, specifically for \(N \ge \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 has 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. 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})\). \(\square\)

Remark 6

1. (a)

   As for \((P_{1})\), except if all the coordinates of \(\overline{d}^{\overline{N}}\) are equal, the solution for \((P_{2})\) will not be unique, because two different coordinates can be permuted to generate a different one.

2. (b)

   When \(\omega _{\delta } > 0\) (see (6)), the hypothesis \(d_{min} > 0\) is necessary for proving the existence of solution for \((P_2)\). In contrast, when \(\omega _{\delta } < 0\), it is easy to show that problem \((P_2),\) with only the lower bound constraints \(d_i \ge 0\), has as solution \((\overline{N},\overline{d})\) with \(\overline{N} = 1\) and \(\overline{d} = \psi _T(1),\) see [9].

3. (c)

   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 unique feasible 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 \ge \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 any solution.

From now on, the restriction will be taken as one of equality. Here, applying the same procedure as for \((P_1)\) in Sect. 2, the objective function will read

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}\).

1. (i)

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

   $$\begin{aligned} (P_{2}^{N,+}) \text { Maximize } \displaystyle \sum _{i=1}^{N} d_{i}, \text { subject to } d \in \mathbb {K}_{2}^{N}, \end{aligned}$$

   where

   $$\begin{aligned} \mathbb {K}_{2}^{N} = \lbrace d \in \mathbb {R}^{N} : E_{T}(N,d) = \gamma _{T}, d_{min} \le d_{i} \le d_{max}, 1 \le i \le N \rbrace . \end{aligned}$$
2. (ii)

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

   $$\begin{aligned} (P_{2}^{N,-}) \text { Minimize } \displaystyle \sum _{i=1}^{N} d_{i}, \text { subject to } d \in \mathbb {K}_{2}^{N}. \end{aligned}$$
3. (iii)

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

Proof

It suffices to note that \(\alpha _{0} - \beta _{0}\alpha _{T}\delta /\beta _{T} = - \beta _{0} \delta \omega _{\delta },\) where \(\omega _{\delta }\) is defined in (6). \(\square\)

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 Sect. 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,

1. (i)

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

2. (ii)

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

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

   with

   $$\begin{aligned} K = \dfrac{N\varphi _{T}(d_{max}) - \gamma _{T}}{\varphi _{T}(d_{max})-\varphi _{T}(d_{min})} \in \mathbb {N} \cup \{0\}, \ \ \text{ or } \end{aligned}$$

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

   with

   $$\begin{aligned} K = \lfloor \dfrac{N\varphi _{T}(d_{max}) - \gamma _{T}}{\varphi _{T}(d_{max})-\varphi _{T}(d_{min})}\rfloor , \end{aligned}$$

   (26)

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

   $$\begin{aligned} \varphi _{T}(d^{*}) = \gamma _{T} - K \varphi _{T}(d_{min}) - (N-K-1) \varphi _{T}(d_{max}). \end{aligned}$$

   (27)

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 (6)), 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. 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\).

1. (a)

   If \(\rho _T \in \mathbb {N}\), \(\rho _T \ge 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}).\)

2. (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:

   1. (i)

      \(\overline{N} = \lceil \rho _T \rceil\) and \(\overline{d}^{\overline{N}}=(d_{min},...,d_{min}),\) or

   2. (ii)

      \(\overline{N} = \lfloor \rho _T \rfloor\) and \(\overline{d}^{\overline{N}}=(\overline{d}_{1},...,\overline{d}_{1})\), with \(\overline{d}_{1}=\psi _T(\overline{N})\).

Proof

It follows the same lines to that of Theorem 7.

\(\underline{\text{Case a).- Assume} \ \rho _T \in \mathbb {N},\, \rho _T\geq2.}\),   

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\)), it is enough to consider the auxiliary function \(\phi _T(N) = N\psi _T(N).\) Once more, it follows easily that \(\phi _T\) is a strictly increasing function. Since we are assuming \(\rho _T \in \mathbb {N}\) and \(\rho _T \ge 2\), then \(\phi _T\) 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= \psi _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

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

Taking into account the expression of \(\overline{d}_1\) and that \(N_2 \varphi _T(d_{min}) = \gamma _T\) (by the definition of \(\rho _T\)), it can be checked that \(H_2'(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_2(\alpha _T/\beta _T)=0\) and from the assumption \(\omega _{\delta } >0\) (see (6)), it follows that \(H_2(\alpha _0/(\beta _0 \delta )) > 0\), which is equivalent to (28).

\(\underline{\text{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= \psi _T(N_1).\) 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})\). \(\square\)

Finally, the case \(\omega _{\delta } < 0\) is studied in the next theorem:

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

1. (a)

   \(\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\};\) otherwise,

2. (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 (26) and \(d^\star\) satisfies (27), both with \(N = \overline{N}.\)

Proof

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 the proof of Theorem 8, 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

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, inequality (29) is equivalent to

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

By construction, we know that

This is equivalent to say that \(H_3(\alpha _T/\beta _T) \ge 0\).

If \(H_3\) is an increasing function, since \(\omega _{\delta } < 0\), we will have

which gives (30). So, taking into account that

let us finish the proof by showing that \(H'_3(x) \ge 0\).

If \(N_2d_{min} > N_1d_{max},\) this is true straightforwardly, because we know that \(N_1d_{max} > K d_{min}+d^*+(N_1-K-1)d_{max}.\)

When \(N_2d_{min} \le N_1d_{max},\) we can argue as in the proof of Theorem 8, taking the point \(\tilde{d}= (d_1,...,d_{N_1}) = (N_2/N_1)(d_{min},...,d_{min}),\) that satisfies the bounds restrictions and

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'_3(x) \ge 0\). \(\square\)

Remark 7

1. (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 was described in [9], while considering nonnegative doses.

2. (b)

   For the hypofractionated case, the single exposure is chosen in [9] 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 difficulty is overcome here and we get the optimal number of dose fractions directly (and their values).

It is possible to show that all the above possibilities mentioned in Theorems 12 and 13 can appear in practice by means of examples. For the reader’s convenience, we have summarized the complete algorithm for the resolution of the problem \((P_2)\) in Table 3.

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)\). We have proved that they essentially depend on the sign of the quantity \(\omega _{\delta } = \alpha _T/\beta _T - \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 (among others [2, 6, 9]). Moreover, this is consistent with some clinical findings as noted in [9].

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 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. The study concerning the derivation of the optimal number of doses N had already been performed for example in [7] in the hyperfractionated case, but (as far as we know) it is completely new for the hypofractionated case.

Let us emphasize that the algorithms (described in Tables 2 and 3) can be implemented quite straightforwardly using any programming language to make them more accessible.