A mathematical optimization model for spatial adjustments of dose distributions in high dose-rate brachytherapy

High dose-rate brachytherapy is a modality of radiation therapy used for cancer treatment, in which the radiation source is placed within the body. The treatment goal is to give a high enough dose to the tumour while sparing nearby healthy tissue and organs (organs-at-risk). The most common criteria for evaluating dose distributions are dosimetric indices. For the tumour, such an index is the portion of the volume that receives at least a specified dose level (e.g. the prescription dose), while for organs-at-risk it is instead the portion of the volume that receives at most a specified dose level. Dosimetric indices are aggregate criteria and do not consider spatial properties of the dose distribution. Further, there are neither any established evaluation criteria for characterizing spatial properties, nor have such properties been studied in the context of mathematical optimization of brachytherapy. Spatial properties are however of clinical relevance and therefore dose plans are sometimes adjusted manually to improve them. We propose an optimization model for reducing the prevalence of contiguous volumes with a too high dose (hot spots) or a too low dose (cold spots) in a tentative dose plan. This model is independent of the process of constructing the tentative plan. We conduct computational experiments with tentative plans obtained both from optimization models and from clinical practice. The objective function considers pairs of dose points and each pair is given a distance-based penalty if the dose is either too high or too low at both dose points. Constraints are included to retain dosimetric indices at acceptable levels. Our model is designed to automate the manual adjustment step in the planning process. In the automatic adjustment step large-scale optimization models are solved. We show reductions of the volumes of the largest hot and cold spots, and the computing times are feasible in clinical practice.

Given a dose distribution that satisfies clinical evaluation criteria, it is common (depending on clinical routines or preferences of the planner) that the plan is inspected visually to find contiguous volumes with a too high dose, called hot spots. If necessary, the dose plan is then adjusted manually, to remove or reduce the hot spots. The manual adjustment process is described in, for example, Maree et al (2019). This task is dependent upon the skills of the planning staff and can be time consuming. Due to the difficulties of this process and the lack of established evaluation criteria that take spatial properties into account, it is relevant to consider developing mathematical optimization models for improving spatial properties of a dose distribution.

Aim
With the aim of automating the adjustment step that is currently performed manually, we propose an optimization model that takes a tentative plan, which is acceptable with respect to clinical evaluation criteria, to improve it with respect to spatial properties. The goal is to reduce the prevalence of both hot spots and cold spots (contiguous volumes with a too low dose) in the PTV and thus make the dose distribution more spatially homogeneous, while retaining the evaluation criteria from the clinical treatment guidelines at acceptable levels.

Outline
The paper is organized as follows. The remainder of this section contains an overview of aggregate and spatial evaluation criteria from the literature. Despite the lack of spatial evaluation criteria in the clinical treatment guidelines, a few have been proposed for BT. We also present such criteria proposed for external beam radiation therapy (EBRT), although our focus is on BT. In section 2 the mathematical optimization model is presented and exemplified, and in section 3 we present the settings for the computational study to evaluate our model. In section 4 we show computational results, to give a broad view of what changes can be seen from the adjustment step. Section 5 contains a discussion and suggestions for future research, and finally in section 6 we conclude the paper.

