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.


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)(8)(9)(10)(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)(17)(18)(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)(21)(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)(21)(22)(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)(28)(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

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-volumehistogram constraint (DVH) term, and minimum DVH term.
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, d denotes the prescribed dose, and e 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, Q, 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: i=1 Ω i .Specifically, the subsets are defined as Equation 2: In summary, the objective function f(d) can be represented as Equation 3 (33): where D represents the dose influence matrix, x represents the spot weights that need to be optimized, D s x denotes the dose distribution at the s-th voxel.The matrices are related as follows:

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: Where 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 nonconvex 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: where L(x) represents the calculated LET distribution based on the spot weights x, and l 0 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 .

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: 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: This formulation introduces z 1 and z 2 as dummy variables.Specifically, z 1 separates the dose and LET objectives, while z 2 is employed to decouple the MMU constraint.The function m(z 2 ) enforces the MMU constraint and is defined as Equation 9: +∞, otherwise : : 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.

The algorithm via ADMM
The augmented Lagrange function of Equation 8 is indeed as Equation 10: where u 1 and u 2 are the scaled dual variables associated with the equality constraints, and r 1 and r 2 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: For the x sub-problem (also noted as the dose sub-problem), its specific form is to solve the following problem (Equation 12): Therefore, updating x comes from the optimal condition of the dose sub-problem taking ∂ x L = 0 (29,34), and its solution form is shown in Equation 13: 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 1 sub-problem, referred to as the LET sub-problem, is to solve the following problem (Equation 14): 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 2 sub-problem (Equation 15): it has analytical solution given by Equation 16(31): where the projection operator S(z, g) is defined as Equation 17: 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 2 sub-problem, ensuring that the MMU constraint is satisfied.

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 3 .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.

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 18in matRad (39): where V 95 represents the volume of the target receiving at least 95% of the prescription dose, V is the target volume, and VD 95 is the total volume enclosed by at least 95% of the prescription isodose line.An ideal CI value is 1, indicating perfect conformity.

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, Frontiers in Oncology frontiersin.organd the information regarding the computational efficiency of different methods is provided in Table 2.

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.

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.

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   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.

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)(42)(43)(44)(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.

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.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.

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 1
The dose and LET evaluation parameters of all plans.

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