Adjustable robust treatment-length optimization in radiation therapy

Technological advancements in radiation therapy (RT) allow the collection of biomarker data on individual patient response during treatment. Although biomarker data remains subject to substantial uncertainties, information extracted from this data may allow the RT plan to be adapted in an informative way. We present a mathematical framework that optimally adapts the treatment-length of an RT plan based on the acquired mid-treatment biomarker information, and also consider the inexact nature of this information. We formulate the adaptive treatment-length optimization problem as a 2-stage problem, where after the first stage we acquire information about model parameters and decisions in stage 2 may depend on this information. Using Adjustable Robust Optimization (ARO) techniques we derive explicit optimal decision rules for the stage-2 decisions and solve the optimization problem. The problem allows for multiple worst-case optimal solutions. To discriminate between these, we introduce the concept of Pareto Adjustable Robust Optimal (PARO) solutions. In extensive numerical experiments based on liver cancer patient data, ARO is benchmarked against several other static and adaptive methods, including robust optimization with a folding horizon. Results show good performance of ARO both if acquired mid-treatment biomarker data is exact and inexact. We also investigate the effect of biomarker acquisition time on the performance of the ARO solution and the benefit of adaptation. Results indicate that a higher level of uncertainty in biomarker information pushes the optimal moment of biomarker acquisition backwards.


Introduction
In radiation therapy (RT), the goal is to deliver a curative amount of dose to the target volume(s), while keeping the dose to all organs-at-risk (OARs) at a minimum.As the radiation beam delivers energy to all tissues that are in its path, the OARs will (often) inevitably receive some dose as well.In practice, for all OARs it is determined how much dose they can safely tolerate before serious complications arise, and the goal is to design a treatment plan that keeps the dose to all OARs within tolerable limits.In order to do so, the treatment plan is optimized spatially and temporally.
Spatial optimization exploits the fact that, by mounting the beam head on a gantry, the tumor can be targeted from various angles.It aims to find the combination of beam angles and weights that gives the best trade-off between tumor dose conformity and healthy tissue sparing.There is a large body of literature on this topic, see for example Shepard et al. (1999), Ehrgott et al. (2008) for reviews.The result of spatial optimization is a dose distribution, which indicates how much dose will be received by each voxel (3-dimensional subvolume) of the tumor and OARs.
Temporal optimization, on the other hand, is concerned with determining the optimal treatment length, or optimal number of treatment fractions.It is based on the concept of fractionation: compared to tumor cells, healthy tissues have better repair capabilities between fractions (Fowler 1989, Withers 1985).At the same time, longer treatments give the tumor the possibility to proliferate, which endangers treatment efficacy.Therefore, there might be an optimal number of treatment fractions that maintains the best balance between tumor proliferation and OAR recovery.Treatments with a higher number of fractions and a lower dose per fraction than the conventional regimen are known as hyperfractionated treatments.Treatments with a lower number of fractions and a higher dose per fraction than conventional are known as hypofractionated treatments.
Recently, spatiotemporal optimization has been gaining some attention (Unkelbach et al. 2013b, Unkelbach and Papp 2015, Kim et al. 2015, Ajdari and Ghate 2016, 2017b, Saberian et al. 2017, Adibi and Salari 2018).In these studies, the spatial and the temporal component are optimized simultaneously.Solving the resulting optimization problems is in general more computationally demanding.
Technological advances in treatment monitoring through imaging and other forms of data acquisition allow for a more accurate assessment of radiation response (Baumann et al. 2016).Biological-based adaptive treatments aim to monitor the treatment, acquire mid-treatment biomarkers, and adapt the remainder of the treatment course accordingly.From a mathematical optimization perspective, this presents an adaptive optimization problem.
Literature on biological-based adaptive treatment optimization is limited.Ghate (2011) and Kim et al. (2012) propose a theoretical stochastic control framework to optimally adapt the beam intensities over a fixed number of fractions, based on hypothetically-observed tumor states.Saberian et al. (2016b) concretize this theoretical framework, using simulated hypoxia status as biomarker.Ajdari et al. (2018) considers a mathematical model where they adaptively determine the optimal number of treatment fractions, in order to minimize the total number of tumor cells remaining (TNTCR) at the end of the treatment course.After each fraction they observe the tumor cell density in each voxel, and use this to re-optimize a spatiotemporal model.In robust optimization terminology, this is known as a folding horizon (FH) approach.
The limited availability and accuracy of required biomarkers poses a primary challenge for adaptive treatments (Baumann et al. 2016).Any biomarker information from data gathered during treatment remains subject to uncertainties, stemming from both measurement errors and the inexactness in the translation of measured data to model parameters.Therefore, it is crucial that any adaptive treatment optimization takes this into account.Ajdari et al. (2019) provide an overview of relevant mathematical (optimization) tools.Robust Optimization (RO) (see, e.g., Ben-Tal et al. (2009), Bertsimas et al. (2011), Gorissen et al. (2015) for an overview) has been the predominant method for dealing with uncertainties in radiation therapy treatment planning.Amongst others, it has been applied for dealing with motion uncertainty (Chan et al. 2006, Bortfeld et al. 2008) and delineation uncertainty (Balvert et al. 2019).
Adjustable Robust Optimization (ARO) (Ben-Tal et al. 2004, Yanıkoğlu et al. 2019) is an extension of RO that takes into account the flow of information over time and exploits the fact that some decisions need to be taken only after the data has (partially) revealed itself.In the standard paradigm, ARO assumes that the revealed data is exact; de Ruiter et al. (2017) introduces ARO methodology for the case when revealed data is not exact but provides only an estimate of the true parameters.This paper considers robust adaptive treatment-length optimization based on mid-treatment acquired biomarker information, taking into consideration the inexact nature of the acquired biomarker data.Our main contributions are: • We develop mathematical tools based on ARO that enable us to (i) optimally adapt the treatment length after acquiring mid-treatment biomarker information, (ii) analyze the influence of biomarker information uncertainty.To the best of our knowledge, this is the first application of ARO to an RT problem.
• We present explicit optimal decision rules for a difficult (non-convex, mixed-integer) yet practically relevant ARO problem.
• We show that there are multiple optimal solutions for the worst-case scenario, and that these perform very differently in other scenarios.To handle this, we introduce the concept of Pareto Adjustable Robustly Optimal (PARO) solutions, a generalization of Pareto Robustly Optimal (PRO) solutions (Iancu and Trichakis 2013) to two-stage robust optimization problems.In case the acquired biomarker data is exact, PARO solutions are obtained.
• We perform a computational study to determine the optimal timing of acquiring biomarker information in case biomarker quality improves over time.Later biomarker acquisition also limits adaptation possibilities, and the optimal balance depends on the improvement rate.
The remainder of this paper is organized as follows.Section 2 introduces the BED model and the nominal treatment-length optimization problem.Section 3 introduces the adaptive treatmentlength optimization problem under the assumption of exact biomarker information and solves this using ARO techniques.Section 4 generalizes this to inexact biomarker information.Section 5 reports and discusses results of the numerical experiments.Finally, Section 6 states several concluding remarks.
Notation.All variables and constants are 1-dimensional (belong to R or N) unless indicated otherwise.In functions, a semicolon (;) is used to separate variables and constant arguments from uncertain parameters.Optimal solutions to optimization problems are indicated with an asterisk ( * ).Properties of optimal solutions to optimization problems have calligraphic font (e.g., ARO) to distinguish them from methods with the same or similar abbreviations.

