Sage Journals: Discover world-class research

Abstract

Background

Spatially fractionated radiotherapy (SFRT) intentionally creates spatially-modulated peak-valley dose patterns for improving the tumor control or/and the sparing of organs at risk, compared to conventional radiotherapy (CONV). Evaluating non-uniform dose distributions using conventional physical dose does not fully account for the biological effectiveness of SFRT. In this study, we adopt the equivalent uniform dose (EUD) as a surrogate metric to evaluate non-uniform dose distributions.

Purpose

This work will develop a SFRT treatment planning method with EUD optimization, which is the first-of-its-kind to the best of our knowledge.

Methods

The SFRT scenario of proton GRID with uniform target dose is considered in this work for dose-only optimization (DO) and joint dose and peak-valley-dose-ratio (PVDR) optimization (JDPO) respectively. In addition to dose and PVDR optimization, SFRT treatment planning also optimizes the EUD. The EUD used in this work is based on cell survival as modeled by the linear-quadratic model. The EUD optimization problem is solved by: (1) iterative convex relaxation to decouple the nonconvex dose-volume-histogram constraint; (2) linearized alternating direction method of multiplier to efficiently handle the nonconvex minimum-monitor-unit constraint and nonlinear EUD objective.

Results

EUD optimization reduced the EUD and increased survival fraction. For example, for a head-and-neck patient, EUD optimization decreased the brainstem EUD from 5.52% to 3.86% for DO, and from 5.78% to 4.82% for JDPO, and increased the brainstem survival fraction from 72.9% to 81.6% for DO, and from 71.5% to 76.6% for JDPO. Furthermore, EUD optimization preserved PVDR from DO or JDPO.

Conclusions

A novel EUD optimization method is proposed for SFRT that can reduce the EUD and increase the survival fraction.

Keywords

spatially fractionated radiotherapy (SFRT)equivalent uniform dose treatment planning inverse optimization

Introduction

Compared to conventional radiotherapy (CONV), spatially fractionated radiation therapy (SFRT) intentionally delivers highly spatially modulated dose to the tumor target or/and the organs-at-risk (OAR), characterized by the peak-valley dose ratio (PVDR), for improving the therapeutic ratio.¹ Evaluating SFRT plan quality using physical dose, such as mean, peak or valley dose, may not fully account for its biological effectiveness.

Equivalent uniform dose (EUD) is a metric that quantifies the uniform dose distribution required to produce the same biological effect for a given nonuniform dose distribution. For example, it is shown in ref.^2,3 that the EUD is a useful quantity to design corrections for various beam parameters and analyze the effects of relative biological effectiveness for different particle types and energies. In ref.,⁴ it is demonstrated that the EUD is a superior predictor of cell survival compared to the physical dose. This motivates the direct optimization of EUD during non-uniform dose distribution, which will be the primary objective of this work.

While the EUD optimization method developed in this work can be applicable to various SFRT scenarios, this work considers a specific SFRT scenario of proton GRID with uniform target dose (UTD)⁵ to demonstrate the efficacy of EUD optimization for OAR sparing. The reasons to consider UTD are: (1) UTD is consistent with CONV and favored by preclinical evidences that the valley dose correlates the best with tumor control;⁶ (2) the purpose of EUD optimization in this work is for OAR sparing, for which UTD conveniently allows the standard plan normalization to targets, so that the OAR sparing can be evaluated under the same plan normalization.

Methods

SFRT Treatment Planning with EUD Optimization

The SFRT treatment planning considered in this work is for UTD proton GRID with dose-only optimization (DO) and joint dose and PVDR optimization (JDPO) respectively. In addition to the optimization of dose and PVDR objectives, EUD is also optimized. The general form of SFRT treatment planning problem with EUD optimization is

min_{x \in R^{N_{x}}} f (d, Ω) + c_{TV} \sum_{a = 1}^{N_{a}} g (d_{a}) + c_{EUD} \sum_{i = 1}^{N_{oar}} E U D (d_{O A R_{i}})

