A Molecular Dynamics Approach for Optimizing Beam Intensities in IMPT Treatment Planning

This paper uses a molecular dynamics (MD) method for intensity-modulated proton therapy treatment planning optimization to overcome the problem of gradient-based optimization methods such as quasi-Newton, which is sensitive to starting conditions and is easily trapped in local minima. We implemented a molecular dynamics (MD) method and a quasi-Newton method for plan optimization. Three types of cancer cases, prostate, head-and-neck and lung cancer, were tested with three starting different initial conditions. Over-all, the MD method consistently resulted in solutions with lower objective function values (OFVs) compared to those from the quasi-Newton algorithm. Furthermore, the MD method converged on the same OFV regardless of its initial starting points used for the prostate cancer case.


Introduction
Radiation therapy uses ionizing radiation to treat malignant tumors by damaging the DNA of cancer cells to block their ability to proliferate. The purpose is to deliver a prescribed dosage of radiation to the targeted tumor while minimizing the dose deposition in healthy tissues [1]. As a new and advanced radiation treatment modality, proton therapy has rapidly gained interest [2] [3]. The deposited dose of a proton beam starts from a low entrance level, its energy increases gradually while increasing depth, then suddenly jumps to a sharp peak known as the Bragg peak. Once the dose deposition reaches a few millimeters beyond this peak, it falls sharply to zero (Figure 1). This shows a proton beam's accuracy-delivering almost no dosage to regions beyond the targeted area, a How to cite this paper: Liao, L., Lim  feat not typically feasible for photon beams commonly used in intensity modulated radiation therapy (IMRT). The most advanced proton therapy delivery technique is the intensity-modulated proton therapy (IMPT) which uses scanning beamlets (pencil beams) whose weights (intensities) can be modulated independently for treatment. Using data from imaging technology such as CT (computed tomography) or MRI (magnetic resonance imaging), the beamlets weight profile can be optimized according to the prescribed treatment dose. The optimization problem to determine the optimal beamlets weight is called a fluence map optimization (FMO). As a result, IMPT can deliver a highly conformal radiation dose to the desired target tumor while sparing the surrounding organs at risk (OAR), especially for structures with complex shapes [4]. Note that unlike IMRT with a single fixed energy level (6 MV, 10 MV or 15 MV) per field (i.e., treatment angle), an IMPT field consists of multiple energy layers for different depths of the tumor and each layer has a set of beamlets to cover the tumor region. As a result, solving for the FMO of IMPT typically takes much longer than IMRT.
The FMO problem for intensity-modulated radiation therapy (IMRT) planning has been extensively studied and has been addressed by various solution strategies. These strategies can be grossly classified into two groups: global optimization (GO) and local optimization (LO). GO approaches include linear programming [5] [6] [7], mixed integer programming [8] [9], simulated annealing [10] and genetic algorithms [11] [12]. These approaches are designed to reach a global optimal solution. However, they all require an excessive amount of time for optimization, which is not practical in clinical treatment planning. In addition, the performance of these approaches depends heavily on the choice of parameters [7]. On the other hand, LO approaches can obtain usable solutions within a clinically acceptable planning time window. Therefore, LO approaches have been commonly used in commercial treatment planning systems such as Journal of Applied Mathematics and Physics Some of the well-known LO methods include gradient-based algorithms [13] [14] [15] [16] [17], local neighborhood search [18] and iterative methods [19].
However, it has been reported in the IMRT treatment planning literature that FMO is highly degenerate and in many popular optimization methods, it can be easily trapped in local minima [20] [21] [22] [23]. As noted by Llacer, Agazaryan [24], the commonly used gradient-based methods often result in different solutions when different starting points are used. The FMO problem has even greater degrees of freedom in IMPT than in IMRT. This can lead to a higher chance of being trapped in a local minimum, as was pointed out by Albertini, Hug [25].
Despite the fact that many solution approaches have been proposed to solve the FMO problem, a method of avoiding local minima while converging to a solution within a practical time limit has rarely been reported. Hou et al. [26] [27] proposed a method to formulate the IMRT FMO problem into a molecular dynamics problem which motivated this study. Molecular dynamics is a powerful computational technique that is often used to simulate the physical movement of atoms and molecules in a many-body system. In Hou's paper, the beamlets in IMRT were considered as virtual atoms. The weight of the beamlets were formulated as the positions of the virtual atoms and the objective function value (OFV) of the FMO problem in IMRT was formulated as the potential energy of the dynamic system. In classical molecular dynamics, because the movement of atoms follows Newton's Law of Motion, the dynamic system will relax to an equilibrium state with the lowest free energy. In this process, the position and velocity of the atoms will change with time. Thus, following the MD formulation, the beamlets weight and virtual velocity will update over time and the OFV of the FMO problem will be minimized. The MD method's feature of virtual velocity differs from traditional gradient algorithms in that it only updates the weight. Furthermore, within the FMO problem, virtual velocity can help atoms overcome the barrier and get out of the local minima, e.g., an atom goes into a local minimum but can continue to move pass the saddle points and reach a lower point with sufficient speed. In addition, the search direction in MD follows dynamic equations which enable MD to converge faster than many global optimization methods. These advantages give this method the potential to overcome the local minima problem in FMO for radiation treatment planning without spending excessive computation time.
In this study, a MD method was implemented for the FMO problem in IMPT and tested on three patient cases with different starting points. The primary objective of this study is to demonstrate the MD method can be a viable alternative for solving the IMPT FMO problem and evaluate whether this method can alleviate the impact of the starting points, a major issue of gradient-based methods.
The remainder of the paper is organized as followed: Section 2 describes the optimization model, the solution algorithms are in Section 3 which includes an introduction to the MD method and its formulation of the IMPT FMO problem Journal of Applied Mathematics and Physics and an existing quasi-Newton method formulation is shown for the purpose of algorithmic performance comparison. The data used in the experiment and the initial configuration setups are listed in Section 4. Last but not least, Section 5 provides the results regarding the convergence properties of the MD method and quasi-Newton algorithm.

