On the interplay between robustness and dynamic planning for adaptive radiation therapy^*

Throughout this work T denotes the total number of fractions t, while N denotes the total number of voxels n. The accumulated dose at fraction t is denoted by ${x}_{t}\in {{\mathbb{R}}}^{N}$ which in our framework is referred to as the system at the sampling instant t. Since it is assumed that the patient does not receive any radiation before the start of the treatment at t = 0, the dose x₀ is set to zero. The radiation plan to be applied at t + 1 is given by ${u}_{t}\in {{\mathbb{R}}}^{M}$, where M refers to the number of bixels. In our framework the plans are referred to as control signals. Dose deposition is modeled by a dose deposition matrix $B\in {{\mathbb{R}}}^{N\times M}$. In a deterministic setting, the accumulated dose would be given by x_t+1 = x_t + Bu_t for t = 0,1, ..., T − 1. However, the focus of our framework is on handling interfractional uncertainties. In this work, the interfractional variations ω are assumed to vary from fraction to fraction and thus, modeled as independent and identically distributed (i.i.d.) random variables. The uncertainties are derived from a discretized normal distribution ${ \mathcal N }(0,{\sigma }^{2})$ and summarized in the set Ω. The probability of each i.i.d. uncertainty scenario ω ∈ Ω to occur is denoted by p_ω. The geometric uncertainty occurring at fraction t is modeled by the stochastic shift matrix $S({\omega }_{t})\in {{\mathbb{R}}}^{M\times M}$. In our framework, the impact of interfractional uncertainties on the fractionated treatment is modeled by the linear, stochastic, discrete-time state-space dynamics

$$\begin{eqnarray}&&{x}_{t+1}={x}_{t}+{BS}({\omega }_{t+1}){u}_{t}\quad \mathrm{for}\ t=0,1,\,\ldots ,\,T-1.\end{eqnarray} \tag{ 1 }$$

Our proposed framework for robust adaptive radiation therapy uses dynamic optimization variables in order to deliver the prescribed dose distribution d_T in the presence of uncertainties. In this study, the goal is to achieve a uniform dose in CTV and zero elsewhere. In order to minimize the gap between the accumulated dose x_T and the prescribed dose d_T, we choose the quadratic penalty ${\parallel {x}_{T}-{d}_{T}\parallel }^{2}$, i.e. $f\left({\parallel {x}_{T}-{d}_{T}\parallel }^{2}\right)$ as a planning objective which is widely used in clinical planning systems, where f refers to the particular robust optimization model. The focus of this study is on the mathematical properties and therefore, we follow Unkelbach and Oelfke [14] and Fredriksson [15] among others [13, 16, 17, 23] who used the quadratic penalty in their mathematical studies. The quadratic penalty has the advantage of being convex and thus, the generated plans are globally optimal solutions, which allows for interpreting the results in a straightforward fashion. A variety of robust adaptive strategies are introduced with the purpose to better understand how to handle interfractional uncertainties. It is the goal of every robust adaptive strategy to generate an individual plan u_t at fraction t which is going to be delivered at fraction t + 1, equivalent to the concept of offline adaptive radiation therapy is applied. The robust adaptive strategies are categorized based on conservativeness, the degrees of freedom of their optimization variables and their timing for treatment planning. In order to clearly distinguish between the time instant of plan delivery t and planning, we introduce the index t_P for the latter. The optimization variables u_t, i.e. plans are either i) time-and-scenario-independent, ii) time-dependent or iii) time-and-scenario-dependent. Time- and scenario-independent optimization variables are used in the non-adaptive robust strategy, which refers to applying the same plan throughout the course of treatment. Time-dependent optimization variables are the basis for the time-varying robust and time-varying robust adaptive strategy. Time- and scenario-dependent optimization variables are used in the time- and scenario-dependent robust adaptive strategy and the robust adaptive strategies based on model-predictive control (MPC), referred to as the MPC-strategies. The use of time- and scenario-dependent optimization variables provides the most degrees of freedom to mitigate the effect of interfractional uncertainties on the accumulated dose, while time- and scenario-independent optimization variables provide the least. The performance of the robust adaptive strategies is analysed in comparison with the corresponding non-adaptive robust strategy. By design, the strategies with time-and-scenario-dependent optimization variables allows for temporal and spatial variations in the fluence profiles, which is the case in the time-and-scenario-dependent robust adaptive and MPC-strategies. The fluence profiles are controlled through constraints in which the daily fraction dose may vary within $(1-\gamma )\tfrac{{d}_{T}}{T}\leqslant {{Bu}}_{t}\leqslant (1+\gamma )\tfrac{{d}_{T}}{T}$, where $\gamma \geqslant 0$ specifies the user-defined lower and upper thresholds for the daily fluence u_t. We introduce these constraints and perform simulations with and without them in order to evaluate the effect on the resulting quality of the solutions. The robust adaptive strategies using these constraints are referred to as stabilized time-and-scenario-dependent and MPC-strategies. As previously discussed, it is the overall goal of all robust strategies introduced in this framework to minimize the quadratic penalty of the gap between the final accumulated dose ${x}_{T}={x}_{0}+{\sum }_{\tau =0}^{T-1}{BS}({\omega }_{\tau +1}){u}_{\tau }$ and the prescription dose d_T, but they all differ in their approach to generate the plans u_t to achieve that goal. Ultimately, it is the objective of our study to analyse the mathematical properties and performance of the proposed robust adaptive strategies to solve

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\min }\limits_{{u}_{\tau }\geqslant 0}f\left({\parallel {x}_{{t}_{p}}+\displaystyle \sum _{\tau ={t}_{P}}^{T-1}{BS}({\omega }_{\tau +1}){u}_{\tau }-{d}_{T}\parallel }^{2}\right),\\ {\rm{s}}.{\rm{t}}.\ {x}_{\tau +1}={x}_{\tau }+{BS}({\omega }_{\tau +1}){u}_{\tau }\\ \mathrm{for}\ \tau ={t}_{P},{t}_{P}+1,\ldots ,T-1;{x}_{{t}_{P}}\ \mathrm{given},\end{array}\end{eqnarray} \tag{ 2 }$$