Evaluation criteria and modelling
A common way to visualize and evaluate the three-dimensional dose distribution is the cumulative dose-volume histogram (DVH). The DVH shows, for each dose level, the portion of a structure that receives at least (or at most) this dose. Each point on the DVH-curve corresponds to a dosimetric index (DI); for the PTV it is the portion of the volume that receives at least a specified radiation dose (for example the prescription dose), while for an OAR it is here instead the portion of the volume that receives at most a specified dose. A DI is denoted V s x , where x is the radiation dose, expressed as a percentage of the prescription dose, and s is the structure (PTV or an OAR). In our presentation, if no structure is specified, V x refers to the PTV. Another form of DI is denoted D s y , where y is a portion of the volume of the structure s; this is the lowest dose that is received by the y% of the volume which receives the highest dose. An alternative to the most commonly used cumulative DVH-curve is the differential DVH-curve, which for each dose level shows the portion of the structure which receives exactly this dose.
The DIs are fundamental in the clinical evaluation of dose plans and they are included in the treatment guidelines for HDR BT, see for example Hoskin et al (2013) for the case of prostate cancer. Dosimetric indices have therefore been incorporated in mathematical optimization models for dose planning. A model of this type for HDR BT is referred to as a dose-volume model (DVM) and it was first formulated as a mixedinteger program in Beliën et al (2009). Dose-volume models have further been studied in Siauw et al (2011), Gorissen et al (2013), Guthier et al (2017), Morén et al (2018a), and Morén et al (2019). The DVMs are nonconvex and computationally hard to solve. Therefore, common approaches have been heuristics, see Beliën et al (2009), Siauw et al (2011), and Deist and Gorissen (2016), and approximations, see Guthier et al (2017) and Morén et al (2018b).
From a mathematical point of view dosimetric indices are similar to a quantity in financial risk management, known as value-at-risk (VaR). A convex risk measure which is related to VaR is the conditional value-at-risk (CVaR). While VaR can be defined as the highest dose that is received by a specified portion (α%) of the dose points which receives the lowest dose, CVaR is instead the mean dose that is received by a specified portion (α%) of the dose points that receives the lowest dose. Hence, in the context of radiation therapy CVaR has been referred to as mean-tail-dose. Both VaR and CVaR can be defined similarly for a specified portion that receives the highest dose. Figure 1 illustrates how VaR and CVaR are derived from the DVH-curve. Rockafellar and Uryasev (2002) showed that CVaR can be formulated as a linear optimization model. In radiation therapy, mean-tail-dose was first used in a model for EBRT in Romeijn et al (2003); it has been used in HDR BT models in Holm et al (2013) and Morén et al (2019).
Several DI-based criteria have been proposed to measure the dose homogeneity, since a homogeneous dose distribution is generally considered as favourable, see Wu et al (1988). Examples of such criteria are the dose homogeneity index (Wu et al 1988), the conformation number (Riet et al 1997), and the conformal index (Baltas et al 1998).
One way of making dose distributions more homogeneous is the dwell time modulation restriction (DTMR) which was introduced in Baltas et al (2009) and further studied in Gorissen et al (2013) and Balvert et al (2015). A DTMR is added to ensure that dwell times at neighbouring positions do not differ too much. The effect is that long dwell times are typically shortened and more positions become active (that is, have positive times). Because of the added restriction the optimal value of the primary objective value in such a model is however worsened to some extent (Balvert et al 2015).
Dosimetric indices, mean-tail-dose, and dose homogeneity measures are evaluation criteria that describe aggregate properties of a dose distribution and do not take spatial properties into account. We refer to such criteria as aggregate evaluation criteria.

