A linear optimization model for high dose rate brachytherapy using a novel distance metric

Purpose. We propose a linear network-based optimization model (LNBM) for high dose rate brachytherapy (HDR-BT) that uses a novel distance metric to measure the discrepancy between the dose delivered and the prescription. Unlike models in the literature, LNBM takes advantage of the adjacency structure of the patients’ voxels by formalizing them into a network. Methods. We apply LNBM to a set of 7 cervical cancer cases treated with HDR-BT. State-of-the-art commercial optimization software solves LNBM to global optimality. The results of LNBM are compared with those of inverse planning by simulated annealing (IPSA) based on tumor coverage, dosimetric indices for the critical organs at risk (OARs), isodose contour plots, and two metrics of homogeneity new to this work (hot-spots volumes and diameters). Results. LNBM produces plans with improved tumor coverage and with improved isodose contour plots and dosimetric indices for OARs that receive highest dose (bladder and rectum in this study) when compared with IPSA. Using new metrics of homogeneity, we also demonstrate that LNBM produces more homogeneous plans on these cases. An analysis of the solutions of LNBM shows that they use a significant part of the voxel network structure, providing evidence that the plans produced are different from those created using traditional penalty approaches and are more directly guided by the geometry of the patients’ anatomy. Conclusions. The proposed linear network-based optimization model efficiently generates more homogeneous high quality treatment plans for HDR-BT.

1. Introduction 1.1. Background Brachytherapy (BT) is a procedure that utilizes a radioactive source to target cancer cells. The most prevalent such treatment is high dose rate BT (HDR-BT). In this procedure, a collection of catheters are inserted inside or are positioned in close proximity to the treatment volume. An after-loader is then used to pass a radioactive source through the catheters, pausing for specified amounts of time at certain points in the catheters that are referred to as dwell positions. An important problem is to decide how long to keep the radioactive source at each dwell location (dwell times) to maximize dose delivered to tumors while limiting dose on surrounding organs at risk (OARs).
Historically, there has been two main planning approaches for this problem: forward planning and inverse planning. In forward planning, dwell times are adjusted manually until a plan is found that is acceptable with respect to various criteria set by the oncologist. Finding a good plan using this approach is time consuming as multiple dwell times adjustments are typically required before a final plan is obtained (De Boeck et al 2014). An example of forward planning for HDR-BT is presented in Anacak et al (1997) for breast cancer.
In inverse planning, multiple clinical restrictions are defined for cancer tissues and nearby OARs (De Boeck et al 2014). A mathematical optimization model is then set up and solved using numerical algorithms to produce the treatment plan. Because it is often impossible to create plans that deliver the desired dose exactly, importance factors are given to the different objectives to meet. A plan is then derived that minimizes the total weighted penalties incurred from deviating from desired values. When compared to forward planning, inverse planning tends to consider anatomical information more comprehensively through the optimization process (Reddy 2013).

Brief literature review
Various models have been proposed for inverse planning. To the best of our knowledge, the most prevalent approach to inverse planning is inverse planning by simulated annealing (IPSA) (Lessard and Pouliot 2001). The underlying model seeks to minimize the total penalty of a plan, where a piecewise-linear penalty function is incurred for deviations of delivered dose to voxels from acceptable upper and lower bounds for each tissue. The weights used in the penalty functions might differ by voxel or tissue (Jacob et al 2008). The authors of Pouliot et al (2004) use the same type of penalties with the difference that they apply them to more than one target tissue. In Alterovitz et al (2006), piecewise-linear functions are also used to penalize dose deviations. However, the main contribution of this study is to formulate model as an exact linear programming problem. Because the penalty functions do not need to be symmetric, the model also allows for overdosing and underdosing penalties to be different. It also allows for different penalties to be applied to different tissues. Penalizing the weighted sum of squared deviations between the desired value and the dose delivered at the different voxels has also been proposed; see Milickovic et al (2002), , . The use of penalties reflects the fact that the practical problem being solved is, in fact, one that is multi-objective. In Ruotsalainen et al (2010), the authors develop a model with multiple objectives that simultaneously seeks to maximize the fraction of tumor volume receiving at least the prescription level while minimizing the maximum dose in target tissue and the average dose in OARs.
Most approaches to inverse planning do not directly optimize the isodose contour plots of the produced plans or their dosimetric indices-the evaluation means most often considered by oncologists. Exceptions exists, such as Siauw et al (2011) that imposes dosimetric indices requirements directly as constraints. To the best of our knowledge, the information that a voxel is geometrically adjacent to another is not considered directly in most existing models. We point to Morén et al (2019) for an exception, where a distance-based penalty is applied in the objective to every pair of dose points when they either receive a dose that is too high or too low. The approach we propose goes further in capturing the relative geometrical positions of the voxels. The premise of this paper is that incorporating such information about the anatomy of patients can help better approximate the evaluation means used by oncologists and produce higher quality treatment plans.