where the current system state ${x}_{{t}_{p}}$ is assumed to be always known and f refers to the specific robust model that is combined with a specific treatment strategy. In particular, we are interested in analysing the performance of combining expected-value-, worst-case- or conditional-value-at-risk (CVaR)-optimization with non-adaptive and adaptive strategies.

2.2.1. Non-adaptive robust strategy

2.2.2. Time-varying robust strategy

In contrast to the non-adaptive robust strategy, the plans in the time-varying robust strategy may vary from fraction to fraction and therefore, the degrees of freedom of the optimization variables are increased to include time-dependency, but information on the current accumulated dose is not taken into account either. Thus, problem (2) is solved before the start of the treatment at t_P = 0 for positive time-dependent optimization variables ${u}_{\tau }\geqslant 0\in {{\mathbb{R}}}^{M}\ \forall \tau$. With the time-varying robust strategy, we want to investigate if adding time-dependency as a degree of freedom to the time-and-scenario-independent optimization variables, used in the non-adaptive robust strategy, may give any advantages. In fact, the move from time- and scenario-independent to time-dependent variables is not necessary, if (i) the final accumulated dose is evaluated in a quadratic cost function, and if (ii) the uncertainties that occur at every fraction are assumed to follow the same distribution. Since conditions (i) and (ii) imply that delivering the mean over all time-dependent plans at every fraction will be at least as good delivering the time-dependent plans at the corresponding fraction. In mathematical terms, the choice of the quadratic cost function in our model gives convexity, while the evaluation of the final accumulated dose affected by i.i.d. uncertainties represents the permutation-invariance. The mathematical evaluation of such a case has been analyzed by Parillo [24], p. 100, as reviewed in the following proposition. (Parrilo [24], p. 100).

Proposition 1 Let $f:\left[U\right]\to {\mathbb{R}}$ be a convex and permutation invariant function of a convex problem, where $\left[U\right]$ refers to the T-tuple $({u}_{0},{u}_{1},\ldots ,{u}_{T-1})$ of vectors ${u}_{t}\in {{\mathbb{R}}}^{N}$ for t = 0,1, ..., T − 1 . Then,

$$\begin{eqnarray*}&&f\left(\left[\bar{U}\right]\right)\leqslant f\left(\left[U\right]\right),\end{eqnarray*}$$