Spatial properties
Despite the lack of spatial criteria in the treatment guidelines, there are good reasons to consider spatial properties of dose distributions. In Kirisits et al (2014) it is mentioned that 'the spatial dose distribution could remain as a key element for outcome'. Furthermore, the role of dosimetric indices in the treatment guidelines is also discussed in the scientific community; in a debate (Njeh et al 2015) regarding whether the DVH-based criteria are appropriate as evaluation criteria, Dr Christopher Njeh argues that DVH-based criteria should be replaced with biologically based criteria, because of the lack of spatial information in the former. Indications that spatial properties are important are also presented in Bijl et al (2003); in a study on rats, they found that one large contiguous hot spot in the cervical spinal cord gave more complications than two separate hot spots.
Although there are no criteria in the clinical treatment guidelines that consider spatial properties, there are a few evaluation criteria in the literature that consider spatial properties: for BT there are the contiguous V 150 (Thomas et al 2007, Poulin et al 2013 and the contiguous V 200 (Thomas et al 2007), both related to the corresponding dosimetric index, and for EBRT there is the large cold spot metric (LCS D ), see Said et al (2017). The LCS D takes values between zero and one, and is defined as where S D denotes the number of dose points that receives a dose that is lower than D Gy and M cold (D, z) is the number of cold spots (contiguous dose points receiving less than D Gy) of size equal to z (number of dose points). A high LCS D value indicates that a large portion of the cold dose points are contiguous. (The measure can be defined analogously for hot spots.) It should be noted that dose distributions in BT are considerably more spatially inhomogeneous than their EBRT counterparts. This is inherent to the approach of internal irradiation. In BT, delivering doses well above the prescription dose to a significant part of the PTV is common. Therefore, values from spatial evaluation criteria are hard to compare between BT and EBRT.
The concept of dose-volume histograms has also been modified and used to capture spatial properties. Examples of this are the z-dependent DVH (Cheng and Das 1999), the spatial DVH (Zhao et al 2010) and the subvolume-DVH (Said et al 2017). First, the z-dependent DVH is defined for a specified dose level, and for each two-dimensional cross-section of the medical image the portion that receives at least the specified dose level is plotted. The curves corresponding to several dose levels are plotted in one graph. Secondly, the spatial DVH shows the differential DVH for the PTV, which is divided into periphery, middle and centre dose points. For each of these PTV volumes, the differential DVH is colour-coded and plotted in the same graph. Finally, the subvolume-DVH shows the differential DVH-curve for the three largest contiguous cold spots in the same graph.

Related work
For dose planning of EBRT, an optimization model that takes spatial properties into account is proposed in Süss et al (2013). Their aim is to improve the spatial properties of a tentative dose plan, constructed with a given optimization model. A critical region of dose points is identified and constraints are added to enforce the doses at these points to be in specified intervals (with possibly different intervals for different dose points). The objective function of the spatial optimization model is a weighted sum with components that are based on the given optimization model. Its aim is to deviate as little as possible from the tentative dose plan, both with respect to the given objectives and constraints, and the received dose at each dose point. Another EBRT optimization model considering spatial properties is proposed in Chao et al (2003). Its aim is to remove cold spots in the centre of the tumour while cold spots on the periphery of the tumour are more acceptable. The rationale for this is to ensure a sufficiently high dose to volumes that are likely to contain cancer cells, rather than preventing the cold dose points from forming a contiguous volume.
From a mathematical perspective, our optimization model is related to the dispersion problem, in which spatial properties are central. The common scenario is here to select a number of elements (with known distances between them) and evaluate the selection according to a distance-based measure. Examples of such problems are the maximum dispersion problem (Fernández et al 2013), which is to select a given number of elements such that the minimum distance is maximized, the maximum mean dispersion problem (Prokopyev et al 2009), which is to select elements such that the mean distance is maximized, and a dispersion problem in which the total distance between the chosen elements is maximized (Glover et al 1995), with a knapsack constraint to limit the number of chosen elements.

Proposed model
This work is a continuation and extension of the results presented in the conference proceedings paper by Morén et al (2018b). The optimization model proposed here takes both hot spots and cold spots into account, contrary to the previous work in which only hot spots were included. The model is also formulated in a more general form. Three optimization models are used to construct tentative dose plans which are used as input to the adjustment step, compared to one model in Morén et al (2018b). These three models are not a part of the proposed model, but are only used for evaluation. We also consider tentative plans that have been constructed manually and used in clinical practice. Furthermore, the results and analyses here are more thorough than in the previous work, both in terms of spatial and aggregate properties of the dose distributions.
Similarly to the study by Süss et al (2013), we propose an adjustment step to improve a tentative dose plan. However, our model is to a large extent different from theirs. They have to explicitly define the critical regions to be improved, including dose bounds to be enforced, and the priority is to improve on spatial properties even if this results in a deterioration of clinical evaluation criteria. In our model, the critical regions do not need to be defined explicitly. Instead, the objective function is constructed to identify and reduce any hot and cold spots, which makes our approach more flexible. The constraints in our model ensure that aggregate criteria are retained at acceptable levels, and thus, our priority is the clinical evaluation criteria with the secondary aim to improve upon spatial properties.

Mathematical model
We propose a mathematical optimization model for spatial dose adjustments. The starting point is that we have a tentative dose plan with acceptable values of aggregate evaluation criteria. The adjustment step is independent of the dose planning method used to construct the tentative plan; the tentative plan could be generated automatically using optimization or be a manually constructed dose plan.
The purpose of the objective function is to a find an adjusted dose plan that is more spatially homogeneous, by reducing the prevalence of contiguous hot spots and cold spots. Our approach considers all pairs of dose points, and is based on two assumptions of how the spatial properties of a dose distribution can be improved: (i) The number of dose points receiving a dose that is too high, or too low, should be reduced.
(ii) For a pair of dose points which both receive doses that are too high, or too low, it is better the farther apart the two dose points are.
The first type of improvement is the preferred one, but it can be impossible to reduce the number of such dose points to zero, or even to a large extent, while retaining aggregate evaluation criteria, and thus it is necessary to also try to redistribute them. The objective function of our model is a sum of penalties where each term considers a pair of dose points in the PTV. Both the dose to each of the two dose points and the distance between them are taken into account. Additionally, there are constraints to ensure that aggregate criteria remain at acceptable levels.

General optimization model
The model in its most general form is presented as: Here, t is the vector of dwell times, X is the feasible set, defined by constraints on aggregate evaluation criteria, i and k are indices for dose points, the set T contains all dose points in the PTV, and D i is the dose received at dose point i ∈ T . The functions f and g are designed to penalize doses which are too high and too low, respectively, and the function h(i, k) is a measure of the distance between dose points i and k. The objective is a weighted sum of two components, with the weight w 1 penalizing the hot spot component and the weight w 2 penalizing the cold spot component. The small example in figure 2 illustrates the effect of the model. Here, there are four dwell positions, illustrated by the white squares placed symmetrically around the centre of two dose distributions. Each dose point is coloured according to the received dose, where a dark colour indicates a low dose and a light colour indicates a high dose. In the dose distribution to the left there is a large contiguous hot spot, which is unfavourable. Such a hot spot is avoided in the dose distribution to the right where there instead is two smaller hot spots, which is preferred, provided that the dose plan is otherwise acceptable. The sum of the dwell times is the same in both figures.

Dose-volume constraints
In our implementation of the optimization model, the constraints are based on dosimetric indices, both for the PTV and for each OAR, see Siauw et al (2011) and Deist and Gorissen (2016). This choice of constraints is based on the fact that they are fundamental in the clinical treatment guidelines. The indices, sets, parameters and variables used in the model are presented in table 1 and the constraints defining the set X are given by (3a)-(3g). The actual decision variables in the model are the dwell times t j , j ∈ J. The dose at the dose point i is calculated as j∈J d ij t j . Constraints (3a) and (3b) enforce a bound on a dosimetric index for the PTV, where (3a) ensures that the indicator variables y i , i ∈ T, take the correct values and (3b) imposes a bound on the value of the dosimetric index. Similarly, constraints (3c) and (3d) enforce a bound on a dosimetric index for an OAR, where (3c) ensures that the indicator variables v s i , i ∈ OAR s , s ∈ S, take the correct values and (3d) impose bounds on the dosimetric indices for the OAR s ∈ S. The objective is to minimize (2).
The values of the parameters ω and τ s , s ∈ S, are based on either the clinical treatment guidelines or the values in the tentative plan, possibly allowing small deviations. The linear penalty model (LPM), see Lessard and Pouliot (2001) and Alterovitz et al (2006), is the most commonly used optimization model in clinical practice, and an alternative way to impose aggregate criteria would be to convert the LPM objective value into one or more constraints. However, we would then not have any direct control on the primary evaluation criteria, the DIs, and thus the dose distribution might change in undesired ways.

Experimental design
The functions f , g and h in (2) can be chosen in various ways. To evaluate our model we have tried two versions of the function f , which gives a penalty if the dose is too high, either where H is the Heaviside step function and b is chosen as 200% of the prescription dose. Similarly, for g(D i ) we implemented the quadratic penalty function (4) and the Heaviside step function (5), but instead formulated to penalize a dose that is too low. In the results to be presented, the same type of function is used for both f and g.
As the distance function h(i, j) we have tried both the standard Euclidean distance and the squared Euclidean distance. The results are similar and we therefore only present them for the squared Euclidean distance.

Patient data set
We have tested our models retrospectively on 13 cases obtained from anonymized patients that have been treated for prostate cancer. The patient data comprised the catheter placements and the anatomical structures (PTV, urethra and rectum) that were outlined according to the practice of the clinic that provided the patient data. Clinical treatment guidelines for prostate cancer (Hoskin et al 2013) also mention the bladder as a structure of interest, but it was not outlined in our data because the medical images are from ultrasound. According to common practice we added artificial, healthy tissue surrounding the PTV. The number of catheters was in the range 14-20 and the number of dwell positions in the range 190-352. We used different sets of dose points for the optimization and for the evaluation of the dose plan; for the optimization the number of dose points was in the range 4369-7939 and they were distributed according to Lahanas et al (2000), and for the evaluation the number of dose points was in the range 51 974-134 509 and distributed uniformly with a volume of 1 mm 3 per dose point.

Tentative plans
As tentative dose plans we consider optimized plans for ten patients and clinically used plans for three patients.
For the former tentative plans, we have tried three optimization models. The first model is the dose-volume model (DVM), see Siauw et al (2011) and Deist and Gorissen (2016), with constraints (3a) and (3c)-(3g), and objective to maximize V 100 . The second is the linear penalty model (LPM), see Lessard and Pouliot (2001) and Alterovitz et al (2006), for which we tuned the parameters with respect to the set of patients. The third is a modification of the DVM (referred to as the MDVM) which has constraints (3a)-(3g), and objective to minimize the dose to the points in the PTV which receive a too low dose. The aim of the third model is to find plans with acceptable values of dosimetric indices but also with spatial properties that are poor, with significant cold spots. The purpose of using such tentative plans is only to test our model on cases where we know that spatial properties with respect to cold spots can be improved.
The clinical plans were constructed manually by an experienced planner. The planning includes time demanding adjustments to achieve spatial homogeneity, as is described in section 1. These three patients were considered to be normal cases with respect to their anatomy.
Parameter values, such as the prescription dose and the dose bounds for OAR are presented in table 2. The specific values in the table are not crucial for the analysis of the computational experiments.

Solution method
Initially we tried to solve the optimization model (2), with functions f and g given by (5) and constraints given by (3a)-(3g), as a mixed-integer program. For this purpose we used a state-of-the-art solver, Gurobi (Gurobi Optimization, Inc., Houston, USA) version 7.0.1, see Gurobi Optimization LLC (2018). This attempt turned out to be computationally intractable and no solution progress could be seen even after days of run time. Therefore an approximation of the mixed-integer program was needed.
The Heaviside step function has previously , Guthier et al 2017 been approximated with the smooth, non-linear sigmoid function, where β is a parameter to control the steepness; this approximation turned out to work well in practice also in our setting. We use values of β equal to 2.0 or 4.0, meaning that the function value increases from just above zero to almost one within intervals of lengths 3 Gy and 1.5 Gy, respectively. The approximation (6) is used in both the constraints for the dosimetric indices and in the objective function, yielding a continuous but non-convex model. The resulting model is however computationally tractable.
For all the results that we present, the objective function is based on the approximation (6) of the Heaviside step function (5), since the results using (6) were better than those obtained using (4). To analyze the trade-off between the hot and cold spot objective function components in (2)  Values used for the prescription dose and dose bounds. The values for the PTV and OAR to the left are used in the optimization models to construct tentative plans. The dose targets and the dose bounds are adopted from Siauw et al (2011) (except the parameter for artificial tissue which is here lower, and close to the value used in Deist and Gorissen (2016)). The values to the right are used for adjusting the clinically used dose plans and are taken from the clinical practice. The OAR (urethra and rectum) parameters for the dose target, the dose bound, and the portion are in There are no standard definitions of hot spots or cold spots. We here define a hot spot to be a contiguous volume which receives more than 200% of the prescription dose (dose points in such a volume are referred to as hot dose points), and a cold spot to be a contiguous volume which receives less than 90% of the prescription dose (cold dose points).

Subvolumes
To evaluate spatial properties of the dose distribution we introduce small, contiguous subvolumes of the PTV, where each subvolume consists of a number of dose points. The subvolumes are constructed to be centred between two adjacent dwell positions, and thus the number of subvolumes is of the order of the number of adjacent dwell positions. The purpose of this construction is to be able to detect contiguous hot spots of the type illustrated in figure 2(a). The light region is a typical example of a subvolume. Such a hot spot should, if possible, be avoided.

Computational results
The results are organized in three parts. First we focus on the adjustments of the tentative dose plans which were constructed using optimization models. We here divide the presentation into one part for spatial properties and one part for aggregate measures describing the dose distributions. In the third part we present results for the adjustments of the clinically used dose plans.
The adjustment step yields large-scale optimization problems with objective functions that have of the order of 10 7 terms (pairs of dose points). The mean computing time for the adjustment step is a few minutes, which is clinically feasible. We focus the presentation on the results for tentative dose plans from DVM and LPM, while the results for MDVM plans are presented only to illustrate reductions in cold spots. All in all the analysis in this section is based on comparisons of approximately 2000 dose plans.

Spatial properties
The aim of the first part of the results is to highlight changes in spatial properties of the dose distribution as a result of the adjustment step, and to show how the values of the weights (w 1 , w 2 ) affect the results. There are no established clinically used criteria for evaluation of spatial properties, and therefore dose plan changes are illustrated in several different ways. We first present visualisations of two-dimensional cross-sections of the PTV, and then comparisons of objective function values, volumes of contiguous hot and cold spots, and average doses to the subvolumes that receive the highest and lowest doses.
In the clinic it is common to identify hot spots by visualising dose plans in two-dimensional cross-sections. Figure 3 shows visual examples of hot spots (to the left in each figure) which are reduced by the adjustment step (to the right). Each dot corresponds to a dose point, for which a dark colour indicates a low dose and a light colour indicates a high dose. In the white areas there are no PTV dose points and they correspond to either dwell positions or the urethra. Correspondingly, in figure 4 cold spots (to the left in each figure), with dose points in blue (dark), are reduced by the adjustment step (to the right). (The dose point density is the same in all of the figures, but the actual sizes of the cross-sections vary.) To quantify the effects on the spatial properties we calculate the volumes of the largest contiguous hot and cold spots in the PTV. For comparison the total hot and cold volumes in the PTV are also calculated. These results are presented in table 4 for tentative plans from the DVM and in table 5 for tentative plans from the LPM. The hot and cold contiguous volumes are calculated as connected components in the uniform three-dimensional grid of dose points (Said et al 2017). The largest contiguous volumes decrease for both hot and cold spots, most noticeable when the weight is only on the objective component for hot and cold spots, respectively. The largest contiguous hot spot shrank from, on average, 2.77 cubic centimetres (cc) to 2.00 cc (28% decrease), and from 2.32 cc to 1.66 cc (28% decrease), compared to the DVM and the LPM tentative plans, respectively. The largest contiguous cold spot shrank from, on average, 1.20 cc to 1.14 cc (5% decrease), and from 1.12 cc to 1.11 cc (1% decrease), compared to the DVM and the LPM tentative plans, respectively. It is worth noting that there are larger decreases of the largest contiguous hot and cold spots than in total hot and cold volumes. Therefore, both of the two effects on spatial improvements, described in the beginning of section 2, can be observed in tables 4 and 5. For some patients there are only very small contiguous cold spots and thus there is no room for improvements Table 3. Weights used to compare the trade-off between objective function terms for hot spots (w 1 ) and cold spots (w 2 ). Each column corresponds to a set of weights, used in (2).  19pp) in the adjustment step. Examples of improved spatial homogeneity with respect to cold spots, for tentative plans from the MDVM, can however be seen in table 6. The volumes of the largest contiguous hot and cold spots are related to measures that were proposed and studied in Said et al (2017), see (1). We therefore compared values for the measure LCS D (and the equivalent formulation for hot spots), but the effect from the adjustment step was less clear. A reason for this is that the large cold spot metric should be used to compare two dose plans with roughly the same number of hot or cold dose points, but in our experiments the total hot and cold volumes were reduced.
The purpose of the objective function in (2) is to quantify spatial properties of dose distributions and we therefore also present the objective values of the adjusted dose distributions compared to those of the tentative plans. Figure 5 shows the trade-off induced from varying the weights (w 1 , w 2 ) according to table 3, with respect to the mean values of the hot and the cold spot component, respectively. Considering each component as a separate objective, the figure shows approximate Pareto fronts with respect to the two objectives. The Pareto fronts shown are not completely smooth because the solutions found are not necessarily global optima. The presented values are normalized such that both objective components of the tentative plan have value one. With equal weights both the hot and cold spot components are improved. It is also clear that there is a trade-off between the hot and cold components. For tentative plans obtained from the DVM and the LPM there are several weights which improve both the hot and cold spot components, while for tentative plans obtained from the MDVM both components are improved for all weights.
In figure 6 the mean dose to the subvolume that receives the highest dose is studied. For each patient and for each set of weights in table 3, the change in this mean dose is plotted. A negative value means that the adjustment step decreases the mean dose compared to that of the tentative plan. The improvement is larger for tentative plans obtained from the DVM than for plans from the LPM. Figure 7 instead shows the changes in the mean dose to the subvolumes which receive the lowest dose, which are small for tentative plans obtained from both the DVM and the LPM.
With tentative plans from the MDVM there are however significant improvements in mean dose to both the subvolumes which receive the highest and the lowest doses. Figure 7(c) shows the increase in the mean dose to the coldest subvolume compared to tentative plans from the MDVM, and for all weights and all patients there is an improvement. Further, the pattern of the decrease in the mean dose to the hottest subvolume is similar to the results presented in figure 7(a).