The BED model
The BED model (Fowler 1989, 2010, Hall and Giaccia 2012) states that the biological effect of an N -fraction dose sequence d = (d 1 , . . ., d N ) (in Gray (Gy)) delivered to a tumor volume is given by which is a model governed by a single parameter, the α/β ratio, which signifies the fractionation sensitivity of the tumor tissue.
Let (d * , N * ) denote the optimal solution to (2).A simple analysis in Mizuta et al. (2012) reveals the following important result: otherwise. (3) In both cases the optimal dose d * is such that (2b) is active.It can be shown that the same result holds if we maximize tumor BED subject to an upper bound on OAR BED (Bortfeld et al. 2015).Similar results have been derived for amongst others heterogeneous dose distributions (Unkelbach et al. 2013a) and multiple normal tissues (Saberian et al. 2015).However, all of these approaches assume the tumor and OAR radiosensitivity parameters τ and ρ are known exactly.There is much research on the α/β ratios for different tumor sites (Tai et al. 2008, Son et al. 2013, Klement 2017) and OAR sites (van Leeuwen et al. 2018), but they remain subject to considerable uncertainties.Badri et al. (2016) takes a stochastic programming approach, assuming a normal distribution for τ .Ajdari and Ghate (2017a) takes a robust optimization approach.However, if biomarker information acquired during treatment provides more accurate information on τ than what was available at the start of the treatment, a static stochastic programming or robust optimization approach may be overly conservative.

ARO: Biomarkers provide exact information
We present an adjustable robust optimization approach that optimally adjusts the remainder of the treatment once biomarker information has provided the true value of parameters τ and ρ.

Problem formulation
In order to establish a meaningful model for the adjustable robust optimization approach, we restrict the dose sequence d = (d 1 , . . ., d N max ) in several ways.We set a minimum number of fractions N min .Furthermore, we assume there is a single moment N 1 where we can adapt the treatment.The dose per fraction in the first N 1 fractions is assumed to be the same, denote this by d 1 .Variable N 2 denotes the number of fractions after adaptation; also these fractions have equal dose, denoted by d 2 .This implies with We additionally set the constraint that d 1 , d 2 ≥ d min , for some predetermined value d min .Lastly, we set a maximum dose per fraction in stage 1, to avoid delivering dosages in stage 1 that severely restrict adaptation possibilities in stage 2. We will later impose some restrictions on the allowed combinations of d min , d max 1 and N max

2
. Figure 1 provides a schematic overview of the situation.We wish to maximize the tumor BED, subject to the constraint that the BED to the OAR is below the tolerance level BED tol , given by i.e., the OAR is known to tolerate a total dose of D Gy if delivered in T fractions under dose shape factor ϕ. The dose shape factor is a parameter characterizing the spatial heterogeneity of a dose distribution, for more details see Perkó et al. (2018).Note that BED tol (ρ) is a function of uncertain parameter ρ.
After delivering N 1 fractions we observe the true (ρ, τ ).At the time that the stage-1 dose per fraction d 1 has to be decided, it is only known that (ρ, τ ) belongs to some uncertainty set Z. Put together, the Exact Data Problem (EDP) reads: The value for the stage-1 dose d 1 has to be decided before the value of (ρ, τ ) is revealed; in ARO this is also commonly referred to as a here-and-now variable or decision.The values stage-2 dose d 2 and stage-2 number of fractions N 2 need to be decided only after (ρ, τ ) is revealed as they may depend on the values of these parameters.Hence, they are written as functions d 2 (ρ, τ ) and N 2 (ρ, τ ) of the uncertain parameters (ρ, τ ).In ARO such variables are also commonly referred to as wait-and-see variables or decisions.In this paper, we will adhere to the terms stage 1 and stage 2, however.
We assume box uncertainty of the form: with 0 < ρ L < ρ U and 0 < τ L < τ U .It is assumed that there is a nominal scenario, e.g., parameter values τ and ρ derived from literature.The first observation we make in ( 6) is that the objective is increasing in parameter τ , so in any worst-case realization it will hold that τ = τ L .This observation has consequences for what uncertainty sets Z need to be considered.Due to (3), one can in general distinguish three cases for uncertainty set Z and parameter σ: Case 1) σρ L ≥ τ L : According to (3), for any realization (with τ = τ L ) it is optimal to deliver the minimum number of fractions in stage 2.
Case 2) σρ U ≤ τ L : According to (3), for any realization (with τ = τ L ) it is optimal to deliver the maximum number of fractions in stage 2.
Case 3) σρ L < τ L < σρ U : In the scenario (ρ L , τ L ), it is optimal to deliver the maximum number of fractions in stage 2 according to (3).In the scenario (ρ U , τ L ), it is optimal to deliver the minimum number of fractions in stage 2 according to (3).
In case 1 and 2, ( 6) is easily solved by plugging in the (worst-case) optimal value for N 2 , and solving the resulting 2-variable optimization problem.Therefore, only Case 3 is of interest and in the remainder of this paper it is assumed that Problem ( 6) is a 2-stage non-convex mixed-integer ARO problem, which are generally hard to solve.Nevertheless, due to the small number of variables the problem can be solved to optimality.

Optimal decision rules and worst-case solution
Before we solve (6), we first need some additional definitions and assumptions.The remaining BED tolerance level of the OAR, if N ′ fractions with dose d ′ have been administered, is given by Subsequently, define the function The value of g can be interpreted as the maximum dose that can be delivered in N ′′ fractions if already N ′ fractions of dose d ′ are (scheduled to be) delivered.It is obtained by solving (6b) for d 1 or d 2 .Functions of this form will be used frequently throughout the remainder of this paper.
The following assumption on the relation between d min , d max 1 and N max 2 makes sure that for a given optimal number of fractions, it is feasible to deliver that number of fractions with minimum dose.