Our contributions
We introduce a novel linear network-based optimization model (LNBM) for inverse planning in HDR-BT. This model seeks to minimize a novel distance metric between the prescription and the plan delivered. This metric accounts for spatial information about geometry of the patient's anatomy. It is computed through a network flow subproblem that estimates the minimum amount of modification needed to transform the prescription into the plan delivered by conceptually transporting radiation dose between nearby voxels. We demonstrate numerically that this metrics yield high-quality plans where homogeneity of dose distribution is promoted.
2. Linear optimization model based on novel distance metric (LNBM) 2.1. Linear optimization models for inverse planning Prescriptions for BT cases typically come in the form of a minimum dose desired on tumor voxels along with maximum doses that can be delivered to the voxels of the OARs. Because of physical limitations of the delivery system, doses delivered often deviate from prescribed values. Plans of high quality are those for which the deviation is small.
Models for inverse planning thus often consider a set of dwell positions D that can be activated and a collection of voxels V r that belong to the tumor or to one of OARs. Associated with each voxel i ä V r is a parameter δ i that represents its prescription level and parameters d ik that represent the amount of dose that it receives from dwell position k ä D per unit of time. The problem is then formulated as follows: In this model, nonnegative variable t k describes the dwell time assigned to dwell position k ä D. When selected dwell times are applied to dwell locations, a nonnegative dose y i (which is variable of the model) can be computed at each voxel-node i ä V r , that we call delivered dose. Constraint (2) computes the physical amount of dose that is delivered at voxel i by summing up all dose contributions coming from the different dwell positions that are activated, assuming linear accumulation. Notation δ is used to represent the vector composed of the prescriptions for all of the voxels and y is used to represent the vector composed of the delivered doses y i to all of the voxels. The function f ( · ) that appears in the objective (1) measures the penalty/deviation associated with delivering a dose y compared to the desired value of δ. In many articles, this function is computed by penalizing the individual deviations between y i and δ i over all of the voxels, so that d d ) for predetermined positive weight parameters w i . The single-variable function f i (z i ) is chosen to be either a convex piecewise-linear function (see Lessard and Pouliot (2001), Alterovitz et al (2006), Holm et al (2012)) or a convex quadratic function (see Milickovic et al (2002), , ). When these penalty functions are not symmetric, they can be used to penalize differently over-and under-dosing in the tissues. The geometrical relationships between positions of voxels are, however, not modeled by this sum. Next, we introduce a different form for f ( · ) that explicitly considers the relative positions of the voxels.

A novel distance metric
If we think of delivered doses y and prescribed doses δ as forming distinct distributions, the problem of measuring how much they differ is similar to that of measuring distances between probability distributions. A common way to compute such distance is to use the Wasserstein metric (Panaretos and Zemel 2019). To describe it intuitively, consider the two discrete probability distributions of figures 1(a) and (b). In our setting, the height of each column would correspond to the amount of dose to the corresponding voxel; figure 1(a) would describe delivered doses y in the computed plan whereas figure 1(b) would represent the doses δ desired in the prescription. Although these two distributions differ from each other, we can convert the first into the second by transferring two units from column/voxel A to column/voxel B, which are located approximately three voxels apart, as depicted in figures 1(c) and (d). The total cost of this transfer, which we refer to as transportation cost, can be defined as the product of the amount transported (3 units) by the distance between the origin and destination of the transport ( + 2 1 2 2 distance units). For distributions that are significantly more different, multiple units might need to be transported from multiple origins to multiple destinations. Among all ways to select such moves, we prefer those whose associated total transportation cost is smallest. We refer to this minimum total transportation cost as the (Wasserstein) distance between the two distributions.
Inspired by the above discussion, we propose to measure the distance between a plan and a prescription by computing the minimum cost of (conceptually) transporting radiation dose between voxels to turn the given plan into the prescription. To model this process, we turn the collection of voxels of the problem into a network. By utilizing a network representation, we are able to create a conceptual structure for transporting radiation dose from one plan to another through 'shipments' between nearby voxels, which we refer to as transportation flows. The search for minimum transportation flows can then be modeled as a network flow optimization model, with the delivered plan representing radiation supply and the prescription plan representing radiation demand; see Ahuja et al (1993) for an introduction to network flows. We emphasize that the proposed movement of radiation from voxel to voxel is purely conceptual and that we are not claiming to physically transport radiation dose inside of a patient. Rather, this conceptual point of view allows us to measure distance between plans and, hence, allows us to seek a treatment plan that is 'closest' in distance to the prescription.

