Optimizing linear energy transfer distribution in intensity-modulated proton therapy using the alternating direction method of multipliers

Purpose: This study develop a novel linear energy transfer (LET) optimization method for intensity-modulated proton therapy (IMPT) with minimum monitor unit (MMU) constraint using the alternating direction method of multipliers (ADMM).

Material and methods: The novel LET optimization method (ADMM-LET) was proposed with (1) the dose objective and the LET objective as the optimization objective and (2) the non-convex MMU threshold as a constraint condition. ADMM was used to solve the optimization problem. In the ADMM-LET framework, the optimization process entails iteratively solving the dose sub-problem and the LET sub-problem, simultaneously ensuring compliance with the MMU constraint. Three representative cases, including brain, liver, and prostate cancer, were utilized to evaluate the performance of the proposed method. The dose and LET distributions from ADMM-LET were compared to those obtained using the published iterative convex relaxation (ICR-LET) method.

Results: The results demonstrate the superiority of ADMM-LET over ICR-LET in terms of LET distribution while achieving a comparable dose distribution. More specifically, for the brain case, the maximum LET (unit: keV/µm) at the optic nerve decreased from 5.45 (ICR-LET) to 1.97 (ADMM-LET). For the liver case, the mean LET (unit: keV/µm) at the clinical target volume increased from 4.98 (ICR-LET) to 5.50 (ADMM-LET). For the prostate case, the mean LET (unit: keV/µm) at the rectum decreased from 2.65 (ICR-LET) to 2.14 (ADMM-LET).

Conclusion: This study establishes ADMM-LET as a new approach for LET optimization with the MMU constraint in IMPT, offering potential improvements in treatment outcomes and biological effects.

1 Introduction

The most significant advantage of proton beam therapy is utilizing its Bragg peak characteristics to spare the dose in the healthy tissue or organs-at-risk (OARs) at the distal end of the target. A previous review has reported that the overall dose of the proton is about 60% lower than that of the photon radiotherapy (1). In the past few years, proton beam therapy has been advanced to intensity-modulated proton therapy (IMPT) based on active scanning technology, where it optimizes each spot’s intensities and position as well as energy layer distribution, providing a much better conformability compared to the conventional passive-scattering technique (2, 3).

However, challenges remain in proton beam therapy, in which the relative biological effectiveness (RBE) of proton beam therapy still used at 1.1 for the last half-century. Because RBE is affected by numerous factors, such as type of particle, dose, organ type, biological micro-environment, clinical endpoint, and linear energy transfer (LET), it is very challenging, if not impossible, to apply the variable RBE model in the routine clinical settings. Even though we know the constant 1.1 might not be precise, it is the most acceptable value at present (4). Fortunately, a monotonic relationship exists between RBE and LET along the proton beam path (5, 6). As a physics parameter, the LET distribution can be well defined, calculated, and optimized in the treatment plan (7–11), which has the potential to be used as a surrogate of the RBE dose calculation or correlated to the clinical end point (12, 13). Thus, to provide IMPT with practical guidance on biological effects, existing studies agree that this can be achieved by modulating the LET distribution. Therefore, designing an optimization algorithm to LET optimization is one of the development prospects of IMPT.

The minimum monitor unit (MMU) threshold in IMPT is influenced by various factors, including the noise level of the monitor chamber, beam current stability, beam shutdown time at a specified monitor unit, and interlock requirements of the primary and redundant monitor chambers (14, 15). Incorporating the MMU constraint in IMPT treatment planning ensures the deliverability of planned spots within the predefined machine limitations, thereby enhancing treatment delivery feasibility and accuracy (16–19). This integration of the MMU constraint in IMPT provides an essential mechanism to guarantee the reliable execution of treatment plans, improving overall treatment quality and patient outcomes.

In current studies focused on incorporating LET optimization into IMPT, two primary methods are commonly employed: direct optimization of biological dose (12, 20–22) and simultaneous optimization of LET and physical dose (9, 13, 23, 24). For instance, Unkelbach et al. (20) proposed a strategy to reduce the distribution of LET at OARs by re-optimizing the initial IMPT treatment plan based on the additional biological dose. However, due to uncertainties associated with conversion coefficients and planning constraints in the biological dose, recent studies have increasingly favored simultaneous optimization of both LET and physical dose (13, 23, 24).

Nevertheless, while these studies have made significant contributions to LET optimization in IMPT, they did not consider the crucial MMU constraint in the optimization (9, 12, 20–23), which substantially impacts the accurate delivery of IMPT. Although Liu et al. (13) and Li et al. (24) consider the deliverability by specifying the MMU through the postprocessing procedure, its plan quality in this way is not as good as considering the MMU constraint in optimization (19).

To address this limitation, Li et al. (25) developed an optimization method based on the iterative convex relaxation (ICR) approach that explicitly incorporates the MMU constraint. The ICR method is often used to assist in solving optimization problems in IMPT, such as to address nonconvexity caused by dose-volume-histogram objective (19), nonlinearity caused by the dose rate constraint (26). However, it is important to note that when dealing with the nonlinearity of LET as it relates to the optimization variable, employing the ICR method, which relies on linear approximation, can introduce numerical errors when solving nonlinear systems. In addition, the ICR approach did not consider to transfer high LET distributions from OARs to tumors. Given the characteristics of LET optimization in IMPT, which involves a non-convex MMU constraint, a linear physical dose objective concerning the optimization variable, and a nonlinear LET objective, a straightforward and effective method is needed to decouple the non-convex constraint and efficiently solve the nonlinear LET objective.