with $\left[\bar{U}\right]$ denoting the T-tuple of the mean $\tfrac{1}{T}{\sum }_{t=0}^{T-1}{u}_{t}$ over all u_t for t = 0,1, ..., T − 1.

For the sake of completeness, the complete proof is given in appendix B. This result also indicates that under the conditions of convexity and permutation invariance, the conventional approach to treatment planning implicitly takes fractionation into account. Therefore, evaluation and analysis of the final dose distribution for the time-varying robust strategy are omitted.

2.2.3. Time-varying robust adaptive strategy

2.2.4. Time-and-scenario-dependent robust adaptive strategy

2.2.5. MPC-strategies

In section 2.2, we introduced various strategies to solve (2) in which f refers to a specific robust model. In this section, we will discuss the nature of the robust models chosen to be combined with the non-adaptive and adaptive strategies. expected-value optimization is the least conservative model, while worst-case-optimization the most conservative; and conditional-value-at-risk (CVaR)-optimization represents a compromise between the two. Thus, $f\left({\parallel {x}_{T}-{d}_{T}\parallel }^{2}\right)$ in (2) takes the following form in

• 1.  
  Expected-value-optimization

  $$\begin{eqnarray}&&{\mathbb{E}}\left[{\parallel {x}_{T}-{d}_{T}\parallel }_{D}^{2}\right],\end{eqnarray} \tag{ 3 }$$
• 2.  
  Worst-case-optimization

  $$\begin{eqnarray}&&\mathop{\max }\limits_{{\tilde{\omega }}_{1}^{T}\in {\tilde{{\rm{\Omega }}}}_{1}^{T}}{\parallel {x}_{T}-{d}_{T}\parallel }_{D}^{2}\end{eqnarray} \tag{ 4 }$$

  and
• 3.  
  Conditional-value-at-risk (CVaR)-optimization

  $$\begin{eqnarray}&&\mathop{\min }\limits_{\lambda }\ \lambda +\displaystyle \frac{1}{\alpha }{\mathbb{E}}\left[{\left({\parallel {x}_{T}-{d}_{T}\parallel }_{D}^{2}-\lambda \right)}_{+}\right],\end{eqnarray} \tag{ 5 }$$

where the shorthand z₊ and λ denote $\max \left\{z,0\right\}$ and an additional optimization variable, respectively. The matrix D in (3), (4) and (5) is referred to as the weighting matrix which will be described in more detail later. In order to avoid discontinuities in the worst-case- and CVaR- optimization models, the objective functions in (4) and (5) are reformulated as

$$\begin{eqnarray}&&\mathop{\min }\limits_{z\geqslant 0}\ \{z:\ z\geqslant {\parallel {x}_{T}-{d}_{T}\parallel }_{D}^{2}\ \forall \ {\tilde{\omega }}_{1}^{T}\in {\tilde{{\rm{\Omega }}}}_{1}^{T}\}\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\min }\limits_{\lambda ,\zeta }\ \left\{\lambda +\displaystyle \frac{1}{\alpha }{p}_{\tilde{\omega }}^{T}\zeta :\ \lambda +{\zeta }_{\tilde{\omega }}\right.\\ \ \ \left.\geqslant \,{\parallel {x}_{T}-{d}_{T}\parallel }_{D}^{2}\ \forall {\tilde{\omega }}_{1}^{T}\in {\tilde{{\rm{\Omega }}}}_{1}^{T},\zeta \geqslant 0\right\},\end{array}\end{eqnarray} \tag{ 7 }$$