Assumption 1. It holds that
The particular form of the upper bound on d max 1 will become clear later.Numerical experiments indicate that results are not very sensitive to the choices of d min and d max 1 .In order to solve (6), we take three steps: Step 1) Eliminate variable d 2 (ρ, τ ) by writing it as a function of d 1 and N 2 (ρ, τ ).
Step 2) Determine the optimal N 2 as a function of (ρ, τ ).
Step 3) Solve the resulting problem of variable d 1 .
In what follows, we give a detailed explanation of each of these steps.
Step 1: eliminate variable d 2 (ρ, τ ) ) denote an optimal solution to (6).In the optimum, constraint (6b) holds with equality because it is the only dose-limiting constraint.Solving for d 2 yields This suggests the following objective function for given (ρ, τ ): where, for given ρ, the value g(0, 0, N 1 ; ρ) is the maximum dose that can be delivered in stage 1 due to the nonnegativity restriction on the stage-2 dose.From Assumption 1 it follows that g(0, 0, N 1 ; ρ) ≥ d max 1 for all (ρ, τ ) ∈ Z.According to Lemma 1 in Appendix B, function f is either convex, concave or constant in d 1 .
Due to the above steps, EDP ( 6) is equivalent to max Step 2: determine the optimal N 2 (ρ, τ ).
Definition 1 (Adjustable robustly feasible).A pair Subsequently, define adjustable robust optimality of a solution.

A stage-1 decision
Although definitions are given for ( 14), due to (12) these readily extend to adjustable robust feasibility/optimality for the original EDP (6).The following theorem, similar to the result (3), states the optimal stage-2 decision rules.
Theorem 1.The decision rule Proof.See Appendix A.1.By ( 12), this implies that the unique optimal stage-2 decision rule for d 2 as a function of d 1 and the uncertain parameters is Clearly, these decision rules are non-linear, and in fact split the uncertainty region in two parts: one where it is optimal to deliver the minimum number of fractions N min 2 in stage 2, and one where it is optimal to deliver the maximum number of fractions N max 2 in stage 2.This suggests splitting the uncertainty set as follows: An illustration is provided in Figure 2.
Step 3: solve the resulting problem of variable d 1 .
It remains to determine the optimal dose per fraction in stage 1.Using ( 16) and ( 17), we can reformulate (14) to a problem of only variable d 1 .Note that uncertain parameter τ plays no other role than to determine the optimal number of stage-2 fractions N * 2 .Uncertain parameter τ appears only in the objective and function f is increasing in τ , so it is sufficient to consider only those observations with τ = τ L .
In order to reformulate (14), we make use of the properties of g and f in Lemma 3 in Appendix B. In particular, Lemma 3b states that function f is either increasing or decreasing in ρ for fixed d 1 .Hence, if we move (14a) to a constraint and split according to (18), for both Z min and Z max it is sufficient to consider the constraint for the highest and lowest value of ρ in the uncertainty set.With τ = τ L , this yields the scenarios (ρ L , τ L ) and ( τ L σ , τ L ) for Z min and ( τ L σ , τ L ) and (ρ U , τ L ) for Z max .Due to Lemma 3a, we find the same candidate worst-case scenarios for constraint (14b).
Figure 3 illustrates a possible instance of (20), displaying constraints (20b)-(20d).In this case, the set of optimal solutions is the union of the two intervals for d 1 where constraint (20d) is active.Dose constraints (20e) may cut off part of these intervals.If due to constraint (20e) both these intervals are infeasible, the optimum is at one of the boundaries for d 1 .

Pareto adjustable robustly optimal solutions
Figure 3 illustrates that it is possible that there are multiple optimal solutions to (20).These solutions are ARO stage-1 solutions to the EDP (6).In general, in case there are multiple ARO solutions these may perform vastly different if a non-worst-case scenario realizes (de Ruiter et al. 2016).Iancu and Trichakis (2013) studies static robust optimization problems with multiple robustly optimal solutions, and introduce the concept of Pareto robustly optimal (PRO) solutions.A robustly optimal solution is called PRO if there is no other robustly feasible solution that is has equal or better objective value for all scenarios in the uncertainty set, while being strictly better for at least one scenario.Non-PRO solutions are dominated by at least one PRO solution and are therefore not desired.The concept of Pareto robust optimality closely resembles the concept of Pareto efficiency in multi-criteria optimization (MCO).In MCO, Pareto efficient solutions can only be improved in one criteria at the cost of a deterioration in another criteria.Only Pareto efficient solutions are of interest, and the overall goal in MCO is to compute this set of solutions (known as the Pareto surface).
We generalize the concept of Pareto robust optimality to 2-stage adjustable robust optimization problems.
We also define the concept PARO for the stage-1 decision d 1 or the stage-2 decision (rule) N 2 (•) only.

A stage-1 decision
If there are multiple ARO solutions, we wish to pick one that is PARO.In Iancu and Trichakis (2013) it is shown for linear optimization that, if we optimize over the robustly optimal solutions for a second criterion (scenario) that is in the relative interior of the uncertainty set, the resulting solution(s) are PRO.In a similar fashion PARO solutions to the current problem can be found.
Let X ARO denote the set of ARO solutions to (20).It turns out that consecutively optimizing over an auxiliary scenario where hyperfractionation is optimal and an auxiliary scenario where hypofractionation is optimal yields a set of PARO solutions.Let (ρ aux-min , τ aux-min ) ∈ int(Z min ), where int(•) is the interior operator.Define the auxiliary optimization problem for the hypofractionation scenario: Denote the set of optimal solutions by X aux-min .Similarly, let (ρ aux-max , τ aux-max ) ∈ int(Z max ).Define the auxiliary optimization problem for the hyperfractionation scenario: Note that it uses X aux-min as input, i.e., we solve the auxiliary problems consecutively.Denote the set of optimal solutions by X PARO .
Theorem 2. All solutions in X PARO are PARO.
Solving ( 22) or ( 23) asks to maximize a strictly convex or strictly concave function over a feasible set consisting of a small number of intervals or points.Hence, these auxiliary problems are easily solved.Note that the second auxiliary problem is only relevant if the first auxiliary problem has multiple optimal solutions.Switching their order, and optimizing ( 22) over the set X aux-max may lead to different solutions, and these are also PARO.The two-step approach is necessary; numerical results show that optimizing over one auxiliary scenario may indeed yield non-PARO solutions.

ARO: Biomarkers provide inexact information
In this section we present an adjustable robust optimization approach to solve a more realistic version of the adaptive treatment-length problem.Because in practice it is impossible to exactly determine the α/β parameters from biomarker data, any values for the α/β parameters obtained during treatment are inexact.This section presents a model that accounts for uncertainty in biomarker observations.

