Optimization of spatiotemporally fractionated radiotherapy treatments with bounds on the achievable benefit

Spatiotemporal fractionation schemes deliver different dose distributions in different fractions in an attempt to minimize BED in healthy tissue and maximize BED in the tumor by hypofractionating parts of the tumor while delivering approximately identical doses to the surrounding tissue. The clinical rationale for spatiotemporal fractionation to optimally exploit the fractionation effect was first proposed in the context of proton radiotherapy (Unkelbach et al 2013b, Unkelbach and Papp 2015), where the potential benefit comes from the fact that the dose in the entrance region of a proton beam is largely independent of the beam's range, which provides some flexibility to modify the dose in the tumor without equally affecting the dose in the entrance region. More recently it has also been shown that spatiotemporal fractionation may provide a therapeutic advantage in arc therapy delivered with conventional photon beams (Unkelbach 2015). In this case, arc therapy plans can be created in such a way that each fraction delivers high single-fraction doses to complementary parts of the target volume while creating a similar dose bath in the surrounding normal tissue. This was demonstrated for fractionated radiosurgery treatments of large arteriovenous malformations (Unkelbach et al 2016).

The optimization of spatiotemporally fractionated treatments is far more challenging than conventional IMRT/IMPT optimization. For most commonly used objective functions, the fluence map optimization problem in conventional IMRT planning is a large-scale convex optimization problem that can be solved to global optimality using well-established gradient-based (e.g. quasi-Newton) optimization methods (Bortfeld 2006). Because spatiotemporal planning is based on BED rather than physical dose and the fluence maps for different fractions are distinct, the optimization models for optimal spatiotemporal fractionation are inherently nonconvex (see section 2.2). It is common to use heuristics and local optimization methods for nonconvex models in RT planning, but the global optimality of the resulting plans is rarely discussed. In a recent work, Ajdari and Ghate (2016) studied a treatment planning problem of similar mathematical structure. They proposed a model predictive control approach which computes a treatment plan iteratively by optimizing the dose distributions of remaining fractions after each fraction, with the assumption that the same dose is delivered in every fraction thereafter. This approach requires only the solution of convex optimization problems; however, it does not promote the type of solutions that we seek in spatiotemporal fractionation: plans that deliver high doses to complementary parts of the tumor in distinct fractions. Furthermore, the model predictive control approach comes with no guarantees that the computed plans are anywhere close to globally optimal.

So far there have been no computational tools to bound the maximum achievable benefit from spatiotemporal fractionation. Prior works (Unkelbach and Papp 2015, Unkelbach et al 2016) use local optimization, and it has been unclear if better optimization algorithms may yield substantially better solutions. Our approach combines the local optimization of the nonconvex treatment planning model for spatiotemporal fractionation with the solution of a convex optimization problem that provides a rigorous bound on the maximum achievable benefit from spatiotemporal fractionation. We approach this by formulating the spatiotemporal treatment planning problem as a nonconvex quadratically constrained quadratic programming (QCQP) problem. We then derive a convex relaxation of the QCQP problem in the form of a semidefinite programming (SDP) problem. The SDP relaxation provides a nontrivial lower bound; in particular, this bound is tighter than what can be achieved by replacing each nonconvex constraint in the QCQP model with its convex relaxation. We test our method on two-dimensional slices of 5 liver tumors that represent a variety of patient geometries. Comparing the quality of the locally optimal solutions against the SDP relaxation bounds, we find that the local optimal solutions computed for these cases are indeed nearly globally optimal. While there is no conceptual difficulty in applying the same approach to three-dimensional cases, and in particular the local optimal solutions can be computed using the same methods and software that we used in this paper, the SDP relaxations of the three-dimensional cases cannot be solved in a reasonable amount of time with available off-the-shelf software.

Prior research has addressed the problem of optimizing fractionation decisions based on the BED model (Mizuta et al 2012a, 2012b, Unkelbach et al 2013a, Keller et al 2013, Gay et al 2013, Saberian et al 2015, Saberian et al 2016). These works aim at maximizing the tumor BED subject to BED constraints to the normal tissue. It was shown that the optimal number of fractions depends not only on the ${\alpha}/{\beta}$ ratios of tumor and normal tissues, but also on the dose distribution. These works, however, all assume uniform fractionation, where the same dose distribution is delivered in all fractions, and only the number of fractions needs to be optimized.