s . t . {\begin{matrix} \begin{matrix} d = D x, \\ d_{O A R_{i}} = D_{O A R_{i}} x, i = 1, 2, \dots, N_{oar}, \end{matrix} \\ d_{a} = I_{a} d, a = 1, 2, \dots, N_{a}, \\ x_{j} \in {0} \cup [G_{m i n}, + \infty), j \leq N_{x} \end{matrix}

(1)

In Equation (1), f is the dose objective that enforces the consistency between optimized dose d and prescribed dose values for the voxels, within the region Ω, that violate the prescribed dose values (from dose-volume-histogram (DVH) based constraints^7,8). D is the dose influence matrix. D_OARi is the influence matrix for the i^th OAR. x∈ℝ^Nx is the beam intensity, and it needs to meet the minimum-monitor-unit (MMU) constraint^9–14 with the MMU threshold G_min.

N_a is the number of beam angles and d_a is the beam-eye-view (BEV) 2D planar dose for the a^th angle that is interpolated from d using the linear interpolation operator I_a. g is the regularization function to maximize the PVDR^5,15 for d_a, ie,

g (d_{a}) = ‖ d_{a, v} ‖_{1} - c_{p} ‖ d_{a, p} ‖_{1} + c_{T} ‖ W d_{a} ‖_{1}

(2)where the ||·||₁ is L₁ norm, defined as the sums of the absolute values. d_a,v and d_a,p are the valley and peak dose of d_a, defined by

d_{a, v} = R_{v} d_{a}

d_{a, p} = R_{p} d_{a}

(3)

R_v and R_p are the restriction operator that maps the dose vector d to the peak and valley regions. W is a discrete total-variation operator that map a 2D dose d_a to a vector field, consisting of vertical and horizontal finite difference, ie,

(W d_{a})_{m, n} = {\begin{matrix} {(d_{a})}_{m + 1, n} - {(d_{a})}_{m, n} \\ {(d_{a})}_{m, n + 1} - {(d_{a})}_{m, n} \end{matrix}

(4)where (·)_m,n denotes the entry of a given matrix with row index m and column index n. c_p and c_T are the weights for the peak and total-variation term. c_TV is the objective weight for g.

The third term EUD(d_OAR,i) is the EUD for the i^th OAR. N_oar is the number of OAR. EUD is used to quantify the uniform dose distribution that would have the same biological effect as a given nonuniform dose distribution. Adopting the linear-quadratic model, the survival fraction (SF) for a given OAR is defined as

S F_{O A R} = \sum_{i = 1}^{N V_{O A R}} ν_{i} e x p (- α d_{O A R, i} - β d_{O A R, i}^{2})

(5)where v_i,

\sum_{i = 1}^{N V_{O A R}} v_{i} = 1

is the fractional volume of dose. α and β are the radio sensitivity parameters. d_OARi is the dose value for the i^th voxel and NV_OAR is the number of the OAR. Then, EUD (with c_EUD be the objective weight) for the OAR is the following¹⁶

E U D (d_{O A R}) = \frac{- α + \sqrt{α^{2} - 4 β \ln (S F_{O A R})}}{2 β}

(6)

Optimization Algorithm

The optimization challenges of Equation (1) come from the nonconvexity of DVH and MMU constraints, the non-smoothness of function g, and the non-linearity of the EUD (Equation (6)). The optimization algorithm for solving Equation (1) consists of (1) the iterative convex relaxation technique^17–20 to decouple the DVH constraint, (2) and then the linearized alternating direction method of multipliers (ADMM)²¹ to handle MMU constraint, g and EUD.

We use quadratic functions for f and rewrite Equation (1) as

\begin{aligned} min_{x \in R^{N_{x}}} \frac{1}{2} ‖ A (Ω) - b (Ω) ‖_{2}^{2} + c_{TV} \sum_{a = 1}^{N_{a}} g (I_{a} D x) \\ + c_{EUD} \sum_{i = 1}^{N_{oar}} E U D (D_{O A R_{i}} x) \\ s . t . x_{j} \in {0} \cup [G_{m i n}, + \infty), j \leq N_{x} \end{aligned}

(7)

The active index Ω depends on a given dose distribution and a specific DVH constraint, and thus depends on the x (see Section S1 in the Supplemental material). For k^th iterate x^k, we update the active index Ω^k ⁺ ¹ using the operator H

Ω^{k + 1} = H (D x^{k})

(8)

With Ω^k ⁺ ¹ fixed, we update the beam intensity x by solving the following problem

\begin{aligned} min_{x \in R^{N_{x}}} \frac{1}{2} ‖ A (Ω^{k + 1}) - b (Ω^{k + 1}) ‖_{2}^{2} + c_{TV} \sum_{a = 1}^{N_{a}} g (I_{a} D x) \\ + c_{EUD} \sum_{i = 1}^{N_{oar}} E U D (D_{O A R_{i}} x) \\ s . t . x_{j} \in {0} \cup [G_{m i n}, + \infty), j \leq N_{x} \end{aligned}

(9)

Using Equation (2) and (3), the function g can be rewritten as

g (I_{a} D x) = ‖ R_{v} I_{a} D x ‖_{1} - c_{p} ‖ R_{p} I_{a} D x ‖_{1} + c_{T} ‖ W I_{a} D x ‖_{1}

(10)

The optimization challenges for Equation (9) come from the L1 objectives with linear operator R_v, R_p, I_a and the MMU constraint. To mitigate the difficulty, we transfer the linear operators and the MMU constraint to their auxiliary variables, ie,

\begin{aligned} min_{x \in R^{N_{x}}} \frac{1}{2} ‖ A (Ω^{k + 1}) - b (Ω^{k + 1}) ‖_{2}^{2} + c_{TV} ‖ z_{v, a} ‖_{1} - c_{p} ‖ z_{p, a} ‖_{1} \\ + c_{T} ‖ z_{T, a} ‖_{1} + c_{EUD} \\ \sum_{i = 1}^{N_{oar}} E U D (D_{O A R_{i}} x) \\ s . t . {\begin{matrix} R_{v} I_{a} D x = z_{v, a}, & a = 1, 2, \dots, N_{a} \\ R_{p} I_{a} D x = z_{p, a}, & a = 1, 2, \dots, N_{a} \\ W I_{a} D x = z_{T, a}, & a = 1, 2, \dots, N_{a} \\ x = z \\ z_{j} \in {0} \cup [G_{m i n}, + \infty), & j \leq N_{x} \end{matrix} \end{aligned}

(11)

Due to the nonsmoothness of objective functions and multiple equality constraints, we use first-order optimization methods to solve Equation (11), eg, ADMM^21–27 in this work.

Firstly, we define the augmented Lagrangian for Equation (11)

\begin{aligned} L (x, z, u, z_{v, 1}, \dots, z_{v, N_{a}}, u_{v, 1}, \dots, u_{v, N_{a}}, z_{p, 1}, \dots, \\ z_{p, N_{a}}, u_{v, 1}, \dots, u_{v, N_{a}}, z_{T, 1}, \dots, z_{T, N_{a}}, u_{T, 1}, \dots, u_{T, N_{a}}) \\ = & \frac{1}{2} ‖ A (Ω^{k + 1}) - b (Ω^{k + 1}) ‖_{2}^{2} + c_{TV} \sum_{i = 1}^{N_{oar}} ‖ z_{v, a} ‖_{1} - c_{p} ‖ z_{p, a} ‖_{1} \\ + c_{T} ‖ z_{T, a} ‖_{1} + c_{EUD} \sum_{i = 1}^{N_{oar}} E U D (D_{O A R_{i}} x) + \frac{ρ_{1}}{2} ‖ x - z \\ + u ‖_{1} + \sum_{a = 1}^{N_{a}} \frac{ρ_{2, a}}{2} ‖ R_{v} I_{a} D x - z_{v, a} + u_{v, a} ‖_{2}^{2} \\ + \sum_{a = 1}^{N_{a}} \frac{ρ_{3, a}}{2} ‖ R_{p} I_{a} D x - z_{p, a} + u_{p, a} ‖_{2}^{2} \\ + \sum_{a = 1}^{N_{a}} \frac{ρ_{4, a}}{2} ‖ W I_{a} D x - z_{T, a} + u_{T, a} ‖_{2}^{2} \end{aligned}

(12)where u_v,i, u_v,i and u_v,i are multipliers (dual variables) for the equality constraints.

Then the ADMM algorithm solves Equation (11) by minimizing the augmented Lagrangian Equation (12) with respect to primal and dual variables alternatively as the follows

x^{k + 1} = a r g m i n_{x} L (x, z^{k}, u^{k}, z_{v, 1}^{k}, \dots, z_{v, N_{a}}^{k}, \dots, u_{T, 1}^{k}, \dots, u_{T, N_{a}}^{k})

(13)

z^{k + 1} = a r g m i n_{z} L (z^{k + 1}, z, u^{k}, z_{v, 1}^{k}, \dots, z_{v, N_{a}}^{k}, \dots, u_{T, 1}^{k}, \dots, u_{T, N_{a}}^{k})

(14)

u^{k + 1} = u^{k} + x^{k + 1} - z^{k + 1}

(15)

\begin{aligned} z_{v, a}^{k + 1} = a r g m i n_{z_{v, a}} L (z^{k + 1}, z^{k + 1}, u^{k + 1}, z_{v, 1}^{k}, \dots, z_{v, a}, \dots, z_{v, N_{a}}^{k}, \dots, \\ u_{T, 1}^{k}, \dots, u_{T, N_{a}}^{k}), \end{aligned}

a = 1, 2, \dots, N_{a}

(16)

u_{v, a}^{k + 1} = u_{v, a}^{k} + R_{v} I_{a} D x^{k + 1} - z_{v, a}^{k + 1}, a = 1, 2, \dots, N_{a}

(17)

z_{p, a}^{k + 1} = a r g m i n_{z_{p, a}} L (z^{k + 1}, z^{k + 1}, u^{k + 1}, \dots, z_{p, a}, \dots, u_{T, 1}^{k}, \dots, u_{T, N_{a}}^{k}), a = 1, 2, \dots, N_{a}

(18)

u_{p, a}^{k + 1} = u_{p, a}^{k} + R_{v} I_{a} D x^{k + 1} - z_{p, a}^{k + 1}, a = 1, 2, \dots, N_{a}

(19)

\begin{aligned} z_{T, a}^{k + 1} = & a r g m i n_{z_{T, a}} L (z^{k + 1}, z^{k + 1}, u^{k + 1}, \dots, z_{T, a}, \dots, \\ u_{T, 1}^{k}, \dots, u_{T, N_{a}}^{k}), a = 1, 2, \dots, N_{a} \end{aligned}

(20)

u_{T, a}^{k + 1} = u_{T, a}^{k} + W I_{a} D x^{k + 1} - z_{T, a}^{k + 1}, a = 1, 2, \dots, N_{a}

(21)

For the first subproblem,

\begin{aligned} x^{k + 1} = & a r g m i n_{x} L (x, z^{k}, u^{k}, z_{v, 1}^{k}, \dots, z_{v, N_{a}}^{k}, \dots, u_{T, 1}^{k}, \dots, u_{T, N_{a}}^{k}) \\ = & a r g m i n_{x} c_{EUD} \sum_{i = 1}^{N_{oar}} E U D (D_{O A R_{i}} x) \\ + \frac{1}{2} ‖ A (Ω^{k + 1}) x - b (Ω^{k + 1}) ‖_{2}^{2} + \frac{ρ_{1}}{2} ‖ x - z^{k} + u_{1}^{k} ‖_{2}^{2} \\ + \sum_{a = 1}^{N_{a}} \frac{ρ_{2, a}}{2} ‖ R_{v} I_{a} D x - z_{v, a}^{k} + u_{v, a} ‖_{2}^{2} \\ + \sum_{a = 1}^{N_{a}} \frac{ρ_{3, a}}{2} ‖ R_{p} I_{a} D x - z_{p, a}^{k} + u_{p, a}^{k} ‖_{2}^{2} \\ + \sum_{a = 1}^{N_{a}} \frac{ρ_{4, a}}{2} ‖ W I_{a} D x - z_{T, a}^{k} + u_{T, a}^{k} ‖_{2}^{2} \end{aligned}

(22)

Let

f_{1} (x) = c_{EUD} \sum_{i = 1}^{N_{oar}} E U D (D_{O A R_{i}} x)

(23)

And

\begin{aligned} f_{2} (x) = & \frac{1}{2} ‖ A (Ω^{k + 1}) x - b (Ω^{k + 1}) ‖_{2}^{2} + \frac{ρ_{1}}{2} ‖ x - z^{k} + u^{k} ‖_{1} \\ + \sum_{a = 1}^{N_{a}} \frac{ρ_{2, a}}{2} ‖ R_{v} I_{a} D x - z_{v, a}^{k} + u_{v, a}^{k} ‖_{2}^{2} \\ + \sum_{a = 1}^{N_{a}} \frac{ρ_{3, a}}{2} ‖ R_{p} I_{a} D x - z_{p, a}^{k} + u_{p, a}^{k} ‖_{2}^{2} \\ + \sum_{a = 1}^{N_{a}} \frac{ρ_{4, a}}{2} ‖ W I_{a} D x - z_{T, a}^{k} + u_{T, a}^{k} ‖_{2}^{2} \end{aligned}

(24)we can approximate x^k ⁺ ¹ by performing one-step proximal gradient descent (PGD)²⁸

x^{k + 1} = P r o x_{γ f_{2}} (x^{k} - γ \nabla f_{1} (x^{k}))

(25)where γ > 0 is the step size.

Using the definition of proximity operator and the first-order optimality condition, it is equivalent to the solution of the following linear system

\begin{aligned} [\begin{matrix} (ρ_{1} + \frac{1}{γ}) I + A (Ω^{k + 1})^{T} A (Ω^{k + 1}) + \sum_{a = 1}^{N_{a}} ρ_{2, a} (R_{v} I_{a} D)^{T} R_{v} I_{a} D + \sum_{a = 1}^{N_{a}} ρ_{3, a} (R_{p} I_{a} D)^{T} R_{p} I_{a} D \\ + \sum_{a = 1}^{N_{a}} ρ_{4, a} (W I_{a} D)^{T} W I_{a} D \end{matrix}] x \\ = A (Ω^{k + 1})^{T} b (Ω^{k + 1}) + ρ_{1} (z^{k} - u^{k}) + \sum_{a = 1}^{N_{a}} ρ_{2, a} (R_{v} I_{a} D)^{T} (z_{v, a}^{k} - u_{v, a}^{k}) + \sum_{a = 1}^{N_{a}} ρ_{3, a} (R_{p} I_{a} D)^{T} (z_{p, a}^{k} - u_{p, a}^{k}) \\ + \sum_{a = 1}^{N_{a}} ρ_{4, a} (W I_{a} D)^{T} (z_{T, a}^{k} - u_{T, a}^{k}) + \frac{x^{k}}{γ} - \nabla f_{1} (x^{k}) \end{aligned}

(26)

Observing the positive definiteness of the system matrix, one can apply conjugate gradient method²⁹ to solve it. The detailed computation of the gradient of f₁ is provided in Section S2 in the Supplemental material. Performing one step PGD for the first subproblem is equivalent to replacing the f₁ by its local linear approximation. In the context of optimization, the whole algorithm is called linearized ADMM.

The z subproblem admits the closed-form solution as the follows

(z^{k + 1})_{i} = {\begin{matrix} 0, i f {(x^{k + 1} + u^{k})}_{i} < \frac{G_{m i n}}{2} \\ m a x ({(x^{k + 1} + u^{k})}_{i}, G_{m i n}), o t h e r w i s e \end{matrix}

(27)where the subscript (·)_i denotes the i^th entry of a given vector.

The $z_{v, a}^{k + 1}$ can be updated by soft thresholding as follows

\begin{aligned} z_{v, a}^{k + 1} = s i g n (R_{v} I_{a} D x^{k + 1} + u_{v, a}^{k}) \\ ⊙ m a x (| R_{v} I_{a} D x^{k + 1} + u_{v, a}^{k} | - \frac{C_{T V}}{ρ_{2, a}}), a = 1, 2, \dots, N_{a} \end{aligned}

(28)where

⊙

is the Hadamard operation.

The $z_{p, a}^{k + 1}$ subproblem is nonconvex, we choose the following local minimal

\begin{aligned} _{p, a}^{k + 1} = s i g n (R_{p} I_{a} D x^{k + 1} + u_{p, a}^{k}) \\ ⊙ m a x (| R_{p} I_{a} D x^{k + 1} + u_{p, a}^{k} | - \frac{C_{T V}}{ρ_{3, a}}), a = 1, 2, \dots, N_{a} \end{aligned}

(29)

Similarly, the z^k ⁺ ¹ _T,a is updated by

z_{T, a}^{k + 1} = s i g n (W I_{a} D x^{k + 1} + u_{T, a}^{k}) ⊙ m a x (| W I_{a} D x^{k + 1} + u_{T, a}^{k} | - \frac{C_{T V} C_{T}}{ρ_{4, a}}), a = 1, 2, \dots, N_{a}

(30)

Material

The efficacy of EUD optimization is validated in comparison with CONV and SFRT (DO and JDPO respectively), by comparing the following five different scenarios: (1) conventional IMPT (CONV); (2) DO (c_TV = c_EUD = 0 in Eq (1)); (3) DO with EUD optimization (DO-EUD, c_TV = 0); (4) JDPO (c_EUD = 0); (5) JDPO with EUD optimization (JDPO-EUD, c_TV ≠ 0, c_EUD ≠ 0), for head-and-neck (HN), lung and prostate patients respectively. For the HN, the target volume is in close proximity to critical structures including the spinal cord, brainstem, and parotid glands. In particular, the spinal cord lies within 3-5 mm of the posterior target boundary, requiring high conformity to maintain OAR sparing while ensuring target coverage. In the lung case, the tumor is located in the left upper lobe and is adjacent to the esophagus, heart, and spinal cord. The esophagus may partially overlap with the high-dose region. For the prostate case, the target lies in immediate proximity to the rectum and bladder. The posterior target edge is within 2-3 mm of the rectal wall, and the bladder is typically 5-10 mm anterior to the target, posing a steep dose gradient challenge. The radiosensitivity parameters for the TCP model were selected based on established radiobiological characteristics of the corresponding tissues. A value of α/β=2 was used for the brainstem to represent late-responding central nervous system tissue, α/β=10 for the esophagus as an acutely responding organ, and α/β=1.5 for the rectum to capture its late-reacting nature in prostate radiotherapy.^30–32

A generic proton pencil beam scanning (PBS) system was modeled, with energy layers ranging approximately from 70 MeV to 230 MeV, consistent with the specifications in matRad. Both treatment planning and optimization were implemented entirely within matRad, using its built-in dose calculation. The dose influence matrix D was generated using with 5 mm lateral spacing on 3 mm³ dose grid. The following clinically-used beam angles were used: (45^o, 135^o, 225^o, 315^o) for the HN, (0^o, 120^o, 240^o) for the lung and (90^o, 270^o) for the prostate.

Clinically-used DVH constraints were adopted, and all plans were normalized to D_95% = 100% in CTV. The conformity index (CI) was evaluated, defined as $V_{100, C T V}^{2}$ /( $V_{C T V} \times V_{100}$ ). ( $V_{100, C T V}$ : CTV volume receiving at least 100% of prescription dose; $V_{C T V}$ : CTV volume; $V_{100}$ : total body volume receiving at least 100% of prescription dose; ideally CI = 1).

In the following result section, the dose quantities are in percentage with respect to the prescription dose; the survival fractions are in percentage with respect to total volume of the structures under consideration; the values of PVDR and CI are unitless.

The proposed optimization algorithm is implemented in MATLAB R2023a using an in-house developed script. All experiments are performed on a desktop equipped with an Intel Core Ultra 9 processor, 64 GB RAM, and the Windows 11 (64-bit) operating system. The augmented Lagrangian parameters $ρ_{1}, ρ_{2, a}, ρ_{3, a}, ρ_{4, a}$ a are selected according to the settings described in ref.⁵ The step size is manually set to $γ = 1 e^{- 3}$ for all experiments. The algorithm terminates when the maximum number of iterations is reached.

Results

EUD and SF

EUD optimization reduced the EUD and increased SF for OAR. Specifically, EUD optimization decreased the brainstem EUD from 12.0% to 8.78% for DO and from 12.4% to 9.01% for JDPO in HN, reduced the esophagus EUD from 2.47% to 2.11% for DO and from 2.54% to 2.31% for JDPO in lung, reduced the rectum EUD from 4.45% to 4.16% for DO and from 4.29% to 4.09% for JDPO in prostate. EUD optimization also increased the brainstem SF from 72.9% to 82.6% for DO and from 71.5 to 81.4% for JDPO in HN, increased the esophagus SF from 60.1% to 65.2% for DO and from 59.0% to 62.3% for JDPO in lung, increased the rectum SF from 49.5% to 53.3% for DO and from 51.6% to 54.0% for JDPO in prostate.

OAR Sparing

In terms of OAR sparing, EUD optimization reduced the physical dose to OAR. For example, in HN, EUD optimization decreased the mean dose to the brainstem from 8.40% and 8.23% to 4.98% for CONV and DO respectively, and from 8.76% to 5.17% for JDPO. In lung, the mean dose to the esophagus was reduced from 6.83% and 8.27% to 4.29% for CONV and DO respectively, and from 7.73% to 5.96% for JDPO. For prostate, the mean dose to the rectum decreased from 20.5% to 18.7% for both CONV and DO, and from 20.7% to 19.2% for JDPO. These improvements are also demonstrated through the dose plots (a-e) and the DVH plot (f) and Figure 1(h) at 135°, Figure 2(h) at 240°, and Figure 3(h) at 90°.

Figure 1.

HN. (a) Dose Plot for CONV; (b) Dose Plot for DO; (c) Dose Plot for DO-EUD; (d) Dose Plot for JDPO; (e) Dose Plot for JDPO-EUD; (f) DVH for Brainstem; (g) BEV Dose Profiles; (h) BEV Dose Profiles of One Slice. The Dose Plot Window is [0%, 120%], with 80% (Yellow), 100% (red) Isodose Line, and Targets (Green) Highlighted in Dose Plots.

Figure 2.

Lung. (a) Dose Plot for CONV; (b) Dose Plot for DO; (c) Dose Plot for DO-EUD; (d) Dose Plot for JDPO. (e) Dose Plot for JDPO-EUD; (f) DVH for Esophagus; (g) BEV Dose Profiles; (h) BEV Dose Profiles of one Slice. the Dose Plot Window is [0%, 120%], with 80% (Yellow), 100% (red) Isodose Line, and Targets (Green) Highlighted in Dose Plots.

Figure 3.

Prostate. (a) Dose Plot for CONV; (b) Dose Plot for DO; (c) Dose Plot for DO-EUD; (d) Dose Plot for JDPO; (e) Dose Plot for JDPO-EUD; (f) DVH for Rectum; (g) BEV Dose Profiles; (h) BEV Dose Profiles of One Slice. The Dose Plot Window is [0%, 120%], with 80% (Yellow), 100% (Red) Isodose Line, and Targets (Green) Highlighted in Dose Plots.

PVDR

EUD optimization generally preserved PVDR for DO or JDPO, by comparing DO and DO-EUD or comparing JDPO and JDPO-EUD in Table 1. The PVDR was also compared by the BEV dose profiles (g) in Figures 1-3.

Table 1.

Plan Parameters.

		CONV	DO	DO-EUD	JDPO	JDPO-EUD
HN	D_max (%)	109.5	109.4	117.3	110.8	116.8
	CI	0.81	0.77	0.74	0.81	0.72
	D_mean−body (%)	4.72	4.78	4.86	4.76	4.75
	D_mean−OAR (%)	8.40	8.23	4.98	8.76	5.17
	EUD_OAR (%)	11.8	12.0	8.78	12.4	9.01
	SF_OAR (%)	73.4	72.9	82.6	71.5	81.4
	PVDR_45°	-	1.29	1.30	1.78	1.46
	PVDR_135°	-	1.84	2.06	2.11	2.30
	PVDR_225°	-	1.54	1.55	2.21	1.92
	PVDR_315°	-	1.46	1.52	2.08	1.75
	Time	11.2	29.2	31.1	41.4	46.5
Lung	D_max (%)	113.3	112.0	118.7	114.1	120.0
	CI	0.88	0.85	0.84	0.83	0.82
	D_mean−body (%)	3.95	4.13	4.01	4.05	4.05
	D_mean−OAR (%)	6.83	8.27	4.29	7.73	5.96
	EUD_OAR (%)	2.48	2.47	2.11	2.54	2.31
	SF_OAR (%)	59.9	60.1	65.2	59.0	62.3
	PVDR_0°	-	1.23	1.23	1.54	1.33
	PVDR_120°	-	1.67	1.66	1.80	1.74
	PVDR_240°	-	1.37	1.42	1.60	1.55
	Time	38.7	88.3	105.8	122.6	128.0
Prostate	D_max (%)	108.0	108.1	117.3	110.8	114.5
	CI	0.84	0.82	0.82	0.82	0.82
	D_mean−body (%)	1.93	1.94	1.81	1.79	1.78
	D_mean−OAR (%)	20.5	20.5	18.7	20.7	19.2
	EUD_OAR (%)	4.42	4.45	4.16	4.29	4.09
	SF_OAR (%)	50.0	49.5	53.3	51.6	54.0
	PVDR_90°	-	1.52	1.54	1.64	1.56
	PVDR_270°	-	1.48	1.48	1.64	1.48
	Time	27.2	55.1	56.4	77.4	75.5

The Dosimetric Quantities from top to Bottom: Max Dose for the Target D_max, Conformity Index (CI), Mean Body Dose D_mean−Body, Mean OAR Dose D_mean−OAR, OAR Equivalent Uniform Dose EUD_OAR, OAR Survival Fraction SF_OAR, and PVDR per Field. Optimization Time (Time) is in Seconds. OAR Refers to Brainstem, Esophagus, and rectum for HN, Lung, and Prostate Cases, Respectively.

Discussion

In this work, we have developed a novel EUD optimization method for SFRT treatment planning, with feasibility demonstrated for UTD proton GRID. To address the unique optimization challenges from EUD optimization, we have developed an efficient optimization algorithm via linearized ADMM.

The results demonstrate that the proposed method effectively reduces the EUD of the OARs in both the DO and JPDO models, leading to a higher survival fraction. However, incorporating the EUD objective slightly compromises target coverage, as indicated by an increased maximum target dose and a decreased CI across all cases. Additional parameters presented in Table A1 (Appendix), including the EUD and tumor control probability (TCP) for the target, also show a modest degradation under EUD optimization. Nevertheless, the extent of this degradation is limited. However, EUD optimization results in lower predicted toxicities. As shown in Table A1, we report dosimetric metrics such as V_x−OAR (the percentage or volume of an OAR receiving at least x Gy), and D_1cc−OAR (the minimum dose delivered to the hottest 1 cm³ of tissue) which serves as an important surrogate for toxicity risk in serial organs. From the results, we see that EUD optimization led to notable reductions in D_1cc for critical structures, such as the brainstem and esophagus. These reductions may have meaningful biological implications, including a reduced risk of radiation-induced complications, such as necrosis or neuropathy, acute esophagitis,^33,34 etc Regarding PVDR, EUD optimization generally yields lower values in the DO and JPDO models. However, this reduction is not substantial, and in certain instances, such as PVDR₁₃₅° for HN in the DO-EUD and JPDO-EUD models, the PVDR is even higher.

Note that the joint optimization of dose, PVDR, and EUD objectives is multi-criteria and may pose optimization tradeoffs.³⁵ In this study, we posed the uplimit in max target dose to be 120% of the prescription dose. However, while violating this uplimit, EUD and SF can be further improved by increasing the EUD regularization parameter c_EUD in Equation (1).

In this work, we use EUD as a metric and optimize it to achieve better OAR sparing. However, the traditional EUD model faces significant limitations when applied to SFRT. EUD assumes a homogeneous dose-response relationship across the tumor and normal tissues, which does not align with the highly heterogeneous dose distributions characteristic of SFRT. This model also neglects spatial dose modulation effects and fails to account for biological phenomena such as the bystander effect,^36–38 and immune activation,³⁹ all of which are prominent in SFRT. Consequently, EUD may underestimate or misrepresent the therapeutic potential and toxicity risks in SFRT, limiting its reliability for outcome prediction and treatment optimization in this context. It is important to emphasize that this study serves as a preliminary effort toward establishing an optimization framework that is extensible to more sophisticated radiobiological objectives. For any specific biological model—whether described by partial or ordinary differential equations (PDE/ODE), or explicitly dependent on dose—our framework remains applicable, provided a mathematical relationship between the model and beam intensity can be established. In such cases, the proposed linearized ADMM algorithm can still be utilized. The primary technical concern lies in the convergence of the method. Although a rigorous convergence analysis is not provided in this work, a degree of smoothness with respect to the beam intensity is generally required for the biological model. When the objective function is nonsmooth or exhibits an unbounded Lipschitz continuous gradient, the algorithm may experience numerical instability. To mitigate this issue, advanced strategies such as Bregman iteration⁴⁰ or Moreau smoothing⁴¹ can be employed to enhance stability and facilitate convergence. In summary, this work opens up a number of future research directions in biologically informed radiotherapy optimization.

In this work, we have only considered the EUD minimization for the purpose of OAR sparing. However, the same methodology can be adapted to maximize the EUD for targets with a heterogeneous dose distribution, such as LATTICE^42,43 and modulated-target-dose GRID,^44–47 for which the sign of the regularization parameter c_EUD should be negative.

The developed linearized ADMM can be generally applicable to other RT problems with nonlinear terms, such as energy layer optimization,^48–52 normal tissue complication probability optimization,⁵³ and biological dose optimization.⁵⁴ The potential benefit of linearized ADMM is to improve both solution accuracy and computational efficiency. In contrast, proximal gradient descent (PGD) method^28,55,56 can be easily trapped in local minimums.

All in all, the proposed EUD-based optimization framework demonstrates strong potential for future clinical translation in SFRT planning with proton GRID delivery. By incorporating EUD minimization, the method enables improved sparing of organs at risk while maintaining clinically acceptable target coverage. While EUD may not fully capture the spatial dose heterogeneity inherent in SFRT, it serves as a practical and computationally tractable surrogate for biological effect modulation. The observed reductions in dosimetric metrics indicate possible clinical advantages in mitigating normal tissue toxicity. Moreover, the computational efficiency of the linearized ADMM algorithm enhances the feasibility of integrating this approach into existing treatment planning workflows. With its flexibility to accommodate more advanced radiobiological models, the framework provides a promising foundation for the development of biologically informed radiotherapy planning strategies.

Conclusion

This work has developed a novel SFRT treatment planning method with EUD optimization and demonstrated the efficacy of EUD optimization for reducing EUD and increasing SF to improve OAR sparing.

Supplemental Material

sj-docx-1-tct-10.1177_15330338251380609 - Supplemental material for Equivalent-Uniform-Dose Optimization for Spatially Fractionated Radiation Therapy

Supplemental material, sj-docx-1-tct-10.1177_15330338251380609 for Equivalent-Uniform-Dose Optimization for Spatially Fractionated Radiation Therapy by Ya-Nan Zhu, PhD, Weijie Zhang, PhD, Yuting Lin, PhD, and Hao Gao, PhD in Technology in Cancer Research & Treatment

Footnotes

Acknowledgment

The authors are very thankful to the valuable comments from reviewers. This research is partially supported by the NIH grants No. R37CA250921, R01CA261964.

ORCID iD

Ya-Nan Zhu

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was partially supported by the National Institutes of Health, (grant number R37CA250921, R01CA261964).

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Supplemental material

Supplemental material for this article is available online.

Appendix

Table A1.

Additional Plan Parameters.

		CONV	DO	DO-EUD	JDPO	JDPO-EUD
HN	EUD (Gy)	2.03	2.04	1.96	2.06	1.97
	TCP(%)	80.7	77.2	74.6	81.0	73.1
	D_max−OAR (%)	80.7	77.2	74.6	81.0	73.1
	V_20−OAR (%)	7.90	7.65	4.13	8.51	4.13
	D_1cc−OAR (%)	48.3	44.7	31.0	45.6	32.6
Lung	EUD (Gy)	2.05	2.08	2.04	2.09	2.10
	TCP(%)	98.7	98.8	98.6	98.8	98.8
	D_max−OAR (%)	44.6	48.5	35.7	39.6	37.1
	V_20−OAR (%)	0.30	2.66	0.10	0.25	0.15
	D_1cc−OAR (%)	24.7	35.6	21.4	28.4	23.7
Prostate	EUD (Gy)	1.83	1.83	1.84	1.78	1.75
	TCP(%)	72.7	72.6	73.2	69.3	66.9
	D_max−OAR (%)	101.5	102.5	104.0	104.0	105.1
	V_43−OAR (%)	3.29	3.35	2.75	2.86	2.81
	D_1cc−OAR (%)	99.2	99.0	99.6	99.5	99.5

The Dosimetric Quantities from top to Bottom: EUD and Tumor Control Probability (TCP) for the Target, max OAR Dose D_max−OAR, the Percentage (or Volume) of a Specified Organ That Receives at Least x Gy of Radiation Dose V_x−OAR, minimum Dose Delivered to the Hottest 1 cm³ of Tissue (D_1cc−OAR). OAR Refers to Brainstem, Esophagus, and rectum for HN, Lung, and Prostate Cases, Respectively.

References

Zhang

Mayr

. Spatially Fractionated, Microbeam and FLASH Radiation Therapy: A physics and Multi-Disciplinary Approach. IOP Publishing; 2023.

Meyer

Stewart

Smith

, et al. Biological and dosimetric characterisation of spatially fractionated proton minibeams. Phys Med Biol. 2017;62(24):9260.

Kraus

Winter

Zhang

, et al. Treatment planning study for microbeam radiotherapy using clinical patient data. Cancers (Basel). 2022;14(3):685.

Ahmed

Bicher

Stewart

, et al. Dosimetric quantities and cell survival for spatially fractionated radiation therapy. Front Phys. 2023;10:1064860.

Zhang

Traneus

Lin

Chen

Gao

. A novel treatment planning method via scissor beams for uniform-target-dose proton GRID with peak-valley-dose-ratio optimization. Med Phys. 2024;51(10):7047-7056.

Fernandez-Palomo

Chang

Prezado

. Should peak dose be used to prescribe spatially fractionated radiation therapy? —a review of preclinical studies. Cancers (Basel). 2022;14(15):3625.

Bortfeld

. Clinically relevant intensity modulation optimization using physical criteria. XII International Conference on the Use of Computers in Radiation Therapy, 1; 1997; Medical Physics Publishing.

Mohan

. Algorithms and functionality of an intensity modulated radiotherapy optimization system. Med Phys. 2000;27(4):701-711.

Zhu

Sahoo

Zhang

, et al. Intensity modulated proton therapy treatment planning using single-field optimization: The impact of monitor unit constraints on plan quality. Med Phys. 2010;37(3):1210-1219.

10.

Lin

Kooy

Craft

Depauw

Flanz

Clasie

. A Greedy reassignment algorithm for the PBS minimum monitor unit constraint. Phys Med Biol. 2016;61(12):4665.

11.

Gao

Clasie

Liu

Lin

. Minimum MU optimization (MMO): An inverse optimization approach for the PBS minimum MU constraint. Phys Med Biol. 2019;64(12):125022.

12.

Lin

Clasie

Liu

McDonald

Langen

Gao

. Minimum-MU and sparse-energy-layer (MMSEL) constrained inverse optimization method for efficiently deliverable PBS plans. Phys Med Biol. 2019;64:205001.

13.

Cai

Chen

Fan

Gao

. Minimum-monitor-unit optimization via a stochastic coordinate descent method. Phys Med Biol. 2022;67:015009.

14.

Zhu

Zhang

Lin

Lominska

Gao

. An orthogonal matching pursuit optimization method for solving minimum-monitor-unit problems: Applications to proton IMPT, ARC and FLASH. Med Phys. 2023;50:4710-4720.

15.

Zhang

Lin

Wang

Chen

Gao

. TVL1-IMPT: Optimization of peak-to-valley dose ratio via joint total-variation and L1 dose regularization for spatially fractionated pencil-beam-scanning proton therapy. Int J Radiat Oncol* Biol* Phys. 2023;115(3):768-778.

16.

Bartzsch

Stewart

Meyer

. Radiobiological modeling and equivalent uniform dose for spatially fractionated radiation therapy. In: Spatially Fractionated, Microbeam and FLASH Radiation Therapy: A Physics and Multi-Disciplinary Approach. IOP Publishing; 2023:14–1.

17.

Gao

. Hybrid proton-photon inverse optimization with uniformity-regularized proton and photon target dose. Phys Med Biol. 2019;64(10):105003.

18.

Gao

Lin

, et al. Simultaneous dose and dose rate optimization (SDDRO) for FLASH proton therapy. Med Phys. 2020;47(12):6388-6395.

19.

Gao

Liu

Lin

, et al. Simultaneous dose and dose rate optimization (SDDRO) of the FLASH effect for pencil-beam-scanning proton therapy. Med Phys. 2022;49(3):2014-2025.

20.

Lin

Rotondo

Gao

. An iterative convex relaxation method for proton LET optimization. Phys Med Biol. 2023;68(5):055002.

21.

Ouyang

Chen

Lan

Pasiliao

Jr . An accelerated linearized alternating direction method of multipliers. SIAM J Imag Sci. 2015;8(1):644-681.

22.

Gabay

Mercier

. A dual algorithm for the solution of nonlinear variational problems via finite element approximation[J]. Comput Math Appl. 1976;2(1):17-40.

23.

Glowinski

Marroco

. On the approximation by finite elements of order one, and resolution, penalisation-duality for a class of nonlinear Dirichlet problems[J]. ESAIM: Math Modell Numer Analy-Math Modell Numer Anal. 1975;9(R2):41-76.

24.

Boyd

Parikh

Chu

Peleato

Eckstein

. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations Trends® Mach Learn. 2011;3(1):1-122.

25.

Goldstein

Osher

. The split Bregman method for L1-regularized problems[J]. SIAM J Imag Sci. 2009;2(2):323-343.

26.

Gao

. Robust fluence map optimization via alternating direction method of multipliers with empirical parameter optimization. Phys Med Biol. 2016;61(7):2838.

27.

Wang

Yin

Zeng

. Global convergence of ADMM in nonconvex nonsmooth optimization. J Sci Comput. 2019;78:29-63.

28.

Combettes

Wajs

. Signal recovery by proximal forward-backward splitting. Multiscale Model Simul. 2005;4(4):1168-1200.

29.

Hestenes

Stiefel

. Methods of Conjugate Gradients for Solving Linear Systems (Vol. 49, No. 1). NBS; 1952.

30.

Fowler

. The linear-quadratic formula and progress in fractionated radiotherapy. Br J Radiol. 1989;62(740):679-694.

31.

Bentzen

Constine

Deasy

, et al. Quantitative analyses of normal tissue effects in the clinic (QUANTEC): An introduction to the scientific issues. Int J Radiat Oncol* Biol* Phys. 2010;76(3):S3-S9.

32.

Joiner

van der KogeJ

. Basic Clinical Radiobiology. CRC press; 2009.

33.

Xue

Goldman

Grimm

, et al. Dose-volume effects on brainstem dose tolerance in radiosurgery. J Neurosurg. 2012;117:189-196.

34.

Werner-Wasik

Yorke

Deasy

Nam

Marks

. Radiation dose-volume effects in the esophagus. Int J Radiat Oncol* Biol* Phys. 2010;76(3):S86-S93.

35.

Breedveld

Craft

Van Haveren

Heijmen

. Multi-criteria optimization and decision-making in radiotherapy. Eur J Oper Res. 2019;277(1):1-19.

36.

Cahoon

Giacometti

Casey

, et al. Investigating spatial fractionation and radiation induced bystander effects: A mathematical modelling approach. Phys Med Biol. 2021;66(22):225007.

37.

Balderson

Kirkby

. Potential implications on TCP for external beam prostate cancer treatment when considering the bystander effect in partial exposure scenarios. Int J Radiat Biol. 2014;90(2):133-141.

38.

Powathil

Munro

Chaplain

Swat

. Bystander effects and their implications for clinical radiation therapy: Insights from multiscale in silico experiments. J Theor Biol. 2016;401:1-14.

39.

Mortezaee

Najafi

. Immune system in cancer radiotherapy: Resistance mechanisms and therapy perspectives. Crit Rev Oncol Hematol. 2021;157:103180.

40.

Hanzely

Richtarik

Xiao

. Accelerated Bregman proximal gradient methods for relatively smooth convex optimization. Comput Optim Appl. 2021;79:405-440.

41.

Rockafellar

. Convex Analysis (Vol. 28). Princeton university press; 1997.

42.

Ahmed

Wright

Gupta

Pollack

Ahmed

. On modern technical approaches of three-dimensional high-dose lattice radiotherapy (LRT). Cureus. 2010;2(3).

43.

Zhang

Lin

Wang

Badkul

Chen

Gao

. Lattice position optimization for LATTICE therapy. Med Phys. 2023;50:7359-7367.

44.

Köhler

. Theorie einer Methode, bisher unmöglich unanwendbar hohe Dosen Röntgenstrahlen in der Tiefe des Gewebes zur therapeutischen Wirksamkeit zu bringen ohne schwere Schädigung des Patienten, zugleich eine Methode des Schutzes gegen Röntgenverbrennung überhaupt. Fortschr Geb Roentgenstr. 1909;14:27-29.

45.

Mohiuddin

Curtis

Grizos

Komarnicky

. Palliative treatment of advanced cancer using multiple nonconfluent pencil beam radiation: A pilot study. Cancer. 1990;66(1):114-118.

46.

Gao

Mohiuddin

Hartsell

Pankuch

. Spatially fractionated (GRID) radiation therapy using proton pencil beam scanning (PBS): A feasibility study. Int J Radiat Oncol* Biol* Phys. 2015;93(3):E562.

47.

Gao

Mohiuddin

Hartsell

Pankuch

. Spatially fractionated (GRID) radiation therapy using proton pencil beam scanning (PBS): Feasibility study and clinical implementation. Med Phys. 2018;45(4):1645-1653.

48.

Ruan

Lyu

Zou

Dong

Sheng

. A novel energy layer optimization framework for spot-scanning proton arc therapy. Med Phys. 2020;47(5):2072-2084.

49.

Zhang

Shen

Lin

Chen

Long

Gao

. Energy layer optimization via energy matrix regularization for proton spot-scanning arc therapy. Med Phys. 2022;49(9):5752-5762.

50.

Wuyckens

Saint-Guillain

Janssens

, et al. Treatment planning in arc proton therapy: Comparison of several optimization problem statements and their corresponding solvers. Comput Biol Med. 2022;148:105609.

51.

Engwall

Battinelli

Wase

, et al. Fast robust optimization of proton PBS arc therapy plans using early energy layer selection and spot assignment. Phys Med Biol. 2022;67(6):065010.

52.

Zhang

Long

Lin

Chen

Gao

. A treatment plan optimization method with direct minimization of number of energy jumps for proton arc therapy. Phys Med Biol. 2023;68(8):085001.

53.

Zaider

Amols

. Practical considerations in using calculated healthy-tissue complication probabilities for treatment-plan optimization. Int J Radiat Oncol* Biol* Phys. 1999;44(2):439-447.

54.

Lin

Shen

Chen

Gao

. Biological optimization for hybrid proton-photon radiotherapy. Phys Med Biol. 2024;69:115040.

55.

Beck

Teboulle

. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J Imag Sci. 2009;2(1):183-202.

56.

Shen

Zhang

Lin

Rotondo

Long

Gao

. Beam angle optimization for proton therapy via group-sparsity based angle generation method. Med Phys. 2023;50:3258-3273.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.06 MB