the max function is substituted, the help variable ζ is introduced and where ${p}_{\tilde{\omega }}={p}_{\omega }\times {p}_{\omega }\times \cdots \times {p}_{\omega }$ refers to the probability vector for the realizations ${\tilde{\omega }}_{1}^{T}$. Since the number of possible realizations of uncertainty scenarios grows exponentially with the number of fractions T, problem reformulations may be considered. In case of minimizing the expected-value of the quadratic penalty (3), the optimal solution can be computed in an exact manner independent of T by exploiting the structure of the quadratic penalty, which is demonstrated in detail in appendix A. In worst-case- and CVaR-optimization, an exact optimal solution can be computed depending on computational power and for a much shorter than conventionally used treatment length T. In cases such as T = 30, the solution could be optimized based on Monte Carlo simulated treatments subjected to interfractional uncertainties ω ∈ Ω. CVaR-optimization can be interpreted as a compromise between expected-value- and worst-case scenario optimization, since its optimal solution minimizes the expected-value of the fraction 0 < α ≤ 1 of the worst scenarios ω ∈ Ω [15]. Thus, the parameter α which affects the upper bound of p_ω controls the grade of conservativeness. The closer α is to 1, the less conservative the solution will be, while $\alpha \leqslant {\min }_{\omega \in {\rm{\Omega }}}{p}_{\omega }$ with p_ω > 0 changes (5) to be equivalent to worst-case-optimization in (4). In this study, α is set to 0.4. The concept of CVaR-optimization was introduced by Rockafellar and Uryasev [28] and it is widely used in finance for risk management. With reference to robust planning, Fredriksson [15] introduced a framework for robust planning using CVaR to linearly adjust the grade of conservativeness.

For the sake of readability, we introduce the Euclidean norm with respect to the weighting matrix D in (3), (4) and (5) in order to simplify the extended quadratic penalty function

$$\begin{eqnarray}&&\left[\displaystyle \sum _{r\in { \mathcal R }}{g}_{r}\displaystyle \sum _{n\in {{ \mathcal N }}_{r}}{\eta }_{n,r}{v}_{n,r}{\left({x}_{T,n}-{d}_{T,n}\right)}^{2}\right],\end{eqnarray} \tag{ 8 }$$

which combines the conventionally used importance weights g_r, the relative volumes v_n,r satisfying ${\sum }_{n\in {{ \mathcal N }}_{r}}{v}_{n,r}=1$ for each region of interest (ROI) $r\in { \mathcal R }$ denoting the set of all regions and the scenario-based voxel weights η_n,r. These scenario-based voxel weights η_n,r are a result of a generalized scenario-based robust optimization approach to account for interfractional uncertainties introduced by Fredriksson and Bokrantz [16], which omits the use of CTV-PTV margins. These voxel weights are derived from modeling an uncertainty scenario ξ ∈ Ξ as whole-body shifts m(ξ) and summarizing the voxels part of the CTV under scenario ξ in the set ${ \mathcal C }^{\prime} (\omega )=\left\{n+m(\xi ):n\in { \mathcal C }\right\}$, with ${ \mathcal C }$ referring to the voxels contained in the CTV under the scenario of no shift occurring. Thus, the scenario based voxel weights are equivalent to the reciprocal number of scenarios ξ under which voxel n is contained in the scenario-related set ${ \mathcal C }^{\prime} (\xi )$ as expressed by

$$\begin{eqnarray}&&{\eta }_{n}=\displaystyle \frac{1}{| \left\{\xi \in {\rm{\Xi }}:n\in { \mathcal C }^{\prime} (\xi )\right\}| }.\end{eqnarray} \tag{ 9 }$$

The operator $| \cdot |$ in (9) refers to the cardinality of the set $\left\{\xi \in {\rm{\Xi }}:n\in { \mathcal C }^{\prime} (\xi )\right\}$ and the scenario-based voxel weights η_n satisfy the relation $1/| {\rm{\Xi }}| \leqslant {\eta }_{n}\leqslant 1$ for every voxel n.