The novelty in the idea of spatiotemporal fractionation lies in the fact that there is a potential advantage of delivering distinct dose distributions in different fractions, purely motivated by the basic fractionation effect as described by the standard BED model. There are several extensions of the BED model that describe higher order biological effects such as incomplete repair of radiation damage between fractions, repopulation of tumors over the course of treatment, accelerated repopulation effects, the effect of chemotherapeutic agents, and reoxygenation of hypoxic tumors (Hall and Giaccia 2012). It was found that some of these models give rise to more complex fractionation schemes, i.e. varying doses per fraction (Wein et al 2000, Yang and Xing 2005, Bertuzzi et al 2013, Bortfeld et al 2015, Salari et al 2015). However, the role of such models to guide fractionation decisions in clinical practice has been limited. Instead, spatiotemporal fractionation as described in this paper is purely based on the basic fractionation effect, whose existence and clinical relevance is undoubted.

Another approach in which different dose distributions may arise in each fraction is adaptive radiotherapy (ART) (Yan et al 1997, Lu et al 2008, Kim et al 2012). This technique uses feedback information, such as changes in the patient anatomy, to modify the treatment plan during the course of a fractionated treatment, e.g. replanning to compensate for tumor shrinkage in lung or head and neck cancer (Wu et al 2009, Sonke and Belderbos 2010). Instead, spatiotemporal fractionation shows that there is a benefit of delivering distinct dose distributions in different fractions, even in the absence of any changes of the patient over time.

The natural generalization of the BED model in (1) to fractionated treatments delivering different doses in different fractions is to define the (cumulative) BED as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ab}{\frac{\alpha}{\beta}} \displaystyle b = \sum_{t=1}^N \Big(d_t + \frac{d_t^2}{(\alpha/\beta)}\Big), \label{eq:BED} \nonumber \end{align} \tag{ 2 }$$

where $d_1, \dots, d_N$ are the doses delivered (in any order) in fractions 1 through N. IMRT planning using BED can be performed analogously to conventional IMRT optimization, using similar objective functions and constraints in the treatment planning optimization model, but substituting BED in place of physical dose (Unkelbach and Papp 2015). We define $x_1, \ldots, x_N$ to be the vectors of beamlet weights delivered in fractions 1 through N, and V to be the set of volume elements (voxels) used for dose calculation during the optimization. Using the cumulative BED from (2), we obtain the nonuniform fractionation problem:

$$\begin{align*} \displaystyle \label{eq:IMRT-BED} \min_{x,d,b}\quad & F(b) \nonumber\\ {\rm s.t.}\quad\; & b_v = \sum_{t=1}^N (d_{vt} + \frac{d_{vt}^2}{(\alpha/\beta)_v}) \qquad\qquad \forall\, {v\in V}\nonumber \\ & Dx_t = d_t \,\, \qquad\qquad\qquad\qquad\qquad t=1,\dots, N \nonumber \\ & x_t \geq 0 \qquad\qquad\qquad\qquad\qquad\,\,\,\,\,\,\,\, t=1,\dots, N,\nonumber \end{align*} \tag{ 3 }$$

where D is the usual dose-influence matrix, and the objective function F represents the desired clinical goals of target coverage, conformity, and organ sparing. Defining I as the index set of the objectives, we write the objective function as the sum

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:obj fun sum} F(b) = \sum_{i\in I} w_iF_i(b), \nonumber \end{align} \tag{ 4 }$$

where each term F_i represents a single clinical objective (such as deviation from a prescribed lower or upper bound on the minimum, maximum, or mean BED of a structure), and each positive weight w_i represents the relative importance of a clinical goal.

BED-based optimization can also be performed analogously for conventional, uniformly fractionated treatments. Throughout this work, we use a uniform reference plan as a benchmark to evaluate the benefit of nonuniform fractionation. We obtain the uniform reference plan by solving (3) with the additional constraint that $x_1=\cdots=x_N$. Eliminating the redundant variables, the uniform reference plan is the optimal solution of the following problem:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ab}{\frac{\alpha}{\beta}} \displaystyle \label{eq:IMRT-BED_uni} \min_{x,d,b}\quad & F(b)\nonumber \\ {\rm s.t.}\quad \; & b_v = Nd_v\left(1+\frac{d_v}{(\alpha/\beta)_v}\right) \qquad \forall\, {v\in V}\nonumber \\ & Dx = d \nonumber \\ & x\geqslant 0.\nonumber \end{align} \tag{ 5 }$$

