Bi-Level Optimization to Enhance Intensity Modulated Radiation Therapy Planning

RT planning can be viewed as an optimization process, in which the fluence maps, represented by vectors of non-negative numbers $x$, define the radiation intensities of individual beamlets. Deposition matrices $D$ translate fluencies $x$ to doses deposited in voxels, so that the doses in voxels are computed as product $Dx$. The RT planning goal is to compute fluencies $x$ that deposit prescribed and homogeneous doses to PTVs and acceptable doses to OARs.

gEUD is a biology-motivated measure to evaluate radiation effects, based on the concept of the uniform radiation dose delivered to a patient organ, that causes the same effect as a nonuniform dose (Niemierko, 1996; Q. Wu and Schmidt-Ullrich, 2002). In the case of gEUD, radiation effects in a PTV or an OAR, both referred to as structure $s$, are evaluated by the following function that aggregates these effects over all voxels belonging to the structure $s$:

where $|M_{s}|$ is the number of voxels of the structure $s$; $d_{j}(x)=D_{j}x$ is the radiation dose deposited in voxel $j$ of structure $s$ by fluence map $x$ and the parameter $a_{s}$ represents the radiation effect on the structure. Its value can be taken from the literature or it can be adjusted individually by trial and error for each patient case.

According to Q. Wu and Schmidt-Ullrich (2002), clinically meaningful RT plans can be obtained by computing the maximum of the following function $F(x,\phi)$, built over the gEUD:

where $EUD_{t}^{0}$ is the prescribed dose for $t$-th PTV, ${EUD_{r}^{0}}$ is the maximum uniform dose at $r$-th OAR; $n_{r}$ and $n_{t}$ express the importance of the prescriptions for the corresponding structure; $\phi$ represents the set of parameters involved in the $F$ definition, i.e., $\phi$ is an instance of parameters $n_{t},n_{r},a_{t},a_{r}$ and $EUD_{r}^{0}$ with $t\in T,r\in R$.

As is, the objective function $F(x,\phi)$ only controls underdosage inside PTVs. However, overdosage control in those volumes is also important. For this purpose it is common to define complementary structures, introduced as "virtual PTVs" on Q. Wu and Schmidt-Ullrich (2002), that are treated as OARs. In this work, for each PTV we have defined a virtual PTV. To lighten optimization costs, we have interrelated the respective parameters, namely, for each PTV $t$, the corresponding virtual PTV $r$ has $EUD^{0}_{r}=EUD^{0}_{t}+1$, $a_{r}=-a_{t}$ and $n_{r}=n_{t}$.

In Choi and Deasy (2002) convexity of function (1) was studied. It was shown that for a range of parameters, function (1) is a convex function of fluence maps ($x$). In Romeijn et al. (2004) the convexity of (2) was analyzed and the conclusion was that it is convex for the range of values of the parameters $a_{s}$ considered in practice. This means that IMRT RT plans can be efficiently obtained by gradient methods, as suggested in Q. Wu and Schmidt-Ullrich (2002).

RT planning based on maximization of function $F(x,\phi)$ goes as follows. At the start, the values of $n_{t},n_{r},a_{t},a_{r},\ t\in T,\ r\in R$, are selected according to suggestions from the literature. Next, in successive planning cycles, the parameters $a_{p}$ and $EUD_{p}^{0},\ p\in P$, are manually tuned by the RT planner. This is a tiresome and time-consuming process, requiring high-expertise planners. In this work, we propose to select the gEUD parameters automatically by multi-objective optimization. In the next section, we present this idea in more detail.

Automated parameter tuning can result in several, even many, RT plans generated in one planning session. We present a multi-objective optimization model which serves that purpose. We also outline a novel approach to solve this model by a hybrid optimization scheme, coupling exact and heuristic (evolutionary) optimization methods.

To maximize function (2), the Gradient Descent method has been custom-coded and the resulting algorithm (GD) is the backbone of the EUDGD module.