Optimization Model Formulation
The main purpose of IMPT is to deliver the prescribed conformal radiation dose to the targeted tumor while sparing normal tissues. To achieve this goal, we define a quadratic objective function to quantify the difference between the pre- to receive radiation doses up to a certain amount. However, once the amount is over a tolerance value, a penalty will be imposed on the voxel according to the degree of deviation from the tolerance.
The dose i D in voxel i can be calculated as: where j ω is the weight or intensity of beamlet j, which is the decision variable of our IMPT optimization model. Notation N denotes the total number of beamlets and ij k is the unit dose contribution of the j th beamlet to the i th voxel; ij k is also known as the dose deposition coefficient. Here the values of ij k are calculated using an in-house dose calculation engine for proton beamlets [31]. Journal of Applied Mathematics and Physics Using the notation described above, the dose-based objective function for our optimization model is: where tumor n p and oar m p denote the penalty weights of the n th tumor site and m th OAR, respectively. These weights are often obtained by trial and error by planners (dosimetrists, physicists, etc.), to find a balance between tumor dose coverage and OAR dose sparing to satisfy the clinical criteria. In this study, the model follows the common practice in the medical community of containing a physical constraint: the beamlet weight cannot be negative, i.e., 0, 1, 2, , Hence, our optimization model for the IMPT FMO problem is:

Solution Algorithms
In the following sections, we briefly describe the general concept of MD followed by the use of MD for radiation therapy. For the purpose of algorithm performance comparison, we describe a quasi-Newton method for FMO, introduced by Lomax [32] and has been well accepted in the medical community.

Molecular Dynamics
MD is a computational technique for many-body system simulation that has been widely applied in the material sciences community. In a classical MD model, the physical movements of particles in the system follow Newton's Laws of Motion [33]. Let j x be the position and j v be the velocity of a particle j.
Then, force j f is the product of mass j m and acceleration j a of the particle: where t is the time of the system and E is the potential energy of the system. The force j f is related to the acceleration and can be expressed as the gradient of the potential energy of the particle.
Based on Newton's Laws of Motion, the position, velocity and acceleration of the particle can be described as functions of time t; Therefore, the continuous motion configuration of the system can be calculated by integrating Newton's Laws of Motion. When the system is under the influence of continuous potential energy, the positions and velocities can be approximated using a Taylor series expansion for a small time step t ∆ , 0 t ∆ > : where a, b and c are the second, third and fourth time derivatives of the coordinates.
This Taylor expansion serves as the basis for the most common integrators used in MD calculations. Many classical methods for integrating equations require information from both current and previous steps to update the system.
This means the information from these steps must be stored in memory and the system cannot self-start at the beginning [34]. To resolve this problem, Swope et al. [35] introduced the Velocity Verlet method which requires information from the previous step only:

a t a t t
Therefore, this approach is selected in our algorithm to update the MD system to solve the IMPT FMO problem.

