Sensitivity analysis of Monte Carlo model of a gantry‐mounted passively scattered proton system

Abstract Purpose This study aimed to present guidance on the correlation between treatment nozzle and proton source parameters, and dose distribution of a passive double scattering compact proton therapy unit, known as Mevion S250. Methods All 24 beam options were modeled using the MCNPX MC code. The calculated physical dose for pristine peak, profiles, and spread out Bragg peak (SOBP) were benchmarked with the measured data. Track‐averaged LET (LETt) and dose‐averaged LET (LETd) distributions were also calculated. For the sensitivity investigations, proton beam line parameters including Average Energy (AE), Energy Spread (ES), Spot Size (SS), Beam Angle (BA), Beam Offset (OA), and Second scatter Offset (SO) from central Axis, and also First Scatter (FS) thickness were simulated in different stages to obtain the uncertainty of the derived results on the physical dose and LET distribution in a water phantom. Results For the physical dose distribution, the MCNPX MC model matched measurements data for all the options to within 2 mm and 2% criterion. The Mevion S250 was found to have a LETt between 0.46 and 8.76 keV.μm–1 and a corresponding LETd between 0.84 and 15.91 keV.μm–1. For all the options, the AE and ES had the greatest effect on the resulting depth of pristine peak and peak‐to‐plateau ratio respectively. BA, OA, and SO significantly decreased the flatness and symmetry of the profiles. The LETs were found to be sensitive to the AE, ES, and SS, especially in the peak region. Conclusions This study revealed the importance of considering detailed beam parameters, and identifying those that resulted in large effects on the physical dose distribution and LETs for a compact proton therapy machine.