Problem formulation
The setup for the ARO problem with inexact data is based on de Ruiter et al. (2017).After N 1 fractions we obtain an estimate (ρ, τ ) for (ρ, τ ), the inverse α/β parameters for the OAR and the tumor.It is still assumed that (8) holds for uncertainty set Z. Furthermore, we assume that (ρ, τ ), (ρ, τ ) ∈ Z (as defined in ( 7)) and that (ρ, τ ) Here r ρ and r τ are parameters that define the accuracy of the estimate/observation (ρ, τ ).This can also be written as (ρ, τ ) ∈ {(ρ, τ )} + Ẑ, which is the Minkowski sum of a singleton and a set.For given observation (ρ, τ ), the new upper and lower bounds for (ρ, τ ) are given by Define The set U contains all possible observation-realization pairs.Lastly, we remove Assumption 1 and impose a different (slightly stricter) assumption on the relation between The inexact data problem (IDP) analogous to ( 6) is given by max Let X(ρ, τ, ρ, τ ) denote the feasible region defined by constraints (28b)-(28e) for fixed (ρ, τ, ρ, τ ).For stage-2 variables d 2 (ρ, τ ) and N 2 (ρ, τ ) it is indicated that they are a function of the observations (ρ, τ ) instead of the uncertain parameters (ρ, τ ).

Optimal decision rules and conservative approximation
Depending on both the observed (ρ, τ ) and the quality of the biomarker observation (i.e., r ρ and r τ ), we may be able to immediately determine the optimal value for N 2 .Therefore, we split the uncertainty set for the observations (ρ, τ ).Define fractions can be optimal in stage 2. Subset Z int ID is the area between the red lines.
and N max 2 fractions in stage 2 may be optimal for the true (ρ, τ ).Before this is formalized in a theorem, we state two definitions.