The EUDGD module receives as inputs for each patient case: the patient data, the deposition matrix, the beamlet geometry, the ROI segmentation, parameters $\phi$ and the number of steps. It returns the optimal fluence $x^{*}(\phi)$. To generate feasible fluence maps, the EUDGD module needs additional configuration parameters. In the current software release these parameters are set to default values, but they can be customized. The default GD parameter values reported in this work are: $20000$ descent steps, $2e^{-7}$ step size, and $0.3$ as the maximum beamlet intensity limit. Lastly, a 3x3 smoothing kernel has been applied to all beamlets after every descent step to facilitate the delivery of fluencies defined by RT plans by radiation equipment.

The number of descent steps required by EUDGD to provide acceptable fluence maps and the size of the auxiliary data structures result in high computational and memory demands. To maximize function (2) in reasonable times, we have applied HPC software and hardware techniques, as described in Moreno et al. (2022).

For fluence maps $x^{*}(\phi^{k})$ delivered by the EUDGD module, the EvalMO module computes the values of functions $f_{0}(x^{*}(\phi^{k})),\dots,f_{|P|}(x^{*}(\phi^{k}))$ and sends them to EvolTuning module.

The EvolTuning consists of a multi-objective evolutionary optimization algorithm that derives parameters $\phi$. It is initialized with collections of $n_{s},a_{s},EUD^{0}_{s}$, $s\in T\cup R$. In order to generate clinically acceptable plans and reduce the search space of the evolutionary process, feasibility ranges for each kind of parameter are to be set accordingly.

For the current implementation, the MOEA/D algorithm has been selected and seeded with parameter values recommended in Zhang and Li (2007), with the exception of the number of generations (epochs) that is set to $50$ and the cardinality of populations set to $150$.

EvolTuning uses the evaluations received from EvalMO to generate new collections of (potentially improved) parameters $\phi$. The process stops when the limit of generations is reached. As the result, $150$ mutually nondominated collections $\phi$ of parameters of function (2), each collection yielding radiation plan in the form of fluence map $x^{*}(x)$, 150 plans in all, constitute an approximation of the set of Pareto efficient collections $\phi$ (in terms of functions $f_{0},\dots,f_{|P|}$).

To avoid overwhelming the RT planner with so many plans, the module Decision Tool reduces the set of mutually nondominated collections $\Phi$. It performs this task as follows:

First, plans that minimize objective functions $f_{0},f_{1},\dots,f_{|P|}$ separately are selected, and they define the extreme points of the respective Pareto front. Next, a limited number of plans, fairly distributed between those extreme points, are added to form the reduced approximation of the Pareto front. Observe that even if $f_{0}$ for some plan takes a positive but small value (meaning that there are some constraint violations), it can be still a clinically viable option as long as it offers plausible values of $f_{1},\dots,f_{|P|}$ at the price of a small constraint violation.

The gEUD-based RT planning system, described in the previous section, consists of multiple modules connected by a messaging system based on the ZeroMQ library (ZeroMQ, 2023).

The EvolTuning module is implemented in Java and it takes advantage of the well-known jMetal framework (Durillo and Nebro, 2011). This framework provides well-crafted implementations of multiple algorithms, allowing the user to select the best one for a given problem. Although in this experimentation we have used the MOEA/D algorithm, this module system allows us to swap the evolutionary algorithm with minimal changes.

For the EUDGD and EvalMO modules, we have developed two versions: a CUDA C (NVIDIA, 2023) version for NVIDIA-based GPUs, and an OpenMP C (OpenMP, 2023) version for multicore CPUs. These two versions can be used interchangeably or even in parallel to fully exploit all the available resources of the given platform.

Finally, the Decision Tool is a Python-based script that is called after the EvolTuning process terminates. All the figures shown in this section are generated using Python tools implemented in this module.

Our experiments have been run on the High-Performance Cluster of the SAL (Supercomputación–Algoritmos) research group, located at the University of Almeria. Two kinds of computing nodes have been used. The first node contains two AMD EPYC 7302 (32 CPU cores), 512 GB of DDR4 RAM, and two NVIDIA Tesla V100 (32 GB). The second node contains two Intel Xeon E5-2620v3 (12 CPU cores), 64 GB of DDR3 RAM, and two NVIDIA Kepler K80 (12 GB).

