Adaptive IMRT using a multiobjective evolutionary algorithm integrated with a diffusion–invasion model of glioblastoma

In the LQ model (Hall 1994), the fraction of cells, S, that survive absorbed dose d (Gy), delivered to a region of tissue in n fractions, is given as S(nd) = exp(–αnd –βGn²d²), where α (Gy⁻¹) and β (Gy⁻²) are intra- and inter-track radiation sensitivity parameters and G is the Lea–Catcheside dose-protraction factor (Sachs et al 1997). For a typical radiation therapy treatment consisting of n well separated fractions, G can be approximated by G ≅ 1/n (Stewart and Li 2007), and the LQ survival model can be re-written as S(nd) = exp(–αnd – βnd²). In this study, we used α = 0.12 Gy⁻¹ and α/β = 3 Gy for normal brain tissue.

Equivalent uniform dose (EUD) is a method for quantifying non-uniform dose distributions. EUD is defined as the absorbed dose, which if uniformly delivered to a tissue region of interest, will produce the same biological effect as delivery of the non-uniform dose distribution to the same tissue region of interest (Niemierko 1997). To compute this we solve iteratively for EUD as shown in equation (1):

$$\begin{equation} \sum\limits_i {c_i v_i \,{\rm e}^{ - \alpha d_i - \beta _i \,d_i^2 } } = \left[ {\sum\limits_i {c_i v_i } } \right]{\rm e}^{ - \alpha {\rm EUD} - \beta {\rm EUD}^2 } ,\end{equation} \tag{ 1 }$$

where c_i, v_i and d_i represent the cell density, volume and dose, respectively, at the spatial location indicated by i.

A patient-specific, proliferation-invasion (PI-RT) mathematical model of GBM (Rockne et al 2009, Rockne et al 2010) describes the change in tumor cell density with respect to time and space by quantifying net rates of proliferation (ρ time⁻¹) and invasion (D area/time). The radiosensitivity parameters α and β in the PI-RT model are the same as in LQ cell survival model. Patient-specific PI-RT parameters, except normal tissue radiosensitivity parameters, were determined from routinely available pre-treatment MRI imaging (Swanson et al 2003, 2008a, Harpold et al 2007, Wang et al 2009b, Szeto et al 2009, Rockne et al 2010). Tumor-specific estimates of α have been shown to correlate linearly with the proliferation parameter ρ and can be determined for any patient a priori to radiation treatment, assuming the ratio α/β is assumed to be fixed at 10 Gy (Rockne et al 2010).

Tumor cell density evolves in both space and time and represents net density per unit volume, not the behavior of individual cells. Proliferation is assumed to follow a logistic growth curve while invasion is modeled as a Fickian diffusion process. The model is of the form:

$$\begin{equation} \frac{{\partial c}}{{\partial t}} = \nabla \cdot ( {D\nabla c} ) + \rho c\left( {1 - \frac{c}{k}} \right ) - Rc\left( {1 - \frac{c}{k}} \right ),\end{equation} \tag{ 2 }$$

where ρ is the rate of proliferation, D is the rate of diffusion, k is the maximum allowed tumor cell concentration, c is the tumor cell distribution as a function of space and time:

$$\begin{equation} c = c({\bf x},{\rm }t) \end{equation} \tag{ 3 }$$

and