Definition 6 (Adjustable robust feasibility
The following theorem states the optimal stage-2 decision rules for a given value of d 1 . Theorem 3. Let d 1 be the stage-1 decision of (28).The decision rules and are ARO solutions for N 2 (•) and d 2 (•) in (28) for the given d 1 .
The worst-case optimal decision rule (31) may yield an intermediate value if (ρ, τ ) ∈ Z int ID .If r ρ and r τ are zero, i.e., we have exact data, then it holds that τL = τ and ρL = ρU = ρ.Hence, the two functions f in the RHS of (31) are equal, and the optimal N * 2 is the one that maximizes the resulting function.One can verify that this does not depend on d 1 .Hence, in case of exact data Theorem 3 is equivalent to Theorem 1.
It turns out that, after plugging in ( 31) and (32), and splitting the uncertainty set according to (29), it is not apparent how to determine the optimal stage-1 decision d * 1 for (28).In Appendix A.4 the following lower bound problem to (28), named the Approximate Inexact Data Problem (AIDP), is derived: Compared to (14) for exact biomarker observations, problem (33) has the added constraint (33e); a piecewise convex-concave function p(d 1 ) defined by (B.68) in Appendix A.4. Lemmas 1, 2 and 5 in Appendix B provide information on the shape and intersection points of constraint functions (33b)-(33e).Consequently, the optimal solution(s) is/are easily obtained.
From an optimal solution (d * 1 , q * ) to AIDP (33) an ARF solution to the original IDP (28) can be constructed by omitting q * and using stage-2 decisions (31) and ( 32).This solution is ARF because AIDP is a conservative approximation of IDP.
Figure 5 illustrates a possible instance of (33), displaying constraints (33b)-(33e).In this case, optimal solutions are locations where (33d) and (33e) are both binding, indicated by the circles.Dose constraints (33f) may cut off (some of) these points.If due to constraint (33f) none of the circles are feasible, the optimum is at one of the boundaries for d 1 .Constraint (33e) is the only conservative constraint in (33).Hence, only if the feasible values for d 1 are such that none of the circles in Figure 5 are feasible and constraint (33e) is binding, it is possible that the optimal value of (33) is strictly worse than that of (28).

Pareto robustly optimal solutions to conservative approximation
Figure 5 also illustrates that it is possible that there are multiple optimal solutions to the AIDP (33).Because the AIDP provides a conservative approximation to IDP (28), optimizing over auxiliary scenario(s), as in Section 3, does not necessarily yield a stage-1 decision d 1 that is PARO to the original IDP.It turns out that a PRO solution to AIDP is obtained from the set of robustly optimal solutions to AIDP if we consecutively optimize for two auxiliary observations such that any worst-case realization is in the interior of set Z min resp.Z max .Two important Figure 5: Illustration of (33).Compared to the case with exact data (Figure 3), the thick black curve (constraint (33e)) is extra.The solid, dashed and dotted lines represent constraints (33b), (33c) and (33d), respectively.Locations where multiple constraints are binding are indicated by circles.remarks are in place here.First, a PRO solution to AIDP need not be a PARO solution to IDP, even if it is worst-case optimal to IDP.Second, the required auxiliary scenarios need not exist; their existence depends on the values of r ρ and r τ .Hence, further details are omitted.

Numerical results
This section presents numerical results of the methods presented in Sections 3 and 4. First, Section 5.1 describes the benchmark methods against which we compare the ARO method for EDP and IDP, and Section 5.2 describes the setup of the numerical experiments.

Benchmark static and folding horizon methods
We analyze the performance of the static and folding horizon nominal method (NOM and NOM-FH), the static and folding horizon robust optimization method (RO and RO-FH) and the adjustable robust optimization method (ARO).In the folding horizon approaches only the stage-1 decisions are implemented, and the model is re-optimized for the second stage once the biomarker information is revealed.
The static method NOM optimizes for the nominal parameter values τ and ρ and disregards any uncertainty and adaptability.This method is the same for both EDP and IDP.In stage 2, NOM-FH solves the nominal problem under the assumption that the obtained biomarker estimate is exact (which is a false assumption for IDP).This method does not guarantee robust feasible solution (feasible for all (ρ, τ ) ∈ Z) nor a robustly optimal solution (RO; optimal for the worst-case realization (ρ, τ ) ∈ Z).
The static method RO optimizes for the worst-case realization of (τ, ρ) in the uncertainty set Z, and disregards adaptability.For EDP the method RO-FH solves the same nominal problem as NOM-FH in stage 2; for IDP it solves a static robust optimization problem in stage 2. The uncertainty set is defined by the accuracy of the biomarker information.RO and RO-FH both guarantee an RO solution.
One may add a folding horizon component to ARO (for either EDP or AIDP).This may improve the results in case a suboptimal stage-2 decision rule is used.However, as shown in Sections 3.2 and 4.2, the used stage-2 decision rules are optimal for any realized scenario (and for given stage-1 decision d 1 in case of inexact information).Hence, adding a folding horizon component would not change results.
Table 1 provides an overview of the guaranteed solution properties of the methods.It is important to note that in case of inexact biomarker data (IDP) the methods RO and RO-FH guarantee an RO solution, whereas ARO guarantees only an ARF solution via solving the approximate problem AIDP.Depending on the approximation quality, the ARF solution may be close or equal to an ARO solution.Next to these five methods, we also report the results for the perfect information optimum (PI).This is the attainable optimum if from the start of the first fraction the true (ρ, τ ) is exactly known.It can be formulated by taking the nominal problem and replacing the nominal parameter values by their true values.While in practice not a viable strategy, PI provides information on the value of perfect information, and allows us to put the performance and differences between the other methods in perspective.

Study setup
The data set contains data from 17 liver patients treated with proton therapy at Massachusetts General Hospital (Boston, USA).The mathematical models in Sections 3 and 4 are based on the assumption that there is a single dose restricting OAR.We assume that the single dose restricting OAR is the normal liver itself.For the models in Sections 3 and 4, an instance is defined by a tuple (σ, ϕ, D, T, N 1 , N min , N max , d min , d max 1 ) and the relevant uncertainty sets.Clinically, the patients were treated with either 5 or 15 fractions.We extend the upper bound to allow for more hyperfractionated treatments, and set N min = 5 and N max = 20.We assume the biomarker acquisition is made once N 1 = 4 fractions have been administered.This implies N min 2 = 1 and N max 2 = 16.The mean OAR dose tolerances were derived for both a 5 and a 15 fraction scheme (this is parameter T in the mathematical models).We pick T = 15, as this will be more realistic for hyperfractionated treatments than the 5 fraction variant.The OAR dose tolerance range corresponding with this scheme is D ∈ [20,27].Instead of directly using the OAR dose distribution data from the 17 patients, we use it to generate a larger and more heterogeneous set of phantom patient parameter tuples (σ, ϕ, D).For details see Appendix C.
We define two uncertainty sets for (ρ, τ ).The α/β ratio for the normal liver in the data set is 4, so we set nominal value ρ = 0.25.The tumor α/β ratio in the data set is 10, but we deviate from this.The reason is that condition (8) will not hold for any reasonable uncertainty set for an α/β ratio of 10.In those cases, the optimal stage-2 fractionation is known a priori.Therefore, we set nominal value τ = 0.25, which corresponds to a nominal tumor α/β ratio of 4. In all numerical experiments we consider the two uncertainty sets of To discriminate between multiple ARO solutions, we follow the procedure detailed in Section 3.3 in case of exact biomarker data.The auxiliary scenarios are sampled uniformly from int(Z min ) and int(Z max ).In case of inexact biomarker data, the procedure discussed in Section 4.3 is followed if the required auxiliary observations exist.If such observations exist, we sample uniformly from Z until we have found two auxiliary observations for which any worst-case realization is in int(Z min ) resp.int(Z max ).If such observations do not exist, the robustly optimal solution to AIDP with lowest stage-1 dose is selected.The method RO (and therefore also RO-FH) may also find multiple robustly optimal solutions.For the obtained set of robustly optimal solutions we again follow the procedure detailed in Section 3.3.It turns out that for RO, the robustly optimal solutions often perform identical in non-worst-case scenarios.We optimize over the auxiliary scenarios consecutively; the first auxiliary scenario is the scenario corresponding to int(Z min ).
The minimum dose per fraction is d min = 1.5 Gy.For ARO the stage-1 upper bound d max 1 is set at the upper bound indicated by Assumptions 1 (for EDP) and 2 (for IDP).For RO, parameter d max 1 is set such that it is feasible to deliver N max 2 fractions with dose d min in stage 2 in all scenarios in Z.For NOM (and NOM-FH) and PI parameter d max 1 is set such that this is feasible under the nominal and true uncertainty scenario, respectively.Numerical results indicate that results are not sensitive to the choice of d min (and the corresponding d max 1 ).First, Section 5.3 presents and discusses the results for the problem with exact biomarker observations (EDP) of Section 3.After that, Section 5.4 presents and discusses the results for the problem with inexact biomarker observations (IDP) of Section 4. Lastly, Section 5.5 again considers the inexact biomarker observation case, and varies parameter N 1 in order to determine the optimal moment of biomarker acquisition.All results are averages over a sample of S POP = 50 patients.All computations were performed on a 3.4Ghz Intel-Core i7 PC with 16GB RAM, using the software package MATLAB R2018b (Mathworks, Natick, MA, US).

Results exact biomarker observations
Table 3 presents the results.First, we find that adding the folding horizon component to the static nominal and robust method improves the 5%-quantile and mean performance.For the nominal method, including FH additionally leads to solutions that have no constraint violations and improves the worst-case performance.
All methods except NOM are worst-case optimal (i.e., RO resp.ARO for the static resp.adjustable methods) in both scenarios.However, their mean performance differs, which implies the existence of multiple worst-case solutions.Indeed, results in Table 3 indicate that these methods have different stage-1 decisions d 1 .As indicated in Section 5.2, RO, RO-FH and ARO optimize over auxiliary scenarios in this case; according to Theorem 2 ARO finds a PARO solution this way.Mean tumor BED for ARO is 2.3% higher than for RO-FH and 10.4% higher than for RO in uncertainty set 1.
Overall, methods that deliver a relatively low dose in stage-1 perform better on average than the methods that deliver a relative high dose.This behavior can be explained as follows.For scenarios where hypofractionation is optimal, it is optimal to deliver a low dose per fraction in the first N 1 = 4 fractions, and to deliver a single high dose fraction in fraction 5.This is the best possible approximation for a single fraction treatment under the current set of constraints.On the other hand, for scenarios where hyperfractionation is optimal, it is best to deliver an intermediate dose per fraction throughout all N 1 + N max 2 = 20 fractions.Such treatments best approximate a uniform fractionation scheme.The delivered dose is likely higher than d min , but not very high as it has to be delivered for a large number of fractions.The optimal stage-1 dose best balances this trade-off.
In Table 3, the mean performance of ARO is slightly better than RO-FH, but equal to NOM-FH.Indeed, their stage-1 decision is equal on average and in stage-2 they possess the same information.The good performance of NOM-FH indicates that, in case of exact biomarker information, ignoring uncertainty in stage 1 does not compromise mean performance nor leads to OAR constraint violations.This observation, together with the fact that NOM has large OAR constraint violations, shows that the value of (exact biomarker) information is high.If a folding horizon step is included, solving the static robust problem in stage 1 even yields worse solutions than when the static nominal problem is solved.
The difference between uncertainty sets 1 and 2 is mostly in a larger spread in performance: For all methods both the worst-case and 5% quantile performance are lower in uncertainty set 2 than in uncertainty set 1, while mean performance is approximately equal.Note that the maximum violation percentage remains unchanged between the two scenarios because the bounds on ρ remain unchanged, and this is the only parameter that influences the OAR constraint violation.Furthermore, it is noteworthy that the average stage-1 decisions remain unchanged for all methods as well.This underlines the limited influence of the uncertainty interval for τ .The analyses in Section 3 (and 4 for inexact data) demonstrate that the main role of τ is to determine the optimal N * 2 (ρ, τ ), and that ρ plays a more important role.

Influence of multiple worst-case optimal solutions
To see the difference in mean performance between the multiple worst-case optimal solutions, we compare PARO solution found by the ARO method to the ARO solution that performs worst in the two auxiliary scenarios.Table 4 shows the results.The worst-performing ARO solution has a considerably higher stage-1 dose.This implies that (for the current parameter settings) delivering a high stage-1 dose does not allow as much adjustment possibilities in stage 2 as a low stage-1 dose, but it does allow for adjustments to reach the worst-case optimum.Relative to the results of Table 3, the difference between the best and worst ARO solution is considerable: the worst-performing ARO solution is worse than the RO-FH solution.To see the influence of the order of auxiliary scenarios for ARO, Table 5 additionally shows the results for ARO with (ρ aux-max , τ aux-max ) as the first auxiliary scenario (denoted ARO aux-max ).For the other methods the order of auxiliary scenarios did not influence results.Results indicate that, in that case the resulting PARO solution delivers 3.60 Gy per fraction in stage 1; this is a considerable difference to the 1.50 Gy per fraction for ARO aux-min .The mean tumor BED is higher for ARO aux-min ; this is probably due to the fact that the majority of scenarios in Z were contained in Z min , for different parameter settings this may be opposite.

Out-of-sample performance
To investigate the out-of-sample performance of the methods, we assume a uniform distribution for (τ, ρ) over a larger set than Z.We can write Z as where (ε ρ , ε τ ) is the maximum deviation from the nominal scenario (ρ, τ ).This allows us to define where c > 0 is a parameter.We assume a uniform distribution over the new set Z c .If c = 1, we have Z c = Z, so we sample exactly from Z.If c > 1, we sample from an interval that is c 2 times as large as Z (c times larger for both τ and ρ).For c = 2 we obtain the results in Table 6.The stage-1 dose d 1 is the same as in Table 3 for all methods except PI, because PI is the only method that is aware that the samples are not taken from uncertainty set Z but from Z 2 .For NOM, both the mean and maximum violation percentage has increased considerably.All other methods were able to deal with the out-of-sample realizations and did not have any OAR constraint violations.
None of the methods (other than PI) is worst-cast optimal.This implies that the worst-case scenario of the experiments is a scenario outside of Z, because RO, RO-FH and ARO are optimal for the worst-case realization in uncertainty set Z. RO-FH has marginally better worst-case performance than ARO.The static methods NOM and RO have poor worst-case performance compared to the adjustable methods, which indicates bad performance of the static methods on scenarios outside of Z. Due to larger sampling space (the area of Z 2 is four times the area of Z) the difference between mean and worst-case performance is much larger than in Table 3 for both scenarios and all methods.
In terms of average performance, ARO and NOM-FH stay close to PI, and the difference of these methods with RO-FH is larger than in Table 3. Mean tumor BED for ARO is 8.0% higher than for RO-FH and 37.2% higher than for RO in uncertainty set 1.All together, the results from Table 3 show that, while all adjustable methods are worstcase optimal and do not lead to constraint violations, some outperform others in non-worst-case scenarios.Results from Table 6 imply that this difference is largest for realizations outside Z.

Results inexact biomarker observations
In the problem with in inexact biomarker information (IDP), we do not obtain the true parameter values (ρ, τ ) after N 1 = 4 fractions, but only an estimate (ρ, τ ).As discussed in Section 4, we specify a new uncertainty set Ẑ such that (ρ, τ ) − (ρ, τ ) ∈ Ẑ.Let DQ ∈ [0, 1] indicate the data quality.Then we set Ẑ such that the width of the new uncertainty intervals for τ and ρ is (1 − DQ) times the width of the original intervals [τ L , τ U ] and [ρ L , ρ U ].That is, DQ • 100% can be interpreted as the percentage by which the uncertainty intervals can be reduced due to the observation.The relation with the accuracy parameter r ρ (or similarly r τ ) is given by Note that even DQ = 0 has some value as the new interval is centered around the observation, which already cuts off part of the original uncertainty set Z. We pick DQ = 2/3, so the obtained information after fraction N 1 reduces the size of the interval by 66.7% around the new observation.Variations for DQ are considered in Section 5.5.For all patients the required auxiliary scenarios for the method of Section 4.3 can be found.Table 7 shows the results.First, we consider the nominal methods (NOM and NOM-FH).While NOM-FH performs better than NOM on all performance measures, it performs considerably worse than in the cases with exact biomarker information.Because any acquired biomarker information is now inexact, stage-2 adjustments are based on inexact data which makes them less valuable.For NOM-FH this leads to (potentially large) OAR constraint violations of up to approximately 9% in terms of OAR BED.Hence, in case of inexact biomarker information, robustness needs to be taken into account from the start of stage 1.
In terms of worst-case performance, RO, RO-FH and ARO are equal to PI in both uncertainty sets.This indicates that, although not theoretically guaranteed, ARO finds an ARO solution in all considered scenarios.The mean performance of RO-FH and ARO is further away from PI than in the case with exact biomarker information.This is as expected, as due to inexact observations the possibility for ARO and RO-FH to make adjustments is less valuable, whereas PI is not influenced by this.The same reasoning explains why RO is closer (but still worse on all performance measures) to RO-FH than in Table 3. ARO is the only method (together with PI) that has a different stage-1 decision than in the case with exact biomarker information (Table 3).This is because it is the only method that takes inexactness of biomarker data into account at the start of stage-1.The average stage-1 dose d 1 differs considerably between ARO and RO-FH, whereas their worst-case performance is equal on average (and equal to PI).This demonstrates the existence of multiple worst-case optimal solutions.Whereas optimizing worst-case optimal solutions for ARO over two auxiliary scenarios does not guarantee a PARO solution (see Section 4.3), results in Table 7 indicate that it does yield solutions that performs slightly better on average than RO-FH (for the current order of auxiliary scenarios).

Optimal moment of biomarker acquisition
We extend the previous experiment by varying the information point N 1 , and assuming a (hypothetical) mathematical relationship between information point N 1 and the data quality parameter DQ.With N max the maximum number of fractions, we consider the following three data quality functions: Hence, DQ 1 , DQ 2 and DQ 3 describe a convex, linear and concave relationship between observation moment and data quality, respectively.Figure 6 shows the graphs of the three functions.
Whether DQ 1 , DQ 2 or DQ 3 is most realistic depends on the specific biomarker(s) that is/are used.For some biomarkers the data quality may greatly improve in the first few fractions, with a diminishing improvement in later fractions, while for others, e.g., imaging biomarkers such as [18]F-FDG and [18]F-FLT, the data quality is poor at the first couple of fractions and only increases substantially in later fractions.In practice, some biomarkers, e.g., radiographic information (due to, for instance, interference from acute inflammation in the lung), may also exhibit a decreasing data quality for high values of N 1 .Such behavior is rare, and as such not considered here.Figure 7 shows the results for the numerical experiments where we varied the information point N 1 from 0 to N max − 1 for data quality functions DQ i (N 1 ), i = 1, 2, 3 and uncertainty sets 1 and 2 (see Table 2).It is important to note that as N 1 increases past N min = 5, this also increases the minimum number of fractions correspondingly.
First, we note that in all data quality functions and both uncertainty sets the method NOM-FH has high mean tumor BED, but at the cost of high OAR tolerance violations.Hence, if robustness is not taken into account at the start of stage 1, even a folding horizon method may overdose the OARs.
Secondly, the performance of RO (the only static method) does not depend on N 1 , i.e., its curves in Figure 7 are approximately horizontal.This is as expected, as it does not use the biomarker data at all.Both ARO and RO-FH perform better than RO due to the use of biomarker information, and ARO performs better than RO-FH.In case of concave data quality (most optimistic case, Figures 7e and 7f), these differences are most clearly observable.This implies that there is value in (i) adapting based on acquired information, (ii) taking into account that we can adapt later on when planning the stage-1 dose.
All figures demonstrate that the moment of biomarker acquisition influences the average performance of all methods (except static RO).The optimal moment of biomarker acquisition for PI is N 1 = 4.This is because the minimum number of fractions is N min = 5.Hence, if hypofractionation is optimal we can deliver one more fraction with high dose, and deliver a low dose in the first four fractions.If hyperfractionation is optimal we can deliver 16 more fractions (and get the total maximum of 20) with low dose.Note that PI does not actually make an observation at N 1 , the choice of N 1 solely determines what non-uniform treatments PI can deliver due to the constant dose per stage restriction.Furthermore, having N 1 > 4 also implies that we force the use of more fractions than N min = 5, which is disadvantageous for those scenarios where hypofractionation is optimal.
Figures 7e and 7f show that in case of concave data quality the optimal moment for biomarker acquisition is between N 1 = 6 and N 1 = 10.Hence, due the inexactness of the biomarker data, it is optimal to postpone the biomarker acquisition (compared to PI).The same effect is observable if the data quality decreases even further.In case of linear data quality (Figures 7c and 7d) the Figure 7: Results for varying the information point N 1 from 1 to N max − 1, for data quality functions DQ i (N 1 ), i = 1, 2, 3 and uncertainty sets 1 and 2 (see Table 2).The OAR BED constraint violation (%) of NOM-FH is measured against the right axis.Note that the mean tumor BED is displayed (measured on the left axis), while the methods maximize the worst-case tumor BED.
optimal moment of biomarker acquisition is approximately between N 1 = 12 and N 1 = 14.In case of convex data quality (Figures 7a and 7b) the optimal moment of biomarker acquisition is between N 1 = 16 and N 1 = 20.As the maximum number of fractions is N max = 20, the optimal moment of observation is only a few fractions before the maximum number of fractions.Observing later would severely limit the possibilities for adaptation.At N 1 = 16, we force the use of at least 17 fractions, also in the scenarios where it is optimal to hypofractionate.Nevertheless, current results indicate that such treatments are still preferred over a hypofractionated treatment where the dose per fraction is based on very uncertain data.
To put the performance of all methods in perspective, we consider the results of PI.The mean tumor BED for PI changes considerably with N 1 .There are two factors that influence this: (i) as N 1 increases the minimum number of fraction increases, which is disadvantageous for those scenarios where hypofractionation is optimal, (ii) as we force a constant dose per stage, the moment of observation also determines what type of nonuniform fractionation schemes we can deliver.These two factors together lead to a difference in mean tumor BED of around 7 or 8 Gy between N 1 = 4 and N 1 = 19.In comparison, the maximum difference between ARO and static RO is at most 9.8 Gy in mean tumor BED (concave data quality increase, uncertainty set 2).
Overall, there is a trade-off between early observation for fast adaptation and good data quality.All plots for ARO and RO-FH in Figure 7 have a relatively flat top, indicating that there is a range of preferred values for N 1 , rather than one clear best choice.Nevertheless, results suggest that rate at which data quality improves has an effect on the optimal moment of biomarker acquisition.

Concluding remarks
In case biomarker data is exact, the current model (EDP) has considerable 'adjustment space', so it is not necessary to take robustness into account from the start.It is in particular noteworthy that NOM-FH performs equal to ARO and better than RO-FH on average.In general this suggests that static RO, even with a folding horizon, can be too conservative for two-stage problems with exact data and sufficient adjustment space.
In case biomarker data is inexact, taking into account robustness from the start of stage 1 is necessary to prevent OAR constraint violations in the current model (IDP).Another consequence of inexactness of acquired biomarker data is that it postpones the optimal moment of biomarker observation.By investigating several types of temporal changes in biomarker data quality, the presented numerical results provide insight in the value of biomarker information quality and the influence of biomarker acquisition time.This can guide the effort in acquiring useful biomarkers.
Next to the good numerical performance of ARO, one of its main benefits is that is provides more insight into the optimal stage-2 decisions.In the current application, the optimal stage-2 decisions can be determined.Furthermore, in case of exact biomarker data ARO guarantees a PARO solution, which is valuable due to the substantial difference in mean performance between different ARO decisions.Although not guaranteed to be PARO, also in case of inexact biomarker data ARO provides good performing solutions.To the best of our knowledge, we are the first to introduce the PARO concept.Further research is needed to develop ways to construct PARO solutions for general (linear) adjustable robust optimization.
The current setting can be extended in several ways.In practice the tumor and OAR α/β values would have to be estimated from actual biomarkers (e.g., imaging, blood-based biomarkers, genotyping), which can be incorporated in the model.Furthermore, the approach can be extended to heterogeneous tumor response (different α/β ratios for different tumor subvolumes).To the best of our knowledge, we provide the first application of ARO in RT.Other RT applications may also benefit from ARO, such as re-optimization to account for organ motion or setup errors, or optimization using the MR-linac.
Combining the result for the first and second derivative, we obtain with unique minimizer g(0, 0, N 1 + N 2 ; ρ); • Function f (d 1 , N 2 ; ρ, τ ) is strictly concave in d 1 if τ < ρσ, with unique maximizer g(0, 0, N 1 + N 2 ; ρ); If τ = ρσ, we can rewrite f (d 1 , N 2 ; τ σ , τ ) to The following lemma provides information on the existence and location of these intersection points.We consider only those values for d 1 where function f (d 1 , N 2 ; ρ, τ ) is finite for all (ρ, τ ) ∈ Z.Let has the following real roots for d 1 on the interval [0, d UB ]: Proof.Both parts of the lemma are proved individually.
Proof.We first prove part (a), after that the result of part (b) is easily obtained.It holds that We distinguish 2 cases: In this case, the delivered dose exceeds the dose that is used to set the BED tolerance, which is only possible if the number of fractions N ′ is strictly larger than the reference number of fractions T .Condition (B.53) clearly holds, so g • ϕD ≥ d ′ N ′ σ.In this case, squaring (B.53) on both sides and simplifying results in  Putting all of the above together yields the required result for g, i.e., Lemma 3a.
The partial derivative of f w.r.t.ρ is given by Hence, the sign of the partial derivative of f w.r.t.ρ is equal to the sign of the partial derivative of g w.r.t.ρ.This yields the result of Lemma 3b.
has exactly the same objective value.Plug in definitions (B.56), rewrite the first equation to eliminate g, and plug this in the second equation to get where to obtain the second line we plugged in the definitions (B.56).The '+' solution to (B.62) returns d ′ 1 = d 1 , and the '-' solution returns Furthermore, using the same argument, Therefore, it holds that d ′ 1 ∈ W (ρ, τ ).Now, d 1 and t(d 1 ; ρ, τ ) both solve (B.61) so their objective values are equal.Finally, this implies that t(d 1 ; ρ, τ ) is a solution that is unequal to d 1 that solves (B.60).It is also the unique solution because if τ = σρ, function f (d 1 , N 2 ; ρ, τ ) is a strictly convex or concave function according to Lemma 1, so f (d 1 , N 2 ; ρ, τ ) = z for some constant z ∈ R has either 0, 1 or 2 solutions.
For the second part of the proof, let we also end up with a contradiction.This completes the second part of the proof.
In the following lemma, let I(•|S) denote the indicator function for a set S:  We use this to split the domain [0, d UB ] as follows:  This implies that if d 1 ∈ [d − 1 (η), d − 1 (η−1)], it is optimal to deliver either η or η−1 fractions.If we deliver η fractions, the restricting worst-case scenario is L , τ L ) and the value f is above K for the scenario (ρ U , τ L ).If we deliver η ′ > η fractions, the value for the scenario (ρ L , τ L ) decreases, while the value for the scenario (ρ U , τ L ) increases even further.Hence, delivering η ′ > η fractions cannot be optimal.Similarly, delivering less than η − 1 fractions cannot be optimal.Therefore, if d 1 ∈ [d − 1 (η), d − 1 (η − 1)] it is optimal to deliver either η or η − 1 fractions.This implies that on the interval Function p(d 1 ) is a piece-wise function, on intervals defined by S min and S max it is convex and concave, respectively.On any interval S η function p(d 1 ) is the maximum of a concave and convex function.