For the experimentation, we have used a Head and Neck IMRT real patient case treated with nine radiation beams. In this case, three PTVs with different prescribing radiation dose deposition levels 66 Gy, 60 Gy and 54 Gy (Gray, symbol Gy, is a unit of radiation dose) are defined. The most important is the PTV with the highest prescribed dose (PTV66) because the highest concentration of tumor cells is there. The other two PTVs are treated prophylactically, so the dose distribution homogeneity in them is less critical. In addition, five OARs are singled out: the spinal cord +3mm, the brainstem +3mm, the left salivary gland, the right salivary gland, and the mandible. Also, the normal tissue is defined as a region of the patient body outside all OARs and all PTVs. The dose deposition model contains 30265 beamlets interacting with 94647 voxels representing the irradiated part of the patient body.

Physical restrictions imposed on RT plans by oncology physicians are converted to the bounds defined in the model presented in Subsection 2.2.1. Table 2 presents the respective bound values for organs considered in the experiments. Those values are repeated also in Table 4 and Table 6. The bound values for PTVs are set by the following rules:

$LB_{t}=0.90\times\textit{prescribed dose for PTV${}_{t}$}$,

$\overline{LB}_{t}=0.98\times\textit{prescribed dose for PTV${}_{t}$}$,

$\overline{UB}_{t}=1.02\times\textit{prescribed dose for PTV${}_{t}$}$,

$UB_{t}=1.10\times\textit{prescribed dose for PTV${}_{t}$},t\in T$.

Explicit constraints on DVH shapes are not a part of model (11) as they are implicitly controlled by gEUD. On the other hand, the DVH shape for PTVs ideally should look like the sequence of three-line segments: horizontal (at 100% level), vertical (at the prescribed deposited dose value), and horizontal again (at 0% level) (cf. Figure 2, Figure 4, Figure 6). DVHs are primary tool to verify ex-post the homogeneity of dose distributions in PTVs. As such, they are always visually inspected by RT planers and oncology physicians for a holistic evaluation. A proxy measure of PTV homogeneity is the dose-volume metric Dx$\%$. For a given PTV, Dx$\%$ is the minimal dose deposited in x$\%$ of the most irradiated PTV voxels.

In our case, metrics D98$\%$ and D2$\%$ are used, namely plans are to satisfy

D98$\%\geq$ 95$\%\times\textit{prescribed dose for PTV${}_{t}$}$, and

D2$\%\leq$ 107$\%\times\textit{prescribed dose for PTV${}_{t}$}$ (i.e., in the latter case: the minimal dose deposited in 2$\%$ of the most irradiated PTV voxels should be less or equal than 107$\%\times EUD^{0}_{t}$).

By this, we have the following levels of Dx$\%$

for PTV66 Gy: D98$\%\geq$ 62.70 Gy, D2$\%\leq$ 70.62 Gy,

for PTV60 Gy: D98$\%\geq$ 57.00 Gy, D2$\%\leq$ 64.20 Gy,

for PTV54 Gy: D98$\%\geq$ 51.30 Gy, D2$\%\leq$ 57.78 Gy.

D2$\%$ and D98$\%$ represent two points on the respective DVHs, thus they give a rough characterization of DVHs.

The planning system PersEUD is run three times, each time with a different aim. The aim of the first run is to derive a number (150 by default) of mutually nondominated plans that compromise constraint violations versus the sum of the average of radiation doses deposited in the left and the right salivary glands (a bi-objective optimization problem). In the second run, the aim is to derive plans that compromise constraint violations versus doses deposited in the spinal cord + 3mm (a bi-objective optimization problem), whereas, in the third run, the aim is to derive plans that compromise constraint violations versus doses deposited in the left and the right salivary gland versus doses deposited in the spinal cord +3mm (a tri-criteria optimization problem).

The derived plans are evaluated by medical physicist experts (below: the Experts) who, in their daily work, make treatment plans in a commercial RT planning system in a series of try-and-correct interactions. Once they are satisfied with the result, they submit the plan to the oncology physician responsible for the case for approval. If the plan is disapproved, the process is repeated. In complex cases, two or three cycles may be needed to secure the final physician approval.

