and analysis of mathematical models
Bj¨orn Mor´ena)and Torbj¨orn Larsson
Department of Mathematics, Link¨oping University, Link¨oping, Sweden
˚Asa Carlsson Tedgrena)
Radiation Physics, Department of Health, Medicine and Caring Sciences, Link¨oping University, Link¨oping, Sweden Medical Radiation Physics and Nuclear Medicine, Karolinska University Hospital, Stockholm, Sweden
Department of Oncology Pathology, Karolinska Institute, Stockholm, Sweden
(Received 9 August 2019; revised 12 November 2020; accepted for publication 22 January 2021; published xx xxxx xxxx)
Treatment planning in high dose-rate brachytherapy has traditionally been conducted with manual forward planning, but inverse planning is today increasingly used in clinical practice. There is a large variety of proposed optimization models and algorithms to model and solve the treatment planning problem. Two major parts of inverse treatment planning for which mathematical optimization can be used are the decisions about catheter placement and dwell time distributions. Both these problems as well as integrated approaches are included in this review. The proposed models include linear penalty models, dose–volume models, mean-tail dose models, quadratic penalty models, radiobiological models, and multiobjective models. The aim of this survey is twofold: (i) to give a broad overview over mathematical optimization models used for treatment planning of brachytherapy and (ii) to pro-vide mathematical analyses and comparisons between models. New technologies for brachytherapy treatments and methods for treatment planning are also discussed. Of particular interest for future research is a thorough comparison between optimization models and algorithms on the same dataset, and clinical validation of proposed optimization approaches with respect to patient outcome. © 2021 The Authors. Medical Physics published by Wiley Periodicals LLC on behalf of American Associa-tion of Physicists in Medicine. [https://doi.org/10.1002/mp.14762]
Key words: catheter placement, dose planning, dwell time optimization, intensity modulated brachytherapy, high dose-rate brachytherapy, mathematical optimization, radiobiology
1. INTRODUCTION
The use of mathematical modeling and optimization is today part of routine clinical treatment planning for intensity-modu-lated radiotherapy (IMRT). It is also increasingly used in high dose-rate brachytherapy (HDR BT), especially for prostate cancer. Furthermore, the use of optimization for treatment planning is an active field of research for both external beam radiotherapy (EBRT) and BT.
The subject of this review is the application of mathemati-cal modeling and optimization of treatment planning of HDR BT. It is written for experts of both medical physics and mathematical optimization, as the topic is translational between these fields of sciences. In this section, we give an overview of the review, with its aim, scope, and delimitations. HDR BT is described in Section 2 along with important con-cepts of dose plan evaluation in Section 3. The intention is to provide a brief introduction for those who are not familiar with the common concepts of BT. The role of mathematical optimization is discussed in Section 4, an overview of mathe-matical models for treatment planning is given in Section 5, and analyses of these are presented in Section 6. The purpose of the review is to provide an overview of the proposed and used optimization models, and to provide a better understand-ing of these and on how the choice of model can affect the resulting treatment plan.
Brachytherapy is used to treat a number of different cancer diagnoses; in the reviewed articles on application and model-ing of treatment plannmodel-ing the most common application is to prostate cancer, studied in approximately 70% of the articles. The second most studied treatment site is gynecological can-cer (can-cervical cancan-cer), studied in less than a third of the arti-cles. Although prostate cancer is the most studied treatment site, optimization models and methods are to a large extent applicable also for other treatment sites and more than one site is investigated in several articles.
The treatment planning problem is fundamentally multi-objective, with aims for each anatomical structure, including tumor or target, healthy tissues, and organs at risk (OAR). These aims are often in conflict as a higher tumor dose corre-lates with higher doses to OAR. From an optimization per-spective, there is not a single valid model for the treatment planning, and clinical evaluation criteria and approximations thereof are included in models in a broad variety of ways.
1.A. Scope and delimitations
The review covers the optimization models that have been proposed for HDR BT treatment planning. Models for both dose planning (dwell time optimization) and catheter place-ment are included. References for low dose-rate (LDR) BT and EBRT are generally not included, except when such
model formulations precede HDR BT models, when they have a close resemblance to HDR BT models, or for model formulations for which there are few HDR BT studies.
The only earlier literature review1of optimization models for HDR BT dose planning covers the time period 1990— 2010. Our review extends the earlier in several ways. First, we cover the time period until June 2020. Second, we include models for catheter placement, which is not considered in the earlier review. Third, we include a section with mathematical analyses of properties and relationships of both models and evaluation criteria. Computational algorithms and software for treatment planning are outside the scope of the review. We expect the reader to be somewhat familiar with mathe-matical notation and modeling.
1.B. Abbreviations
In TableI, we introduce the abbreviations that are used for optimization models and evaluation criteria.
2. A BRIEF INTRODUCTION TO BRACHYTHERAPY We provide a brief introduction to brachytherapy for those who are not familiar with it. Readers who are already familiar with BT can skip this section.
In BT, the radiation is delivered by small, sealed radiation sources, of millimeter dimensions, inserted close to or within a target volume (which can be a tumor or an adjuvantly trea-ted volume). BT is used alone or as a boost to EBRT. The placement of the radiation source within the body is in con-trast to EBRT, in which the radiation is delivered from the outside using linear accelerators.
A single source is used in HDR BT, most often of the iso-tope 192Ir. The source is guided through the body using
interstitial preinserted catheters, that is, hollow needles, and/ or intracavitary anatomy shaped applicators, which define the possible locations for the source and are inserted before treat-ment. For simplicity, in our presentation they are all referred to as catheters.
For dose delivery, each catheter is discretized into a num-ber of dwell positions in which the source can dwell. After the catheters have been implanted and the treatment planning is finished, the catheters are connected to a remote after-loader. At delivery, the radiation source steps through each catheter, and at each possible position the source either dwells for a predetermined time, referred to as a dwell time, to irradiate the surrounding tissue, or it passes by to the next position (i.e., the dwell time is zero).
A BT treatment is planned with respect to the specific anatomy of each individual patient. The first step is to acquire images of the treatment site using imaging modalities such as ultrasound, computed tomography, or magnetic resonance imaging. Medical images are generally acquired during (or in close conjunction to) the treatment session, as they must incorporate the inserted catheters, which may alter the anat-omy. The images are used to manually contour and define the anatomical structures of interest, including both the target (or targets) and OAR in its proximity.
There are several interrelated volume definitions regarding the target (or targets),2the two most important being the clini-cal target volume (CTV) and the planning target volume (PTV). The CTV is the volume that should be irradiated, while the PTV contains the CTV and an extra margin because of uncertainties due to movements (e.g., breathing) or other reasons. PTV margins are in general smaller in BT than in EBRT, and the PTV can even be equal to the CTV, because in BT organ motion is less problematic as the radiation source moves together with the irradiated target.
In addition to OAR, an artificial structure of healthy tissue surrounding the PTV can be added, for which restrictions of radiation doses are imposed. The reason for this is to avoid unreasonably high dwell times in the outskirts of the PTV, where distances to the delineated OAR are large. In a prostate cancer study, where no artificial structure of healthy tissue is added, it is described to result in clinically undesired proper-ties of the dose distribution.3
The aim is to irradiate each target with at least a specified prescription dose. The level of the prescription dose depends on the treatment site and the treatment schedule, for example the number of fractions. In our presentation, due to simplic-ity, the target will generally be referred to as the PTV, in sin-gular.
The dose planning, that is, to decide the dwell times, can be conducted either with forward or inverse planning. For-ward planning is a trial-and-error process in which the plan-ner constructs a dose plan manually, which is then evaluated, and if it is not good enough the planner adjusts the dose dis-tribution. This process is repeated until the plan is found good enough for approval, which can be time-consuming. The forward planning is conducted either through manual adjustments of individual dwell times or with graphical aid TABLEI. Overview of the abbreviations used for optimization models and
evaluation criteria.
Abbreviation Meaning Section
BFGS Broyden–Fletcher–Goldfarb–Shanno 5.D
CN Conformation number 3.C
COIN Conformal index 3.C
CVaR Conditional value-at-risk 5.C
DHI Dose homogeneity index 3.C
DI Dosimetric index 3.B
DTMR Dwell time modulation restriction 5.A
DVM Dose–volume model 5.B
EUD Equivalent uniform dose 3.D
gEUD Generalized equivalent uniform dose 3.D
LP Linear programming 5.A
LPM Linear penalty model 5.A
MIP Mixed-integer program 5.B
MTDM Mean-tail-dose model 6.E
NTCP Normal tissue complication probability 3.D
QPM Quadratic penalty model 5.D
from a treatment planning system (TPS). Manual treatment planning for prostate cancer is reported to take up to an hour,4 or 20–35 min,5 and for gynecologic cancer 30–60 min.6
In inverse planning, the starting point is the treatment aims and in the planning step, a dose plan is constructed to achieve or to be as close as possible to these goals, modeled exactly or approximately. Due to the computational complexity of inverse planning, optimization models and algorithms are used to construct the dose distribution.
During the treatment process, the patient is commonly under some form of anesthesia (spinal in general) as the catheters are inserted invasively for most treatment sites. For prostate and cervical cancer the treatment planning is often conducted during the treatment session (“on-line”) in order to plan the treatment using up-to-date anatomical informa-tion. The overall time for the treatment planning is hence an important aspect, and it is beneficial to be able to conduct the planning fast.
3. CLINICAL TREATMENT PLAN EVALUATION The main decisions in treatment planning are the catheter placement (how many to use and where to insert them) and the dwell times. When these have been decided, the dose dis-tribution can be calculated and evaluated.
3.A. Dose calculation and dose points
To calculate the dose, the dose-rate contribution, per sec-ond, from each dwell position is needed. The dose calcula-tions are today generally based on the AAPM TG43 formalism,7–9by which the composition of the patient is con-sidered to be equivalent to water. In reality, patients are of course not water equivalent, but the water approximation of TG43 is generally considered to be good enough in the case of the HDR isotope 192Ir, see, for example, Ref. [10] for details on tissue and applicator heterogeneities in BT. The distribution around the single source is precalculated and scaled with respect to the daily air kerma strength of the radi-ation source. With informradi-ation on source air kerma strength and the dwell times, the total received dose at any point in the contoured structures can be calculated by adding the dose contributions from all dwell positions, which is the dwell time times the dose-rate contribution scaled with the strength of the radiation source.
In the evaluation of dose plans, entities referred to as dose points are used. Dose points constitute a discretization of the contoured structures and each dose point represents a small volume of tissue. One way to generate dose points is accord-ing to Lahanas et al.11
3.B. Dose–volume histogram
An essential concept in clinical treatment guidelines and dose plan evaluation is the dose–volume histogram (DVH),
which plots portions of a structure against dose levels. DVH curves can be presented in a differential version, in which the portion that receives a certain dose level (in Gy) is shown for each dose level, or more commonly in a cumulative version, in which the portion that receives at least a dose is shown for each dose level. In the following, the expression DVH curve refers to the cumulative DVH curve (unless otherwise speci-fied).
Each point on the DVH curve corresponds to a dosimetric index (DI), which can be formulated in two ways. These are denoted Vsxor Dsy, respectively, where x is a dose level, s is a structure, either PTV or part of the PTV, or an OAR, and y is a portion of the structure or a volume (commonly in cc). The DI Vsx is the portion of the volume of the structure that receives at least the specific dose level x, while the DI Dsy is the lowest dose that is received by either a portion of the structure or a volume that receives the highest dose. If no structure is specified for Vx, we refer to the PTV. The actual
values of DIs may depend on factors like image slice thick-ness and the TPS implementation of volume interpolation.12
3.C. Dose homogeneity and conformity
Conformity measures estimate how conformal the pre-scription dose is to the target volume, some of them also include doses to OAR. Homogeneity measures instead con-sider how homogeneous the dose to the target volume is.
A measure that is related to the DVH curve is the dose homogeneity index (DHI),13which is defined for the PTV as
DHI¼V PTV 100 V PTV 150 VPTV 100 , (1)
where 100 and 150 refer to percentages of the prescription dose. This measure takes a value between zero and one. The value of one is the ideal case when the full volume receives at least the prescribed dose and no portion of the volume receives more than 150%.
The conformation number (CN)14 is a conformity mea-sure. In addition to the dose to the PTV, it also takes the por-tion of the OAR that receives a dose that is higher than the prescription dose into account. The portion of the PTV that receives a dose that is too high is, however, not included. The CN is calculated as
CN¼ VPTV100 PTVref Totref
, (2)
where PTVref corresponds to VPTV100 but is expressed as a
vol-ume, and Totref is the total volume (with both PTV and OAR
included) which receives at least the prescription dose. The conformal index (COIN)15 is an extension of the CN which in addition has an extra term which is specific for OAR. The conformal index is defined as
COIN¼ VPTV 100 PTVref Totref Y s∈SO ð1 Vs 100Þ, (3)
where SOis the set of OAR. Both CN and COIN takes a value
between zero and one, with one being the best possible value in the two measures; this corresponds, for both measures, to the whole PTV receiving the prescription dose and no portion of any OAR receiving a dose that is higher.
See Refs. [16] and [17] for reviews of conformity and homogeneity measures.
3.D. Radiobiological concepts
In contrast to the DIs and the homogeneity and conformity measures, which are based solely on the physical dose, radio-biological indices are meant to explicitly capture the radiobi-ological treatment effect of a dose distribution. One such index is the tumor control probability (TCP) which is an esti-mate of the probability of local tumor control.
The corresponding radiobiological index for an OAR is the normal tissue complication probability (NTCP), which estimates the probability that complications occur in an OAR. A summary of how NTCP is used is given in Ref. [18]. Another radiobiological concept is the equivalent uniform dose (EUD).19It captures the radiobiological treatment effect on a tumor with the purpose of comparing the treatment effects from inhomogeneous dose distributions. The EUD value is meant to correspond to the homogeneous dose distri-bution that gives the same treatment effect as the evaluated inhomogeneous dose distribution. The EUD was developed with the generalized EUD (gEUD),20which can be used also for OAR. One way to define the gEUD is given in Ref. [21] as gEUD¼ 1 N∑ N i¼1Doseai 1=a , (4)
where N is the number of dose points in the structure, Dosei
is the dose to dose point i, and a is a parameter for radiosensi-tivity, which depends on the type of structure considered. Compared to TCP and NTCP the gEUD has fewer parame-ters, in this formulation only one, but this parameter is sup-posed to capture both characteristics of the treatment protocol and the type of structure.
Finally, a concept that considers both tumor control proba-bility and risk for complications in OAR is p+.22It estimates the probability of tumor control without any complications.
3.E. Clinical treatment guidelines
There are established guidelines for the different steps in the treatment planning. HDR BT guidelines for prostate cancer23–25 include planning aims for the PTV and OAR, which are typically expressed in terms of DIs. For example, the DI VPTV
100 should be at least 90% with the recommended
goal that it should be above 95%. The guidelines24also rec-ommend reporting a value that is equivalent to the DHI, see Eq. (1).
Corresponding clinical BT treatment guidelines are estab-lished also for other treatment sites; see for example the
guidelines for vaginal cancer,26cervical cancer,27breast can-cer,28and head and neck cancer.29
4. USING OPTIMIZATION IN TREATMENT PLANNING
A general formulation of the treatment planning problem is to decide“how to deliver a high enough dose to the tumor while sparing OAR as much as possible.” This formulation is in itself vague, with expressions such as“high enough” and “sparing OAR as much as possible,” and since these aims are in conflict, compromises are necessary. There is however no objectively best compromise between these aims. Further-more, even though there are goals for the clinical evaluation criteria in the treatment guidelines, it is hard to a priori decide the best treatment plan for a specific patient, because the anatomy and the catheter placement determine which dose distributions are achievable; hence, the compromise between aims are patient dependent. Moreover, there is also uncer-tainty in for example the spread of the tumor cells and in patient radiosensitivity.
The treatment planning problem can therefore be seen as a soft optimization problem. From a modeling point of view, there is no obvious objective function and this also becomes apparent from the broad variety of dose planning models pre-sented in the next section. The lack of a standard problem for-mulation makes treatment planning different from many other mathematical optimization problems, which commonly have well-defined objective functions and constraints that define the properties of feasible and optimal solutions.
The softness of treatment planning is also apparent from the fact that feasibility of a dose plan is not well defined and the presented models consider feasibility in different ways. Furthermore, even though there are specified goals for the clinical evaluation criteria, there are cases which are extra dif-ficult and for which not all goals are attainable. For such a case, the treatment goals must be relaxed to allow any dose plan at all. The other way around, it can also occur that it is possible to find a dose plan with considerably better values than the goal values.
Nevertheless, despite the inherent difficulties in using optimization for HDR BT treatment planning, all the pro-posed models are capable of delivering dose plans that are clinically relevant, possibly after including a trial-and-error process of parameter tuning or manual adjustments of the resulting dose plan.
The optimization models for dose planning that are pre-sented in Sections 5.1–5.7 have the dwell times as the main decision variables and the objective function is based on eval-uation of the received dose at each dose point. In the models for dose planning in Section 5, the catheter placement is pre-determined but these decisions can also be made in combina-tion with the dwell times. When the dose planning is combined with the catheter placement the objective consid-ered is often to use as few catheters as possible. Such models are discussed in Section 5.10.
From a mathematical modeling point of view, there is a large difference between HDR BT and LDR BT, in which seeds are implanted permanently or for a prolonged time in the body. In HDR BT, the dwell times are modeled with con-tinuous variables and the catheters with binary variables, while in LDR BT the seed placements are modeled with bin-ary variables. There are also large differences between treat-ment planning optimization of EBRT and HDR BT. In EBRT, there are both more degrees of freedom and con-straints regarding the dose delivery, and these characteristics put limitations on which modeling techniques that are practi-cally feasible. For HDR BT there are more optimization mod-els which are possible to use in practice.
4.A. Conflicting aims
The main treatment goal is to find a dose plan which is a good trade-off between what is attainable in terms of PTV coverage and in sparing of OAR. Trade-offs are not only pre-sent between different structures but can also be relevant within a structure, for example the trade-off between deliver-ing a dose that is sufficiently high but not too high to a struc-ture, or trade-offs between substructures with different priorities within the PTV (due to for, e.g., the cancer spread). In addition to these main goals, dose homogeneity is also of interest and included in treatment guidelines.
The treatment planning problem is hence fundamentally multiobjective. The problem is, however, handled in a num-ber of ways which do not fully consider it as multiobjective. The single-objective models in the literature are primarily based on one of the following reformulations of multiobjec-tive models: (a) to consider a weighted sum of the objecmultiobjec-tive terms, (b) to put hard constraints on all but one of the objec-tive terms, (c) to minimize the penalty for the worst objecobjec-tive term, as a pessimistic approach, or (d) to consider the objec-tives in a priority sequence, optimizing one objective at the time, and then restating it as a constraint.
A property in common for these reformulations is that the result is a single solution (dose distribution) while the result from a true multiobjective approach typically is a set of solu-tions. Multiobjective approaches are computationally more demanding but have the advantage of allowing the planner to get a better understanding of the possible and necessary trade-offs. Multiobjective approaches to dose planning are discussed in Section 5.6; some of these approaches are based on single-objective models which are solved repeatedly. The literature review on IMRT dose planning found in Ref. [30] is focused on multiobjective formulations and decision-mak-ing.
The multiobjective aspect makes it difficult in general to compare optimization algorithms and their results, and also with clinical plans. Unless one plan is better in all evaluation criteria, a reason that one plan is preferred over another may be because of different priorities. Moreover, to compare the results between multiobjective models, the quality of a set of dose plans should be assessed and compared, which is not straightforward.
The components of the optimization models presented in our review can be divided into three categories depending on how closely related they are to the evaluation criteria in the clinical treatment guidelines. There are components that: (a) are explicitly based on an evaluation criterion, (b) approxi-mates or in some other way capture a clinical evaluation crite-rion, and (c) are not based directly on any clinical evaluation criterion; examples of components of different types are given in the following sections. The different ways of considering clinical evaluation criteria in the modeling are typically related to a trade-off between computing time and solution quality, since modeling approaches that are directly related to clinical evaluation criteria tend to be more computationally expensive.
4.B. Why use optimization?
Manual planning is a time-consuming task of trial-and-er-ror and the usage of optimization models has the potential of saving time, which is beneficial for the patient. Also treatment plan quality can be improved with the use of optimiza-tion.6,31–36 In one study,34 planners were presented with a manually constructed dose plan and five dose plans which were constructed with an optimization model. It turned out that they preferred an optimized dose plan in 53 out of 54 cases. Nevertheless, to fully utilize the potential of optimiza-tion models for dose planning it is important to further develop the dose planning models to include all relevant aspects that are considered in manual planning as well as radiobiological criteria. Further, manual planning requires much experience5 and the use of optimization can make the planning process less dependent of individual and varying skills.
The importance of, and also need for, optimization methods becomes more apparent the more degrees of free-dom there is in the treatment planning. An example of this is the emerging technique of intensity-modulated brachytherapy (IMBT) in which rotational shields are used. The number of possible dwell times are then highly increased, making optimization a necessary tool. A case that can easily be handled with optimization is when it is especially important that a part of the PTV receives a dose that is high enough, or when it should receive a higher dose than other of its parts.
One reason for further developing optimization models and algorithms is the limited amount of computing time available for the treatment planning. Since the treatment plan-ning is often conducted when the patient is under anesthesia, any decrease in the planning time is beneficial. There are no hard limits on how much computing time an optimization algorithm is allowed to use, but finding a dose plan within a few minutes is a reasonable aim. Included in the allowed computing time should also be the time needed for tweaking parameter values (or trying class solutions) and manual adjustments, both which depend on the optimization model. Furthermore, even if an improvement in computing times might not be of direct clinical importance, the increased speed could be invested in solving an improved model, for example with an increased number of dose points, or
allowing the solver to run for a longer time, possibly finding better solutions.
Convexity is an important concept when analyzing opti-mization models and for assessing their computational diffi-culty. For convex models, each local optimum is also a global optimum, hence global optimality is determined by only local information and it is therefore easier to find a global optimal solution and to validate its optimality. Convex models are generally considered to be easy to solve, while non-convex models are sometimes very hard to solve. An important cate-gory of non-convex models are those that contain binary vari-ables. Examples of such models are dose planning models for LDR BT, and HDR BT models that include catheter place-ment decisions. See Appendix A for formal definitions and results regarding convexity.
5. MATHEMATICAL MODELS
In this section, we present an overview of the proposed mathematical optimization models for HDR BT dose plan-ning. Not only prostate cancer is the most frequently studied treatment site but also cervical cancer is commonly studied.
TableIIintroduces the notations in common for the pre-sented models. Subscript T refers to the PTV and O to OAR. The set S contains all structures, including both the PTV and OAR. Optimization models that are formulated for a single PTV and a single prescription dose can be generalized to sev-eral PTVs and prescription doses, for example by considering a weighted sum objective for the PTVs. For the purpose of generality, some optimization models contain objectives which strive to irradiate all dose points with at least a speci-fied level of radiation (Ls); for an OAR s∈S such an objective
might not be relevant and can be modeled by setting Ls¼ 0.
5.A. Linear penalty models
The optimization model most used in clinical practice, and implemented in several commercially available TPSs, gives a penalty for each dose point at which the dose is outside a spec-ified interval. This penalty increases linearly with the devia-tion from the interval. In HDR BT such a method, called IPSA (inverse planning simulated annealing), was first proposed,37 in which the model was solved heuristically with simulated annealing (see e.g., Ref. [38] for a description of simulated annealing). The same model was formulated as a linear pro-gram,39which made it possible to find a globally optimal solu-tion. This linear programming (LP) model is referred to as the linear penalty model (LPM). Another linear penalty model was also proposed;40instead of giving a penalty for each dose point it minimizes the maximum deviation from a specified dose level.
5.A.1. Optimization model
TablesIIandIIIintroduce the notation used for the LPM. The solid line in Fig.1shows the linear dose–penalty rela-tionship for one dose point. There is no penalty when the
dose is within the specified interval (that is, between Ls and
Us), and outside this interval the penalty increases linearly.
The linear penalty model is mathematically formulated as fol-lows. min ∑s∈Swl s∑i∈Psx l iþ ∑s∈Sw u s∑i∈Psx u i (5a) subject to ∑j∈Jdijtj≥ Ls xli i∈Ps, s∈S (5b) ∑j∈Jdijtj≤ Usþ xui i∈Ps, s∈S (5c) xli≥ 0 i∈Ps, s∈S (5d) xu i≥ 0 i∈Ps, s∈S (5e) tj≥ 0 j∈J (5f)
Constraints (5b) and (5d) in combination with the objec-tive function ensure that the penalty variables for the lower dose bounds take the correct values, while constraints (5c) TABLEII. Indices, sets, parameters, and variables in common for the opti-mization models.
Indices
i Index for a dose point j Index for a dwell position s Index for a structure Sets
S Set of structures, including PTV and OAR
SO Set of OAR, including an artificial structure of healthy tissue,
SO⊂S
Ps Set of dose points in structure s∈S
P Set of all dose points, P¼ ∪s∈SPs
PT Set of dose points in the PTV
J Set of dwell positions Parameters
dij Dose-rate contribution from dwell position j∈J to dose point
i∈P
Ls Prescription dose or a lower dose bound for structure s∈S
Us Upper dose bound for structure s∈S
Ms Maximum allowed dose for structure s∈S
ws Non-negative penalty parameter for structure s∈S
Variables
Dosei Radiation dose received at dose point i
tj Dwell time for dwell position j∈J
TABLEIII. Notation used for the linear penalty model.
Parameters wl
s Non-negative penalty for dose being too low in structure s∈S
wu
s Non-negative penalty for dose being too high in structure s∈S
Variables xl
i Penalty variable for dose being too low at dose point i∈Ps, s∈S
xu
and (5e) in combination with the objective function ensure that the penalty variables for the upper dose bounds take the correct values. In this model, any combination of non-nega-tive dwell times, tj, j∈J, define a feasible solution.
The LPM is connected to clinical practice through the use of prescription doses, but it does not relate to the clinical evaluation criteria in any way, and neither the penalty param-eters nor the objective function value in itself have direct clin-ical interpretations. The model can thus be regarded as an ad hoc model. Still, practice has shown that it can be useful, because after proper tuning of the penalty parameters, dose plans that are clinically adequate can be obtained within an acceptable time. However, the parameter tuning is typically a manual task and can itself be time consuming.
Because the objective function in the LPM has no clinical meaning in itself, criteria such as dosimetric indices are instead used to evaluate and compare dose plans. Solutions from IPSA and solutions from the LP formulation were com-pared;39 even though the objective function value was improved significantly by using an LP solver (as the LP solu-tion is a globally optimal solusolu-tion) the DIs were not signifi-cantly different. There are examples that two dose plans with similar DVH curves can differ a lot in terms of linear penalty values,41for example, the objective values of two similar plans is shown to differ by a factor of 12. To further strengthen this comparison, the correlation between the objective function value of the LPM and the DIs was studied and shown to be poor.42
5.A.2. Properties and extensions
It has been observed that IPSA produces dose plans with less homogeneous dwell times than manual planning.43Two approaches, to limit the maximum value for the dwell times and to add an extra artificial structure around the catheters, were tried to improve the homogeneity, and they both achieved more homogeneous dwell times, although not as homogeneous as manual planning.43It was shown that inho-mogeneous dwell times is a property of the LPM and the sim-plex method.44(See, e.g. Ref. [45] for an introduction to LP and the simplex method, which is the standard method for
solving LPs.) An optimal solution to an LP solved with the simplex method is always an extreme point (a vertex) in the polytope defined by the linear constraints. Hence is such an optimal solution to the LPM an extreme point. For any extreme point in the LPM there is a maximum number of dwell positions that can be active, that is positions with a pos-itive dwell time, and it is therefore not surprising that solu-tions to the LPM from the simplex method tend to use few active dwell positions compared to manual planning. It is shown44 that an upper bound on the number of active posi-tions is the number of dose points which receives a dose that is either at their upper or lower bound. A remedy by using piecewise linear penalties is proposed,44 which give increased penalties for dose deviations farther from the speci-fied intervals. Such a model is shown to produce dose plans with more homogeneous dwell times and shorter maximum dwell times compared to the LPM.44An example of a convex piecewise linear penalty function can be seen in Fig. 1, where the dotted line corresponds to a penalty function with six seg-ments. It should also be noted that few active dwell positions is not an inherent property of IPSA because a solution pro-duced by IPSA is not necessarily an extreme point.
An extension of the LPM that increases the degrees of freedom is to allow penalty parameters in a structure to be set differently for individual dose points. It was suggested for IMRT46 to be able to put a higher penalty on points in the middle of the PTV, but this has to our knowledge not yet been clinically evaluated. Because of the increased flexibility, the parameter tuning for such a model becomes more complex than for LPM.
5.A.3. Dwell time modulation restriction
The concept of dwell time modulation restriction (DTMR)47 was introduced to make the dwell times more homogeneously distributed. However, their DTMR formula-tion was not mathematically defined. Compared to the model without DTMR, the studied model47gives similar results in terms of PTV coverage (VPTV100), but a reduction in dose to ure-thra, and also more homogeneous dwell times. Three choices of DTMR are suggested and evaluated,48
tj1 ≤ ð1 þ γÞtj2, j1∈J, j2∈Γj1, (6a) tj1 tj2 ≤ θ, j1∈J, j2∈Γj1, (6b) 1 2∑j1∈J∑j2∈Γj1 tj1 tj2 2 ≤ ρ, (6c)
whereΓj denotes the set of dwell positions adjacent to j in
the same catheter, andγ, θ, and ρ are parameters. These for-mulations extend both an LPM and a model which explicitly includes constraints on dosimetric indices (see the next sec-tion). For each of the constraints (6a)–(6c), there is a signifi-cant deterioration of VPTV
100. A component aiming at reducing
the variance in dwell times was added49 to IPSA and its effects were studied,50showing a reduction in PTV coverage of 1.5%–2.5%.
FIG. 1. The solid line is the linear penalty function for one dose point for
each feasible dose level. The dotted line is the piecewise linear penalty func-tion, here with six segments (the horizontal segment included). The dashed function is a piecewise quadratic penalty function.
Adding a constraint to a model imposes a restriction and the objective value of the restricted model can only become worse than that of the original model (or possi-bly remain the same). Results from DTMR studies,47,48,50 are inconclusive about the size of this effect, but there is an evident trade-off between imposing a more restrictive DTMR and keeping the PTV coverage at a high level. Both Refs. [47] and [50] have a parameter between zero and one to control the DTMR, where zero is the least restrictive choice, and this parameter value is suggested47 to be in the range 0.1–0.2. A DTMR is implemented in the treatment planning system SagiPlan but it is sug-gested51 that the DTMR parameter should be zero, corre-sponding to no restriction.
Another way of imposing a DTMR is proposed,52 in which they add an objective term penalizing the two-norm of the dwell times. It is referred to as a Tikhonov regularization. Also with this approach there is a trade-off between the DTMR term and the (quadratic) objective term which penal-izes dose deviations.
5.B. Dose-volume models
Since dosimetric indices are the primary evaluation criteria for dose plans in the clinical treatment guidelines it is a very reasonable approach to explicitly include them in the dose planning model. By doing so we obtain a dose–volume model (DVM). The first optimization model in HDR BT to explicitly include a DI was proposed in Ref. [53], in which a linear penalty model is extended with a constraint on a DI for one OAR. A DVM was, however, studied earlier for EBRT.54 Such models for HDR BT were further developed and studied.41,55–58
The DI is related to concepts in other fields of research; the counterpart in finance is called value-at-risk (VaR),59and concepts similar to DIs also appear in the field of chance con-straints.60
5.B.1. Numerical example
To illustrate how the DIs are calculated a small numerical example of a dose distribution in the PTV is given in TableIV. With a prescription dose of 9 Gy, the portion of the PTV that receives at least the prescription dose, that is VPTV100, is 80%, and the lowest dose to the 10% and 80% of the dose points that receives the highest dose, that is DPTV
10 and D PTV 80 ,
equals 18 and 9 Gy, respectively.
5.B.2. Optimization model
TablesIIandVintroduce the notation used for the DVM. The following model is adapted from Refs. [55] and [56] (but without the DTMR constraint from the latter reference).
max 1 jPTj∑i∈PT yi (7a) subject to ∑j∈Jdijtj≥ LTyi i∈PT (7b) ∑j∈Jdijtj≤ Ms ðMs UsÞvsi i∈Ps, s∈SO (7c) ∑i∈Psv s i≥ τsjPsj s∈SO (7d) yi∈ f0,1g i∈PT (7e) vsi∈ f0,1g i∈Ps, s∈SO (7f) tj≥ 0 j∈J (7g)
The notation j j is used for the cardinality of a set. The objective function value is the VPTV100 value and the set of con-straints (7b) ensures that each binary indicator variable yi
takes the correct value, that is, one if the dose is sufficiently high and zero otherwise. Therefore, to calculate the value of the DI VPTV100 it is enough to count the indicator variables which take the value one, and normalize the objective func-tion with the number of dose points in the PTV to get a value between zero and one. Similarly, the set of constraints (7d) ensures that the value Vs
100 is at leastτ
s for OAR s∈S O. The
set of constraints (7c) ensures that each binary indicator vari-able vs
itakes the value one only if the dose is sufficiently low.
These constraints, referred to as “Big-M” constraints, also impose a hard upper bound, Ms, on the dose received in each dose point. This upper bound can be chosen to be large enough to not cut off any solutions of interest or take the value of a clinically relevant upper bound, if available.
5.B.3. Properties and solutions methods
The DVM is a mixed-integer programming (MIP) model. For an introduction to mixed-integer programming, see, for example, Ref. [61]. Its non-convexities in the field of IMRT were studied,62showing that there can be multiple local min-ima which must be handled (the definition of a local mini-mum is however not clear). The existence of local minima has been further studied.63–65 An equivalence relationship between solutions to the LPM and solutions to the linear pro-gramming relaxation (i.e., instead of being either zero or one, the binary variables can take any value between zero and one) of the DVM has been presented.57
Solving the model directly with a general-purpose MIP solver has proved to be hard53,55,56,58 and due to time con-straints the MIP solver has to be stopped prematurely, before proving optimality. The computational difficulties of this model are not surprising because models with indicator vari-ables and constraints are known to have weak linear
TABLEIV. Small numerical example of a dose distribution.
Dose point i 1 2 3 4 5 6 7 8 9 10
programming relaxations and to be hard to solve.66However, it is argued that Big-M formulations remain the most practi-cal way for solving DI problems (which they refer to as VaR problems).67 A common alternative approach is to use the meta-heuristic simulated annealing to find good results within a short time. This has been studied in Refs. [53,55,
56].
Another approach to handle DVMs is to approximate the Heaviside step functions (corresponding to the binary indica-tor variables) with smooth nonlinear sigmoid func-tions,6,35,68–70
fðDoseiÞ ¼ 0:5ð1 þ tanh β Doseð ð i LTÞÞÞ, (8)
where β is a parameter to control the steepness of the func-tion; see Fig. 2 for an illustration. Good results from this approach compared to manual planning, LPM and IPIP55are shown,6within only very short computing times, in the order of seconds. However, because of the approximations there are no guarantees that the constraints on the DIs are satisfied. Also, it is neither known if the solution is the global optimum nor are any optimistic bounds on the optimal objective value known (as is the case when a MIP solver is used and stopped prematurely). Approximations of the Heaviside step functions have been studied also in IMRT.71They were approximated with concave functions, and in addition to constraints on DIs the model had an objective function which aimed at finding a homogeneous dose distribution. The LP relaxation of the DVM is solved in Ref. [72].
A clear advantage of using a MIP solver to solve an exact formulation is that an optimistic bound on the optimal objec-tive value is available during the solution process. Combining the optimistic bound with the pessimistic bound from the best found feasible solution gives an interval for the optimal VPTV100 value. This is useful because the planning process includes decisions on the trade-offs between multiple criteria, and to make an informed decision it is relevant to know what is attainable for each evaluation criterion. However, with com-puting time as a limiting factor, improvements of the solution methods are necessary to fully capitalize on this potential.
A DVM-based approach which includes uncertainty in the target volume delineation was considered.72 Instead of the standard margin around the CTV to construct the PTV they
considered scenarios corresponding to different delineations of the PTV. This is motivated by the observations73that the extra margins of the CTV lead to dose escalation. The approach72 is formulated as an optimization model which maximizes target coverage for the worst-case scenario.
5.C. Mean-tail-dose models
A measure that is related to the DVH and the dosimetric indices is the conditional value-at-risk (CVaR), which was introduced as a financial risk measure.74In the financial con-text, CVaR can be used to limit the mean loss in a specified portion of the scenarios with the worst outcome. In radiother-apy CVaR has also been referred to as mean-tail-dose (MTD), which we will use in the following. In the context of radiotherapy the MTD either quantifies the mean dose to the α% of the structure that receives the lowest dose or the mean dose to theβ% of the structure that receives the highest dose. To distinguish between them, LCV aRs
αdenotes the MTD for
the portion α of the structure s∈S that receives the lowest dose and UCV aRs
β denotes the MTD for the portionβ of the
structure s∈S that receives the highest dose. It has been shown that MTD measures can be modeled with a linear for-mulation,74either to maximize (or lower bound) LCV aRsαor to minimize (or upper bound) UCV aRsβ. An overview of applications where CVaR have been used can be found in Ref. [75].
Figure3 shows how the LCV aRsα and the UCV aRsβ mea-sures are related to the differential DVH curve. This curve shows, for each level of radiation dose, the portion of the vol-ume that receives this particular dose level. The mean dose in the grey area is LCV aRs
α, and is calculated as LCV aRs α¼ 1 αjPsj∑i∈Ps:Dosei≤ D s 1αDosei: (9)
The mean dose in the black area is UCV aRsβ, and is calcu-lated as UCV aRs β¼ 1 βjPsj∑i∈Ps:Dosei≥ D s βDosei: (10)
The MTD measures can be regarded and used as a convex approximation of dosimetric indices. Even though the MTD
FIG. 2. Shows the Heaviside step function (the solid blue line) and the
smooth nonlinear sigmoid approximation (the dashed line).
TABLEV. Notation used for the dose–volume model.
Parameters
τs Portion of a structure, here used for OAR s∈S O
Variables
yi Binary indicator variable for PTV dose points, i∈PT, that takes the
value one
if the dose is high enough, that is at least LT, and zero otherwise
vs
i Binary indicator variable for OAR dose points, i∈Ps, s∈SO, that
takes
the value one if the dose is low enough, that is at most UT, and zero
is related to the DVH curve and its meaning is tangible, in contrast to the objective of the LPM, it is only an approxima-tion of goals in the clinical treatment guidelines.
5.C.1. Numerical example
To continue the numerical example in Section 5.2.1, one can also calculate the mean-tail-doses. From the dose distri-bution in Table IV we get LCV aRPTV20 ¼5þ62 ¼ 5:5 Gy (i.e., the mean dose to the 20% dose points that receive the lowest dose is 5.5 Gy) and UCV aRPTV10 ¼181¼ 18 Gy.
5.C.2. Optimization model
An optimization model for dose planning of IMRT with MTD constraints was first formulated by Romeijn et al.76. For HDR BT, optimization models with MTD constraints have been proposed.58,77,78
The following model, based on the model formulation,46 shows an example of a MTD model with constraints that impose upper bounds on UCV aRsβs, s∈S, and lower bounds
on LCV aRsαs, s∈S. The objective function is based on the
LPM (presented in Section 5.1), but with piecewise linear penalties. The notation is given in TablesIIandVI.
min ∑s∈S∑q∈Ql sw l sq∑i∈Psx l iqþ ∑s∈S∑q∈Qu sw u sq∑i∈Psx u iq (11a) subject to ∑j∈Jdijtj ≥ Lsq x l iq i∈Ps, q∈Qls, s∈S (11b) ∑j∈Jdijtj ≤ Usqþ xuiq i∈Ps, q∈Qus, s∈S (11c) ζs þ 1 βsjP Tj∑i∈PTξ s i ≤ Usβs s∈S (11d) ζs 1 αsjP Tj∑i∈PT ξs i ≥ L s αs s∈S (11e) ξs i ≥ ∑j∈Jdijtj ζ s i∈Ps, s∈S (11f) ξs i ≥ ζ s ∑ j∈Jdijtj i∈Ps, s∈S (11g) ξs i ≥ 0 i∈Ps, s∈S (11h) ξs i ≥ 0 i∈Ps, s∈S (11i) ∑j∈Jdijtj ≤ Ms i∈Ps, s∈S (11j) xl iq ≥ 0 i∈P,q∈Q l s, s∈S (11k) xu iq ≥ 0 i∈P,q∈Q u s, s∈S (11l) tj ≥ 0 j∈J (11m)
The objective function (11a) and constraints (11b), (11c), (11k) and (11l) are similar to those presented and explained for the LPM in Section 5.1. The difference is that the penal-ties are here piecewise linear with a number of segments; see Fig.1 for an illustration of a piecewise linear penalty func-tion. Constraints (11j) constitute upper dose bounds for all dose point. The linear formulation of UCV aRs
βs is given in
constraints (11d), (11f) and (11h), as presented by Rockafellar and Uryasev79. The auxiliary variables ζs take the values of the dosimetric indices Ds
βs. The variables ξ
s
i are only used to
define the maximum of the two right-hand-side values in con-straints (11f) and (11h). The linear formulation of LCV aRsαs is
similarly defined with constraints (11e), (11g) and (11i). The feasible region defined by constraints (11d)–(11i) is a convex set, as shown in Ref. [74]. Furthermore, this model is a restriction of the LPM, as the feasible set is made smaller by adding constraints (11d)–(11i).
FIG. 3. Example of a differential DVH curve. The mean dose in the grey area is the LCVaR value and the mean dose in the black area is the UCVaR value.
TABLEVI. Notation used in the mean-tail-dose model.
Indices
q Index for a segment in the piecewise linear penalty function Sets
Ql
s Set of penalty segments for dose being too low in structure s∈S
Qu
s Set of penalty segments for dose being too high in structure s∈S
Parameters
αs A portion of the structure s∈S
βs
A portion of the structure s∈S Ls
αs Lower bound on the value of LCV aRsαs, s∈S Us
βs Upper bound on the value of UCV aRsβs, s∈S wl
sq Non-negative penalty for segment q∈Qlsfor the dose being too
low in structure s∈S wu
sq Non-negative penalty for segment q∈Q u
sfor the dose being too
high in structure s∈S Ls
q Lower bound of segment q∈K l
sfor structure s∈S
Us
q Upper bound of segment q∈K u sfor structure s∈S Variables ζs Auxiliary variable ζs Auxiliary variable ξs
i Auxiliary variable which is used to capture the maximum of two
expressions ξs
i Auxiliary variable which is used to capture the maximum of two
expressions xl
iq Penalty variable for segment q∈Q l
sfor dose being too low at the
dose point i∈Ps, s∈S
xu
iq Penalty variable for segment q∈Qusfor dose being too high at the
The latter MTD constraints for LCVaR and UCVaR are combined with an objective function to maximize an approxi-mation of the DHI,77and the resulting model is formulated and solved as an LP. A DVM is extended with a MTD com-ponent to maximize the LCVaR of the PTV,58 giving a weighted sum objective. The purpose of this extension is to avoid tumor cold spots. (The extended model was also solved faster than the DVM.) A variant of LCVaR and UCVaR con-straints where a number of dose points, which are expected to receive a very high or very low dose, are removed has also been suggested.78
5.D. Quadratic penalty models
In HDR BT it has traditionally been more common with linear models, such as the LPM, while in IMRT it has been more common with quadratic penalty models (QPM). There are however such models also for HDR BT dose plan-ning,80,81 which were developed to the HIPO82 (hybrid inverse treatment planning optimization) algorithm. During the last decade, several QPMs have been proposed.52,83–88 Similarly to the LPM, except for the prescription doses, the objective of the QPM does not directly correspond to the clinical evaluation criteria. In Section 5.4.1 we introduce for-mulations for the QPM.
5.D.1. Optimization models
We first present a generic formulation of a QPM, which is formulated in for example Ref. [83]. For each dose point, there is a penalty if the dose deviates from a specified dose. The notation used in the following models is given in TableII. min ∑s∈Sws∑i∈PsðDosei L sÞ2 (12a) subject to ∑j∈Jdijtj¼ Dosei i∈Ps, s∈S (12b) tj≥ 0 j∈J (12c)
The purpose of the model is to penalize doses being higher or lower than specified values. The objective function can instead be formulated as a piecewise quadratic penalty function, giving penalties if the doses are outside specified intervals, see for example Ref. [86]. Such a function is illus-trated in Fig.1.
The alternative formulation (13) is constructed by remov-ing the non-negativity constraints (12c) and replacing the dwell time variables tjwith t2j.
min ∑s∈Sws∑i∈Psð∑j∈Jdijx 2 j LsÞ 2 (13) Here, x2
j is used to denote the dwell time tj. This gives an
unconstrained model, which was solved81,88 with the quasi-Newton method L-BFGS (limited memory Broyden– Fletcher–Goldfarb–Shanno).89For an introduction to nonlin-ear programming in general, see Ref. [90].
In the studies,84,86,87 an extra component is added to the objective function penalizing heterogeneous dwell times. Both the one-norm and the two-norm of the differ-ence between adjacent dwell positions as regularization are evaluated.84 The former formulation is used in Ref. [86], and the latter formulation in Ref. [87]. The purpose of these regularizations is the same as that of the DTMRs (presented in Section 5.1.3), to achieve a more homogeneous distribution of the dwell times. A regular-ization term is added,52 referred to as Tikhonov regular-ization, which is the two-norm of the dwell times. Here, the focus is rather on reducing long dwell times than achieving more homogeneous dwell times.
The model (12a)–(12c) is convex while the model (13) is non-convex; the details are given in Appendix B.
5.E. Radiobiological models
Instead of taking the detour of evaluating a dose plan based on physical dose and its coupling to the dosimetric indices, the anticipated radiobiological effects of the dose distribution can be explicitly incorporated in an optimiza-tion model. Examples of such attempts are optimizaoptimiza-tion models based on EUD and on TCP. The use of these radiobiological models has drawbacks however, a major one being the difficulties of estimating generally reliable values of parameters. Both the use of EUD models and TCP models are discussed,91 although with focus on EBRT.
5.E.1. EUD models
The radiobiological index gEUD is used in two HDR BT studies. A nonlinear optimization model for prostate cancer treatment based on gEUD is proposed,92with the aim of spar-ing OAR. A more homogeneous dose distribution with respect to COIN is also obtained. The index gEUD is also used in an optimization model for cervical cancer.93In both these HDR BT studies, the objective function used is the gEUD formulation21proposed in the context of IMRT. Equa-tion (14) shows this formulation using the notation given in TablesIIandVII.
max 1 1þ EUDT0 EUDT wT Y s∈SO 1 1þ EUDs EUDs 0 ws (14)
The aim with the tumor component is that the EUD value should be above the specified value, while for the OAR the EUD values should be below the specified values.
Furthermore, it is observed that overdosage of the tumor has only a small impact on the gEUD, while it is very sensi-tive to underdosing of the tumor.93
There has been a larger interest in optimization models based on gEUD in IMRT, possibly because nonlinear models are more common in IMRT, while linear models such as the LPM and the DVM are more common in HDR BT.
5.E.2. TCP models
The radiobiological index TCP is included in a bi-objec-tive model for cervical cancer.94The two objective function components are: (a) TCP, based on the formulation in Zaider and Minerbo95, and (b) a weighted sum of binary indicator variables (for dose points), each of which equals one if the dose is sufficiently high or sufficiently low, depending on the type of structure. This corresponds to a sum of dosimetric indices. This formulation is both nonlinear and non-convex. The primary goals in the solution process in Ref. [94] are the DIs and the model is solved as a MIP, using branch-and-cut, without the TCP component. Whenever an integer solution is found, a local search heuristic is however initiated to improve its TCP value.
In the study35on prostate cancer they consider the radiobi-ological objective LTCP (logarithmic TCP),96which with the notation given in TablesIIandVII, is formulated as
LTCP¼ 1 jPsj∑ j Psj i¼1eαðDoseiL sÞ , (15)
where s is the urethra. LTCP is included as one objective using α ¼ 0:5 in their multiobjective approach, in order to minimize the dose to the urethra. It is used as an approxima-tion of Durethra
1% .
Moreover, in an LDR BT study,97 estimates of TCP and NTCP are included. Because the OAR have different radiobi-ological sensitivity, small changes in the DVH curve for the rectum had a big impact on the NTCP, while a large change in the DVH curve for the urethra gave only a small change in the NTCP.
Regarding underdosing the tumor, it was shown by Tom´e and Fowler98that cold volumes reduce the TCP value signifi-cantly even when these constitute only 1% of the tumor. Related to this is an IMRT study on a version of EUD,99 denoted tail EUD, focused on increasing the dose to the cold-est part of the tumor. Hence, the purpose of the tail EUD is the same as that of LCVaR. An LCVaR component is added to the DVM,58also with the purpose to increase the dose to the coldest part of the tumor.
5.F. Multiobjective models
Because the evaluation of treatment plans takes multiple criteria into account, a frequent modeling approach to dose planning is multiobjective optimization models. In HDR BT, there are several studies on such models.3,33,80,81,88,100–106For
all of them, the final result is a set of dose plans that satisfies the Pareto optimality conditions, and there is no decision on which the best plan is.
With the goal to minimize a number of objective functions fiðxÞ, i ¼ 1,...,n, a dose plan x is said to dominate a dose plan y if the following two criteria are satisfied.
1. fiðxÞ ≤ fiðyÞ, i ¼ 1,...,n,
2. fkðxÞ< fkðyÞ holds for some k ∈f1,...,ng.
The Pareto (optimality) front is the set of all non-domi-nated dose plans, or in other words, solutions on a Pareto front cannot be improved in one objective without sacrificing another.
Figure 4 shows an example of a solution space, where each solution is plotted according to the values of two objec-tives, f1 and f2, which both should be minimized. The solu-tions that are plotted with filled circles belong to the Pareto front. For these points it is not possible to improve the value of one objective without making the other worse.
An advantage of presenting a set of Pareto optimal solu-tions instead of a single solution is that the trade-offs can be made clear for the planner, for example as shown in Fig.4. For some Pareto solutions, a large gain in one objective value is possible at the of cost of only a small deterioration in the other, while for some Pareto solutions this is not the case. Such an illustration is also useful for choosing a treatment plan, as it is important to know what objective values are at all attainable.
5.F.1. Single-objective approaches
A composite objective function which is a weighted sum of multiple quadratic objectives is considered.80,101 They solve a large number of nonlinear optimization problems with different weights, each of them giving a candidate to the Par-eto front. Also in Ref. [81] is the Pareto front obtained by varying the objective weights.
Moreover, in the studies33,88,105, single-objective models are considered and solved repeatedly to generate a good approximation of the Pareto front. The model33,105is based
TABLEVII. Notation used in the models with radiobiological indices.
Parameters EUDs
0 Is a structure specific parameter for s∈S, estimated from
dose-response data
α Parameter related to cell survival Variables
EUDs The value of EUD for structure s∈S, as calculated by
expression (4)
FIG. 4. Example of a Pareto front. All points (circles) correspond to a
solu-tion plotted according to two objective values, f1and f2, which both should
on linear penalties, and is solved with IPSA. For updating the objective weights, they calculated upper and lower estimates of the Pareto front, and used this information to get the next solution. The model88is based on quadratic penalties.
In one study on cervical cancer36, multiple objectives are considered in sequence, giving a number of single-objective problems with an increasing number of constraints. The order of the objectives is tuned to generate acceptable trade-offs.
5.F.2. Multiobjective approaches
Simulated annealing, which is popular in HDR BT dose planning, is a metaheuristic which considers a single solu-tion at each iterasolu-tion. In contrast to simulated annealing and the previously described approaches for composite objective functions, genetic algorithms belong to the class of population-based metaheuristics, which means that there is a set of solutions that is updated at each iteration.38 Genetic algorithms are the most commonly studied meta-heuristic methods in HDR BT that allow for dealing with the objectives in a truly multiobjective manner,3,100–104 since with a fitness criterion that is based on domination there is no need to combine all objectives into a single value. Furthermore, in the case of non-convex objectives, a weighted-sum method may not be able to find all solutions on the Pareto front, which is a disadvantage compared to dominance based methods, see Ref. [107], p. 73. With convex objectives, a weighted-sum formulation can be seen as equivalent to a formulation in which the objectives are constrained.108
With multiple objectives, it is not easy to visualize the set of possible Pareto solutions and the trade-offs that have to be made. The trade-offs are mainly illustrated by plotting the Pareto front with respect to two objectives at a time.80,101
Another approach to visualize the trade-offs is considered in the studies34,102,103,106(which are based on the evolutionary algorithm presented by Luong, Poutr´e, and Bosman109). Here, the multiobjective problem is reduced to a bi-objective prob-lem by partitioning the objectives into two categories, either related to PTV coverage or to OAR sparing, and for each cate-gory the objective value is defined as the worst objective value among those in the category, in comparison to a target value for each objective. Dosimetric indices for high doses to the PTV (VPTV
200 and V PTV
150) are considered as hard constraints.
The objective related to coverage, least coverage index (LCI), is formulated as
LCI¼ minfVprostate100 95%,Vvesicles80 95%g, (16) where 95% is the target value in the clinical treatment proto-col. The objective related to OAR sparing, least sparing index (LSI), is formulated as
LSI¼ minfδðDbladder 1cc Þ,δðD bladder 2cc Þ,δðD rectum 1cc Þ,δðD rectum 2cc Þ,δðD urethra 0:1cc Þg, (17) where Ds
yis the lowest dose that is received by the y cc of the
structure s that receives the highest dose; see Section 3.2. The functionδ is defined as
δ ¼ Ds,max y D
s
y, (18)
where Ds,maxy is a threshold value from the clinical treatment protocol for structure s and for volume y. The value of LCI is a portion and the value of LSI is in Gy. (The formulations (16)–(18) are taken from Ref. [103].) For both LCI and LSI a positive value indicates that the clinical requirements are sat-isfied, and thus it is easy to verify for which solutions, on the Pareto front, the requirements are satisfied.
The bi-objective formulation described above allows for not only an easier visualization but also reduces the possibili-ties for the planner to select a dose plan, because a dose plan that is not on the Pareto front might be the preferred one. An extension to this approach is presented3,104, but in addition they also consider the number of catheters as an objective to minimize (see Section 5.10).
An evolutionary algorithm is used to automate the parame-ter tuning for IPSA and HIPO.110They conclude that more sat-isfactory plans are found with HIPO than with IPSA. Furthermore, class solutions for different prioritization between target and OAR are also constructed. The results from IPSA and HIPO are also compared with the ones from the bi-objective formulation34and a slight advantage to the latter one is seen.
5.F.3. Interactive approaches
A multiobjective approach different from the ones pre-sented above is proposed.111Instead of providing the planner with a set of dose plans to choose among, the planner is active during an iterative and interactive optimization pro-cess. In each iteration the planner has five options for each objective:
1. to improve it,
2. to improve it up to a specified level,
3. to allow for a deterioration to a specified level, 4. to keep it fixed, or,
5. to not bound it at all.
The result from an iteration can either be a single dose plan or a set of dose plans based on convex combinations of pairs of previous dose plans. This approach is based on the algorithm NIMBUS.112
The approach PNaV78 (pareto navigation and visualiza-tion) is based on a CVaR constrained model and a library of plans is created by varying model parameters. Furthermore, to generate new plans based on the plans in the library, an additional optimization problem is solved according to the practitioner’s wishes. A visualization tool is also provided to illustrate the necessary trade-offs. The method was tested on prostate, cervical and breast cancer.
5.G. Modeling of spatial properties
The optimization models and evaluation criteria that have been presented so far evaluate the dose distribution on an
aggregate level only, and how the doses are spatially dis-tributed are not in any way considered. For example, dosimet-ric indices only count the number of dose points in which the dose is high enough or low enough, and no consideration is given to how these dose points are spatially distributed in the volume.
An optimization model for IMRT that takes also the spa-tial dose distribution into account is proposed,113 while for HDR BT such a model is proposed.69,70 Both approaches take a tentative dose plan and tries to improve upon spatial properties.
The aim in Ref. [113] is to improve a tentative dose plan by removing“critical spots,” which are either regions in the PTV where the radiation doses are too low (“cold spots”) or regions in OAR where the radiation doses are too high (“hot spots”). In their model, there are constraints to adjust the doses within the critical spots while the objective strives to maintain key properties of the given plan, including original objective func-tion values, feasibility and the dose at each dose point.
The aim in Ref. [70] is to reduce PTV hot spots and cold spots in a tentative dose plan. Constraints are maintaining the dosimetric indices from the tentative plan and the purpose of the objective function is to improve the plan’s spatial proper-ties.
5.G.1. Optimization models
The complete models70,113 are as follows, with the nota-tion introduced in TablesIIandVIII.
mint ∑b∈IfðfbðtÞ f 0 bÞ 2 þ∑c∈Icπcmax 0, gð cðtÞÞ þ∑i∈Pri Dosei dose0i 2 (19a) subject to Li≤ Dosei≤ Ui i∈Is (19b) ∑j∈Jdijtj¼ Dosei i∈P (19c) t∈X (19d)
Formulation (19a)–(19d) is the model given in Ref. [113]. As this study is conducted for IMRT, the variables t can be seen as treatment plan parameters. The purpose of the objec-tive (19a) is to maintain values from the tentative plan of the objectives, constraints, and doses at dose points. The set of constraints (19b) ensures that the doses are adjusted within the critical spots to respect the dose bounds Liand Ui.
Formulation (20a)–(20c) below is the model in Ref. [70]. mint w∑i,l∈PT:i≠l
fðDoseiÞf ðDoselÞ
hi,l þ ð1 wÞ∑i,l∈PT:i≠l
gðDoseiÞgðDoselÞ
hi,l
(20a) ∑j∈Jdijtj¼ Dosei i∈P (20b)
t∈X (20c)
Here, X is based on constraints on dosimetric indices for both PTV and for OAR; see Section 5.2 for the mathematical formulation of these constraints. Functions f and g are based on the sigmoid expression (8), giving penalties for doses being too high and too low, respectively. Parameter values for the DIs are based on the tentative dose plan or clinical treat-ment guidelines, possibly allowing a small deviation.
5.H. Intensity-modulated brachytherapy
A recently proposed treatment technique is the IMBT, in which rotating shields are used to create a dose distribution that is not radially symmetric around the source.114The rotat-ing shields introduce an extra degree of freedom in the plan-ning. IMBT techniques are also known as dynamically modulated brachytherapy or direction-modulated brachyther-apy (DMBT) or rotating shield brachytherbrachyther-apy (RSBT).
These techniques have the potential to significantly spare OAR,115,116with doses to urethra that are 29–44% lower with IMBT than with conventional HDR BT.115The most common radiation source in conventional HDR BT is192Ir. Due to the shield thickness required to shield 192Ir, alternatives have been suggested for IMBT, such as an electronic source117or other radiation isotopes.118When for example153Gd is used, TABLEVIII. Notation used in the spatial optimization models.
Indices
b Index for an objective function c Index for a constraint i Index for a dose point l Index for a dose point Sets
If Set of objective functions
Ic Set of constraints
Is Set of dose points in the critical spot
X Set of feasible dwell times (or treatment plan parameters) Functions
fb Objective function b∈If, in original model formulation
gc Constraint c∈Ic, in original model formulation
f Non-negative penalty function, giving a penalty for doses being too high
g Non-negative penalty function, giving a penalty for doses being too low
Parameters dose0
i Received dose at dose point i∈P, in the tentative dose plan
f0b Objective value of function fb, b∈If, in the tentative dose plan
hil Distance measure between dose points i and l
πc Non-negative penalty parameter for infeasibility in constraint gj,
c∈Ic
ri Non-negative penalty parameter for dose deviation at dose point
i∈P
Li Lower dose bound for dose point i∈Is
Ui Upper dose bound for dose point i∈Is
w Non-negative penalty parameter for weighting objectives, w∈½0,1