Figure 1 shows an example of a reference plan and a nonuniform 5-fraction treatment of a patient harboring a large liver tumor. The reference plan (figure 1(c)) was computed by solving the BED-based treatment plan optimization model (5). Figure 1(a) shows the 5 dose distributions of a nonuniformly fractionated plan. Both plans were computed by optimizing the same objectives described below in section 3 to represent identical clinical goals. The reference plan yields a conventional treatment that irradiates the tumor to the prescribed dose in each fraction. The nonuniformly fractionated plan delivers high single-fraction doses to parts of the tumor while delivering a similar low-dose bath to the surrounding tissue in each fraction. The high single-fraction doses in the tumor allow for a reduction in the total physical dose delivered (figure 1(b)). Since the dose to the surrounding normal tissue is approximately uniformly fractionated, this also yields a reduction of BED in the normal tissue.

Zoom In Zoom Out Reset image size

It is well-known that conventional fluence map optimization is a convex optimization problem when the clinical goals are modeled using a convex objective function F of the physical dose. Despite the apparent nonconvexity introduced by the quadratic equality constraints in the uniform fractionation problem (5), this model is convex when the physical doses d are restricted to clinically relevant values and F is a piecewise quadratic penalty function similar to the dose-based objective functions (Unkelbach and Papp 2015). This is the case for the objectives used in this paper. Hence, the uniform reference plan can be computed using the gradient-based local optimization methods commonly used in IMRT plan optimization.

The nonuniform fractionation problem (3), however, is nonconvex ⁴ . This nonconvexity is an inherent characteristic of the problem that cannot be eliminated by reformulating the model, as the following argument shows: by definition, permuting the fractions in an optimal treatment plan leads to another optimal plan, yet the average of these $N!$ treatment plans will be a plan with identical fractions, which in general will not be optimal. Thus, gradient-based optimization methods used to solve the uniform model can only yield locally optimal solutions for the nonuniform model. Figure 2 illustrates this problem, showing that substantially different locally optimal spatiotemporal plans may exist for the same case.

Zoom In Zoom Out Reset image size

A drawback of formulation (3)–(4) is that the optimal solution does not realize the maximum benefit of spatiotemporal fractionation over conventional fractionation in any given objective. Instead, the benefit is distributed among the terms of the objective function; in other words, the benefit is realized as a combination of smaller improvements with respect to the different clinical goals. In the liver cases we study in section 4, the primary gain is expected to be in lowering BED to the liver without compromising target coverage and conformity. This suggests an alternative formulation in which we minimize only one of the clinical objectives of (4), say F₁, while constraining the solution to be at least as good as the uniform reference plan with respect to each other objective F_i ($i\in I, i\neq 1$).

Let $b^*$ be the BED distribution delivered by the uniform reference plan, which is obtained by solving (5), and let F₁ be the penalty for the mean BED in the liver exceeding zero. We modify (3) to include the constraints that the solution must be as good as the uniform reference plan with respect to all objectives besides F₁, and we obtain the constrained nonuniform spatiotemporal fractionation problem:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ab}{\frac{\alpha}{\beta}} \displaystyle \label{eq:IMRT-BED_constr} \min_{x,d,b}\quad & F_1(b) \nonumber \\ {\rm s.t.}\quad \; & F_i(b) \leqslant F_i(b^*) \quad\quad\quad\quad\quad\quad\quad\quad i\in I,\; i\neq 1 \nonumber \\ & b_v = \sum_{t=1}^N (d_{vt} + \frac{d_{vt}^2}{(\alpha/\beta)_v}) \quad\quad\quad\quad \forall\,v\in V\nonumber \\ & Dx_t = d_t \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad t=1,\dots, N \nonumber \\ & x_t \geqslant 0 \quad \quad\quad\quad\quad\quad \quad\quad\quad\quad\quad t=1,\dots, N.\nonumber \end{align} \tag{ 6 }$$

The first set of constraints ensures that the improvement in the objective F₁ is not at the cost of sacrificing the other clinical objectives; the computed spatiotemporal plan is either preferable or identical to the uniform reference plan with respect to every objective.