In the whole process, they deal with one RT plan at a time. The three runs of PersEUD produce 450 plans, fully automatic. For proof-of-concept analyses, the Experts selected one plan from those produced in each PersEUD run and pre-selected by the Decision Tool module.

To fulfill the aim of the first PersEUD run, the second objective $f_{1}$ is the sum of the average radiation doses deposited in both glands, $|P|=1$, and $f_{0}$ is defined by formula (8).

From the resulting RT plans, Experts select one that does not violate any constraint ($f_{0}(x)=0$). This plan is denoted as Plan 1 and it is presented in Table 2. Table 3 presents the parameters of gEUD measures and of the function $F(x,\phi)$ for this plan, Figure 2 presents its dose-volume histograms, and Figure 3 the dose distributions for this plan on two exemplary cross sections of the patient irradiated part (head and neck).

As seen in Table 2, this plan provides adequate dose distributions in PTVs, and, as intended, preserves both salivary glands, and it does this surprisingly well. The average doses are 14.26 Gy for the right and 12.56 Gy for the left salivary gland. These values are well below the imposed bound of 26 Gy. All other OARs with no protection priority are also below their acceptable maximal doses. The most notable is the maximum dose in the spinal cord +3mm, equal to 41.43 Gy, while the bound is 50 Gy, and the maximum dose in the brainstem +3mm, equal to 28.67 Gy, while the bound is 60.00 Gy.

The actual average dose deposited in PTV66 (the most important PTV) is 66 Gy, as prescribed, and D98$\%\geq$ 62.70 Gy and D2$\%\leq 70.62Gy$.

The actual average dose deposited in PTV60 is 60.04 Gy and that is almost the perfect match as this value is within $\pm$2$\%$ range from the target value. However, D98$\%$ is significantly unmet (52.06 Gy vs 57.00 Gy expected) due to the protection of the salivary gland. D2$\%$ is also exceeded but by a clinically insignificant value.

The actual average dose deposited in PTV54 is 53.47 Gy and that is also almost the perfect match as this value is within $\pm$2$\%$ range from the target value. Moreover, D2$\%\leq 57.78$Gy, but D98$\%\not\geq 51.30$.

As mentioned earlier, in regions treated prophylactically (in the considered case PTV 60 and PTV54) slight violations of D98$\%$ and D2$\%$ related constraints are acceptable.

To fulfill the aim of the second PersEUD run, the second objective $f_{1}$ is the average radiation dose deposited in the spinal cord +3mm, $|P|=1$, and $f_{0}$ is defined by formula (8).

From the resulting RT plans, the Experts select one that does not violate any constraint ($f_{0}(x)=0$). This plan is denoted as Plan 2 and is presented in Table 4. Table 5 presents parameters of gEUD and of function $F(x,\phi)$ for this plan, Figure 4 presents its dose-volume histograms, and Figure 5 the dose distributions for this plan on two exemplary cross sections of the patient irradiated part (head and neck).

In Plan 2 the maximum dose deposited in the spinal cord +3mm is 11.02 Gy, noticeably lower than in the Plan 1. Although protection of the salivary glands in this plan is not a priority, the average dose is 22.19 Gy and 22.08 Gy in the left and right salivary gland, respectively, which is far below the dose tolerance for these organs. This is at the price of a slight violation of the maximal dose allowed in the mandible (70 Gy) where the actual dose deposited is 71.13 Gy.

The actual average dose deposited in PTV66 is 66 Gy, as prescribed. However, D98$\%\not\geq$ 62.70 Gy (violation by 1.83 Gy), and D2$\%\not\leq$ 70.62 Gy (violation by 1.04 Gy).

The actual average dose deposited in PTV60 is 59.68 Gy and that is almost the perfect match as this value is within $\pm$2$\%$ range from the target value. Moreover, D98$\%\geq$ 57.00 Gy and D2$\%\leq$ 64.20 Gy.

The actual average dose deposited in PTV54 is 54.30 Gy and that is also almost the perfect match as this value is within $\pm$2$\%$ range from the target value. However, D98$\%\not\geq$ 51.30 Gy, but D2$\%\leq$ 57.78 Gy.

As mentioned earlier, in regions treated prophylactically (i.e. PTV 60 and PTV54) slight violations of D98$\%$ and D2$\%$ related constraints are acceptable.