On the other hand, beside the proton physical dose distribution, its relative biological effectiveness (RBE) should also be taken into account. [3][4][5] However, according to recent publications, 3,6 there is a significant variability in the RBE of protons as a function of depth or particle energy in the beam. In radiation dosimetry, linear energy transfer (LET) is one of the fundamental variables employed to derive the RBE. 7 According to the recently published AAPM TG-256, voxel-by-voxel dose-averaged LET can be employed as a valuable tool for biologically optimized treatment planning even without knowing dose-and tissue endpoint-specific RBE values accurately. 4 Therefore, it is important to provide accurate proton LET distributions with rigorous sensitivity analysis, in addition to the physical dose, for clinical applications.
Monte Carlo (MC) simulation, as a gold standard tool in simulating complex radiation transport, 8 plays an increasingly important role in proton therapy. [9][10][11] Moreover, MC calculated LET values can be efficiently employed in the optimization of proton treatment planning systems (TPS) for clinical applications. 12 Considering these advantages, MC simulation "can be an alternative or complimentary source of dosimetric data for developing, configuring, and validating analytical dose algorithms in clinical TPS". 1 In order to derive physical dose and LET distributions of a proton therapy unit by MC simulations, all major mechanical components of the treatment nozzle should be modeled in detail. 9,13,14 However, even detailed simulation of all machine components cannot account for deviations, especially for the radiation source, from factory specifications. 11,15,16 Furthermore, the source information provided by the manufacturer is often limited to spot size and nominal energy, and customization of the MC model is required to match its results with the measured data.
The large number of adjustable parameters in a clinical proton therapy system (e.g., average energy, energy spread, spot size, beam offset from central axis, etc.) demands a thorough sensitivity analysis that provides important characteristics that are difficult or impossible to measure. In addition, the routine quality assurance processes can be significantly facilitated by correlating the adjustable simulation parameters with measured dose distributions.
In this study, physical dose and LET distributions of a passive double scattering compact proton therapy unit, known as Mevion The aim of this work is to present guidance on the correlation between treatment nozzle and proton source parameters, and physical dose distribution to the following: for researchers modeling clinical proton beam systems, and clinical medical physicists tasked with physically tuning their passive double scattering compact proton therapy unit to bring beam parameters to within clinically acceptable levels. Moreover, this work proposes to create a reference library to troubleshoot of the machine installed and commissioned in the S. Lee Kling Proton Therapy Center at Barnes-Jewish Hospital in St. Louis, MO, USA.
To the best of our knowledge, this is the first reported simulation and sensitivity analysis of the Mevion S250 using MCNPX. The Mevion S250 machine is particularly noteworthy because of its unique beam characterization, which is due to the lack of energy selection and beam transportation systems, and mounting interest in single-room proton unit. 17 The Mevion S250 has 24 different beam options divided into large, deep, and small groups. 2 Each beam configuration is generated using a unique arrangement and combination of different beam line components ( Fig. 1). In order to acquire physical dose distribution data for each configuration, three sets of measurements are taken: (a) pristine peak, (b) lateral profiles in air, and (c) spread out Bragg peak (SOBP).
To obtain pristine peaks and SOBPs, a parallel-plate chamber (PPC05, IBA Dosimetry) was used to measure percent depth-dose curves in a 3D scanning tank (Blue phantom, IBA Dosimetry America, Bartlett, TN, USA) at nominal source-to-surface distance (SSD) of 200 cm with radiation isocenter placed on the water surface. Lateral profiles in air were measured using a diode Edge detector (Sun Nuclear Corporation, Melbourne, FL, USA). The general guideline for acquiring beam data for photon machines has been described in the report of AAPM Task Group 106. 18 For the MC simulations, MCNPX was used in this work. 19 Computations were performed using the facilities at the Washington University, Center for High Performance Computing. The simulation was started by using the manufacturer's specifications for all dimensions and materials of each beam component. Then, average energy, energy spread (FWHM), spot size, and First Scatter (FS) thickness were tuned to match the measured data (experimental results).
For SOBP simulations, a user defined beam current modulation (BCM) sequence was defined considering the rotation angle of the range modulator wheel. 10,20 According to the method described by Polf et al. individual pristine Bragg peaks were created and weighted to form a uniform and flat SOBP with the desired modulation. 10 For the pristine peak and lateral profiles, for each dose point, MCNPX derived data were compared with the measurements by calculating the local difference. In addition, for the penumbra region of the profiles, distance between the 80% and 20% dose levels was compared. SOBP was evaluated by comparing the simulated results and measured data, on the SOBP width as defined by the proximal 95% to the distal 90% dose, beam range, and the depth of distal 20% dose. The beam range was defined as the depth of 90% dose (D90%) on distal fall-off.
The simulations in water and air were performed with 2.0 × 10 9 and 7.5 × 10 8 histories respectively. Generally, for the so-called "good practice" in MC simulations, enough history should be calculated to BARADARAN-GHAHFAROKHI ET AL.
| 27 ensure that MC results have at least 1σ (k = 1, 67% confidence index) <1% statistical uncertainty at depths of interest in the water and air. 19 The mesh tally detector was used due to its functionality for proton dose calculations. 19 The photons, electrons, protons, and positrons were suppressed for simulations with a cutoff energy of 990 eV, 57.3 keV, 5 keV, and 56.6 keV respectively.
For comparison purposes, the results were then compared with the derived values from TOPAS (version 2.0) simulations. 20

2.A.2 | LET
After tuning the treatment nozzle for all the 24 options, both track-averaged LET t and dose-averaged LET d , were calculated according to the method by Guan et al. 7 It should be noted that, LET t calculated as the arithmetic mean value of the fluence spectrum, matches the definition by the ICRU, and LET d is a quantity that accounts for both physical dose and LET, to predict biological effects. 7,19 To calculate both LET t and LET d , the detector cells were modeled as spheres on the central axis of the beam in water phantom.
The MCNPX LET special tally was employed to record flux over the cells as a function of stopping power instead of energy. 19 Using this tally, the recorded values in the energy bins are interpreted as stopping power values (units of MeV/cm).