The alternating direction method of multipliers (ADMM) is a simple but powerful algorithm well suited to distributed convex and some nonconvex optimization problems (27–29). ADMM has played a significant role in the inverse optimization of IMPT’s treatment planning (15, 30, 26, 31, 32). It enables solving the inverse optimization problem with the MMU constraint in IMPT while considering the dose optimization objective and other optimization objectives simultaneously.

To summarize, this study aims to investigate the potential of LET optimization with the MMU constraint in IMPT using the alternating direction method of multipliers (ADMM-LET). This approach allows for the escalation of the LET distribution in the target while mitigating the LET distribution in the OARs.

2 Methods and materials

2.1 Formalize the dose objective

The dose objectives are integrated into the ADMM-LET framework. In order to facilitate the algorithm’s presentation, we initially express the original dose objectives formally. These objectives encompass several terms: the squared deviation term, squared over dosage term, squared under dosage term, maximum dose-volume-histogram constraint (DVH) term, and minimum DVH term.

$\begin{array}{ll}\begin{array}{l}f\left(d\right)=p_1\sum _{i\in {\text{Ω}}_{\text{ROI},1}}{\left(d_i-{\widehat{d}}_1\right)}^2\\ +p_2\sum _{i\in {\text{Ω}}_{\text{ROI},2}}\text{Θ}\left(d_i-{\widehat{d}}_2\right){\left(d_i-{\widehat{d}}_2\right)}^2\\ +p_3\sum _{i\in {\text{Ω}}_{\text{ROI},3}}\text{Θ}\left({\widehat{d}}_3-d_i\right){\left(d_i-{\widehat{d}}_3\right)}^2\\ +p_4\sum _{i\in {\text{Ω}}_{\text{ROI},4}}\text{Θ}\left(d_i-{\widehat{d}}_4\right)\text{Θ}\left({\tilde{d}}_{\text{max}}-d_i\right){\left(d_i-{\widehat{d}}_4\right)}^2\\ +p_5\sum _{i\in {\text{Ω}}_{\text{ROI},5}}\text{Θ}\left({\widehat{d}}_5-d_i\right)\text{Θ}\left(d_i-{\tilde{d}}_{\text{min}}\right){\left(d_i-{\widehat{d}}_5\right)}^2\end{array}& (1)\end{array}$

In Equation 1, f(d) represents the dose fidelity term in the sum of squares. d_i represents the dose of the i-th voxel, $\widehat{d}$ denotes the prescribed dose, and $\tilde{d}$ represents the dose at the prescribed volume. The term is calculated over the region of interest (ROI), corresponding to the voxel index Ω. The weights assigned to each term in the objective function are denoted by p. The Heaviside function, Θ, is used as a step function to enforce certain conditions based on the dose values.

The objective function consists of the squared deviation term, which is convex, as the first term. However, the second to fifth terms introduce non-convexity due to the presence of the Heaviside function. To address this non-convexity, the active set, denoted as Ω_d, needs to be determined. The active set Ω_d is obtained as the union of five subsets: ${\Omega }_d= {\cup }_{i=1}^5{\Omega }_i$. Specifically, the subsets are defined as Equation 2:

$\begin{array}{ll}\begin{array}{l}{\text{Ω}}_1={\text{Ω}}_{\text{ROI},1}\hfill \\ {\text{Ω}}_2=\{k\in {\text{Ω}}_{\text{ROI},2}|d_k\ge {\widehat{d}}_2\}\hfill \\ {\text{Ω}}_3=\{k\in {\text{Ω}}_{\text{ROI},3}|d_k\le {\widehat{d}}_3\}\hfill \\ {\text{Ω}}_4=\{k\in {\text{Ω}}_{\text{ROI},4}|d_k\ge {\widehat{d}}_4\wedge d_k\le {\tilde{d}}_{\text{max}}\}\hfill \\ {\text{Ω}}_5=\{k\in {\text{Ω}}_{\text{ROI},5}|d_k\le {\widehat{d}}_5\wedge d_k\ge {\tilde{d}}_{\text{min}}\}\hfill \end{array}& (2)\end{array}$

In summary, the objective function f(d) can be represented as Equation 3 (33):

$\begin{array}{ll}f(d)=\sum _{s\in {\text{Ω}}_d}{p}_s{\Vert d_s-{\widehat{d}}_s\Vert }_2^2=\sum _{s\in {\text{Ω}}_d}{p}_s{\Vert D_{s}x-{\widehat{d}}_s\Vert }_2^2={\Vert Ax-b\Vert }_{{\text{Ω}}_d}^2.& (3)\end{array}$

where D represents the dose influence matrix, x represents the spot weights that need to be optimized, D_sx denotes the dose distribution at the s-th voxel. The matrices are related as follows: $A^{T}A={\displaystyle {\sum }_s{D}_s^T{p}_s{D}_s}$ and $A^{T}b={\displaystyle {\sum }_s{D}_s^T{p}_s{\widehat{d}}_s}$.

2.2 The LET objective

In order to enhance the LET at the target region and reduce the LET at the OARs, this study employed the squared under LET term for the target and the squared over LET term for the OARs (24), the form of the LET objective is shown in Equation 4:

$\begin{array}{ll}\begin{array}{l}ɡ(LET)=w_{let}\sum _{i\in {\text{Ω}}_{\text{CTV}}}\text{Θ} \left(L\widehat{E}{T}_{\text{CTV}}-LE{T}_i\right){ \left(LE{T}_i-L\widehat{E}{T}_{\text{CTV}}\right)}^2\\ +w_{let}\sum _{i\in {\text{Ω}}_{\text{OARs}}}\text{Θ} \left(LE{T}_i-L\widehat{E}{T}_{\text{OARs}}\right) {\left(LE{T}_i-L\widehat{E}{T}_{\text{OARs}}\right)}^2\end{array}& (4)\end{array}$

Where

$\begin{array}{ll}LE{T}_i=\frac{1}{d_i}\sum _{j=1}^{N_s}(D_{ij}{L}_{ij}{x}_j),i=1,\cdots ,N_v,& (5)\end{array}$

LET_i, defined as Equation 5, represents the LET distribution at the i-th voxel, L the LET influence matrix, N_s the number of spots, N_v the number of voxels. D_ij and L_ij represent the dose and LET delivered by the j-th spot located at the i-th voxel, respectively.

Similar to the dose objective, the LET objective is also non-convex due to the presence of the Heaviside function. The active set Ω_l can be determined as the same set as Ω_d. Thus, the LET objective can be expressed as Equation 6:

$\begin{array}{ll}g(L(x))=w_{let}{\Vert L(x)-l_0\Vert }_{{\text{Ω}}_l}^2.& (6)\end{array}$

where L(x) represents the calculated LET distribution based on the spot weights x, and l₀ denotes the desired LET value. This formulation captures the deviation between the calculated LET distribution and the desired LET value within the active set Ω_l.

2.3 The LET optimization problem

Overall, the LET optimization problem involves optimizing the spot weights x by considering both dose and LET objectives. Simultaneously, the problem takes into account the MMU constraint. The objective is to find the optimal distribution of spot weights that satisfies the dose requirements, enhances LET in the target region, and minimizes LET in the OARs while adhering to the MMU constraint. The LET optimization problem considered in this study is defined as Equation 7:

$\begin{array}{ll}\begin{array}{c}\underset{x}{\min } f(d)+ɡ(LET),\\ \text{s}.\text{t}.\{\begin{array}{l}{d}_i=\sum _{j=1}^{N_s} (D_{ij}{x}_j),i=1,\cdots ,N_v,\hfill \\ LE{T}_i=\frac{1}{d_i}\sum _{j=1}^{N_s} (D_{ij}{L}_{ij}{x}_j),\hfill \\ x\in \left\{0\right\}{\cup }^{}[ɡ,+\infty ),\hfill \end{array}\end{array}& (7)\end{array}$

where g represents the MMU threshold of the proton therapy system.

To facilitate the demonstration of the ADMM algorithm for solving this problem, we first express it in standard form as Equation 8:

$\begin{array}{ll}\begin{array}{c}\underset{x,z_1,z_2}{\min } f (x)+ɡ (z_1)+m (z_2),\\ \text{s}.\text{t}.\{\begin{array}{l}x-z_1=0,\hfill \\ x-z_2=0,\hfill \end{array}\end{array}& (8)\end{array}$

This formulation introduces z₁ and z₂ as dummy variables. Specifically, z₁ separates the dose and LET objectives, while z₂ is employed to decouple the MMU constraint. The function m(z₂) enforces the MMU constraint and is defined as Equation 9:

$\begin{array}{ll}m(z_2)=\{\begin{array}{l}0,\text{\hspace{1em}if\hspace{1em}}{z}_2\in \left\{0\right\}{\cup }^{ }[ɡ,+\infty ),\hfill \\ +\infty ,\text{\hspace{1em}otherwise}.\hfill \end{array}.& (9)\end{array}$

By incorporating the dummy variables and introducing the appropriate constraints, the LET optimization problem is transformed into a standard form that can be effectively solved using the ADMM algorithm.

2.4 The algorithm via ADMM

The augmented Lagrange function of Equation 8 is indeed as Equation 10:

$\begin{array}{ll}\begin{array}{ll}L(x,z_1,u_1,z_2,u_2)=\hfill & f(x)+g(z_1)+m(z_2)\hfill \\ +{\rho }_1\Vert x-z_1+u_1{\Vert }^2+{\rho }_2\Vert x-z_2+u_2{\Vert }^2,\hfill \end{array}& (10)\end{array}$

where u₁ and u₂ are the scaled dual variables associated with the equality constraints, and ρ₁ and ρ₂ are the penalty parameters controlling the strength of the penalty terms.

The ADMM iteration loop, which aims to solve the optimization problem defined by Equation 8, can be summarized as Equation 11:

$\begin{array}{ll}\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{x}^{k+1}=\underset{x}{\arg \min }L\left(x,z_1^k,u_1^k,z_2^k,u_2^k\right),\\ z_1^{k+1}=\underset{z_1}{\arg \min }L\left(x^{k+1},z_1,u_1^k,z_2^k,u_2^k\right),\end{array}\\ u_1^{k+1}=u_1^k+x^{k+1}-z_1^{k+1},\\ z_2^{k+1}=\underset{z_2}{\arg \min }L\left(x^{k+1},z_1^{k+1},u_1^{k+1},z_2,u_2^k\right),\end{array}\\ u_2^{k+1}=u_2^k+x^{k+1}-z_2^{k+1}.\end{array}& (11)\end{array}$

For the x sub-problem (also noted as the dose sub-problem), its specific form is to solve the following problem (Equation 12):

$\begin{array}{ll}{x}^{k+1}=\underset{x}{\arg \min }{\Vert Ax-b\Vert }^2+{\rho }_1{\Vert x-z_1^k+u_1^k\Vert }^2+{\rho }_2{\Vert x-z_2^k+u_2^k\Vert }^2,& (12)\end{array}$

Therefore, updating x comes from the optimal condition of the dose sub-problem taking ∂_xL = 0 (29, 34), and its solution form is shown in Equation 13:

$\begin{array}{ll}{x}^{k+1}={[A^{T}A+\left({\rho }_1+{\rho }_2\right)I]}^{-1}\left[{\rho }_1\left(z_1^k-u_1^k\right)+{\rho }_2\left(z_2^k-u_2^k\right)+A^{T}b\right]& (13)\end{array}$

The dose sub-problem is formulated as a differentiable least squares problem, requiring the solution of a system of linear equations. In this problem, the coefficient matrix is symmetric and positive definite. To solve this problem efficiently, numerical algebra techniques such as the preconditioned conjugate gradient method (PCG) (35) can be utilized.

The specific form of the z₁ sub-problem, referred to as the LET sub-problem, is to solve the following problem (Equation 14):

$\begin{array}{ll}{z}_1^{k+1}=\underset{z_1}{\arg \min }g(z_1)+{\rho }_1\Vert x^{k+1}-z_1+u_1^k{\Vert }^2,& (14)\end{array}$

It presents a nonlinear least squares problem, which differs from the dose sub-problem. The most suitable method is the gradient descent method to address this sub-problem within the LET optimization framework. Three main types of solvers are commonly used: the gradient algorithm, the quasi-Newton algorithm, and the trust region algorithm. For example, the Barzilar-Borwein (BB) method (36), the limited memory Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS) (37) and the Levenberg-Marquardt-Fletcher method (LMF) (38).

For the z₂ sub-problem (Equation 15):

$\begin{array}{ll}{z}_2^{k+1}=\underset{z_2}{\arg \min }m(z_2)+{\rho }_2\Vert x^{k+1}-z_2+u_2^k{\Vert }^2,& (15)\end{array}$

it has analytical solution given by Equation 16 (31):

$\begin{array}{ll}{z}_2^{k+1}=S(x^{k+1}+u_2^k,ɡ),& (16)\end{array}$

where the projection operator S(z, g) is defined as Equation 17:

$\begin{array}{ll}S(z,ɡ)=\{\begin{array}{c}0,z<ɡ/2\\ \max (z,ɡ),\text{\hspace{0.17em}if\hspace{0.17em}}z\ge ɡ/2\end{array}.& (17)\end{array}$

In other words, if the value z is less than half of the threshold g, the projected value is set to 0. Otherwise, if z is greater than or equal to g/2, the projected value is the maximum of z and g. Applying this projection operator can efficiently solve the z₂ sub-problem, ensuring that the MMU constraint is satisfied.

2.5 Materials

Three representative cases were selected for testing the proposed method, including brain (50 Gy in 25 fractions), liver (60 Gy in 30 fractions), and prostate (60 Gy in 30 fractions). To generate the dose influence matrix and LET influence matrix, an open-source treatment planning platform called matRad (39) was utilized. The beam spot lateral spacing was set to 5 mm, and the dose grid had a resolution of 3 mm³. It is worth noting that the MMU threshold, denoted as g, was set to 1.1 MU, which corresponds to the limit for the Varian ProBeam system (32).

A fair comparison was performed between the new method (ADMM-LET) and a state-of-the-art method (ICR) (25). The ‘ICR’ refers to just considering the dose optimization, while the ‘ICR-LET’ denotes the LET optimization using the ICR method. Four proton beams (45°, 135°, 225°, 315°) were used for brain and liver, and two proton beams (90°, 270°) were used for prostate. The beam angle setting was the same as that of the published ICR method for better comparison of method. Besides, it is worth noting that all plans, including ‘ICR’, ‘ICR-LET’, and ‘ADMM-LET’, utilized the same dose objective. Furthermore, the LET objective was consistent for the ‘ICR-LET’ and ‘ADMM-LET’ plans, ensuring a fair and meaningful comparison between the methods.

The ADMM-LET results presented in this study utilize the BB method as the solver for the LET sub-problem. The details on determining the LET sub-problem solver, including the BB, L-BGFS, and LMF method, are provided in Supplementary Material Section 1.

2.6 Evaluation plan

The dose distribution was evaluated using the conformity index (CI) and dose-volume histogram (DVH), and the LET-volume histogram (LVH) was evaluated to assess the LET distribution. The CI measures the conformity of target coverage and is defined as Equation 18 in matRad (39):

$\begin{array}{ll}CI=\frac{V_{95}^2}{V\times V{D}_{95}},& (18)\end{array}$

where V₉₅ represents the volume of the target receiving at least 95% of the prescription dose, V is the target volume, and VD₉₅ is the total volume enclosed by at least 95% of the prescription isodose line. An ideal CI value is 1, indicating perfect conformity.

3 Results

Figures 1–3 present the dose and LET distributions and their corresponding volume histograms for three representative examples. The assessment for all plans is summarized in Table 1. In addition, the convergence comparison of solutions between the ICR-LET plan and the ADMM-LET plan is presented in Figure 4, and the information regarding the computational efficiency of different methods is provided in Table 2.

Figure 1

Figure 1 The brain case. A representative slice of dose distribution for (A) ICR plan; (B) ICR-LET plan; (C) ADMM-LET plan (the display window is [0, 120%] of prescription dose). A representative slice of LET distribution for (D) ICR plan; (E) ICR-LET plan; (F) ADMM-LET plan (the display window is [0, 6] keV/µm). (G, H) Comparison of DVH, LVH between ICR (dashed line), ICR-LET (dotted line), and ADMM-LET (solid line). D, Dose; L, LET.

Figure 2

Figure 2 The liver case. A representative slice of dose distribution for (A) ICR plan; (B) ICR-LET plan; (C) ADMM-LET plan (the display window is [0, 120%] of prescription dose). A representative slice of LET distribution for (D) ICR plan; (E) ICR-LET plan; (F) ADMM-LET plan (the display window is [0, 7] keV/µm). (G-H) Comparison of DVH, LVH between ICR (dashed line), ICR-LET (dotted line), and ADMM-LET (solid line). D, Dose; L, LET.

Figure 3

Figure 3 The prostate case. A representative slice of dose distribution for (A) ICR plan; (B) ICR-LET plan; (C) ADMM-LET plan (the display window is [0, 120%] of prescription dose). A representative slice of LET distribution for (D) ICR plan; (E) ICR-LET plan; (F) ADMM-LET plan (the display window is [0, 7] keV/µm). (G-H) Comparison of DVH, LVH between ICR (dashed line), ICR-LET (dotted line), and ADMM-LET (solid line). D, Dose; L, LET.

Table 1

Table 1 The dose and LET evaluation parameters of all plans.

Figure 4

Figure 4 Convergence comparison of solutions between the ICR-LET plan with the ADMM-LET plan. (A) Brain case, (B) Liver case, and (C) Prostate case.

Table 2

Table 2 The dimension information of the optimization problem and the computational efficiency of different methods.

3.1 The brain case

The ICR, ICR-LET, and ADMM-LET plans demonstrate similar dose distributions. The ICR and ADMM-LET plans achieve a CI of 0.91, and the CI of the ICR-LET plan is 0.92. Representative dose distribution slices for these plans are shown in Figures 1A–C, while Figure 1G compares the corresponding DVH. The LET distribution, as depicted in Figures 1D–F, indicates improved results for the ICR-LET and ADMM-LET plans compared to the ICR plan. The ADMM-LET plan exhibits the most favorable LET distribution among the three plans. Specifically, the mean LET (unit: keV/µm) within the CTV increased from 3.76 (ICR) to 4.61 (ICR-LET) and 4.83 (ADMM-LET), respectively. The maximum LET (unit: keV/µm) at the optic chiasm decreased from 5.13 (ICR) to 4.11 (ICR-LET) and 3.97 (ADMM-LET), respectively.

3.2 The liver case

For the liver case, both the ICR-LET and ADMM-LET plans exhibit similar dose distributions. The CI of the ADMM-LET plan (0.81) and the ICR-LET plan (0.88) is slightly lower than that of the ICR plan (0.89). Figures 2A–C display representative dose distribution slices for the three plans, and Figure 2G compares the corresponding DVH. As depicted in Figures 2D–F, the LET distribution shows improved results for the ICR-LET and ADMM-LET plans compared to the ICR plan. Notably, the ADMM-LET plan exhibits the most favorable LET distribution among the three plans. Specifically, the mean LET (unit: keV/µm) within the CTV increased from 3.55 (ICR) to 4.98 (ICR-LET) and 5.50 (ADMM-LET), respectively.

3.3 The prostate case

The three plans yield similar dose distributions, with respective CI values of 0.72, 0.75, and 0.73. Figures 3A–C display representative dose distribution slices for these plans, and Figure 3G compares the corresponding DVH. Representative slices of the LET distribution in Figures 3D–F indicate improved LET distributions for both the ICR-LET and ADMM-LET plans compared to the ICR plan. Notably, the ADMM-LET plan demonstrates the most favorable LET distribution among the three plans. Specifically, the mean LET (unit: keV/µm) within the CTV increased from 3.24 (ICR) to 3.32 for both the ICR-LET and ADMM-LET plans. The mean LET (unit: keV/µm) at the rectum decreased from 2.87 (ICR) to 2.65 (ICR-LET) and 2.14 (ADMM-LET), respectively. Moreover, the mean LET (unit: keV/µm) at the bladder decreased from 1.56 (ICR) to 1.51 (ICR-LET) and 1.40 (ADMM-LET), respectively.

4 Discussion

The topic of LET optimization in proton beam therapy has garnered significant interest within the particle therapy community (8–10, 12, 13, 20–23, 25). Among the various optimization frameworks available, ADMM has emerged as a prominent approach for addressing this challenge and has been successfully employed in numerous inverse optimization studies (15, 30, 26, 31, 32). The results of this study demonstrate that ADMM-LET has the potential to effectively regulate the distribution of LET while preserving the dose distribution. By leveraging the ADMM framework, the proposed method balances dose objectives, LET objectives, and the MMU constraint. ADMM-LET enables the generation of treatment plans that optimize the dose distribution and control the LET distribution, leading to improved treatment outcomes.

Some existing studies only spare the LET distribution in OARs (12, 25, 40). However, it is also important to transfer high LET distribution from OARs to the tumors (21, 23, 24). Therefore, the optimization objective of LET in this study is not only to reduce the high LET distribution in OARs but also to transfer the high LET distribution from OARs to the tumors.

LET optimization in IMPT is a bi-objective optimization problem, where the goal is to find a trade-off between dose distribution and LET distribution. The solution to the problem requires how the bi-objective trade-off is considered. ADMM-LET achieves this trade-off by adjusting the weight value of the objectives. This trade-off method is not specific to ADMM-LET and can be applied to other optimization algorithms. In Supplementary Material Section 2, the test results demonstrate that optimizing the LET distribution more can lead to a compromise in the dose distribution. This trade-off between the objectives is inherent due to the physical characteristics of proton beams. Through the ADMM-LET platform, we could find a suitable weight value that can ensure the dose distribution and modulate the LET distribution well.

On the other hand, this study’s results and phenomena are based on a fixed number of gantry angles in the IMPT. The degree of freedom will significantly increase with the new development of Spot-scanning Proton Arc therapy (SPArc). Recent publications indicated a better ability of LET modulation with more beam angles (41–45). In SPArc therapy, the beam angle or arc trajectory selection and its associated spot and energy layer placement or optimization play a critical role in LET optimization (24, 41). Furthermore, to make the LET-integrated SPArc therapy more challenging, the treatment delivery time plays a key role in the degree of freedom. The number of spots, the energy layer, and the proton gantry’s mechanical constraint determine the total delivery time. Thus, this study serves as a starting point toward the advantage of SPArc LET optimization.

A recent study has highlighted that while ADMM effectively solves the optimization problem in IMPT with small MMU thresholds, it may face challenges when dealing with larger MMU thresholds (19). In such cases, alternative methods are required to tackle the LET optimization problem. One promising approach for addressing the optimization problem with larger MMU thresholds is the stochastic coordinate descent method (19). However, the specific implementation and utilization of this method for LET optimization in IMPT still needs further investigation. Besides, the ADMM-LET method proposed in this study currently relies on the open-source dose calculation and optimization toolkit [matRad (39)] for proton therapy treatment planning. This is still a long way from a commercial treatment planning system (TPS). Therefore, this is a direction that this study needs to strive for in the future.

5 Conclusion

This work proposed an ADMM-based approach for solving the LET optimization problem with a non-convex MMU constraint in IMPT. The results demonstrated that the ADMM framework could balance the dose and LET objectives while considering the MMU constraint, which has the potential advantage of modulating the LET distribution in IMPT treatment plans.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving humans were approved by Ethical Committees of Union Hospital, Tongji Medical College, Huazhong University of Science and Technology. The studies were conducted in accordance with the local legislation and institutional requirements. The ethics committee/institutional review board waived the requirement of written informed consent for participation from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements.

Author contributions

QF: Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft. XZ: Investigation, Software, Visualization, Writing – original draft. RD: Investigation, Software, Visualization, Writing – original draft. YQ: Investigation, Software, Visualization, Writing – original draft. LZ: Methodology, Validation, Writing – original draft. XL: Methodology, Validation, Writing – original draft. XD: Conceptualization, Supervision, Validation, Writing – review & editing. GL: Conceptualization, Data curation, Funding acquisition, Methodology, Project administration, Supervision, Writing – review & editing. SD: Conceptualization, Data curation, Funding acquisition, Methodology, Project administration, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study is supported by the National Natural Science Foundation of China (No. 12005072), the Open Research Funding of Hubei Key Laboratory of Precision Radiation Oncology (No. jzfs003, No. jzfs014).

Acknowledgments

The authors would like to thank the editor and two reviewers for their constructive comments and suggestions that have helped improve this manuscript’s quality.

Conflict of interest

XD reports honorarium from Ion Beam Application IBA and Elekta outside the work presented here. XD received industry research funding from Elekta and RadioMed outside the work presented here.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2024.1328147/full#supplementary-material

References

1. DeLaney TF. Proton therapy in the clinic. In: Meyer J, Dawson L, Kavanagh B, Purdy J, Timmerman R, editors. IMRT IGRT SBRT- Advances in the Treatment Planning and Delivery of Radiotherapy, vol. 43. Basel, Switzerland: Frontiers of Radiation Therapy and Oncology (2011). p. 465–85.

Google Scholar

2. Lomax A. Intensity modulation methods for proton radiotherapy. Phys Med Biol (1999) 44:185–205. doi: 10.1088/0031-9155/44/1/014

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lomax AJ, Boehringer T, Coray A, Egger E, Goitein G, Grossmann M, et al. Intensity modulated proton therapy: A clinical example. Med Phys (2001) 28:317–24. doi: 10.1118/1.1350587

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Paganetti H, Niemierko A, Ancukiewicz M, Gerweck LE, Goitein M, Loeffler JS, et al. Relative biological effectiveness (RBE) values for proton beam therapy. Int J Radiat OncologyBiologyPhysics (2002) 53:407–21. doi: 10.1016/S0360-3016(02)02754-2

CrossRef Full Text | Google Scholar

5. Wilkens JJ, Oelfke U. A phenomenological model for the relative biological effectiveness in therapeutic proton beams. Phys Med Biol (2004) 49:2811–25. doi: 10.1088/0031-9155/49/13/004

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Carabe A, Moteabbed M, Depauw N, Schuemann J, Paganetti H. Range uncertainty in proton therapy due to variable biological effectiveness. Phys Med Biol (2012) 57:1159–72. doi: 10.1088/0031-9155/57/5/1159

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Frese MC, Wilkens JJ, Huber PE, Jensen AD, Oelfke U, Taheri-Kadkhoda Z. Application of constant Vs. variable relative biological effectiveness in treatment planning of intensitymodulated proton therapy. Int J Radiat OncologyBiologyPhysics (2011) 79:80–8. doi: 10.1016/j.ijrobp.2009.10.022

CrossRef Full Text | Google Scholar

8. Grassberger C, Trofimov A, Lomax A, Paganetti H. Variations in linear energy transfer within clinical proton therapy fields and the potential for biological treatment planning. Int J Radiat OncologyBiologyPhysics (2011) 80:1559–66. doi: 10.1016/j.ijrobp.2010.10.027

CrossRef Full Text | Google Scholar

9. Giantsoudi D, Grassberger C, Craft D, Niemierko A, Trofimov A, Paganetti H. Linear energy transfer-guided optimization in intensity modulated proton therapy: feasibility study and clinical potential. Int J Radiat OncologyBiologyPhysics (2013) 87:216–22. doi: 10.1016/j.ijrobp.2013.05.013

CrossRef Full Text | Google Scholar

10. Deng W, Yang Y, Liu C, Bues M, Mohan R, Wong WW, et al. A critical review of LET-based intensity-modulated proton therapy plan evaluation and optimization for head and neck cancer management. Int J Particle Ther (2021) 8:36–49. doi: 10.14338/IJPT-20-00049.1

CrossRef Full Text | Google Scholar

11. Feng H, Shan J, Anderson JD, Wong WW, Schild SE, Foote RL, et al. Per-voxel constraints to minimize hot spots in linear energy transfer-guided robust optimization for base of skull head and neck cancer patients in IMPT. Med Phys (2022) 49:632–47. doi: 10.1002/mp.15384

PubMed Abstract | CrossRef Full Text | Google Scholar

12. An Y, Shan J, Patel SH, Wong W, Schild SE, Ding X, et al. Robust intensity-modulated proton therapy to reduce high linear energy transfer in organs at risk. Med Phys (2017) 44:6138–47. doi: 10.1002/mp.12610

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Liu C, Patel SH, Shan J, Schild SE, Vargas CE, Wong WW, et al. Robust optimization for intensity modulated proton therapy to redistribute high linear energy transfer from nearby critical organs to tumors in head and neck cancer. Int J Radiat OncologyBiologyPhysics (2020) 107:181–193. doi: 10.1016/j.ijrobp.2020.01.013

CrossRef Full Text | Google Scholar

14. Clasie BM, Kooy HM, Flanz JB. PBS machine interlocks using EWMA. Phys Med Biol (2016) 61:400–12. doi: 10.1088/0031-9155/61/1/400

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Gao H, Clasie B, Liu T, Lin Y. Minimum MU optimization (MMO): an inverse optimization approach for the PBS minimum MU constraint. Phys Med Biol (2019) 64:125022. doi: 10.1088/1361-6560/ab2133

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zhu XR, Sahoo N, Zhang X, Robertson D, Li H, Choi S, 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:1210–9. doi: 10.1118/1.3314073

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Cao W, Lim G, Li X, Li Y, Zhu XR, Zhang X. Incorporating deliverable monitor unit constraints into spot intensity optimization in intensity-modulated proton therapy treatment planning. Phys Med Biol (2013) 58:5113–25. doi: 10.1088/0031-9155/58/15/5113

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Shan J, An Y, Bues M, Schild SE, Liu W. Robust optimization in IMPT using quadratic objective functions to account for the minimum MU constraint. Med Phys (2018) 45:460–9. doi: 10.1002/mp.12677

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Cai J-F, Chen RC, Fan J, Gao H. Minimum-monitor-unit optimization via a stochastic coordinate descent method. Phys Med Biol (2022) 67:015009. doi: 10.1088/1361-6560/ac4212

CrossRef Full Text | Google Scholar

20. Unkelbach J, Botas P, Giantsoudi D, Gorissen BL, Paganetti H. Reoptimization of intensity modulated proton therapy plans based on linear energy transfer. Int J Radiat OncologyBiologyPhysics (2016) 96:1097–106. doi: 10.1016/j.ijrobp.2016.08.038

CrossRef Full Text | Google Scholar

21. Bai X, Lim G, Grosshans D, Mohan R, Cao W. A biological effect-guided optimization approach using beam distal-edge avoidance for intensity-modulated proton therapy. Med Phys (2020) 47:3816–25. doi: 10.1002/mp.14335

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Gu W, Ruan D, Zou W, Dong L, Sheng K. Linear energy transfer weighted beam orientation optimization for intensity-modulated proton therapy. Med Phys (2021) 48:57–70. doi: 10.1002/mp.14329

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Cao W, Khabazian A, Yepes PP, Lim G, Poenisch F, Grosshans DR, et al. Linear energy transfer incorporated intensity modulated proton therapy optimization. Phys Med Biol (2017) 63:015013. doi: 10.1088/1361-6560/aa9a2e

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Li X, Ding X, Zheng W, Liu G, Janssens G, Souris K, et al. Linear energy transfer incorporated spot-scanning proton arc therapy optimization: A feasibility study. Front Oncol (2021) 11:698537. doi: 10.3389/fonc.2021.698537

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Li W, Lin Y, Li H, Rotondo R, Gao H. An iterative convex relaxation method for proton LET optimization. Phys Med Biol (2023) 68:055002. doi: 10.1088/1361-6560/acb88d

CrossRef Full Text | Google Scholar

26. Gao H, Lin B, Lin Y, Fu S, Langen K, Liu T, et al. Simultaneous dose and dose rate optimization (SDDRO) for FLASH proton therapy. Med Phys (2020) 47:6388–95. doi: 10.1002/mp.14531

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Glowinski R, Marroco A. FrenchSur l’approximation, par elements finis d’ordre un, et la resolution, par penalisation-dualite, d’une classe de problemes de dirichlet non lineaires. Rev Fr Autom Inf Rech Oper (1975) 9:41–76. doi: 10.1051/M2AN/197509R200411

CrossRef Full Text | Google Scholar

28. Gabay D, Mercier B. A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Comput Mathematics Appl (1976) 2:17–40. doi: 10.1016/0898-1221(76)90003-1

CrossRef Full Text | Google Scholar

29. Boyd S, Parikh N, Chu E, Peleato B, EcKstein J. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations Trends Mach Learn (2010) 3:1–122. doi: 10.1561/2200000016

CrossRef Full Text | Google Scholar

30. Gao H. Hybrid proton-photon inverse optimization with uniformity-regularized proton and photon target dose. Phys Med Biol (2019) 64:105003. doi: 10.1088/1361-6560/ab18c7

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Gao H, Liu J, Lin Y, Gan GN, Pratx G, Wang F, et al. Simultaneous dose and dose rate optimization (SDDRO) of the FLASH effect for pencil-beam-scanning proton therapy. Med Phys (2022) 49:2014–25. doi: 10.1002/mp.15356

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Zhang W, Li W, Lin Y, Wang F, Chen RC, Gao H. 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 OncologyBiologyPhysics (2023) 115:768–78. doi: 10.1016/j.ijrobp.2022.09.064

CrossRef Full Text | Google Scholar

33. Zhang G, Long Y, Lin Y, Chen RC, Gao H. A treatment plan optimization method with direct minimization of number of energy jumps for proton arc therapy. Phys Med Biol (2023) 68:085001. doi: 10.1088/1361-6560/acc4a7

CrossRef Full Text | Google Scholar

34. Gao H. Robust fluence map optimization via alternating direction method of multipliers with empirical parameter optimization. Phys Med Biol (2016) 61:2838–50. doi: 10.1088/0031-9155/61/7/2838

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Barrett R, Berry M, Chan TF, Demmel J, Donato J, Dongarra J, et al. Templates for the solution of linear systems: building blocks for iterative methods, 2nd Edition. (1994) Philadelphia, PA: SIAM.

Google Scholar

36. Barzilai J, Borwein JM. Two-point step size gradient methods. IMA J Numerical Anal (1988) 8:141–8. doi: 10.1093/imanum/8.1.141

CrossRef Full Text | Google Scholar

37. Liu D, Nocedal J. On the limited memory BFGS method for large-scale optimization. Math Programming (1989) 45:503–28. doi: 10.1007/BF01589116

CrossRef Full Text | Google Scholar

38. Nocedal J, Wright SJ. Numerical Optimization. New York, NY: Springer (2006). doi: 10.1007/978-0-387-40065-5

CrossRef Full Text | Google Scholar

39. Wieser H-P, Cisternas E, Wahl N, Ulrich S, Stadler A, Mescher H, et al. Development of the open-source dose calculation and optimization toolkit matRad. Med Phys (2017) 44:2556–68. doi: 10.1002/mp.12251

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bai X, Lim G, Wieser H-P, Bangert M, Grosshans D, Mohan R, et al. Robust optimization to reduce the impact of biological effect variation from physical uncertainties in intensity-modulated proton therapy. Phys Med Biol (2019) 64:025004. doi: 10.1088/1361-6560/aaf5e9

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ding X, Li X, Zhang JM, Kabolizadeh P, Stevens C, Yan D. Spot-scanning proton arc (SPArc) therapy: The first robust and delivery-efficient spot-scanning proton arc therapy. Int J Radiat OncologyBiologyPhysics (2016) 96:1107–16. doi: 10.1016/j.ijrobp.2016.08.049

CrossRef Full Text | Google Scholar

42. Li X, Liu G, Janssens G, De Wilde O, Bossier V, Lerot X, et al. The first prototype of spot-scanning proton arc treatment delivery. Radiotherapy Oncol (2019) 137:130–6. doi: 10.1016/j.radonc.2019.04.032

CrossRef Full Text | Google Scholar

43. Liu G, Li X, Zhao L, Zheng W, Qin A, Zhang S, et al. A novel energy sequence optimization algorithm for efficient spot-scanning proton arc (SPArc) treatment delivery. Acta Oncol (2020) 59:1178–85. doi: 10.1080/0284186X.2020.1765415

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Liu G, Li X, Qin A, Zhou J, Zheng W, Zhao L, et al. Is proton beam therapy ready for single fraction spine sbrs? - a feasibility study to use spot-scanning proton arc (SPArc) therapy to improve the robustness and dosimetric plan quality. Acta Oncol (2021) 60:653–7. doi: 10.1080/0284186X.2021.1892183

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Liu G, Zhao L, Li X, Zhang S, Dai S, Lu X, et al. A novel ultra-high dose rate proton therapy technology: spot-scanning proton arc therapy FLASH (SPLASH). Int J Radiat OncologyBiologyPhysics (2023) 117(3):730–7. doi: 10.1016/j.ijrobp.2023.05.012

CrossRef Full Text | Google Scholar