To fulfill the aim of the third PersEUD run, the second objective function $f_{1}$ is the average of the radiation doses deposited in the left and the right salivary glands, the third objective function the is the maximal radiation dose deposited in the spinal cord +3mm, $|P|=2$, and $f_{0}$ is defined by formula (8).

From the resulting RT plans, the Experts select one that does not violate any constraint ($f(x)=0$). This plan is denoted as Plan 3 and is presented in Table 6. Table 7 presents parameters of gEUD and of function $F(x,\phi)$ for this plan, Figure 6 presents its dose-volume histograms, and Figure 7 the dose distributions for this plan on two exemplary cross sections of the 3D model of the patient irradiated part (head and neck).

In Plan 3, the average dose deposited in the right and the left salivary gland is 13.94 Gy and 13.14 Gy, respectively. That is much lower than the maximal admissible value of 26 Gy, while the maximal dose deposited in the spinal cord +3mm is 24.06 Gy (the maximal acceptable value is 50 Gy). Although Plan 1 and Plan 2 offer slightly lower doses in respective priority protected OARs (the salivary glands in the first case, and the spinal cord +3mm in the second case), Plan 3 well balances the doses delivered to both salivary glands and to the spinal cord +3mm.

The actual average dose deposited in PTV66 is 66 Gy, as prescribed, and D98$\%\geq$ 62.70 Gy, and D2$\%\leq$ 70.62 Gy. The actual average dose deposited in PTV60 is 60.2768 Gy, within $\pm$2$\%$ range from the target value. Moreover, D98$\%\not\geq$ 57.00 Gy and D2$\%\leq$ 64.20 Gy. The actual average dose deposited in PTV54 is 55.34 Gy and that is also almost the perfect match as this value is within $\pm$2.5$\%$ range from the target value. However, D98$\%\not\geq$ 51.30 Gy, but D2$\%\leq$ 57.78 Gy.

As mentioned earlier, in regions treated prophylactically (i.e. PTV 60 and PTV54) slight violations of D98$\%$ and D2$\%$ related constraints are acceptable.

Table 8 provides a comparison of the three plans. At a glance, the Experts disapproved Plan 2 due to its underperformance in PTV66, the most important PTV in the considered case. The explanation for the disapproval is that such an underperformance cannot be outweighed by the excellent performance in the spinal cord +3mm and the brainstem +3mm. From the remaining two, the Experts prefer Plan 3 because it offers low radiation doses deposited in both salivary glands, the spinal cord +3mm, and brainstem +3mm.

To verify the clinical viability of Plan 3 (below: P3), we compare it with two plans prepared by the Expert in Eclipse^TM Treatment Planning system (below: Eclipse) from Varian Medical Systems, namely:

• R1

  Starting plan: an RT plan generated automatically by Eclipse at the start of the planning process.

• R2

  Blind-expert plan: an RT plan prepared by the Expert without seeing P3.

All plans have to satisfy the same constraints (specified in Table 6). The comparisons are by DVH and, for illustration, by one representative cross section of the patient model.

The first comparison (Figure 8) is between P3 and R1. As displayed, P3 reduces the doses delivered to all OARs, with comparable or better dose distributions in PTVs.

When comparing P3 with R2 (Figure 9), R2 provides a better dose distribution in PTVs, however at the expense of higher dose deposited in OARs, especially in both salivary glands and the spinal cord +3mm. Notably, in P3 the doses deposited in the salivary glands are 20$\%$ lower than in R2. And in P3 the dose deposited in PTV66 is within the specified bounds.

Next, we compare with Guided-expert plan (R3), a plan prepared by the Expert using P3 as the Eclipse starting plan, with R2 (Figure 10). R3 provides a similar to R2 dose distribution in PTVs and a significant reduction in the doses deposited in OARs. This shows that plans produced by the PersEUD system can be used to enhance RT plans produced by commercial systems currently in use.

As preliminary conclusions from our work, it can be assumed that using plans from PersEUD as starting points in RT planning systems routinely used in an oncology hospital should improve RT plan quality and reduce RT plan preparation times. This claim has to be validated in more clinical cases.