Network structure
A crucial component of the procedure described above is the specification of the different ways in which dose can be transported between voxels. If we allowed for direct transport from each voxel to each other voxel, we would obtain a bipartite network that has as many nodes in its left and right partitions as there are voxels in the plan; see figure 2 for a graphical depiction where only 6 voxels and their connections are represented. Because for practical cases in inverse planning it is common to consider upwards of 10 000 voxels, this would lead to a large network with approximately 100 000 000 connections between left and right partitions. The resulting model, although simple in structure, is too large to be solved in reasonable time using commercial software.
To circumvent this difficulty, we propose to model the transportation of radiation using networks that have fewer arcs and take into account the relative geometrical positions of the different voxels and the fact that only certain voxels (those corresponding to OARs and tumors) within a patient are considered in optimization models. Depending on whether adjacency is considered only between voxels of a CT-scan slice, or between voxels present in successive CT-scan slices, we obtain two different networks.
The first, which we refer to as 3D-network, is built as follows: • Nodes: we create a node for each voxel in each slice. We refer to these nodes as voxel-nodes. Voxel-nodes can either be part of OARs or the target. We also create a collection of artificial voxels, the number of which is at the discretion of the user, and whose purpose we describe later. We denote the voxel-nodes by V r and the artificial nodes by V a .
• Arcs: we create directed connections (arcs) between certain nodes of the network. We include an arc between two voxel-nodes if their associated voxels touch (i.e. they have a face in common). A network with only these arcs could be disconnected when OARs and targets do not touch. Similarly, a restriction of this network to a slice could also be disconnected; see figure 2(a). To ensure that the network is sufficiently connected, we proceed as follows. For each voxel-node v inside of the tumor and for each OAR, we consider the six principal directions d (west, east, north, south, up, and down) and determine the closest voxel w d in the OAR that can be reached from v by moving in the direction d. If such voxel exists, we include an arc between v and w d . It follows that each tumor voxel v has no more than six arcs that connect it to each OAR. We do not include the reverse arcs since we do not wish to allow excess dose at OARs to serve a way to mitigate a deficit of dose in the tumor. We refer to all arcs created above as regular arcs. Finally, we create a pair of arcs with opposite directions between every voxel-node and every artificial node. We refer to the arcs incoming to an artificial node as overdosing artificial arcs and refer to the arcs outgoing from an artificial node as underdosing artificial arcs. We denote the set of regular arcs as A r and the set of artificial arcs as A a .
• Arc capacities: we associate a transportation capacity with each of the artificial arcs. This capacity can be selected to be any nonnegative real number or can also be selected to be infinite. We denote the capacity of arc (i, j) ä A a by u ij .
• Arc costs: we associate a nonnegative unit transportation cost to regular and artificial arcs. For regular arcs, this transportation cost is selected as a function of the distance between the two corresponding voxels and the type of voxels they are (OAR or target). Making this transportation cost a function of the distance has the advantage of capturing how far apart voxels are when they are not adjacent. Making this transportation cost a function of the type of voxels has the advantage of allowing penalties to be tailored to different tissues. We denote the transportation cost of arc (i, j) ä A as C i,j . We provide details and intuition for setting transportation costs in section 3.1.
An advantage of the definition of the network presented above is that it is relatively sparse in terms of the number of arcs that it creates and nodes that it uses. Further, its structure captures the relative position of the voxels with respect to each other within each slice and within nearby slices. A graphical depiction of the part of the network corresponding to a slice is given in figure 3. In this figure, we see that voxel-nodes in the upper boundary of the sigmoid and the rectum are connected with all of the tumor voxels located to their north. Although only a subset of these arcs are necessary to connect the tumor to OARs, including them all appears to be beneficial when solving our models. We also mention that, as is common in the literature, we only consider voxels that belong either to the tumor or to OARs, as this keeps models smaller. A second approach to create a network is described next. In this network, which we refer to as 2D-network, transportation along regular arcs of the network can take place between voxels if they belong to the same slice of the patient. The 2D-network is built in the same way as the 3D-network with the exception that regular arcs between slices are removed.
We refer to a network N created in any of the above ways as N(V, A) where V is the collection of its nodes (V = V r ∪ V a ) and A is the collection of its arcs (A = A r ∪ A a ).
We next comment on how this network is used to compute distance between plans. To each voxel-node of V, we associate a supply and a demand. We set the demand of node i to be the prescription δ i at voxel i. We let the supply at each voxel-node to be the delivered dose y i . Whatever discrepancy exists between the supply and demand at each voxel-node must be transported across the arcs of network N(V, A), respecting their capacities, to even out supply and demand. Because traditional network flow models require that flow is conserved, this transportation problem would typically prove infeasible if the total supply did not exactly match the total demand. This is the reason for the existence of artificial nodes that can be used, at will, to eliminate the discrepancies between supply and demand that cannot be handled by the rest of the network. The transportation cost of artificial arcs involving voxel-node i can then be understood as penalties for discrepancies observed between dose delivered and prescription at voxel i. More generally, arc transportation costs act as a way to measure the amount by which delivered plan and prescription differ.
2.4. Linear optimization models for inverse planning with novel distance metric We next modify Models (1)-(4), replacing the function f ( · ) used in its objective with the distance metric we just described. Using the network notation introduced in section 2.3, this linear network-based optimization model (LNBM), which takes the same form independent of whether N(V, A) is a 3D-or a 2D-network, is as follows: This model has variable x ij to model the conceptual transport of radiation dose between the delivered plan and the prescription over the arcs (i ,j) of the network. The cost associated with this transport measures the distance between the plans. It computes transportation costs along the regular arcs and along the artificial overdosing/ underdosing arcs of the network. Hence, objective function (5) seeks to minimize distance between plans by minimizing the total transportation cost associated with converting the delivered plan into the prescription. The left-hand-side of (7) computes the difference between the prescription and the dose delivered at voxel i. When this difference is positive, voxel i is underdosed, and the rest of the network (including artificial nodes) must provide for the deficit of dose through the arcs that are incoming to voxel i. When this difference is negative, voxel i is overdosed, and the rest of the network (including artificial nodes) must help diffuse the excess of dose through the arcs that are outgoing from voxel i. Constraint (8) limits the amount of dose that can be transported on artificial arcs. Constraints (9)-(11) require that dwell times, doses delivered, and doses transported are nonnegative. If we ensure that, for each voxel-node, there is at least one outgoing and one incoming artificial arc with infinite capacity, then LNBM always has an optimal solution because it is feasible and is not unbounded. By selecting the transportation costs on regular arcs to be very large compared to those of artificial arcs, the model can be shown to reduce to (1)-(4) in which the penalty functions f i ( · ) described in section 2.1 are piecewise-linear convex functions with a number of pieces related to the number of artificial arcs; see Supplementary Material A for a discussion. When transportation costs on regular arcs are chosen smaller, LNBM differs from penalty models in that it incorporates spatial information about the voxels through the network representation. Since the model definition then implies that it is more cost-effective to redistribute dose locally rather than over longer distances, the network model should favor dose distribution homogeneity.