| SENSITIVITY ANALYSIS
The sensitivity of the model to changes in machine proton source parameters was analyzed by varying the Average Energy (AE), Energy Spread (ES), Spot Size (SS), Beam Angle (BA), and Offset from central Axis (OA). FS thickness was the only treatment nozzle parameter used for the sensitivity investigations. In this step, all 24 options of Mevion S250 were simulated in different stages to obtain the uncertainty of the derived results on the depth of the pristine peak, shape, and symmetry of the resulting dose profiles, as well as LET t and LET d distributions. For the depth of pristine peak, the dose distribution was evaluated using the distal 90% (D90%). Because, each SOBP was created by superimposing single Bragg peaks, the sensitivity study of the pristine peaks also reflects uncertainties associated with SOBPs. 10,13,15 The flatness and symmetry of the profiles were analyzed using the method proposed by Prusator et al. 20 as follow: where, D min and D max are, respectively, the minimum and maximum doses within the central axis of the beam to the 80% dose levels.
where, LD integral and RD integral are the integral doses of the left and right side of the radiation field respectively.
The stages included in the sensitivity study were the AE (±9%),

| RESULTS
With the number of simulated histories used, the uncertainty associated with statistical (quantum) uncertainty in the MC-calculated results in air was less than 0.50% at all distances in the transverse plane. For the large group, the statistical uncertainty for calculations in water was less than 0.61% for the depths of 25.0 cm (D90% of the deepest Bragg Peak). Whereas, for the deep and small groups, it was 0.94% and 0.69% for D90% of 31.9 and 20.0 cm, respectively.
This statistical uncertainty made it feasible to investigate noticeable effects on the physical dose distribution due to slight changes of the sensitivity study parameters. In other words, high precision was obtained in the results of the simulations, which was due to simulating large numbers of histories.   Altering the proton beam energy to 109% showed a strong effect on the small group profile flatness (up to 8.6%), whereas results on the symmetry of the profiles were not significant (less than 2.2%) and remained within the statistical uncertainty of the MC calculations. Figure 5, demonstrates the variations of the flatness and symmetry of lateral profiles for the small group due to 9% increase in the AE. Figure 6, shows the absolute difference between the baseline and sensitivity derived width of the pristine peaks and peak-to-plateau ratios (peak-to-plateau Sensitivity /peak-to-plateau Baseline ) of the studied groups (large, deep, and small). Due to 20% increase in the ES, the maximum increase in the width of the pristine peaks were 4.1 and 2.9 mm for the large and deep groups respectively, whereas, the maximum decrease in the peak-to-plateau ratio were 4.5% and 1.6%, respectively. For the small group, the maximum differences for the width of the pristine peaks and peak-to-plateau ratios were 4.4 mm and 9.5%, respectively.

4.B.3 | Spot size
It was revealed that changes of the SS resulted in no significant effect on the pristine peaks D90% (<1 mm). Figure 7 gives the changes of distal width and peak-to-plateau ratio of the pristine peaks due to uncertainty associated with SS in large, deep, and small groups. It was observed that, considering up to 3 mm uncertainty of the SS, there was no significant effect on the distal width (<1 mm) and peak-to-plateau ratio (<1%) owing to the beam spot size change.
In our study, the distal width of a Brag peak was defined as the difference between the 10% and 90% dose on the distal fall-off.
Noticeable effect on the distal width (about 2 mm) was seen for SS of 10 mm, especially for the small group (Fig. 7).
Therefore, the beam SS have little influence on pristine peaks and the effects might be too subtle to produce significant changes to the distal width of the dose distribution.  On the other hand, a 0.6°deviation from normal was found to change the flatness by 2.7% for small group and 1.8% for deep group ( Table 2). The maximum effect was seen for the small group  significantly decreased flatness and symmetry of the profiles, especially for small option (31.0% and 11.4%). Figure 9 gives variations of the flatness and symmetry of the large group profiles due to uncertainty in the OA. Similar results were observed for the SO. In other words, considering the statistical uncertainty of the MC results, OA was indistinguishable from SO. Figure 10 shows changes in the flatness and symmetry of the profiles due to uncertainty associated with SO from the central axis.

4.B.4 | Beam angle
Moreover, 3-15 mm OA resulted in change of D90% (up to 3 mm). According to our results, small groups was found to be more sensitive on the OA and SO small changes compared to the other treatment nozzle configurations.

4.B.6 | First Scatter thickness
In this work, it was hypnotized that changes in the thickness of the lead FS have a much more prominent effect on the depth of D90%.
The maximum thickness of the FS is 8.370 mm (large group-option 12) to 1.322 mm (deep group-option 17) depending on the treatment nozzle configuration. Therefore, options 1 and 13 with corresponding FS thickness of 6.167 mm and 1.322 mm, respectively, were selected for the sensitivity analysis.
For the small group, there is no FS in the treatment nozzle configuration; therefore, we were unable to investigate the sensitivity of the derived dose distribution due to uncertainty in FS thickness.