MD Method for FMO
In IMPT, the optimization problem can be formulated as a dynamic system with N virtual atoms [26]. Each beamlet weight ( j ω ) is assumed to be the position ( j x ) of a virtual atom j in 1-D dimension. The objective function (F) can be considered as the potential energy (E) of the system. As a result, the dynamic equations for virtual atom j can be expressed as: In this study, we followed the approach of Hou et al. [26], in which the mass of the virtual atom j equals the summation of the unit dose contribution of all voxels influenced by the jth beamlet, written as Hence, we calculate the updating beamlet weights by Equation (12) and we can combine Equation (12) with Equation (4) to calculate the trajectory of the OFVs.
In physics, temperature is used to specify the thermodynamic state of a system. In the MD system, temperature (T) is related to the kinetic energy of the system and can be calculated as: Journal of Applied Mathematics and Physics  (13) where k is the Boltzman constant and N is the total number of particles in the system.
The MD system will converge to an equilibrium state with the lowest free energy. Note that free energy consists of kinetic energy and potential energy. Therefore, the potential energy equals the free energy only when kinetic energy is zero, i.e., the temperature of such system is zero. Thus, the objective function (potential energy) of our FMO model is minimized when the system reaches an equilibrium state with zero system temperature. However, in physics, the kinetic energy and potential energy of a dynamic system follow the law of energy conservation. Although kinetic and potential energy will interchange continuously, the total energy will remain unchanged. This can create an issue of convergence to a specific point because the atoms can still carry significant speed when they reach that point. As a result, atoms may move away from an optimal point with the lowest potential energy. In order to solve this problem, a "friction" to the system (i.e., a damping factor to the MD system) is added to slow the movements of atoms as they get closer to an optimal point. The following damping function is applied in our algorithm: As we mentioned above, their speed may cause the atoms to pass the optimal point and create an issue of convergence. On the other hand, when the virtual atoms become trapped in local minima, a proper speed may help them continue to move and get out of those local minimum points. Using this feature, we employ temperature scaling to adjust the velocities to help the atoms move out from the local minimum points. From the updating function, Equation (12), we define the scaling function as: where 0 T is the initial temperature and d T is the desired temperature.
An important physical constraint of IMPT planning is that the beamlet weight

Quasi-Newton Algorithm
The quasi-Newton algorithm is often used for solving the FMO problem in radiation treatment planning; it has been implemented in Eclipse, a leading radiation mentioned that different gradient-based algorithms for FMO have similar convergence properties and the quasi-Newton method appears to be the fastest one.
So, in this study, we implemented a quasi-Newton method which was proposed by Lomax [32] for the algorithm comparison purpose.
For a standard Newton method, the updating function of a beamlet weight at each iteration is defined as:  [15]. Another limitation of many existing approaches in the clinic is that algorithm implementations are of-Journal of Applied Mathematics and Physics ten based on an unconstrained minimization problem even though a radiation treatment plan requires that all beamlets be nonnegative. Wu and Mohan [14] described a method in which the constraint is checked at the end of each iteration and for any beamlet weights with negative values, they are reset to zero. Liu et al. [30] used 2 j ω instead of j ω in the objective function. In our study, we applied the damping factor introduced by Lomax [32] in the form of: Therefore, the iterative updating function is: Similar to the MD method, the quasi-Newton algorithm stops when either one of the following two conditions is met: 1) a preset number of iterations is reached or 2) there is no improvement in the OFV after a certain number of consecutive iterations.

Numerical Experiments
All algorithms were implemented in C++ and all experiments were performed on a 64-bit Linux workstation with 128 GB of memory and quad Intel Xeon

Patient Data
The three clinical cancer cases from The University of Texas MD Anderson Cancer Center selected for this study was a prostate cancer, head-and-neck cancer and a lung cancer case (see Figure 2). The corresponding beam angles for each case are marked by arrows F1, F2 and F3, respectively.
The beam angles, number of beamlets in each beam, volumes of interest and number of voxels of each volume for each case are listed in Table 1. The prostate case involved a medium-sized tumor that required only a simple treatment plan in which parallel-opposed fields were used. The head-and-neck case had a small Journal of Applied Mathematics and Physics target size but some very subtle OARs such as the optical chiasm. The lung case represents a large sized tumor case. It contains twice as much target volume than the prostate case and twice as large in total volume compared to the head-and-neck case. Because of this, more beamlets were required to cover the tumor for the lung case.
The planned doses and penalty weights for the corresponding VOIs in a dose-based objective function for the three IMPT cases are listed in Table 2. The same penalty was applied to different initial conditions in each case. Because the highest priority was to satisfy the tumor coverage and dose uniformity requirements, the penalty for the target was high. Meanwhile, the target dose to the OARs is set to 0 Gy, which means we wish to minimize the dose on OARs as low as possible.
The parameter values of each algorithm were assigned the same values for all