In order to focus on the mathematical properties of the proposed framework, a one-dimensional patient phantom, as illustrated in figure 1, is considered. This one-dimensional model schematically represents a slice of a two-dimensional phantom or an intersection of a sagittal and transversal cut of a three-dimensional patient-geometry. The one-dimensional phantom geometry is discretized into 40 voxels. The phantom contains one clinical target volume (CTV) and two organs at risk (OARs) which are asymmetrically located around the CTV. The CTV is located between -1.2 and 1.2 cm, while the organs at risk are located in the intervals [ − 2, −2.2] cm and [2,3] cm. The CTV is an extension of the visible tumor and accounts for the microscopic spread of cancer cells. As a consequence, dose coverage of the CTV during the whole course of treatment has to be maintained despite the presence of interfractional variations in order to guarantee a successful treatment outcome. According to our model, the CTV is irradiated by a perpendicular oriented field, while the absorbed dose in each voxel is modeled by Gaussian functions at a spacing of 1.5 mm with a standard deviation of 3 mm, which is represented by the dose-deposition matrix B. As a consequence of the idealized phantom geometry, the number of bixels M and their positions are assumed to be identical with the number of voxels N and their location. Thus, the dimensions of the plans, the shift-matrix and dose-deposition matrix are simplified to ${u}_{t}\in {{\mathbb{R}}}^{M}$, $S(\omega )\in {{\mathbb{R}}}^{N\times M}$ and $B\in {{\mathbb{R}}}^{M\times M}$, respectively. In consideration of reducing computational effort, interfractional geometric uncertainties are modeled as rigid whole-body shifts such that the dose received by each voxel can be computed and tracked in a straightforward manner. Moreover, the link between parameter changes in robust adaptive strategies and features in the accumulated dose may be clearer than for two- or three-dimensional geometries. Organ-deformations, however important, are not considered in the proposed framework, such that finding a direct relation between parameter changes in the model and the resulting physical dose may not be compromised.

Zoom In Zoom Out Reset image size

The computational study is performed in MATLAB version 9.1 using the IBM optimization solver CPLEX in the studio version 12.6.3. As previously discussed, we assume that the one-dimensional phantom geometry is subjected to rigid whole-body shifts, which can be interpreted as movement in the anterior-posterior direction. These random shifts ω ∈ Ω are derived from the discretized normal distribution ${ \mathcal N }(0,0.25)$, which is in accordance with practical studies on interfractional prostate motion [29–31]. The discretized scenarios ω ∈ Ω are obtained by adapting the normal distribution to the voxel grid according to its mean and standard deviation [17]. The set Ω = {− 2, −1, 0, 1, 2} which is included in the robust optimization models contains the uncertainty scenarios measured in voxel shifts with the corresponding normalized probabilities p_ω = {0.092 4, 0.241 4, 0.332 4, 0.241 4, 0.092 4}, Throughout all optimization models and treatment strategies the optimization weight g_CTV for the CTV is set to 100, while the weight for the right OAR and left OAR is set to 10 and the external is set to 1. The study of the proposed framework is conducted in two stages.

First, the feasibility study is carried out in order to study the trade-off between mathematical optimality and computational feasibility of our proposed framework and the feasibility of adaptive strategies in hypofractionated radiotherapy. The trade-off is evaluated by comparing the objective function values at the optimal solution of every strategy. Among the presented strategies, the non-adaptive robust strategy is considered the lower benchmark and the time-and-scenario-dependent strategy the upper benchmark, since it takes into account uncertainties over the whole treatment horizon. In the mathematical study, the impact of the choice of optimization variables, i.e.: time-dependent or time-and-scenario-dependent, on the optimal solution is examined. For the sake of computational tractability, the total number of fractions T is set to five and the prediction horizon K is set to one, two and three fractions. As a consequence, the evaluation of the time-varying robust adaptive and MPC-strategies relative to the robust non-adaptive and time-and-scenario-dependent robust adaptive strategy is carried out in a fair manner.

Second, in the conventional fractionation length study we study the performance of our proposed strategies for a population of 100 treatments, each consisting of 30 fractions with interfractional uncertainties. These uncertainty scenarios are generated from the same normal distribution ${ \mathcal N }(0,0.25)$ as those accounted for during the optimization process, but truncated in a different manner in order to test our strategies for a larger set of uncertainty scenarios {− 3, −2, −1, 0, 1, 2, 3}. For the sake of a fair comparison, the robust non-adaptive, time-varying robust adaptive and MPC-strategies are evaluated for the same treatment population. The plans for the robust non-adaptive and time-varying robust adaptive strategy combined with worst-case- and CVaR-optimization for T = 30 are computed based on Monte Carlo simulations, as explained previously. In order to investigate the statistical accuracy and confirm the predictive value of our study, the framework is evaluated for two additional populations with 150 and 200 treatments. From all three populations, the total sum of squares of voxel dose deviations from the prescribed dose and root mean square errors are obtained and compared with each other. Since these values are of the same order of magnitude, we conclude that the population of 100 treatment is sufficiently large.