4.B.7 | LET
For the small group due to 9% increase in the AE, in the peak region  Figure 12, illustrates changes of the LET t and LET d distributions due to 3% and 5% increase in the AE for the large group.
As ES becomes larger (up to 20%), for all the groups, the absolute maximum of LET d and LET t become lower (up to 1.67 keV/μm −1 and 1.23 keV/μm −1 , respectively), and less steep at the end of the range (Fig. 13). Increasing the ES resulted in decrease in LET t and LET d for the peak regions, especially for the small group (up to 9.2% and 9.   6. Absolute difference between the baseline and sensitivity derived width of the pristine peaks (a) and peak-to-plateau ratios (b) due to change in the ES for the studied groups (large, deep, and small). Peak-to-plateau ratio was derived based on the ratio of the peak-to-plateaus of sensitivity results to the baseline values (peak-toplateau Sensitivity /peak-to-plateau Baseline ).

| DISCUSSION
There is no literature available on the sensitivity analysis of physical dose and LET distributions of Mevion S250 as a passive double scattering compact proton therapy unit.
The MCNPX benchmarked physical dose with the measured data showed up to 2% or 1 mm discrepancy (Table 1; Fig. 2  The increase in distal width and decrease in peak-toplateau ratio of the pristine peaks due to uncertainty associated with SS in large, deep, and small groups. Increase in distal width was calculated based on the absolute difference of the baseline and sensitivity derived results. Peak-to-plateau ratio was derived based on the ratio of the peak-to-plateaus of sensitivity results to the baseline values (peak-to-plateau Sensitivity /peak-to-plateau Baseline ). Large Option Deep Option Small Option F I G . 8. Variations of the peak-to-plateau ratio due to changes of the incident angle of the proton beam. Peak-to-plateau ratio was derived based on the ratio of the peak-to-plateaus of sensitivity results to the baseline values (peak-to-plateau Sensitivity /peak-toplateau Baseline ).
T A B L E 2 Variations of the flatness and symmetry of the profiles due to changes of the incident angle of the proton beam.  of the MCNPX MC code lead to results that are within to the clinically used quality assurance criteria. 24 Results of the MCNPX simulation were in reasonably good agreement with the TOPAS simulations (<2% of the dose points compared and <1 mm for the distances) ( Table 1). Small differences between the MCNPX and TOPAS, especially for the lateral beam profiles, may be due to the shape of the aperture cut-out, which was relevant for penumbra region and was more difficult to model accu- On the other hand, it has been reported that, the beam SS may influence the peak-to-plateau ratio of pristine Bragg peaks, and accordingly the SOBPs uniformity. 15 Some publications have shown that the beam spot size adjustments can be based on the steepness of the distal fall-off and the peak/plateau ratio, both being quite sensitive to this parameter. 10,13 For the Mevion S250 proton system, we observed up to 2 mm increase in the distal width of the pristine peaks due 7 mm increase in the SS, especially for the small group ( Fig. 7).

BA (°)
We found that BA, OA, and SO significantly decreased the flatness and symmetry of the profiles (Figs. 8-10). Up to 2.0°tilt of proton beam resulted in nonsignificant effect on the D90% (<1 mm), while variations in the peak-to-plateau ratio were significant (Fig. 8).
For the Mevion S250, the slope of the lexan and lead layers of second scatter are greater near its central axis. 2  OA required to induce a 2% change in symmetry was roughly 25% greater for the large option than for the small configuration (Table 2).
Moreover, flatness of the profile was more sensitive to SO and OA changes near the central axis than in the periphery of the field, especially for small options. For deep group, the 5 mm OA from the central axis can significantly change the symmetry of the profile (up to 7.5%), while for small group lesser changes of OA (4 mm) was needed to produce a same effect (Figs. 8 and 9). This results in another source of uncertainty that can be included in a reference library guiding clinical medical physicists and engineers to troubleshoot and repair the machine, and also to tune the beam parameters to within clinically acceptable levels. Therefore, the proton beam line and second scatter must be taken into account as a system to evaluate changes in the shape of the profiles. While for a proton system already in clinical use, determining reasonable tolerances in all the moving parts in the beamline can be difficult, some   measurements, i.e. beam profile due to second scatter offset, during acceptance testing of a new system may be helpful for later Monte Carlo commissioning work.
In addition to the uncertainties studied in this work, physics constants of various materials used in the construction of the treatment nozzle may produce another source of uncertainty. 16 Previous studies have shown a significant sensitivity in MC calculated dose distribution on the variations in properties of materials used in passive scattering proton therapy treatment nozzles. 15 Based on the sensitivity study performed on IBA (Louvain-la-Neuve, Belgium), the authors have reported that slight changes in density of the materials of field shaping parts clearly influenced the range and uniformity of dose distribution. 15 These results can be used to improve quality assurance procedure or speed up commissioning process, especially for the commissioning of MC models of clinical passively scattered proton beams.
Commissioning a MC model of a passively scattered proton can be a more rigorous and difficult process than it is for more standard treatment planning softwares. Correlating the many adjustable simulation parameters of the nuzzle with measurable dose distributions can notably facilitate the commissioning process. Moreover, to speed up the quality assurance of this system our results are an effective means of relating nuzzle parameters to clinical measurements.
The RBE of proton therapy is a function of dose, tissue endpoint, and energy deposition characteristics. 4 In this regard, the LET can be used to parameterize the latter for proton beams, by taking into account range uncertainties, for a given dose and biological endpoint. 3,6 It is known that, the LET distribution in a proton beam depends on the range, 25 therefore, uncertainty associated with range may affect LET distribution. In this study, it was hypothesized that changing the beam parameters including; AE, ES, and SS, alters the LET t and LET d distribution for a given proton beam. Generally, it is known that LET t and LET d are less sensitive in the plateau region, however, both LETs are highly sensitive to variations in energy near the Bragg peak when proton energy become low. 3,7 It has been stated that, for the passive-scattering unmodulated monoenergetic proton beams of 250 MeV at the Proton Therapy Center at Houston In our study, considering the sensitivity study of the LET, it was found that in the peak region, depending on the option, both LET t and LET d were sensitive to changes of AE, ES, and SS ( Fig. 12-14).
However, compared to the LET d which showed higher sensitivity in the plateau region, LET t showed higher sensitivity in the peak (up to 1.8%).
According to the previous reports for biological dose calculations, LET d is more appropriate than the LET t at therapeutically relevant dose levels. 3 On the other hand, the proton biology experiments have shown the role of LET in the plateau region for determining cell kill is small. 7 Therefore, in line with previous recommendation of Guan et al. 3 , we recommend the use of LET t in the dose plateau region due to its characteristics of continuous increase along beam path and lower sensitivity to beam uncertainties. It means that a spatially variant switch between the use of LET t and LET d to quantify the LET is recommended for biological studies.
Recently, the idea of adaptive treatment planning by LET painting has started its development in the framework of TPS. 4 Based on this idea, for a passively scattered proton treatment plan, optimization algorithm can attempt to minimize the volume of normal tissues exposed to high LET d , resulting in reducing radiation-induced toxicity.

| CONCLUSION S
This study presents a detailed sensitivity analysis of most important but often poorly specified beam parameters required for simulating a gantry-mounted passively scattered proton system. Our results revealed the importance of these parameters specially those resulted in large effects on the physical dose distribution and/or LETs, i.e., average proton beam energy, initial energy spread, spot size, and offset from the central axis. The findings can be used as a useful tool when quality assurance of this system. Moreover, the sensitivity analysis can also be used to aid machine design for determining reasonable tolerances in all the moving parts in the beamline. The simulation results from the sensitivity analysis can be utilized to construct a reference library to guide troubleshooting and repairing for the machine as well.

ACKNOWLEDG MENT
Monte Carlo computations were performed using the facilities provided by the Center for High Performance Computing at Washington University School of Medicine, partially provided through grant NCRR1S10RR022984-01A1.