The convex relaxation of the constrained nonuniform model (6) can be derived analogously to how we obtained the relaxation (9) of the model (3). For simplicity of notation, we assume (consistently with the application to liver tumors) that the mean dose objective $F_1(b)=(\frac{1}{|V_1|}\sum_{v\in V_1} b_v - m_1^{hi})_{+}^2$ is the primary objective. Then the convex optimization model bounding the minimum value of $F_1(b)$ from below is the following:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:SDPrelaxation-constrained} \min_{x,X,p,q,r} \quad & r_1 \nonumber \\ {\rm s.t.} \quad \;\;\; & \sum_{v\in V_i} p_{iv}^2 \leqslant F_i(b^*) \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \forall i \in I^+\nonumber \\ & \sum_{v\in V_i} q_{iv}^2 \leqslant F_i(b^*) \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \forall i \in I^-\nonumber \\ & r_i \leqslant F_i(b^*) \,\,\,\, \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \!\!\! \forall i \in I^m,\; i \neq 1\nonumber \\ & p_{iv} \geqslant N \left\langle \left(\begin{array}{@{}cc@{}} 1 & x^T \nonumber \\ x & X \end{array}\right), \left(\begin{array}{@{}cc@{}} -\frac{b_{iv}^{hi}}{N} & \frac{e_v^TD}{2} \nonumber \\ \frac{D^Te_v}{2} & C_v\end{array}\right)\right\rangle \,\,\quad\quad\quad \quad\quad\quad\quad\quad \forall i \in I^+,\; \forall v\in V_i\nonumber \\ &p_{iv} \geqslant 0 \,\,\,\,\, \quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \forall i \in I^+,\; \forall v\in V_i\nonumber \\ &q_{iv} \geqslant -N \left\langle \left(\begin{array}{@{}cc@{}} 1 & x^T \nonumber \\ x & X \end{array}\right), \left(\begin{array}{@{}cc@{}} -\frac{b_{iv}^{lo}}{N} & \frac{e_v^TD}{2} \nonumber \\ \frac{D^Te_v}{2} & C_v\end{array}\right)\right\rangle \quad\quad\quad\quad\,\, \quad\quad\quad \forall i \in I^-,\; \forall v\in V_i\nonumber \\ &q_{iv} \geqslant 0 \!\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \forall i \in I^-,\; \forall v\in V_i \nonumber \\ & r_i \geqslant \frac{N}{|V_i|} \sum_{v\in V_i} \left\langle \left(\begin{array}{@{}cc@{}} 1 & x^T \nonumber \\ x & X \end{array}\right), \left(\begin{array}{@{}cc@{}} -\frac{m_i^{hi}}{N} & \frac{e_v^TD}{2} \nonumber \\ \frac{D^Te_v}{2} & C_v\end{array}\right)\right\rangle \!\quad \quad\quad\quad\quad\quad \forall i\in I^m \nonumber \\ & r_i \geqslant 0 \!\!\!\!\!\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \forall i\in I^m \nonumber \\ & \left(\begin{array}{@{}cc@{}} 1 & x^T \nonumber \\ x & X \end{array}\right) \succcurlyeq 0,\;\; x \geqslant 0,\;\; X \geqslant 0.\nonumber \end{align} \tag{ 10 }$$

Although the optimization problems (9) and (10) are convex, their solution is not straightforward using the gradient-based algorithms commonly used in treatment planning optimization, such as the projected gradient descent method and its variants. This is because computing the violation of the linear matrix inequality constraints $\newcommand{\e}{{\rm e}} \left(\begin{array}{@{}cc@{}} 1 & x^T \\ x & X \end{array}\right) \succcurlyeq 0$ and projecting on the set of points satisfying this inequality require an expensive matrix factorization in each step of the method. Optimization problems involving this type of constraints are called semidefinite programs (Vandenberghe and Boyd 1996, Boyd and Vandenberghe 2004). They have been extensively studied in the numerical optimization literature, and several reliable and efficient algorithms are available for their solution.

Formally, a semidefinite program is a convex optimization problem with a linear objective function and with constraints that are all either linear inequalities or linear matrix inequalities of the form $A(x)\succcurlyeq 0$, where $A(\cdot)$ is an affine function that maps the optimization variables x to the space of real symmetric matrices. The problems (9) and (10) have convex quadratic objectives and constraints, rather than only linear ones, but it can be shown that convex quadratic inequalities can be equivalently written as linear matrix inequalities (Vandenberghe and Boyd 1996); hence their inclusion is without loss of generality.

Large-scale semidefinite programs are routinely solved using Newton-type methods implemented in widely available optimization software such as SeDuMi (Sturm 1999) and MOSEK (MOSEK ApS 2017). In our experiments we used MOSEK to solve the semidefinite programs (9) and (10).

4.1.1. The BED prescription and $\newcommand{\ab}{{\alpha}/{\beta}} \ab$ ratios