In the feasibility study, the objective function value for every robust strategy is evaluated at its respective optimal solution and summarized in table 1 for a mathematical and objective analysis of the proposed strategies. We quantify to which extent the objective function values of the various robust adaptive strategies decrease in comparison with the corresponding non-adaptive robust strategy as the degrees of freedom in the optimization variables increase. Thus, the comparison is performed among the strategies and not among the robust models. Moreover, we investigate how well the MPC-strategies approximate the time-and-scenario-dependent strategy. The objective function values of the time-and-scenario-dependent strategy in combination with CVaR-optimization could not be computed because of computational issues related to limited memory, which underlines the importance of approximation methods, e.g. MPC. Overall, the evaluation of objective function values indicates a number of trends.

Table 1.  Evaluation of the final doses x₅ in all models and strategies for $| {\rm{\Omega }}| =5$ and number of voxels $N=40$.

	non-adaptive	time-varying adaptive	MPC1	MPC2	MPC3	time- & scenario dependent
${\mathbb{E}}$	0.208 6	0.186 8	0.199 1	0.179 4	0.174 3	0.172 9
Stabilized ${\mathbb{E}}$	n/a	n/a	0.209 7	0.201 3	0.200 3	0.200 2
Worst-case	0.483 2	0.341 3	0.343 1	0.317 9	0.305 8	0.297 7
Stabilized worst-case	n/a	n/a	0.347 9	0.327 8	0.319 8	0.316 7
CVaR	0.255 2	0.237 1	0.257 5	0.230 2	0.219 5	−a
stabilized CVaR	n/a	n/a	0.260 6	0.244 5	0.241 1	−a

Note.

^aThe optimal solution of the time-and-scenario-dependent strategy combined with CVaR-optimization could not be computed due to computational limitations.

First, the objective function values of the robust adaptive strategies are lower than those of the respective non-adaptive strategy, which confirms the property of the non-adaptive robust strategy as the lower benchmark and the time-and-scenario-dependent strategy as the upper benchmark in the context of mathematical accuracy. Second, the adaptive strategies which consider uncertainties over a prediction horizon larger than one fraction, i.e. the time-varying robust adaptive, MPC2- and MPC3-strategies achieve lower objective function values than the MPC1-strategy, which indicates a likely benefit of taking into account possible uncertainties over a prediction horizon K > 1. The decrease of objective function value in the time-and-scenario-dependent and MPC-strategies stems from optimizing the plans u₀ and u_t with respect to possible realizations of uncertainty scenarios over the horizon T or K. However, there is a risk of overcompensation relative to the prescribed fraction dose $\tfrac{{d}_{T}}{T}$. Therefore, we introduce stabilizing constraints to keep the fraction dose within the interval $(1-\gamma )\tfrac{{d}_{T}}{T}\leqslant {{Bu}}_{\tau }\leqslant (1+\gamma )\tfrac{{d}_{T}}{T}$, where $\gamma \geqslant 0$ specifies the user-defined lower and upper thresholds for u_t. In our framework, γ is set to 0.05. The nature of overcompensation is discussed later in this section. The trend of decreasing objective function values continues for the stabilized time-and-scenario-dependent- and MPC-strategies as illustrated in table 1, even though it is less pronounced since the stabilizing constraints decrease the feasible region. However, in case of combining CVaR optimization with the MPC1-strategy this trend may not be as strong as in the other robust models, since its objective function value is slightly larger than that of the non-adaptive robust strategy. Third, the combination of worst-case-optimization with the time-and-scenario-dependent- and MPC-strategy seems to lead to the largest improvement in objective function value from the non-adaptive strategy. Thus, worst-case-optimization may be most compatible with time-and-scenario-dependent variables. It should be pointed out that the MPC-strategies are approximations of the time-and-scenario-dependent strategy and that the objective function values are computed from all $| {\rm{\Omega }}{| }^{T}$ attainable final dose distributions and therefore, small differences in objective function values should be interpreted in relative terms

