Objective. Delineating and planning with respect to regions suspected to contain microscopic tumor cells is an inherently uncertain task in radiotherapy. The recently proposed clinical target distribution (CTD) is an alternative to the conventional clinical target volume (CTV), with initial promise. Previously, using the CTD in planning has primarily been evaluated in comparison to a conventionally defined CTV. We propose to compare the CTD approach against CTV margins of various sizes, dependent on the threshold at which the tumor infiltration probability is considered relevant. Approach. First, a theoretical framework is presented, concerned with optimizing the trade-off between the probability of sufficient target coverage and the penalties associated with high dose. From this framework we derive conventional CTV-based planning and contrast it with the CTD approach. The approaches are contextualized further by comparison with established methods for managing geometric uncertainties. Second, for both one- and three-dimensional phantoms, we compare a set of CTD plans created by varying the target objective function weight against a set of plans created by varying both the target weight and the CTV margin size. Main results. The results show that CTD-based planning gives slightly inefficient trade-offs between the evaluation criteria for a case in which near-minimum target dose is the highest priority. However, in a case when sparing a proximal organ at risk is critical, the CTD is better at maintaining sufficiently high dose toward the center of the target. Significance. We conclude that CTD-based planning is a computationally efficient method for planning with respect to delineation uncertainties, but that the inevitable effects on the dose distribution should not be disregarded.
1.Introduction
Accurate delineation of regions of interest (ROIs) is a vital part of radiation therapy treatment planning. Conventionally, delineation is performed by an experienced radiation oncologist in a manual and time-consuming process. Yet, it is sometimes referred to as the weakest link in the radiotherapy chain, and many studies suggest that there is considerable inter-observer variability between clinicians, especially with regard to the clinical target volume (CTV) (Weiss and Hess 2003, Unkelbach et al 2020, Fiorino et al 2020). The CTV is defined in ICRU report 50 (International Commission on Radiation Units and Measurements 1993) as the volume suspected to contain microscopic tumour infiltration with clinically relevant probability. A CTV delineated by a clinician will thus depend largely on two factors: the perceived probability distribution of the tumorous volume and the threshold at which probability of tumor presence is considered relevant. The investigations in this work address the ambiguity associated with the second factor. We will thus assume a known probability distribution over the potential target shapes and assess the merit of a recently proposed approach of moving away from the threshold-based, binary definition of the CTV.
This approach, which accounts for CTV delineation uncertainties explicitly in planning, is to use what is known as the clinical target distribution (CTD). The CTD is a distribution over the potentially tumorous voxels that for each voxel specifies a probability of tumor presence. In treatment planning, one may then use the CTD as voxel-wise weights in the optimization functions, to give higher priority to high-risk voxels. This approach was proposed by Shusharina et al (Shusharina et al 2018) and has since then been explored further by Ferjancic et al (Ferjančič et al 2021) and applied in a robust optimization context by Buti et al (Buti et al 2021). The idea of including voxel-wise probabilities of ROI occupation in optimization was proposed by Baum et al (Baum et al 2006) for managing overlapping margins in prostate treatments. Unkelbach and Oelfke then demonstrated that the approach was equivalent to minimizing the expected value of certain objective functions with respect to the ROI-delineation uncertainty (Unkelbach and Oelfke 2006). Ideally, CTD-weighted optimization will assign dose even to low-probability regions if there is little conflict with sparing organs at risk (OARs), while balancing OAR sparing and target coverage based on the probabilities in regions where the objectives conflict. Compared to more advanced approaches, e.g. the tumor control probability maximization by Bortfeld et al (Bortfeld et al 2021), this method has the advantage that the scaling of voxel weights in optimization preserves convexity and does not introduce any additional computational complexity.
In the present paper, we investigate the implications of CTD optimization compared to using some optimally chosen margin with respect to the underlying tumor infiltration model, the dose deliverability conditions, and the evaluation criteria. The comparison is based on the trade-offs between a target coverage criterion and penalties associated with dose to healthy tissue. The primary target coverage criterion considered is formulated as the probability of (almost) all parts of the target receiving (almost) the prescribed dose. We have suspected that the low voxel-weight toward the edge of the CTD would result in plans with sub-optimal trade-offs between this criterion and the conflicting objectives, and that there rather exists some margin which is more efficient in the described sense. For cases when a proximal, critical OAR does not allow satisfactory values of the primary target coverage criterion, we consider additional criteria. In addition, we view the methods in the light of previously developed frameworks for managing geometric uncertainties, to better understand and compare the methods.
2.Method
2.1.Notation
Any treatment planning problem in the present paper is considered as an optimization problem based on a discretization of the treated geometry into a grid of m voxels. The set of ROIs, denoted by R, is partitioned into the set of OARs and the set of targets, denoted by O and T, respectively. The voxel index set of any ROI r, r ∈ R, is then denoted by 
or more specifically by 
or 
if the type of the ROI is known. Any voxel index 
has a corresponding relative volume vr,i which is the ratio between the volume of voxel i that belongs to r, and the volume of r. It follows that 
. For the purposes of the present paper, it is useful to also define the absolute volumes 
for which it holds that 
where Vr is the volume of r.
2.2.Preliminaries
The dose-at-volume (DaV) will be used to evaluate target coverage. For an ROI r, it is defined as the greatest dose that is received by at least the fraction v of its volume
where the dose vector 
is a function of the optimization variables 
and 
is the feasible set, dependent on the delivery technique. DaV is non-convex and inherently hard to optimize, and will in this work thus only be used as a target coverage evaluation criterion. Instead, the optimization is performed with respect to presumably correlated surrogate functions. In the present paper, the standard notation D100v is used for DaVv when presenting results. To limit the dose to r, we employ maximum dose functions of the form
where 
denotes a reference dose level ideally not to be exceeded and (·)+ is the positive part operator. Function (2) is convex and thus suitable for optimization. Since the focus in this work is on implications on target coverage, maximum dose functions are not only used in optimization but also for evaluation of overdosage. In what follows, analogous notation is used for dose-promoting functions in target volumes.
2.3.Problem setting
We consider phantoms with a known GTV and OAR. The CTD extends away from the edge of the GTV and toward the OAR, and is the region at risk of microscopic tumor infiltration. The external ROI is defined to comprise the full phantom volume. Both 1D and 3D setups are considered to first display the basic properties of the methods in the simplest setting, and then show their implications in a setting which is more realistic. The 1D setup and a 2D slice of the 3D setup are illustrated in figures 1(a) and (b), respectively. To avoid ambiguity when assessing resulting plans, we do not account for overlap between the conflicting regions. Instead, we assume that the OAR acts as an impenetrable barrier for the tumor infiltration, and should ideally not receive any dose.
2.3.1.Tumor expansion model
To model microscopic tumor spread we consider an isotropic infiltration model away from the GTV, proposed in (Buti et al 2021). Therein, the tumor spread can be characterized by a single random variable S with probability density function ρ(s). The probability of the tumor to have spread beyond a certain distance s from the GTV is then given by
2.3.2.Discretization into scenarios
For computational purposes, the distribution of S is discretized into a scenario set 
. For distributions with negligible tail, this discretization involves setting a cut-off value to be considered as the worst case, and renormalizing the scenario probabilities, which are then denoted by 
. In the following, each scenario s maps to a voxel index set 
and the corresponding relative volumes. In 3D, this mapping is defined by the ROI algebra functionality in RayStation (RaySearch Laboratories AB, Stockholm).
2.4.Establishing and managing conflicting objectives
Given the unknown extent of the true target tS , one must define appropriate criteria by which to evaluate a given treatment plan. In the present paper, we assume that any tumorous voxel should ideally receive almost the full prescribed dose, by considering the target coverage metric
which is the probability that the DaV is at least 
. Since our concern is the minimum or near-minimum dose, we typically set both γ and v close to or equal to 1.
There are typically a number of criteria of minimal dose to healthy tissue which stand in conflict with maximizing the target coverage criterion (4). Thus, we add an objective that is a weighted sum of maximum dose functions on form (2), resulting in a bi-objective problem. Because DaV is intractable in optimization, we consider a special case of objective (4) with γ = v = 1 (the probability that the minimum target dose is at least 
). The simplicity of the infiltration model implies that a set of Pareto efficient solutions of the bi-objective problem can be obtained by solving a set of simpler subproblems (see appendix), using the 
-constraint method with a set 
of acceptable values of the target coverage metric (4):
Here, 
is the set of voxels that require dose for j to be a lower bound of objective (4). Although 
is random, each 
is cheaply computed in advance by computing the infiltration distance corresponding to j and generating the resulting target expansion. For computational purposes, we replace the linear inequalities (5b ) by a single function inequality
In summary, we have replaced the true problem of maximizing objective (4) with a surrogate problem (5), in return for tractability. The -constraint methods ensures that the set of solutions to problem (5) are Pareto efficient with respect to the bi-objective problem for the special case of objective (4). The set of associated solutions will act as a reference for the methods presented next.
2.5.Enumeration of CTV margins
Although the -constraint method can be used to bound the target coverage probability in our simple model, it may not be most desirable to do so. First, it is possible that maximizing function (4) with v = γ = 1 requires too much dose, and that striving for a more lenient, clinically relevant goal (e.g. 
) could reduce integral dose and overdosage to OARs. Second, our empirical experience shows that finding a feasible solution to a problem with such a nonlinear constraint (or equivalently a large number of linear constraints) may be computationally inefficient in practice. Instead, we solve a relaxation of problem (5) by moving the left hand side of (3), scaled by an importance weight wt > 0, to the objective. With this relaxation we arrive at conventional margin-based planning where the CTV margin is determined by the probability level j
2.6.The CTD as voxel-wise weights in the objective
We now arrive at the CTD approach which consists of a formulation which is similar to problem (4) but with the first term in the objective replaced by wt fCTD(d(x)), where
in which each 
denotes the probability of voxel i being tumorous. Each pi is calculated from equation (3) using voxel-wise distances to the surface of the GTV.
2.7.The methods in a robust optimization context
Using the scenario discretization in section 2.3.2, the tumor infiltration problem can be seen through the robust optimization framework developed for geometric uncertainties by Fredriksson (Fredriksson 2012). In this framework, margin-based planning is equivalent to planning with respect to a single scenario, or to worst-case optimization if the tail of the probability distribution which falls outside the margin is negligible. Since a larger margin represents a scenario with a target voxel index set that contains those of smaller margins, the largest considered margin will always make up the worst case with respect to minimum dose functions. Varying the size of the margin is then equivalent to varying the length of the neglected distribution tail. Then, a discretization made so that the largest target expansion scenario maps to 
implies that
In contrast, using the CTD objective is analogous to the target presence probability approach in Baum et al (Baum et al 2006). Similarly to what was noted by Unkelbach and Oelfke (Unkelbach and Oelfke 2006), the CTD objective (5) is approximately equal to the expected value of the minimum dose function over the scenario set, making the CTD approach analogous also to the expected-value minimization in Fredriksson (Fredriksson 2012). The approximation is increasingly accurate with the number of voxels m and the number of scenarios 
2.8.Numerical cases
Each method was implemented on a 1D and a 3D phantom, respectively. For the 1D case, a phantom was modeled as the interval [0, L], like in figure 1(a), with L = 18 cm. The GTV and an OAR were represented as the sub-intervals [0, L/3] and [2L/3, L], respectively. The CTD was characterized as an expansion of the GTV by S cm, with S uniformly distributed in the interval [0, L/3]. Dose delivery was modelled with equally spaced Gaussian shaped spots with standard deviation σ = 0.5 cm and inter spot distance 1.5σ. The leftmost spot was centered at l = −1. For each approach, the model was discretized into m = 9000 voxels using MATLAB (2021a, The MathWorks, Inc., Natick, Massachusetts, United States), and solved with the sequential quadratic programming (SQP) solver SNOPT (7.7 Stanford Business Software, Stanford, California). In 3D, a model of a water phantom comprising m = 168 × 177 × 177 voxels, each of size 1 × 1 × 1 mm, was created in a research version of the treatment planning system RayStation. The GTV and an OAR were modelled as two spheres, both of radius 1.5 cm and 0.3 cm apart. The CTD was then characterized by an expansion of the GTV by S from a truncated normal distribution between 0 and 2 cm, with mean 0.6 cm and standard deviation 0.4 cm, followed by a subtraction of the overlap with the OAR. The phantom setup is shown in figure 1(b). To limit bias toward a specific modality, results were generated for both intensity-modulated proton therapy (IMPT), with two laterally opposed beams, and tomotherapy. The problems were solved using 200 iterations of RayStation's internal SQP algorithm, which was sufficient for convergence to a resulting dose distribution. For tomotherapy, the resulting optimization problems were non-convex due to a filtering step: leaf opening times were either set to 0 or bounded from below at iteration 100, depending on their current value. Spot weight filtering was not performed for IMPT, resulting in convex problems.
In section 2.4 we considered a special case of (4) in which γ and v are both close to 1. In this case, the functions used to penalize high dose and their corresponding weights were chosen as in Case A in table 1, There, maximum dose objectives were used to limit the doses to the external and OAR. When evaluating target coverage, the metrics of interest were 
in 1D and 
for γ = 0.95, 0.97, 0.99 in 3D. To address cases when such strict criteria may be physically incompatible with strict limitations on dose to the OAR, the implications of parameter values chosen as in Case B were also investigated. There, maximum dose constraints were used on the external and OAR to bound the dose from above. In both cases, the target metrics for the resulting plans were evaluated by considering 1000 samples of S, for each of which DaV values were calculated. The currently considered target metric was then compared to the weighted sum of maximum dose functions used as dose-limiting objectives, normalized by the maximum across the plans currently under comparison.
Table 1.Parameter values of the dose-limiting maximum dose functions.
| Case | r | wr |  | |
|---|---|---|---|---|
| 1D | 3D | |||
| External | 100 | 100 | 100 | |
| A | External | 0.01 | 0 | 0 |
| OAR | 1 | 0 | 0 | |
| External | ∞ | 101 | 105 | |
| B | External | 0.01 | 0 | 0 |
| OAR | ∞ | 50 | 75 | |
| OAR | 1 | 0 | 0 | |
3.Results
Parts of the results are presented as trade-off curves in some evaluation criteria space, in which a normalized dose-limiting objective (5a ) (to be minimized) is plotted against the target coverage metric (4) (to be maximized) for various values of v and γ. Each curve represents plans for a set of target weights wt . A plan is dominated in the evaluation criteria space if there is another plan under comparison which is not worse in terms of any of the criteria, and strictly better in at least one. Plans which are not dominated are efficient. The normalization of the dose-limiting objective was done by division with the maximum value across the compared plans.
The remaining results consist of comparison of selected dose distributions. For each combination of modality and case, a qualitative selection of similar plans was made to illustrate the characteristics of doses produced by the respective methods. The selection was made by considering the trade-off curves in the evaluation criteria space, and comparing plans in close proximity to each other. In figures 3, 4, 7 and 8, plans selected for comparison are accentuated with an enlarged marker where necessary.
3.1.1D results
Figure 2(a) shows the trade-off curves in the 1D experiment. Notably, for high values of wt , results are similar between the methods, in the sense that the trade-off curves are close in the criteria space. This is emphasized further by the similarity of the dose profiles between the CTD plan and the margin-based plans in figure 2(b). For low wt , the doses from the CTD plan are clearly not efficient due to the early and slow dose fall-off, which can also be seen from the dose profiles in figure 2(c). Going forward, we mainly consider results for values of wt which are empirically determined to result in CTD plans with satisfactory target coverage, replicating the approach in Buti et al (Buti et al 2021). We believe that this is similar to what could be expected from a manual planner in a clinical planning process. Figure 2(d) shows the 1D dose profiles for Case B. Clearly, larger margins suffer more from variation in dose close to the conflict region, while smaller margins result in more hom*ogeneous dose. The dose of the CTD plan is a compromise between the extremes, and better at maintaining high dose toward the center of the target compared to the doses from using large margins.
3.2.Case A: near minimum target dose as first priority
Figures 3 and 4 show the trade-offs of the plans from each of the methods considering a set of various target weights wt , for IMPT and tomotherapy, respectively. For the highest dose threshold (γ = 0.99), the plans from the CTD approach are typically dominated by the margin-based plans, but the difference seems to decrease with increasing wt . In addition, it can be noticed that the plans from the CTD approach are increasingly efficient for more lenient target coverage criteria (γ = 0.97, 0.95), sometimes even dominating the most efficient margin-based plans. Notably, the plans based on formulation (5) are never efficient with respect to the evaluation criteria.
Figures 5 and 6 show example doses and corresponding dose-volume histograms (DVHs) from plans for each of the methods for IMPT and tomotherapy, respectively. The plans were chosen based on proximity with respect to the evaluation criteria from figures 3 and 4. For both modalities, the CTD-based plan achieves better sparing of the OAR for comparable target coverage. However, it seems to come at the cost of increased total dose, as seen especially from the DVH curve of the external for tomotherapy in figure 6(d). For IMPT, the CTD approach gives a noticeably higher dose near the edge of the CTD. For tomotherapy, the dose differences are larger for certain beam angles near the conflict region.
3.3.Case B: strict oar constraint as first priority
Given the constraint on maximum OAR dose, physical limitations on dose conformity make it physically impossible for the target coverage criterion (4) to approach 1 given the values of v and γ in Case A. In addition, the probability of meeting a certain DaV threshold may be insufficient to score a plan in terms of target coverage, given that a satisfactory dose level in the target may vary between high and low conflict regions. Thus, figures 7 and 8 show trade-offs for Case B in terms of the sample 10th percentile and mean target D98, in addition to criteria (4) with v = 0.98, γ = 0.95. Here, wt = 100 was kept constant, as changing wt showed little empirical effect on the trade-offs, likely due to the dominance of the maximum dose constraints constraints compared to the remaining dose-limiting terms in the objective. For both modalities, the CTD plan dominates the larger margins for the three presented target coverage criteria. For the margin-based plans, increasing the margin beyond a certain size often has a detrimental effect on the trade-off.
As for Case A, figures 9 and 10 show doses and DVHs from plans for each of the methods for Case B. For IMPT, the regions in which the doses differ are similar to as in Case A, but the differences are larger in magnitude. The result is a reduced OAR-sparing effect for the CTD approach, as well as an additional increase in the dose difference near the edge of the CTD. The opposite effect is seen for tomotherapy, for which the OAR-sparing effect of the CTD approach is increased, while the total dose increase is decreased.
4.Discussion
In the present work, we have investigated the efficacy of CTD optimization in terms of creating favorable trade-offs between target coverage and overdosage, with respect to the specified criteria. For comparison, conventional margin-based plans were created by enumeration over various CTV margin sizes. The methods were applied to two phantom geometries to investigate the difference in dose between the resulting plans. In this setting, CTD optimization was shown to be slightly inefficient in the case of the strict near-perfect target coverage criteria from Case A. In Case B, considering a critical, proximal OAR, the plan from the CTD approach dominated those resulting from large CTV margins.
The results for Case A agree with the hypothesis that CTD-based doses are somewhat less conformal, given that low-weight voxels at the edge of the CTD drive the optimal solutions toward giving some but not full dose to the outer regions of the target. Another, related effect of using the CTD approach is that of additional total dose, when compared to a margin-based plan selected to be close in the evaluation criteria space. This result is seen across both Case A and B for both modalities, although in varying magnitude. It should be noted that the increase is not only in the aforementioned regions near the edge of the CTD, but also in low-dose regions, as seen from the DVH curves for the external in figures 5(d), 6(d), 9(d), and 10(d). On the other hand, the CTD approach consistently outperforms the margin-based plans when only the dose to the OAR is considered. Whether this is true in general, or dependent on that the weights to dose-limiting functions to the OAR and the external have been held fixed during the experiments, remains open for investigation.
When introducing CTD-weighted optimization, Shusharina et al (Shusharina et al 2018) considered a constraint on mean OAR dose. We think that for the phantom geometries considered in this work, mean and maximum dose constraints have similar effects on the achievable target dose distributions. It is therefore understandable that the results for Case B agree well with those in (Shusharina et al 2018).
Here, we present a plausible explanation of the relative advantage of the CTD approach for Case B. For IMPT, using the CTV0.95 margin results in visible underdosage of a small part of the GTV closest to the OAR, as seen in figure 9(a). In a sense, the effect is similar to the corresponding case in 1D, illustrated in figure 2(d), in which solutions with respect to large margins achieve a little extra dose to the very edge of the margin by allowing a cold spot closer to the center. We believe that this behavior of the dose is largely a result of the quadratic nature of the dose-promoting minimum dose functions, in which minimization of large deviations from the prescription dose is prioritized. Under these conditions, CTD optimization mitigates this effect by prioritizing central voxels when assigning them greater relative weight in the optimization problem. Although the corresponding result is not seen for tomotherapy in the slice shown in figure 10, the probabilities of satisfying 
in figures 7 and 8 suggest that this prioritization of central regions of the target is seen across both modalities.
In the present work, we have disregarded the difficulty in determining the optimal combination of weight and margin. Figure 2(a) shows that in 1D, varying only wt in CTD optimization has similar effect as varying the margin size and wt in the enumeration approach. The trade-offs for Case A suggest that this result also applies for more lenient target coverage criteria in 3D, which is an indication that the CTD approach can produce competitive plans with significantly less computational complexity. The optimal margin is likely to vary with the treated geometry, tumor infiltration model, and delivery technique, making it difficult to determine in advance. This difficulty may increase further in even more complex planning situations such as in robust optimization with respect to setup and range uncertainties. In this sense, CTD optimization may be considered as a simple and practical alternative to finding a desirable combination of both margin size and objective weights, since it only requires tuning of the weight on the target.
The results above apply to the considered isotropic tumor infiltration model. For this model, the information contained in the CTD also fully defines the probability distribution over the sample space of potential target volumes. In a more general setting for target delineation uncertainty, the CTD corresponds to each voxel's marginal probability of tumor presence, which specifies the probability of being tumorous without conditioning on the state of any other voxel. In this sense, it is similar to the concept of spatial probability maps (SPMs) put forward by van Rooij et al (van Rooij et al 2021). In a general setting, marginal probabilites like those from SPMs provide only partial information about the full probability space. Going forward, it may be relevant to investigate the implications of expected value-based optimization (like the CTD approach) with respect to target delineation uncertainty in such a setting.
Data availability statement
The data cannot be made publicly available upon publication because they are owned by a third party and the terms of use prevent public distribution. The data that support the findings of this study are available upon reasonable request from the authors.
: Appendix. Multi criteria optimization
Consider the bi-objective problem
A point 
is Pareto efficient for problem (8) if 
such that fi (x) ≤ fi (x*), i = 1, 2, with the inequality satisfied strictly for at least one of the objectives.
The -constraint method consists of, for a given value of , solving
It can be shown that a solution x* is Pareto efficient with respect to problem (8) if and only if it is a solution to problem (9), with = f2(x*). Thus, by fixing and solving problem (9) for x*, the value of f2(x*) can be controlled while ensuring Pareto efficiency. For a proof and more information on MCO and the -constraint method, the reader is referred to Miettinen (Miettinen 1998).