C Details study setup
The data set contains data from 17 liver patients treated with proton therapy at Massachusetts General Hospital (Boston, USA).All patients were treated with either 5 or 15 fractions.The data is a subset of the data set used in Perkó et al. (2018), for more details on patient selection we refer to that paper.We assume that the target is the gross tumor volume (GTV), and the single dose restricting OAR is the healthy liver itself.
We generate a set of phantom patient parameters (tuples (σ, ϕ, D)) of size S POP = 50.A larger and more heterogeneous phantom patient sample allows for a better comparison between ARO and the benchmark methods.Suppose the OAR consists of n voxels, with s i the dose sparing factor of voxel i.Let s denote the relative mean OAR dose, i.e., the mean OAR dose sparing factor: which is a measure for the dose heterogeneity, for details see Perkó et al. (2018).That paper fits the relative mean OAR dose to the OAR dose shape factor using a two-term power series model:

Figure 2 :
Figure 2: Split of uncertainty set Z according to (18).The circles indicate the locations of the candidate worst-case scenarios for (14).
Figure 4: The uncertainty set Z (solid lines) for (ρ, τ ) is split into Z min ID , Z int ID , Z max ID , according to (29).Subset Z int is the area between the dotted and dash-dotted curves.If (ρ, τ ) ∈ Z int both N min 2 and N max 2