Methods and materials
We use Institutional Review Board (IRB) approved data, CT scans and contours for 7 sequential patients with cervical cancer obtained from the Department of Radiation Oncology at the University of Minnesota. The corresponding anonymized DICOM files are accessed through MATLAB R2019a. For these cases, there are three OARs that must be protected from radiation (rectum, bladder, and sigmoid) and one target. For the patients in our study, tumor sizes average 83.275 cm 3 and range from 70.940 to 117.838 cm 3 . The images are voxelized through MATLAB codes we developed. To create the data our optimization model requires, we chose voxels of dimension 3 mm × 3 mm × 3 mm to keep models tractable. This leads to a number of voxels ranging from 12, 424 to 14, 435. To evaluate the quality of the plans developed (e.g. to compute DVHs and dosimetric indices), we use a finer voxelization with voxels of dimension 1 mm × 1 mm × 1 mm. The total number of voxels (including

Parameters selection
For the ensuing computational tests, we use networks with a single artificial node. We set the capacities of all artificial arcs to be infinite. This intuitively implies that the artificial arcs model a v-shaped penalty function with two slopes. Setting the transportation costs of regular arcs is important to complete the definition of the penalty function and to tailor it to the application. In our experimentation, we have found that the transportation cost structure used does not need to be complicated in order to obtain high quality plans. In particular, we have been successful at obtaining good quality plans using the following transportation cost structure that relies on only five parameters.
To promote homogeneity of the tumor dose, we discourage the transportation flow of radiation inside of the tumor by setting the corresponding transportation costs to be a large number M. Similarly, to prevent overdosing of the OARs and underdosing of the tumor, we discourage the transportation of flow of radiation from OAR voxels to tumor voxels by setting the corresponding transportation costs to be M. For the other regular arcs, we select the transportation cost to be a multiple C of the distance L ij between the corresponding voxels. We create transportation costs that depend on the distance between voxels to limit the ability of the model to reconcile dose deviations by assigning them to remote voxels. We set the transportation cost of artificial arcs to one of four constant values based on their type. The artificial overdosing arcs from OARs and from the tumor are assigned transportation costs of CO max and CT max , respectively. The artificial underdosing arcs from OARs and from the tumor are assigned a transportation cost of 0 and CT min , respectively. Table 1 summarizes this transportation cost structure. In this table, we use V n to denote the collection of voxels corresponding to OARs and V t to denote the collection of voxels corresponding to the tumor. We have V r = V t ∪ V n . In our computation, we use the following numerical values: = CO 1000 max , = CT 10 max , = CT 1000 min , and C = 1. We then choose M = 10 000, so that it is larger than all other transportation costs. Even though these constants are common across patients, the penalties imposed are individualized to each patient since they depend on the network structure of OARs and tumor voxels, which is unique to each patient.
We show in section 4 that this choice of parameters consistently produces plans of similar or higher quality than those produced using IPSA. The ease of creating such high-quality plans without fine-tuning is a notable advantage of LNBM.
While the above discussion provides one way to select model parameters, there are many alternative choices. For example, adjusting transportation costs for a specific substructure could increase (or decrease) the penalty

For arc
for deviating from the prescribed dose for that structure relative to others. Moreover, the introduction of artificial nodes, combined with changes in arc capacity, can create more intricate piecewise-linear penalty functions with varying slopes for different substructures, thereby offering even more precise control over dose deviations in particular substructures.

Method evaluation
We compare the plans produced by our model with those obtained using IPSA. The version of IPSA we consider is that which is used in practice by the medical physicists who created the patients' plans at the University of Minnesota. Its penalty parameters are selected to be the default parameters used by the physicists to obtain the original treatment plans they generate for the oncologist. These plans might not be the ones ultimately used for the patients as the final plans are obtained after a set of interactions between the oncologist and the physicist that would be difficult to replicate in this study. Further, to ensure consistency in the comparison of results, we compute the dosimetric indices of the IPSA plans directly from their dwell times using our MATLAB codes. The differences between the indices we compute and those generated by IPSA are small and we attribute them to the difference in voxelization used in the two codes. We acknowledge that, by adjusting the penalties in IPSA, it might be possible to produce different (possibly better) plans. However, we point out that the same can be said for LNBM. Given this situation, one could think of comparing the best possible plans obtained with IPSA and LNBM, respectively. This, however, is a difficult task since the problem has many objectives. We therefore believe that comparing a version of IPSA whose default parameters have been tuned at the clinic over many years for the particular type of cancers we study and a naive implementation of LNBM provides a reasonable assessment of the potential of our model. We evaluate plans using three approaches. The first relies on dosimetric indices, the second on isodose contour plots, and the third, which is new to this work, aims to measure the prevalence of hot-spots in the plan.

Dosimetric indices
In our evaluation, we use the dosimetric indices VT 100 , DT 90 , DB 2cc , DB 0.1cc , DR 2cc , DR 0.1cc , DS 2cc , DS 0.1cc , where T stands for Target, B for Bladder, R for Rectum, and S for Sigmoid, since they are the dosimetric metrics used in the practical evaluation of the cases in our data set, consistent with Pötter et al (2018). These measures are extracted from dose-volume histograms (DVHs) (Drzymala et al 1991). For q ä (0, 100), D q is the dose d (measured as a fraction of the prescription dose) such that q% of the structure receives a dose higher than or equal to d, and (1 − q)% of the structure receives a dose lower than or equal to d. Similarly, for x 0, D xcc is the dose d (measured as a fraction of the prescription dose) such that x cubic centimeters of the structure receives a dose higher than or equal to d, and the rest of the structure volume receives a dose lower than or equal to d. For a dose d (measured as a fraction of the prescription dose), V d is the fraction of the volume receiving at least d% of the prescribed dose.
Hot-spots evaluation Hot-spots, which we think of as connected collections of voxels receiving much higher doses than desired, are common in HDR-BT as dose accumulates at every dwell location where the source stays sufficiently long. This feature makes it difficult, or even impossible, to obtain homogeneous plans. We hypothesize that, although LNBM does not contain explicit constraints to do so, it produces plans that are more homogeneous and have less pronounced hot-spots compared to IPSA.
We introduce next measures to evaluate the preponderance of hot-spots in treatment plans. We say that a voxel is an α-hot-voxel if the dose it receives is at least a multiple α (which we call hot-spot factor) of the average dose δ in the tissue. We use average dose in the tissue instead of prescription to keep the measure relevant in plans where the tumor is under-dosed. In our computation, we use two values of α, namely α = 1.5 and α = 2. When two α-hot-voxels v b and v e are connected, in that there exists a collection of α-hot-voxels v 0 , K, v k with v 0 = v b and v k = v e such that voxels v i and v i+1 are geometrically adjacent for i = 0, K, k − 1, then we say that they belong to the same hot-spot. Hot-spots partition the set of α-hot-voxels into (non-connected) clusters. We define the volume of a hot-spot to be the sum of the volumes of all the voxels it contains. We define the radius of a hot-spot to be the maximum distance between two points in the voxels it contains. A cumulative frequency distribution plot for the volume of hot-spots, CF V , is a graph that indicates, for a volume v measured in cc, the number of hot-spots with volume smaller than or equal to v. A cumulative frequency distribution plot for the diameter of hot-spots, CF D , is a graph that indicates, for a diameter d, the number of hot-spots with diameter smaller than or equal to d. Figure 4 displays CF V and CF D for a treatment plan whose hot-spots are depicted in a 3D plot that can be found in figure B1 in Supplementary Material B. In figure 4(a) we can see that 73 hot-spots throughout the tumor regions have a volume of 1.5cc or less. We can also see in figure 4(b) that 70 hot-spots have a diameter of 20 mm or less. This illustrates that only two hot-spots are relatively large for this plan.
Hot-spots can be easily derived from information about the proximity of voxels and the dose they receive. We first inspect all of the voxels in the tumor, as this is the structure for which we want to monitor hot-spots. For each α-hot-voxel i, we create a set S i that contains voxel i together with all of its direct neighbors (up, down, west, east, north, and south) that are also α-hot-voxels. We then sequentially merge all such sets S i and S j that have a nonempty intersection. Once this step is completed, we have obtained a collection of sets of voxels that are such that (i) each pair of voxels within a set is connected by a path of adjacent α-hot-voxels, and (ii) there is no path of adjacent α-hot-voxels between voxels belonging to different sets. Hence, each of these sets corresponds to a hotspot. The volume of the hot-spot is then computed directly from the number of voxels it contains. The diameter of the hot-spot is computed by considering each pair of voxels it contains.

Results
We compare the quality of LNBM plans with those obtained using IPSA. Comparing plans can be difficult when multiple dosimetric indices are computed and these indices are all significantly different. Comparison, however, becomes clearer if the plans exhibit the same value for one specific dosimetric index. For this reason, we compare plans after they are normalized around one specific dosimetric index, which we choose to be DT 90 in our study. Normalization around DT 90 can be achieved by scaling the dwell times of one plan so that the value of the DT 90 after scaling is close to the value it takes in the reference plan. After normalization, the quality of the two treatment plans can be compared based on the other metrics. The maximum normalization discrepancy across all cases in our study is 0.3%. This discrepancy is calculated as

Comparison between LNBM with 3D-and 2D-networks
First, we investigate whether using 3D-networks or 2D-networks in LNBM provide plans that significantly differ from each other. In principle, using a 3D-network could yield better plans as it allows for a wider range of options for the model to transform delivered dose into the prescription. We use LNBM to obtain treatment plans for all of the patients in our data set. We create the networks so that they contain a single artificial node. We set the transportation costs for all arcs as described in section 3.1. Once plans are obtained, we normalize the results with respect to DT 90 . The computed dosimetric indices are reported in tables 2 and 3. On average, we see that VT 100 and DT 90 change by −0.17% and −0.1%, respectively, when going from a 3D-to a 2D-network whereas the metrics for Bladder (DB 2cc and DB 0.1cc ), Rectum (DR 2cc and DR 0.1cc ) and Sigmoid (DS 2cc and DS 0.1cc ) change by (−0.18% and −0.33%), (0.09% and 0.76%), and (−0.05% and −0.49%), respectively. We conclude that the dosimetric indices obtained with both approaches are very similar since the maximum difference is always less than 1%. Further, when superimposing dwell times prescribed by the two plans on their dwell positions, plans appear quite close in the dwell positions they select and the dwell times they use; see figure C1 in supplementary material C.1 for details. Finally, we observe that both model variants produce solutions reasonably fast. The time required by CPLEX to optimize LNBM averages approximately 40 seconds on our seven test instances. Detailed statistics can be found in table C1 in Supplementary Material C.2 for both 3D-and 2D-networks. These times do not include the time required for data processing, including voxelization and computation of the dose parameters d ik . On average, the model produces solutions faster when a 2D-network is used. Because the plans they produce are relatively similar, we therefore recommend the use of LNBM with a 2D-network and we focus the rest of our computation on 2D-networks.

Comparison between LNBM with 2D-network and IPSA
In the ensuing subsections, we compare various aspects of the plans obtained with LNBM (using a 2D-network) and of those obtained with IPSA. Before we do so, we observe that both methods generate short treatment plans, as measured by the total activated dwell times and that IPSA generates slightly shorter treatment plans on average; see table 4.

Dosimetric indices
In table 5, we present dosimetric indices and their mean (Mean), minimum (Min), maximum (Max) and standard deviation (SD) for all patients analyzed in this study for both LNBM and IPSA plans. In the plans produced by IPSA, bladder and rectum are the two OARs that receive the highest doses whereas the sigmoid dose tend to be very low. Plans produced by LNBM tend to be more homogeneous and reduce the Max values on most OARs. Since SD values for all dosimetric indices are smaller for OARs in LNBM compared to IPSA, we also infer that LNBM produces plans more similar across patients. Further, tumor tissues receive more dose in LNBM as Min values for tumor tissues are higher in LNBM. Hence, LNBM achieves lower dose on the bladder    and rectum at the expense of an increase in the dose delivered to the sigmoid, which remains low. Figure 5 gives a visual comparison of the dosimetric indices for the plans obtained with LNBM and IPSA. In table 6, we summarize the performance of the two models based on the mean difference of their dosimetric indices (using IPSA as the base method). On average, LNBM improves the tumor coverage VT 100 by 1.1%, DB 2cc by 4.96% , DR 2cc by 5.91%, DB 0.1cc by 19.96%, and DR 0.1cc by 19.43%. Evaluating the two models using the statistical test used in (Holm and Larsson 2013) on different dosimetric indices, we see that the difference in metrics for LNBM and IPSA are significantly different (p-value >0.05).

Isodose contour plots
We present isodose contour plots in figures 6-7 for selected slices on the top and middle sections of Patient 3. Plots for the other patients can be found in supplementary material D and are qualitatively similar. Four isodose contour lines are depicted for each slice. The outermost line connects voxels receiving 50% of the prescription level, the second outermost line connects voxels at the prescription level, the third outermost line connects voxels receiving 150% of the prescription level. The innermost line, when it exists, connects voxels receiving 200% of the prescription level. From these figures, we observe that, for the middle layers (those represented in figure 6) the curves produced by both methods are similar. For the top layers (those represented in figure 7), LNBM is more successful at delivering the prescription to the tumor, since the prescription level curve contains the tumor in the second row whereas it does not in the first row.
A careful observation of the figures also shows a difference in the hot-spots present in the different plans. In particular, the areas of the tumor that receive a dose of more than 150% or 200% of the prescription seem to be more widespread in the plans generated by IPSA than in the plans generated by LNBM for the middle layers; see the first row in figure 6. In contrast, the areas receiving higher doses in the top layers tend to be larger in the plans generated by LNBM than those generated by IPSA, partly because LNBM has better tumor coverage in these regions; see figure 7. Similar behavior is observed for the other patients. The presence of many disconnected hotspots is, however, a characteristic likely to make the plans generated by IPSA less desirable than those generated by LNBM as it corresponds to very inhomogeneous plans.  Figure 8 presents CF V and CF D plots for two of the patients (Patients 1 and 2) in our data set when the hot-spotfactor α is set to 1.5. Plots for Patient 3 are similar to those of Patient 1 and plots for Patients 4, 5, 6 and 7 are similar to those of Patient 2. They can be found in supplementary material E. Figure 9 presents CF V and CF D plots for three of the patients (Patients 1, 2 and 6) in our data set when the hot-spot-factor α is set to 2. Plots for Patient 3 are similar to those of Patient 1. Plots for Patients 4 and 5 are similar to those of Patient 2, and plots for Patient 7 are similar to those of Patient 6. They can be found in Supplementary Material E. Dashed curves correspond to plans obtained with IPSA. Solid curves correspond to plans obtained with LNBM. In these figures, we prefer to see curves that do not rise as high, as they correspond to plans with fewer hot-spots. We also prefer curves that do not reach as far to the right as they correspond to plans that have smaller hot-spot volumes and hot-spot diameters. We see that LNBM most often reduces the maximum diameter and the total volume of hot-spots when compared with IPSA. For most of the patients, LNBM tends to generate smaller hot-spots as solid curves, with limited exceptions, are shorter than dashed ones in both CF V and CF D plots. Further, for most of the patients, LNBM tends to generate fewer hot-spots as solid curves lay mostly below the dashed ones in both CF V and CF D plots. Even in the few cases where the solid curve is above the dashed one, we see that LNBM produces much smaller hot-spots whereas IPSA creates fewer but much more voluminous hot-spots.

Network analysis
In section 2.4, we have mentioned that, when transportation costs for the regular arcs of the network are large, LNBM reduces to traditional penalty models; see Supplementary Material A for a more detailed discussion. In this section, we analyze how much of the network structure (not including artificial arcs) is used in the plans that we derive. Analyzing the magnitude of transportation flows assigned to variables x ij in solutions to LNBM is an indication of how much they can differ from those obtained with simpler penalty models.
In figure 10 we represent the network transmission variables x ij in LNBM for two slices of Patient 3; the pictures for other layers and other patients show similar patterns. In this figure, the thickness of the arcs of the network is proportional to the value of the transmitted dose/transportation flow over that specific arc. If no dose is transmitted on an arc, then this arc is not represented. In the plan that LNBM produces, almost every layer takes advantage of network arcs by transmitting dose from hot regions to the cold tissues. Figure 11(a) depicts the collection of all regular arcs and voxel-nodes in the network of Patient 3. A careful inspection of this picture reveals that this network does not have arcs connecting consecutive horizontal slices. This is because we use a 2D-network for this case. Figure 11(b) shows the arcs of the network that are used in the treatment plan generated by LNBM for that patient. We observe that LNBM uses many of its arcs. This in turn suggests that these plans would have been difficult to obtain using single-voxel-based penalty models. An analysis of the differences between dwell times and positions selected by LNBM and those selected by IPSA can be found in Supplementary Material F.

Conclusion
We introduce a new model for inverse planning in HDR-BT. This model seeks to minimize the distance between the plan delivered and the prescription using a concept for distance that has not been used in this setting. In particular, the distance is measured by the change in voxel doses required to turn the delivered plan into the prescription. Unlike traditional penalty models, LNBM considers the relative positions of the voxels. LNBM has several advantages and contributes to the literature as follows: • LNBM is versatile and its parameters can be tuned to yield other existing penalty-based approaches. As a result, it can produce plans generated by these models.
• The geometrical interpretation of the parameters of LNBM makes it easy to tune. In particular, high quality LNBM plans can be obtained without fine-tuning of the parameters.
• On a test set of seven clinical cases of HDR-BT for cervical cancer patients, the solutions obtained with LNBM tend to outperform those obtained using IPSA on a large array of measures that are commonly used to evaluate these plans (dose-volume histograms, tumor coverage, dosimetric indices for tumor regions and critical OARs, K) and isodose contour plots.
• LNBM produces plans that are more homogeneous when evaluated through DVH-inspired graphs we introduced to capture the preponderance and magnitude of hot-spots in plans. We believe that this is due to the way we measure distance between plans, which takes into account the relative position of voxels.
• LNBM is tractable and sufficiently fast to be used in practice, as globally optimized plans can be easily obtained using commercial solvers such as CPLEX.
We conclude by mentioning that, even though this paper focuses exclusively on HDR-BT, we envision that the distance metric between plans we introduced could also be applied to other radiation planning problems.