IMPT Starting Conditions
In this study, the IMPT plans of each of three cases were obtained from three different initial conditions (Figure 3). These initial conditions have been described by Albertini et al. [25] and described briefly below: 1) Forward wedge (FW). All beamlet weights are set the same creating a wedge-shaped dose that has a high dose at the proximal edge and a low dose at the distal edge (Figure 3(a)).
2) Inverse wedge (IW). The beamlet weights are set to distal tracking creating an inverse wedge-shaped dose that has a very low dose to the proximal edge and a high dose to the distal edge (Figure 3(b)).
3) Spread-out Bragg peak (SOBP). The beamlet weights are arranged to deliver a flat dose on the targeted area (Figure 3(c)).

Results
Choosing the prostate cancer case as a representative, the plots of OFV as a function of the number of iterations using the two algorithms are shown in Figure 4.
For each of the three different initial conditions, the MD method consistently reached a lower OFV than the quasi-Newton method. Furthermore, in the same iteration, the MD method converged faster in all instances. Going beyond the iteration did not make much change in OFV. Figure 4 also shows several bumps   We make several observations from this comparison study.
Observation 1: In all cases, the MD method outperformed the quasi-Newton method in terms of OFV.
Observation 2: For each problem instance, the MD method consistently converged on the same objective value regardless of the starting conditions while Overall, the MD method required longer computation times before it reached one of the stopping criteria. For the prostate and head-and-neck cases, the quasi-Newton method had a greater time variation among the initial conditions while the computation times for the MD method differed slightly depending on the initial conditions. For the lung case, the quasi-Newton method reached convergence points about 25% faster than the MD method in three initial conditions.
In Table 3 The dose-volume histograms (DVHs) of the scanning target volume (STV) and femoral heads for the prostate case, the DVHs of the clinical target volume (CTV), brain, brainstem, and optic chiasm for the head-and-neck case and the DVHs of the planning target volume (PTV), spinal cord and esophagus for the lung case are shown in Figure 6. These DVHs show the results using the MD Abbreviations: NA means the stopping criteria were met before the OFV was reached. Figure 6. The dose-volume histograms (DVHs) of three cases using the molecular dynamics method or the quasi-Newton method starting from three initial conditions: forward wedge (FW), inverse wedge (IW) and spread-out Bragg peak (SOBP). Journal of Applied Mathematics and Physics method starting from different initial conditions were almost identical. In contrast, the results using the quasi-Newton method varied from case to case. The higher OFVs may reflect worse target coverage and OARs sparing in DVHs.
The dose distributions in the prostate cancer case using the MD and quasi-Newton methods are shown in Figure 7. The result doses were normalized to 95% of the target volume received at 100% of the prescribed dose for all plans. For the dose distribution, the MD method yielded the same result regardless of the initial conditions used, whereas the results from the quasi-Newton method varied from one case to the next. The quasi-Newton method for the head-and-neck case was the most sensitive to different initial conditions while the lung case was the least sensitive. Furthermore, the dose distribution from the MD method was clinically better (i.e., more conformal) than the quasi-Newton method.

Conclusion
Our MD method was developed to overcome the local entrapment issue observed in many gradient-based algorithms that are extensively used in radiotherapy planning systems in clinics. A previously reported quasi-Newton method for fluence map optimization was implemented as a benchmark for our proposed MD method. The performance of the MD method was compared to the quasi-Newton method using three representative clinical cancer cases. For each method, three popular initial conditions were used for each case. We found that the MD method consistently produced solutions that were the same or within a negligible margin of error regardless of the initial conditions used. The quasi-Newton method, however, was sensitive to its initial conditions by converging on noticeably different solutions. Compared to the quasi-Newton method, the MD method requires longer total computation time to reach the stopping points Figure 7. Intensity-modulated proton therapy dose distributions from the molecular dynamics (MD) method and quasi-Newton method for the prostate cancer case. A-C: Dose distributions using the MD method starting from forward wedge (FW), inverse wedge (IW) and spread-out Bragg peak (SOBP), respectively. D-F: Dose distributions using the quasi-Newton method starting from FW, IW and SOBP, respectively. The orange, red and magenta isodose lines indicate relative doses of 50%, 100% and 110%, respectively. Journal of Applied Mathematics and Physics especially in large tumor sized case. However, to reach the same OFVs, fewer number of iterations were required for the MD method. This indicates that the MD method can not only solve the FMO problem effectively, but also has the potential to obtain a clinically feasible solution in a short time. Future work can include the implementation of the MD method using parallel computation in a multi-core desktop environment.