Figure 6 :
Figure6: The biomarker data quality is a function of the number of delivered fractions N 1 after which it is acquired.We consider three functions DQ i (N 1 ), i = 1, 2, 3.

Figure 10 :
Figure10: The graph of the fit (C.84), which relates the mean OAR dose sparing factor to the shape for the liver patient population (solid line), along with 10% bounds (dashed and dotted lines).
s) −b + c,(C.84)where the fit to the data yields a = 0.339, b = 1.304 and c = 0.940, with an adjusted R 2 of 0.987.Figure10displays the graph of (C.84), along with 10% bounds.The set of phantom patient parameters is generated based on (C.84).The mean OAR dose sparing factors in the data set Mizuta et al. (2012) to the OAR in fraction t is given by σd t .For notational convenience, let τ (for tumor) and ρ (for risk) denote the inverse α/β ratio of the tumor and OAR volume, respectively.Mizuta et al. (2012)consider the problem of minimizing OAR BED subject to a lower bound BED pres T on tumor BED.The number of fractions N is restricted to be at most N max .The problem reads min The generalized dose-sparing factor σ denotes the fraction of mean tumor dose that the OAR receives on average.This enables expressing the OAR dose in terms of Schematic overview of the model of Section 3.There are 3 variables: d 1 , d 2 and N 2 .First, we deliver N 1 fractions of dose d 1 per fraction.After this, we observe (ρ, τ ).Subsequently, we deliver N 2 fractions of dose d 2 per fraction.