In terms of evaluating the resulting treatment quality throughout a hypofractionated treatment, we gain the following insights. First, the non-adaptive robust and time-varying robust adaptive strategy generate quite uniform dose profiles, as illustrated in figures 2(a)–(d). The x-axis represents the one-dimensional phantom where the location of the CTV and OARs is indicated through vertical separation lines. The dose levels of the realized accumulated dose distributions and plans are given along the y-axis, ranging from 0 to 1.2 relative to the prescriptions dose level 1. The most likely accumulated dose realization in every fraction is indicated in black, while the realizations within one standard deviation are visualized in dark grey. The remaining less likely realizations are shown in light grey. Second, the plans obtained from the time-varying robust adaptive strategy in combination with expected-value-optimization become less uniform as the treatment progresses and the time horizon (T − t) shrinks, as demonstrated in figure 2(b). These plans are optimized in an exact manner by taking advantage of the structure of the expected quadratic penalty, as demonstrated in appendix A, which leads to a linear and quadratic dependency on the shrinking time horizon in the Hessian $(T-t){\mathbb{E}}\left[{S}^{T}{B}^{T}{DBS}\right]\,+(T-t)(T-t-1)$ ${\mathbb{E}}{\left[S\right]}^{T}{B}^{T}{DB}{\mathbb{E}}[S]$ of (3). Thus, as the plans for the last fractions are computed, the contribution of the squared mean value of $(T-t)(T-t-1){\mathbb{E}}{\left[S\right]}^{T}{B}^{T}{DB}{\mathbb{E}}[S]$ decreases and eventually disappears for t = T − 1. Therefore, the plans become less uniform because the term $(T-t){\mathbb{E}}\left[{S}^{T}{B}^{T}{DBS}\right]$ squares the impact of the uncertainties modelled by the shift matrix S. Since this is a result of minimizing the expected quadratic penalty, the dose profiles obtained from the time-varying adaptive strategy combined with worst-case- and CVaR-optimization remain uniform, as illustrated in figures 2(c) and (d), respectively. CVaR-optimization is a compromise between expected-value- and worst-case-optimization in terms of conservativeness which is reflected in the shape of the resulting dose profiles. Third, the plans obtained from the time-and-scenario-dependent robust adaptive and the MPC-strategies are optimized in order to account for the possible realizations of uncertainty scenarios over the horizon T or K, which may cause overcompensation as mentioned earlier. The order of magnitude of overcompensation seems to correlate with the planning horizon T and K. We observe that the longer the planning horizon, the larger the overcompensation will be at the edges of the CTV which are most likely to be affected by interfractional uncertainties, as illustrated in figures 3(a) and (b) for the combination of expected-value-optimization with time-and-scenario-dependent strategy and the MPC-strategy with K = 3, respectively. In expected-value optimization, the overcompensation takes shape in the form of overdose at the edges of the CTV. In worst-case-optimization, the impact of future uncertainties is handled by a flat dose-fall-off at the edges in order to reduce the risk to the nearby OARs in the event of the largest uncertainty, as shown in figure 3(c). Furthermore, worst-case-optimization results in a smaller spread of the various realized accumulated dose distributions x_t+1. As a result of the stabilizing constraints, the accumulated doses and plans are more uniform, as illustrated in figures 4(a)–(c). Moreover, the spread of dose realizations around the most likely outcome is decreased. In the non-adaptive robust and time-varying robust adaptive strategy, the stabilizing constraints are not imposed since these strategies generate plans which do not exceed the prescribed fraction dose across the whole CTV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the conventional fractionation length study, we evaluate the interplay between robustness and non-adaptive and adaptive strategies for a conventionally long treatment with 30 fractions. The performance of the non-adaptive robust, time-varying robust adaptive, MPC1-, MPC2- and MPC3-strategy combined with expected-value-, worst-case- and CVaR-optimization is measured by evaluating their final dose distribution x_T for a population of 100 treatments with predictable interfractional uncertainties.