$$\begin{equation} R = \left\{ \begin{array}{@{}ll@{}} 0,&{\rm between \; daily\; fractions} \\ 1 - S(\alpha ,\beta ,d({\bf x},t)),&{\rm otherwise} \\ \end{array} \right.\end{equation} \tag{ 4 }$$

where d(x, t) is the radiation dose distribution and S is the estimated surviving fraction using the LQ cell survival model with the parameters α and β.

The GBM model uses a 'tip of the iceberg' view of clinical imaging, where the boundaries of the enhancing regions on T1Gd and T2 MRIs are associated with isosurfaces of constant density. These isosurfaces are used to infer a gradient of tumor cells, characterized by an invisibility index, the ratio of model parameters D/ρ (Rockne et al 2010, Szeto et al 2009, Wang et al 2009b, Harpold et al 2007). The GBM model is a form of the well-studied Fisher's equation, with solutions that approach a traveling wave with a constant velocity of $2\sqrt{D}$ ρ (Fisher 1937). This is in agreement with observed linear velocities on serial MRI of untreated low-grade gliomas followed over as many as 15 years (Mandonnet et al 2003) as well as untreated GBM (Swanson et al 2002). With these two equations, we can solve for the patient-specific model parameters D and ρ algebraically as in (Rockne et al 2010, Szeto et al 2009, Wang et al 2009b, Harpold et al 2007).

Results presented in the work are based on data obtained from two patients. The GBM model parameters used for patients were derived from MRI images taken before and after standard radiation therapy. The values used for patient 1 with relatively radio-resistant glioblastoma and for patient 2 with relatively radio-sensitive disease are given in table 1 (patients 6 and 5 in Rockne et al 2010, respectively). Results were calculated for the standard of care plan that was generated for clinical use that delivered 61.2 Gy in 1.8 Gy daily fractions to the CTV over 7 weeks. In order to better quantify improvement from radiation therapy, the GBM model was used to predict what the tumor size (T2 MRI visible radius) would be for patient 1 and patient 2 over time if no radiation therapy was used. Quantitative results of tumor control are reported in terms of 'Treatment Gain' which is the difference of the time it takes for this baseline to reach a 4 cm T2 MRI visible radius versus the time it would take using various forms of radiation treatment as calculated by the model. Normal tissue EUD is used to quantify risk to normal tissue for these simplified example patients.

A MOEA was used for IMRT optimization (Holdsworth et al 2010, Holdsworth et al 2011, Holdsworth et al 2012). The MOEA itself consists of two interacting algorithms each with its own set of objectives: (1) a deterministic IMRT optimization algorithm that generates IMRT plans by minimizing a penalty function that is a weighted sum of clinical objectives and (2) an evolutionary algorithm that selects superior (non-dominated) plans based on a separate set of decision criteria. The population of superior plans is used to generate new formulations of the penalty function which can be minimized to produce new IMRT plans. Using this process of natural selection, the MOEA will produce sets of IMRT plans that approach the Pareto Front and represent the best possible tradeoff amongst decision criteria.

Deterministic multiobjective IMRT algorithms minimize fitness criteria, and then use those same criteria to judge plan quality in the selection process; however, this requires that all fitness criteria are of a mathematical form that lends itself to efficient optimization. The MOEA is a stochastic multiobjective algorithm that uses a set of penalty objectives for individual IMRT optimization, and uses a separate set of decision criteria to judge plan quality in the selection process of the evolutionary algorithm. While the penalty objectives in our algorithm are quadratic, the decision criteria have no mathematical restrictions and can be designed solely to emulate to most appropriate clinical goals. Another advantage is that the set of decision criteria can be kept small for more efficient multiobjective optimization (Holdsworth et al 2011), while the set of penalty objectives can be expanded for greater flexibility in the search for more optimal IMRT plans (Holdsworth et al 2012).

The IMRT optimization utilized in our MOEA minimizes a fairly standard penalty function that is a weighted sum of quadratic penalty objectives:

$$\begin{equation} \sum \sum H_i ({D_{{\rm ref}, i}} - {d_{i,j}} )\;W_i ({D_{{\rm ref},i}} - {d_{{i},{j}}})^2 + F_s \end{equation} \tag{ 5 }$$

where d_i, j is dose in the jth voxel of the ith structure, W_i is the weight of the ith structure, D_ref,i is the reference dose parameter that defines the ith structure's penalty objective, F_s is a smoothing matrix that ensures the fluence can be translated into a deliverable beam, and H_i (D_{ref, i} – d_i, j) is a Heaviside function that sets the weight to 0 in OAR voxels when D_{ref, i} > d_i, j and is 1 otherwise. Each iteration, the evolutionary algorithm generates sets of W_i and D_ref, based on recombination and mutation from existing non-dominated plans that formulate a new penalty function. This penalty function is minimized generating a new IMRT plan. For this work, the quality of the plan is judged by three decision criteria described below, and all plans are compared and any dominated plans are eliminated prior to beginning the next iteration of the MOEA.

This stochastic algorithm has the advantages of being able to incorporate any form of decision criteria and of not requiring a priori determination of the importance of each criterion. The algorithm returns Pareto optimal plans with respect to the decision criteria utilized in the multiobjective optimization. In this work, three decision criteria were used. Two criteria were used to evaluate all plans: (1) minimizing the maximum dose per fraction to any voxel and (2) minimizing normal tissue EUD. The third decision criterion was varied in order to see the effects of altering target optimization objectives: (3a) tumor cell kill, (3b) tumor survival after 1 week of treatment or (3c) tumor survival 11 weeks after 1 week of treatment. For criteria (3b) and (3c), the GBM model was used to predict tumor growth from that stage of the optimization process. These criteria were used in the weekly adaptive optimization process. Each fraction was scaled such that the EUD to normal brain was less than or equal to that of the standard-of-care plan and such that the maximum dose per fraction to any voxel in the scaled dose distribution was also limited by a given maximum dose (3 Gy, 8 Gy, or no limit).

In this work, the optimal distribution for each week of treatment is generated in 300 iterations, and 7 weeks of treatment are generated for each scenario. Optimization of all of the distributions throughout the 34 fractions to be delivered to one patient takes anywhere from 6 min to 10 h. When maximum cell kill is used as target decision criteria and there is a low allowed maximum dose, the adaptive optimization is much faster. If decision criteria for the target is calculated using the model to simulate 12 weeks of response for each iteration and if no maximum dose is used, the optimization takes much longer. The optimal distribution to be used in the GBM model for each week of treatment is selected based on a weighted sum of decision criteria. The GBM model then simulates of growth, diffusion and radiation damage for five days using this radiation dose distribution and two weekend days without radiation. The resulting tumor cell distribution is used for the following week's optimization.

For ease of model integration and computation, simulations were carried out in a simplified three-dimensional geometry where the tumor cell density was spherically symmetric and the structures in the optimization consisted of concentric shells. Each shell had its own proportion of tumor cells and genes (objectives and weights used in the penalty function) to be optimized by the MOEA. In this work, the simulations included a few sets of simple decision criteria, but other objectives such as dose to specific structures, DVH information or any quantifiable metrics specific to a given clinical situation can be used to either judge or scale dose distributions.

Adaptive therapy was implemented by performing weekly MOEA IMRT optimizations using the most up-to-date predictions of tumor cell distributions obtained from the GBM model. Each week of treatment is simulated with the GBM model using the generated dose distribution to determine the tumor cell density used in the next round of optimization. Figure 1 demonstrates the iterative dialogue between the two algorithms.

Zoom In Zoom Out Reset image size

The results of the optimizations were compared using four metrics:

• (1)  
  tumor radius as would be visible on T2 MRI;
• (2)  
  dose–volume histograms of total normal tissue EUD outside of the T2 enhancing region;
• (3)  
  dose–volume histograms of total tumor tissue EUD;
• (4)  
  treatment gain, the delay in days for the T2 radius to reach 4 cm compared to the untreated control.

In each simulation, three decision criteria were used to judge plan quality. Two of these criteria are used in every simulation while the third criterion was varied to test different models of tumor control. Results are reported in terms of the T2 MRI radius over time as predicted by the GBM model, treatment gain, and normal tissue EUD.

The predicted radius of the T2 MRI region over time is plotted for the standard plan and for the results of five different MOEA optimizations in figure 2. Three of the optimized plans were generated using minimum tumor cell survival after 1 week of treatment as the tumor control decision criteria with three different maximum dose/fraction restrictions: 3 Gy, 8 Gy, and no restriction. The other two MOEA optimizations shown use an 8 Gy/fraction maximum dose restriction using either maximum tumor cell kill per fraction or minimum tumor survival 12 weeks later for tumor control decision criteria. Similar results are shown in figure 3 for the second example patient.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A summary of results are in table 2 below. Normal tissue EUD is calculated by solving equation (1) iteratively. Treatment gain is the added number of days for the T2 MRI radii to reach 4 cm using RT versus no treatment.

Note that for some simulations, the radius plot shows a steep jump in radius before resuming linear growth. This does not represent change in overall growth rate, but instead occurs when a large region of the tumor cell population achieves the minimum T2 cell density still visible on an MRI image. This reflects the complex interplay between a diffuse distribution of gliomas cells and the ability to reach a local threshold of detection to be visible on MRI.

For both patients, optimized plans with the maximum dose per fraction limited to 3 Gy resulted in response slightly worse than the standard plan for one-fourth the dose to normal tissue. For patient 1 that had lower parameter values for the net proliferation rate, ρ, and radiosensitivity, the treatment gain was lower for all RT plans, and results were similar for all optimized plans with maximum dose per fraction > 3 Gy. For patient 2 that had larger values for the net proliferation rate, ρ and radiosensitivity, tumor response was highly dependent on both the maximum dose and the tumor control decision criterion implemented. When no maximum dose was used, the treatment gain improved to 157 days; however, the high dose in small volumes of tissue would be prohibitive for practical use in the clinical. The MOEA optimized treatment scheme that used tumor survival 12 weeks post treatment demonstrated the best simulated performance for this patient with a treatment gain of 153 days. All results were simulated by the GBM model for two simplified example patients.

Cumulative dose volume histograms (DVHs) are used to demonstrate how different tumor decision criteria and maximum dose restrictions affect dose distributions. In figure 4, DVHs using total dose delivered in 34 fractions to the tumor cell population for the standard plan compared to the five same MOEA optimized dose distributions as in figures 2 and 3 for the second patient are displayed.

Zoom In Zoom Out Reset image size

The maximum dose restriction directly affects the dose delivered to the bulk of the tumor; however, the predicted tumor control of the plan that used 11 weeks after 1 week of treatment cell survival as decision criteria had a much lower dose level and performance that was comparable to the MOEA plan that did not have a maximum dose restriction according to the GBM model. This result implies that the location of tumor cells targeted for radiation therapy is as important as dose level to the tumor core when it comes to predicted long term tumor control. There is wide patient-to-patient variation in the extent of diffusely invaded tumor cells peripheral to the imaging abnormality (Harpold et al 2007, Wang et al 2009a, Szeto et al 2009, Rockne et al 2010) that supports the need for patient-individualized treatment optimization to carefully incorporate this diffuse disease extent.

Without setting a maximum dose, the optimization can reach doses over small volumes that are far higher than current clinical practice. By forcing the maximum dose of optimizations to be lower, distributions will be more acceptable for clinical use, but tumor control may be sacrificed. As previously discussed (Holdsworth et al 2011), the specific goals or decision criteria used during multiobjective optimization play an important role in the final dose distributions and the quality of the radiation therapy with regard to what is best for the patient.

In figure 5, DVHs of the normal cell population for the optimized dose distributions for the second patient are compared to the DVHs resulting from the conventional clinical radiation therapy. Minimum tumor cell survival after 1 week of treatment is used as decision criteria in the MOEA for three optimizations using different restrictions on maximum dose. DVHs for normal tissue are also displayed for two other MOEA optimizations using maximum tumor cell kill and minimum tumor cell survival 12 weeks after radiation therapy that are also displayed (both with an 8 Gy maximum dose restriction).

Zoom In Zoom Out Reset image size