Table 1 :
Guaranteed solution properties of the five methods.

Table 2 :
Uncertainty sets Z used in the experiments.

Table 3 :
Results for experiments with exact biomarker observations and uniform sampling over Z.All results are averages over a sample of S POP = 50 patients.Per patient all reported statistics are computed over a sample of size S α/β = 200 of uncertain parameters τ and ρ from a uniform distribution over Z.All methods optimize for worst-case tumor BED, which is displayed in bold.

Table 4 :
Comparison between the best (PARO) and worst performing ARO solutions.Per patient all reported statistics are computed over a sample of size S α/β = 200 of uncertain parameters τ and ρ from a uniform distribution over Z.

Table 5 :
Results for ARO for both orders of auxiliary scenarios.Method ARO aux-min optimizes first for (ρ aux-min , τ aux-min ) and ARO aux-max optimizes first for (ρ aux-max , τ aux-max ).Results are displayed for uncertainty set 1, uncertainty set 2 is similar.OAR constraint violations are zero in all cases.

Table 6 :
Results for experiments with exact biomarker observations and uniform sampling over Z c , with c = 2.All results are averages over a sample of S POP = 50 patients.Per patient all reported statistics are computed over a sample of size S α/β = 200 of uncertain parameters τ and ρ from a uniform distribution over Z c .All methods optimize for worst-case tumor BED w.r.t.Z; displayed worst-case tumor BED is w.r.t.Z 2 .

Table 7 :
Results for experiments with inexact biomarker observations and data quality DQ = 2/3.All results are averages over a sample of S POP = 50 patients.Per patient all reported statistics are computed over a sample of size S α/β = 200 of uncertain parameters τ and ρ from a uniform distribution over Z.All methods optimize for worst-case tumor BED, which is displayed in bold.
are the roots of this concave parabola if they are finite.The smaller root, d − 1 (N ′′ ), is finite if and only if N ′′ ≤ T .The larger root, d + 1 (N ′′ ) is finite if and only if N ′ ≤ T .