Their performance in terms of CTV-coverage under uncertainty is illustrated by dose-probability histograms in figure 5. The dose-probability histograms show the probability that the final accumulated dose distribution x_T provide CTV-coverage in accordance with the ICRU guidelines [32] which recommend that 99% of the CTV should receive at least 95% of prescription dose. The probability that 99% of the CTV will receive at least a certain dose level or above is plotted on the y-axis. In order to facilitate the analysis of the performance of the various robust non-adaptive and adaptive strategies, the recommended dose level of 95% of the prescription dose is displayed in figure 5. According to the results illustrated in figures 5(a)–(c), the majority of strategies and robust models appear to provide sufficient target coverage as given by the ICRU-guidelines with similarly high probability. This result indicates that robust optimization in combination with the conventional non-adaptive strategy ensures CTV-coverage, if the uncertainties occurring throughout the course of treatment are predicted correctly and are included in the robust optimization problem. Among the investigated robust models, worst-case-optimization may be the most favourable model to be combined with the non-adaptive treatment approach, as demonstrated in figure 5(b).

Zoom In Zoom Out Reset image size

Concerning the interplay between robust models and the adaptive strategies, expected-value optimization seems to be a good match for the time-varying robust adaptive strategy, which achieves target coverage with the highest probability compared to the MPC-strategies, as shown in figure 5(a). According to our evaluation of CTV-coverage probability for the various robust strategies combined with expected-value optimization in figure 5(a), the MPC-strategies which are based on time-and-scenario optimization variables may not be as a good match for expected-value optimization as the time-varying robust adaptive strategy. However, it should be noted that the MPC-strategies are approximative methods in contrast to the time-varying adaptive strategy, which in combination with expected-value optimization gives exact solutions. Since the expected-value over the anticipated scenarios gives little weight to extreme scenarios with a low probability, the use of time-and-scenario-dependent optimization variables may not be suitable. On the other hand, time-and-scenario-dependent optimization variables appear to be good combination for the more conservative optimization models as shown in figures 5(b) and (c) according to which the MPC-strategies provide target coverage with the highest probability. In the linear formulations of our worst-case- and CVaR-optimization problems (7) and (6), every realization of uncertainty scenarios is represented by a corresponding constraint. Time-and-scenario-dependent optimization variables seem to exploit this feature and thus, the additional computational effort may be worthwhile in the case of combining MPC-strategies with conservative robust optimization models. Based on the observation that all MPC-strategies show very similar trends in terms of CTV-coverage, the most suitable prediction horizon cannot be clearly identified. The similarity of these trends among the MPC-strategies may be caused by their rather short prediction horizons relative to the treatment length of 30 fractions.

Concerning the risk of high dose exposure in the OAR, the accumulated dose in the OARs is evaluated by evaluating the dose levels of the 90th percentile over the 100 treatments. We evaluate the 90th percentile in order to visualize the order of magnitude of an upper dose threshold for 90% of the treatments with the robust optimization models which according to our framework are most likely to provide sufficient CTV-coverage. Since the right OAR is closer to the CTV, it may likely receive a higher dose than the left OAR in the event of interfractional uncertainties. Since expected-value- and worst-case-optimization appear to provide the highest probability of CTV-coverage, we choose to draw attention to their impact on the right OAR. The voxels of the right OAR which are in proximity to the CTV receive the largest amount of dose as shown in figures 6(a) and (b), which illustrate the  90th percentile for each voxel in the right OAR after 100 treatments. In case of expected-value optimization, the time-varying robust adaptive strategy leads to the least dose exposure to the right OAR among the adaptive strategies, but its dose levels are comparable to those of non-adaptive strategy. Since the combination of the time-varying robust adaptive strategy with expected-value optimization gives the highest probability of CTV-coverage, we may conclude that a robust adaptive strategy based on expected-value optimization may be most effective with time-and-scenario-independent optimization variables. In case of worst-case-optimization, all strategies lead to similar upper dose levels in the right OAR, indicating that robust adaptive strategies may not increase the risk of high dose exposure to the OARs. Moreover, the upper dose thresholds of the accumulated dose in the right OAR are higher than in combination with expected-value optimization. However, these differences between non-adaptive and adaptive strategies and the robust models are rather small since the simulated treatments follow the same distribution as those included in the robust models.

Zoom In Zoom Out Reset image size

Overall, we conclude that if the treatment is supposed to be adaptive with a high adaptation frequency throughout the course of treatment , it may be the most beneficial to combine expected-value optimization with time-and-scenario-independent optimization variables. If however, the treatment is supposed to be executed in a non-adaptive manner or only adaptive in rare cases, a more conservative model in combination with time-and-scenario-dependent optimization variables may be recommended.