4.1.2. Computing the uniformly fractionated reference plan

4.1.3. Computing the spatiotemporally fractionated plans

Figures 1 and 3 show a graphical comparison of the uniform reference plan with locally optimal nonuniformly fractionated solutions for Cases 1 and 2; the graphical results for the remaining three cases can be found in the supplement (stacks.iop.org/PMB/63/015036/mmedia). Table 1 summarizes the computed bounds and benefits observed in all five cases. In each case, the nonuniform plan is constrained to maintain the same target coverage as the uniform reference plan, yet the mean liver BED in the nonuniform plans is substantially lower than in the uniform plans (see table 1). Case 1 displays the smallest reduction in the mean liver BED with an approximately 13% improvement, and Case 2 has the largest reduction of approximately 34% in the mean liver BED. Note that these improvements are achieved without sacrificing any other clinical objective, as by definition, the spatiotemporal plans computed by solving (6) are at least as good as the uniform reference plan with respect to every objective. In particular, the spatiotemporal plans are as conformal as the corresponding uniform reference plans.

Zoom In Zoom Out Reset image size

The effectiveness of a spatiotemporal treatment can be quantified by its equivalent dose DEQ5, defined via

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\ab}{\frac{\alpha}{\beta}} \displaystyle {\rm BED} = \sum_{t=1}^5 \left(d_t + \frac{d_t^2}{\alpha/\beta} \right) = {\rm DEQ5} + \frac{({\rm DEQ5})^2/5}{\alpha/\beta}, \nonumber \end{align*}$$

which yields

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\ab}{\frac{\alpha}{\beta}} \displaystyle {\rm DEQ5} = 5 \left( -\frac{\alpha/\beta}{2} + \sqrt{\left( \frac{\alpha/\beta}{2} \right)^2 + \frac{\alpha/\beta}{5} \sum_{t=1}^5 \left(d_t + \frac{d_t^2}{\alpha/\beta} \right)}\, \right). \nonumber \end{align*}$$

The DEQ5 is the cumulative dose distribution of a uniformly fractionated 5-fraction treatment that achieves the same BED as the given spatiotemporal treatment. (Note the difference between DEQ5 and the also commonly used equieffective dose, EQX, which is the total physical dose that needs to be delivered in a uniform treatment with a dose per fraction of X Gy to achieve the same BED.) Using DEQ5 has the advantage that the spatiotemporal plan can be directly compared to the physical dose distribution of the uniform reference plan. The comparison is shown in panels (c)–(e) of figures 1 and 3 for Cases 1 and 2, and in the Supplement for the remaining cases. By computing the difference between the two plans, we observe that the nonuniformly fractionated plans maintain the same tumor BED as the uniformly fractionated plans while reducing the mean dose and BED in healthy liver tissue. The same can also be demonstrated by comparing the DVH curves of the physical dose and DEQ5 of the two plans. For Case 1, these curves are included in the supplement (figure 7). The DVH curves show that there is substantial dose reduction in the liver (both the DEQ5 and physical dose curves shift to the left), and that this is a consequence of lower physical dose (but nearly identical DEQ5) delivered to the tumor.

4.2.1. Bounds on the maximum achievable liver BED reduction

Because the nonuniform fractionation problem is nonconvex, the nonuniformly fractionated solutions presented in this paper are only certified to be locally optimal. To gauge the quality of these solutions, for each case we computed a lower bound on the minimum liver BED, as described in section 3.1. These lower bounds, also shown in table 1, are close to the mean liver BED values achieved by the locally optimal nonuniform ly fractionated solutions. For example, in Case 1, spatiotemporal fractionation reduced the mean liver BED to 75.9 Gy, which is a reduction of 8.7 Gy from the uniform reference plan. The lower bound from the SDP relaxation is 73.4 Gy, which means that no treatment plan can achieve a reduction of more than 11.2 Gy from the uniform reference plan. It is important to note that the value of 73.4 Gy is only a lower bound, and we are not able to guarantee the existence of a nonuniform ly fractionated solution that achieves this value. The mean liver BED value for the true globally optimal plan may be anywhere between 75.9 and 73.4 Gy.