Aggregate measures
In this section, we show results for aggregate evaluation criteria, including dosimetric indices. Figure 8 illustrates the trade-off between V 200 and mean-tail-dose. The latter is calculated for the 1% of the PTV which receives the lowest dose. Each data point is based on mean values for the ten patients, and results are shown for the tentative dose plans and for the adjusted plans obtained with the objective weights in table 3. The mean value of the tentative plans is shown as a red cross in each plot, while the blue, larger asterisk corresponds to the adjusted dose plans from equal weights on hot and cold spots. The mean-tail-dose is not improved much, neither with tentative plans from the DVM nor with plans from the LPM. For most weights, the value of V 200 is however improved significantly in the adjusted plans as compared to the tentative plans.  . Shows the dose to a two-dimensional cross-section of the PTV for one patient. Here the focus is on cold spots. The dose distributions to the left are tentative plans from the MDVM and to the right are plans obtained from the adjustment step. Table 4. Shows the largest contiguous hot and cold spot (in cc), and the total hot and cold volume (in cc) for the tentative plan DVM and for three sets of weights (number 1, 8 and 15 in    Finally, in tables 7 and 8 we present results for aggregate evaluation criteria and results related to dwell times, in terms of mean values for the ten patients. The mean-tail-dose is calculated as previously stated. The effect from varying the weights is clearly seen in table 7, although the changes are small in some of the evaluation criteria. A large weight on the hot spot component tends to decrease V 200 but increase V 150 , while a large weight on the cold spot component tends to decrease V urethra 120 . The table only shows mean values but these results are consistent for the ten patients. For example, the standard deviations for the differences between the adjusted dose plans and the tentative plans for LPM and DVM were, for all weights, less than 0.5% , 1.5%, and 0.06 Gy for V 100 , V 200 , and  (2), obtained from various weights (w 1 , w 2 ), see table 3, as mean values for the ten patients. The tentative plan is the red cross, with normalized objective value equal to one for both hot and cold spots, and the blue asterisk is the adjusted dose plan with equal weights on hot and cold spots. mean-tail-dose, respectively. Furthermore, for all of these ten patients, for tentative dose plans obtained from the LPM and DVM, and for weights (w 1 , w 2 ) = (1, 0), the values of V 100 and V 200 decreased, which shows consistency of the adjustment step. For tentative plans obtained from the MDVM the standard deviations were slightly higher due to larger variations in potential for improvements, but the results were still consistent.   An increase in the number of active dwell positions is seen in the adjusted dose plans, see table 8, for all weights and for all tentative plans. Such an increase can be considered as beneficial, because it typically results in dose plans which are more robust to for example errors in the catheter placement (Balvert et al 2015). As the sum of the dwell times does not change much, the adjusted dose plans have shorter average dwell times. Also the      (w 1 ,w 2 ) = (1,0), PTV 2: (w 1 ,w 2 ) = (1,0), rectum 2: (w 1 ,w 2 ) = (1,0), urethra 3: (w 1 ,w 2 ) = (1,1), PTV 3: (w 1 ,w 2 ) = (1,1), rectum 3: (w 1 ,w 2 ) = (1,1), urethra 4: (w 1 ,w 2 ) = (0,1), PTV 4: (w 1 ,w 2 ) = (0,1), rectum 4: (w 1 ,w 2 ) = (0,1), urethra (c) Figure 9. DVHs for the three patients. The solid, dashed, dotted and dashed-dotted curves correspond to the clinical plan (referred to as 1), and the adjusted plans with weights equal to (1, 0) (referred to as 2), (1, 1) (3), and (0, 1) (4), respectively. Blue curves (rightmost) correspond to PTV, red curves (leftmost) to rectum, and black curves (in the middle) to urethra. The notations D r 85 and D u 95 are used for the DI for rectum and urethra, respectively. The vertical lines are the dose targets in table 2.
longest dwell time is generally shorter in the adjusted dose plans, and the standard deviation of the dwell times is smaller. To conclude, the dwell times become more homogeneous after the adjustment step.

Clinical dose plans
In the remainder of the section we show results for the adjustment step for clinically used plans. Figure 9 shows DVHs for the three patients. A significant increase of the value of V 100 , which is the main treatment goal, is seen in both figures 9(a) and (c) for the weights (w 1 , w 2 ) = (0, 1) and (w 1 , w 2 ) = (1, 1), respectively. This comes at the cost of a higher dose to the urethra, specially for the third patient. The choice of equal weights on hot and cold spots seems to give a good trade-off between increasing V 100 , when this is possible, and reducing the overdosage to the PTV (V 150 and V 200 ).
Extensive manual adjustments had been made for these clinical dose plan and thus there is only limited room for improvement of spatial properties. In particular, both the largest hot spots and the largest cold spots are very small for all patients. Therefore the first effect, to reduce the number of dose points receiving a dose that is too high or too low (described in the beginning of section 2), is here the dominating one. There are in general improvements of the average doses to the subvolumes that receive the highest and the lowest dose, although the improvements are small. For weights (w 1 , w 2 ) = (1, 1) there is an average decrease in the mean dose to the hottest subvolume of 12% and an average increase in the mean dose to the coldest subvolume of 9%. Also the values of both objective components are improved for weights (w 1 , w 2 ) = (1, 1), with a decrease of on average 90% and 58% for the hot and cold spot component, respectively. The initial objective values are however also small. Further, results similar to those presented in table 8 are seen for the clinically used dose plans, that is, there is a similar pattern in the increase of the number of active dwell positions.
Finally, figure 10 shows visually apparent changes in the dose distribution.

Discussion and future research
The proposed optimization model adjusts a tentative dose plan in order to improve its spatial properties, in terms of contiguous volumes in the PTV in which the radiation dose is too high or too low. The tentative plan can be generated in any way but should have aggregate evaluation criteria, such as dosimetric indices, at acceptable levels. For tentative plans obtained from the optimization models DVM, LPM and MDVM, larger improvements in spatial properties are seen for hot spots than for cold spots. One explanation for this is that there are only small cold spots in half of the studied tentative plans, which is in turn because the models used to generate the plans have a higher priority on giving sufficiently high radiation dose to the PTV than to avoid giving it a too high dose. However, when model MDVM is used to intentionally generate tentative plans with significant underdosage, cold spots are also reduced significantly, which verifies that our optimization model is able to reduce cold spots as well. This capability of reducing cold volumes can also be observed in the adjustment of the clinical plans. Results regarding the dwell time distribution are presented in table 8. It shows an increase in the number of active dwell positions, which is generally considered to be beneficial. The tentative plans obtained from the DVM have fewer active dwell positions and longer maximum dwell times than the plans from the LPM, although the total dwell times are very similar. This is somewhat surprising as the LPM is known to use relatively few dwell positions (Holm et al 2012).
Our spatial adjustment model is non-convex and there are no guarantees to find (or prove) a global optimum. Thus, our results give lower bounds on the potential for improvements of spatial properties. Further, a solution to the model with the sigmoid approximation (6) of the Heaviside step function is not guaranteed to satisfy the dosimetric index constraints, and in a few cases the DI constraint for the urethra was not satisfied. A dosimetric index constraint can also be violated because the sets of dose points used for dwell time optimization and for dose plan evaluation are different. Even though the dosimetric index constraints are sometimes violated, for the two reasons mentioned above, the average dose to the urethra was only affected to a small extent. There are even examples of a decrease in the average dose to the portion of the urethra which receives the highest dose, while the dosimetric index became worse.
The focus in this study has not been the solution method for (2) and there is room for improvements in both computing times and solution quality with a tailored implementation. Computing times of a few minutes are short enough for clinical practice. However, shorter computing times might be advantageous, since it would allow comparing a number of scenarios (obtained by varying the objective weights).
The clinical dose plans which we have studied had been manually adjusted to improve upon spatial properties. There is therefore only a limited room for improvements. Nevertheless, we show in section 4.3 that some spatial properties can still be improved in the adjustment step. For future research, it would be of interest to consider a set of dose plans for which further manual adjustments of spatial properties would be necessary (or desired) to attain clinically acceptable dose plans. This would allow for a comparison of manual and automatic adjustments.
There is a lack of well-defined criteria for evaluating spatial properties of a dose distribution and we therefore analyse spatial properties in several ways, such as the largest contiguous volumes, the worst (hottest and coldest) subvolumes, and the objective values. We have also shown illustrations of reductions in hot and cold spots, see figures 3-4. For future research, it is important to develop well-defined criteria for spatial properties. In addition, there is a need for clinical studies on how the spatial properties of the dose distribution affect the treatment outcome.

Conclusions
We have proposed a general optimization model for adjusting and improving spatial properties of a tentative dose plan, while retaining the levels of the aggregate criteria, specifically the dosimetric indices, obtained from the tentative plan. Although our primary goal is to improve upon spatial properties, the results in section 4.3 show that aggregate criteria can also be improved. The optimization model is based on the assumptions that spatial properties can be improved by reducing the total volume where the dose is too high or too low, and by spreading out dose points where the dose is too high or too low. Our computational experiments show that both these effects are achieved; hence, the model works as intended. Furthermore, we show computing times for the adjustment step of a few minutes.