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  . As a new and advanced radiation treatment modality, proton therapy has rapidly gained interest   . 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
Figure 1. A comparison of dose deposition: proton (high energy) versus proton (low energy) versus photon.
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  . 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    , mixed integer programming   , simulated annealing  and genetic algorithms   . 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  . 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 Eclipse [Varian Associates, Palo Alto, CA] and Pinnacle [Philips, Milpitas, CA]. Some of the well-known LO methods include gradient-based algorithms      , local neighborhood search  and iterative methods  . 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     . As noted by Llacer, Agazaryan  , 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  .
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.   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 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.
2. 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 prescribed dose ( , where n is the index of each organ of interest) and the actual dose ( , where i is the voxel index) delivered to the patient. Although different types of objective functions are reported in the literature   , dose-based quadratic objective functions are commonly used in the medical physics community      for optimizing beam intensities in radiation therapy and is the base model for this paper.
A dose-based objective function F is composed of two parts: for the nth target, , and for the mth OAR, . Because the primary goal of treatment planning is to obtain an actual radiation dose profile that is identical or nearly identical to the prescribed dose level on the target, can be defined as the deviation of the resulting actual dose on voxel i from the target prescription dose :
where is the total numbers of voxels in the nth target. Similarly, can be defined as
where is the total number of voxels in the mth OAR, is the specified tolerance dose for the organ, and is defined as . Note, we introduced the step function because healthy organs are often allowed 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 in voxel i can be calculated as:
where 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 is the unit dose contribution of the jth beamlet to the ith voxel; is also known as the dose deposition coefficient. Here the values of are calculated using an in-house dose calculation engine for proton beamlets  .
Using the notation described above, the dose-based objective function for our optimization model is:
where and denote the penalty weights of the nth tumor site and mth 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., . Hence, our optimization model for the IMPT FMO problem is:
3. 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  and has been well accepted in the medical community.
3.1. 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  . Let be the position and be the velocity of a particle j. Then, force is the product of mass and acceleration of the particle:
where t is the time of the system and E is the potential energy of the system. The force 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 , :
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  . To resolve this problem, Swope et al. (1982)  introduced the Velocity Verlet method which requires information from the previous step only:
Therefore, this approach is selected in our algorithm to update the MD system to solve the IMPT FMO problem.
3.2. MD Method for FMO
In IMPT, the optimization problem can be formulated as a dynamic system with N virtual atoms  . Each beamlet weight ( ) is assumed to be the position ( ) 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.  , 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 , where is the total number of voxels influenced by beamlet j. Using the velocity Verlet method mentioned in Section 3.1, the dynamic updating equations for the IMPT FMO problem is 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:
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 is the initial temperature and is the desired temperature.
An important physical constraint of IMPT planning is that the beamlet weight cannot be negative. Hou and Wang  suggested a barrier potential with an infinite height at to impose this constraint. The virtual atoms are reflected by changing the sign of their velocity each time they try to pass the barrier: , when .
The proposed MD method to optimize the IMPT FMO problem is outlined below. The algorithm stops when either of these following conditions are met: 1) a certain number of iterations is reached or 2) there is no change (smaller than , i.e. 10−4) in the objective value for a certain number of consecutive iterations.
In our implementation, if a new solution increases the OFV to a value larger than
3.3. 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 treatment planning system        . Zhang, et al. (2004) 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  for the algorithm comparison purpose.
For a standard Newton method, the updating function of a beamlet weight at each iteration is defined as:
Because calculating the Hessian matrix and its inverse can be computationally expensive due to the large size of the beamlet vector, the matrix is often approximated by its diagonal matrix   . Another limitation of many existing approaches in the clinic is that algorithm implementations are often based on an unconstrained minimization problem even though a radiation treatment plan requires that all beamlets be nonnegative. Wu and Mohan  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.  used instead of in the objective function. In our study, we applied the damping factor introduced by Lomax  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.
4. 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 E5649 2.53 GHz processor.
Both MD method and quasi-Newton method were used to solve the FMO problems for plan optimization. Different types of cancer patient cases were tested to represent varying complexity and tumor size for the study. There were three different point generation methods used to create different initial conditions for each case. The final objective function value and the objective function values during the optimization process were recorded for result analysis.
4.1. 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 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
Figure 2. The three clinical cancer cases selected for the study and their corresponding field directions. (a) Prostate cancer case; (b) Head-and-neck cancer case; (c) Lung cancer case; the tumors are contoured in red.
Table 1. The intensity-modulated proton therapy beam angles, number of beamlets in each beam, VOIs and number of voxels of each VOI for the three cancer cases.
Abbreviations: VOI: volume of interest, STV: scanning target volume, CTV: clinical target volume and PTV: planning target volume.
Table 2. Dose-based objective function parameters used for optimizing the intensity-modulated proton therapy plans.
cancer cases. The stopping criteria were either: 1) the number of iterations reached 10,000 or 2) there was no change in the OFV for 20 consecutive iterations, where .
4.2. 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.  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)).
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
Figure 3. The initial beamlet weights for a single intensity-modulated proton therapy field from three initial conditions: (a) Forward wedge; (b) Inverse wedge; (c) Spread-out Bragg peak (SOBP).
Figure 4. The OFVs of different methods as functions of the iteration steps for the intensity-modulated proton therapy processes. The green arrows show bumps in the MD optimization process, which do not occur when the quasi-Newton method is used.
along the optimization process of MD. This indicates the OFVs were not monotonically decreased over the iteration which was not the case for the quasi-Newton method where the OFV decreased smoothly and monotonically over time. Similar results were also observed in the other two cases.
Figure 5 compares OFVs and the computation times obtained from different initial conditions for the three cases using the MD and quasi-Newton algorithms. 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
Figure 5. The objective function values (OFVs) and computation times for the molecular dynamics (MD) and quasi-Newton methods for intensity-modulated proton therapy planning for the three cancer cases starting from forward wedge (FW), inverse wedge (IW) and spread-out Bragg peak (SOBP) initial conditions.
OFVs obtained by the quasi-Newton method varied substantially when different starting conditions were utilized.
Observation 3: The IW condition for the quasi-Newton method yielded the worst OFVs.
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 number of iterations required to reach the same or final OFVs for two algorithms are listed. To reach the final OFV, the number of iterations was influenced by the starting points for both algorithms. For example, the MD method took 3750, 3334, 3210 iterations to reach final points in the FW, IW, SOBP condition of the prostate case respectively. Compare to the quasi-Newton method, which took 3921, 5210, 1952 iterations to converge to final results. However, for all three cases, fewer iterations were needed by the MD method to reach the same OFVs compared to the quasi-Newton method. For example, in the prostate of FW condition, to reach the objective function value 2500, the MD method only took 213 iterations compare to the quasi-Newton method, which took 1880 iterations. This leads us to believe the MD method could be used for FMO to obtain a reasonably good treatment plan in a relative short computation time.
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
Table 3. The number of iterations required to reach the same and final OFVs using MD method and quasi-Newton method with different starting conditions.
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).
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.
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.
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.
 Cao, W., et al. (2013) Incorporating Deliverable Monitor Unit Constraints into Spot Intensity Optimization in Intensity-Modulated Proton Therapy Treatment Planning. Physics in Medicine & Biology, 58, 5113.
 Zhang, X., et al. (2010) Intensity-Modulated Proton Therapy Reduces the Dose to Normal Tissue Compared with Intensity-Modulated Radiation Therapy or Passive Scattering Proton Therapy and Enables Individualized Radical Radiotherapy for Extensive Stage IIIB Non-Small-Cell Lung Cancer: A Virtual Clinical Study. International Journal of Radiation Oncology Biology Physics, 77, 357-366.
 Wu, X., Zhu, Y. and Luo, L. (2000) Linear Programming Based on Neural Networks for Radiotherapy Treatment Planning. Physics in Medicine and Biology, 45, 719.
 Langer, M., et al. (1996) A Comparison of Mixed Integer Programming and Fast Simulated Annealing for Optimizing Beam Weights in Radiation Therapy. Medical Physics, 23, 957-964.
 Bortfeld, T., et al. (1990) Methods of Image Reconstruction from Projections Applied to Conformation Radiotherapy. Physics in Medicine and Biology, 35, 1423-1434.
 Pflugfelder, D., et al. (2008) A Comparison of Three Optimization Algorithms for Intensity Modulated Radiation Therapy. Zeitschrift für Medizinische Physik, 18, 111-119.
 Lomax, A., Pedroni, E., Schaffner, B., Scheib, S., Schneider, U. and Tourovsky, A. (1996) 3D Treatment Planning for Conformal Proton Therapy by Spot Scanning. In: Proceedings of 19th LH Gray Conference, BIR, London, 67-71.
 Lim, G.J. and Cao, W. (2012) A Two-Phase Method for Selecting IMRT Treatment Beam Angles: Branch-and-Prune and Local Neighborhood Search. European Journal of Operational Research, 217, 609-618.
 Soukup, M., et al. (2009) Study of Robustness of IMPT and IMRT for Prostate Cancer against Organ Movement. International Journal of Radiation Oncology Biology Physics, 75, 941-949.
 Albertini, F., Hug, E. and Lomax, A. (2010) The Influence of the Optimization Starting Conditions on the Robustness of Intensity-Modulated Proton Therapy Plans. Physics in Medicine and Biology, 55, 2863-2878.
 Hou, Q., et al. (2003) An Optimization Algorithm for Intensity Modulated Radiotherapy—The Simulated Dynamics with Dose-Volume Constraints. Medical Physics, 30, 61-68.
 Lim, G.J., Choi, J. and Mohan, R. (2008) Iterative Solution Methods for Beam Angle and Fluence Map Optimization in Intensity Modulated Radiation Therapy Planning. OR Spectrum, 30, 289-309.
 Romeijn, H.E., Dempsey, J.F. and Li, J.G. (2004) A Unifying Framework for Multi-Criteria Fluence Map Optimization Models. Physics in Medicine and Biology, 49, 1991-2013.
 Li, Y., Zhang, X. and Mohan, R. (2011) An Efficient Dose Calculation Strategy for Intensity Modulated Proton Therapy. Physics in Medicine and Biology, 56, N71.
 Swope, W.C., et al. (1982) A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. The Journal of Chemical Physics, 76, 637-649.