We compared the mean liver BED reduction in the locally optimal nonuniform ly fractionated plan with our bound on the maximum possible reduction by computing the fraction of the gap between the mean liver BED in the uniform reference plan and the lower bound from the SDP relaxation that the nonuniform ly fractionated plan closes. Formally, if m_ref and m_sp are the mean liver BEDs achieved by the uniform reference plan and the locally optimal spatiotemporal plan, respectively, and m_SDP is the lower bound obtained from the SDP relaxation, then the fraction of the gap closed by the spatiotemporal plan is the ratio

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:gap} \frac{m_{\rm ref}-m_{\rm sp}}{m_{\rm ref}-m_{\rm SDP}}. \nonumber \end{align} \tag{ 11 }$$

In particular, a ratio of $100\%$ is the best possible value, which is attained if the spatiotemporal plan matches the SDP lower bound, and $0\%$ is the worst possible value, attained when the uniform plan cannot be improved upon using spatiotemporal fractionation. Note that this ratio being close to $100\%$ does not only speak directly for the quality of the spatiotemporal plan; it also demonstrates the sharpness of the lower bound. For example, in Case 1, the nonuniform plan's reduction of 8.7 Gy is 78% of the upper bound on the maximum possible BED reduction, which is 11.2 Gy. The nonuniformly fractionated plans closed 78–97% of the gap between the mean liver BED of the uniform reference plan and our lower bound on the lowest achievable mean liver BED.

4.2.2. The optimal number of fractions

The results indicate that spatiotemporal treatments can achieve substantial reductions in both BED and physical dose to the liver and other normal tissue. The approximately 15–20% mean liver BED reduction is consistent with the benefit of spatiotemporal fractionation observed in previous work for brain lesions (Unkelbach et al 2016). Our work further shows that local optimization techniques provide high-quality plans that are close to realizing the maximum potential normal tissue dose reduction.

The results in this paper are limited to the optimization of treatment plans for two-dimensional slices of clinical liver cases in order to demonstrate the concepts introduced in the paper, but in a practical setting the computations would have to be carried out for the entire patient. There is no additional difficulty in computing locally optimal nonuniform solutions for three-dimensional cases with the same algorithms that we used in our study (Unkelbach et al 2016). Therefore, the main challenge in extending this study to three-dimensional cases is the computational complexity of solving the SDP relaxations in section 3. In the nonuniform fractionation optimization problem there is at least one variable for each beamlet in the fluence maps and at least one constraint for every voxel associated with a quadratic penalty function. In the SDP relaxation, with the introduction of the matrix variables X, the number of variables is roughly squared. SDP problems of this size cannot be solved for three-dimensional cases in a reasonable amount of time with the available off-the-shelf software. Even with the two-dimensional case, the solver took up to 7 hours to solve the SDP relaxation. It is an interesting direction for future research to devise customized algorithms for the solution of these large-scale SDP problems that arise for the three-dimensional cases. For now, the results obtained for the two-dimensional slices suggest that the locally optimal solutions computed for the nonuniform fractionation problem are close to global optimality.

As mentioned in section 4, the treatment plans in this work have a higher number of beams than a conventional IMRT plan to approximate a VMAT plan, where we expect nonuniform fractionation to display the most benefit. The nonuniformly fractionated plans that are most effective in lowering BED and total physical dose are those in which the tumor is hypofractionated while the surrounding tissue is fractionated. In other words, these plans deliver a high single-fraction dose to parts of the tumor during each fraction and a consistent lower dose to the liver and other healthy tissue throughout all of the fractions. VMAT treatments are particularly suitable for delivering such a nonuniformly fractionated plan, thanks to their characteristic low-dose bath delivered to the healthy tissue.

Further research is needed to investigate the impact of various sources of uncertainty on spatiotemporally fractionated treatments. We anticipate that the effect of small changes in the radiobiological parameters is negligible, as they only slightly affect what fractionated treatments are isoeffective. Interfractional patient motion and soft tissue deformation is a greater concern for spatiotemporal treatments in which each fraction delivers precariously aligned dose distributions with sharp dose gradients in the interior of the target structure. Future research will consider including additional safeguards such as bounds on the dose gradients as well as stochastic and robust optimization methods to explicitly incorporate uncertainty into the planning. We also note that spatiotemporal fractionation can be applied to treatments with fewer distinct fractions that are delivered a number of times each; for example, a 15-fraction scheme can be designed (entirely analogously to our plans) that consist of only 3 distinct fractions that are delivered 5 times each. The repeated delivery of identical fractions may mitigate somewhat the effects of random interfraction motion. Finally, other technologies, currently under development, may also make the safe delivery of spatiotemporal liver SBRT treatments feasible in the future. These include MLC and couch tracking, and image-guided radiotherapy using an MR-linac.