Virtual particle Monte Carlo: A new concept to avoid simulating secondary particles in proton therapy dose calculation

Abstract

Background: In proton therapy dose calculation, Monte Carlo (MC) simulations are superior in accuracy but more time consuming, compared to analytical calculations. Graphic processing units (GPUs) are effective in accelerating MC simulations but may suffer thread divergence and racing condition in GPU threads that degrades the computing performance due to the generation of secondary particles during nuclear reactions.

Purpose: A novel concept of virtual particle (VP) MC (VPMC) is proposed to avoid simulating secondary particles in GPU-accelerated proton MC dose calculation and take full advantage of the computing power of GPU.

Methods: Neutrons and gamma rays were ignored as escaping from the human body; doses of electrons, heavy ions, and nuclear fragments were locally deposited; the tracks of deuterons were converted into tracks of protons. These particles, together with primary and secondary protons, are considered to be the realistic particles. Histories of primary and secondary protons were replaced by histories of multiple VPs. Each VP corresponded to one proton (either primary or secondary). A continuous-slowing-down-approximation model, an ionization model, and a large angle scattering event model corresponding to nuclear interactions were developed for VPs by generating probability distribution functions (PDFs) based on simulation results of realistic particles using MCsquare. For efficient calculations, these PDFs were stored in the Compute Unified Device Architecture textures. VPMC was benchmarked with TOPAS and MCsquare in phantoms and with MCsquare in 13 representative patient geometries. Comparisons between the VPMC calculated dose and dose measured in water during patient-specific quality assurance (PSQA) of the selected 13 patients were also carried out. Gamma analysis was used to compare the doses derived from different methods and calculation efficiencies were also compared.

Results: Integrated depth dose and lateral dose profiles in both homogeneous and inhomogeneous phantoms all matched well among VPMC, TOPAS, and MCsquare calculations. The 3D-3D gamma passing rates with a criterion of 2%/2 mm and a threshold of 10% was 98.49% between MCsquare and TOPAS and 98.31% between VPMC and TOPAS in homogeneous phantoms, and 99.18% between MCsquare and TOPAS and 98.49% between VPMC and TOPAS in inhomogeneous phantoms, respectively. In patient geometries, the 3D-3D gamma passing rates with 2%/2 mm/10% between dose distributions from VPMC and MCsquare were 98.56 ± 1.09% in patient geometries. The 2D-3D gamma analysis with 3%/2 mm/10% between the VPMC calculated dose distributions and the 2D measured planar dose distributions during PSQA was 98.91 ± 0.88%. VPMC calculation was highly efficient and took 2.84 ± 2.44 s to finish for the selected 13 patients running on four NVIDIA Ampere GPUs in patient geometries.

Conclusion: VPMC was found to achieve high accuracy and efficiency in proton therapy dose calculation.

Keywords: GPU acceleration; Monte Carlo; intensity-modulated proton therapy; real-time adaptive treatment planning; secondary particles.

Figure 1.

Diagram of converting the track histories of a realistic proton and its secondaries in a conventional Monte Carlo simulation into two virtual particles (a) Track history of a realistic proton. (b) the doses of electrons, heavy ions, and nuclear fragments were locally deposited. (c) neutrons and gamma rays were considered to escape from human bodies and thus ignored. (d) Track history of two virtual particles to have the same path from the starting point to R2. These two VPs are independent from each other. This figure is only used to demonstrate the concept of VP.

Figure 2.

Workflow of VPMC. The upper box shows the loading process of the pre-calculated databases of the CSDA model in range shifter, the CSDA model in patient geometry, and the LAE model. The lower box shows the simulation process of VPs. The steps in the green box indicates the simulations in range shifter, while the steps in the blue box indicates the simulations in patient geometry.

Figure 3.

Illustration of the homogenous phantom (a) and inhomogeneous phantom. The beam enters the phantoms at the center of the y-z plane. The larger white boxes are composed of the water (HU = 0), the smaller grey boxes are composed of the blocking material (HU = 1000).

Figure 4.

The IDD curves (a)(e), log-scale lateral profiles at the Bragg peak (c)(g) of the dose distribution generated by TOPAS (magenta), MCsquare (blue), and VPMC (orange), and the corresponding IDD (b)(f) and lateral profile (d)(h) differences between MCsquare and TOPAS (blue) and the differences between VPMC and TOPAS (orange), in a homogeneous phantom (top row) and in an inhomogeneous phantom (bottom row). RBE = 1.1.

Figure 5.

Comparison of dose profiles on typical transverse planes of a prostate patient without range shifter between VPMC and MCsquare. (a) Dose map from MCsquare, (b) Dose map from VPMC, (c) Absolute dose difference map between the MCsquare calculated dose and the VPMC calculated dose. The white arrow indicates the position and direction of the dose profiles in (c). The red curve in (d) is the dose profile from MCsquare, while the blue curve is the one from VPMC. RBE = 1.1.

Figure 6.

Comparison of dose profiles on typical transverse planes of a chest wall patient with ERS between VPMC and MCsquare. (a) Dose map from MCsquare, (b) Dose map from VPMC, (c) Absolute dose difference map between the MCsquare calculated dose and the VPMC calculated dose. The white arrows indicate the position and directions of the dose profiles in (d) dose profile in X direction and (e) dose profile in Y direction. The red curves in (d) and (e) are the dose profile from MCsquare, while the blue curves are the ones from VPMC. RBE = 1.1.

Figure 7.

Comparison of dose profiles on typical transverse planes of a H&N patient with RS between VPMC and MCsquare. (a) Dose map from MCsquare, (b) Dose map from VPMC, (c) Absolute dose difference map between the MCsquare calculated dose and the VPMC calculated dose. The white arrows indicate the position and directions of the dose profile in (d) dose profile in X direction and (e) dose profile in Y direction. The red curves in (d) and (e) are the dose profile from MCsquare, while the blue curves are the ones from VPMC. RBE = 1.1.

Figure 8.

Comparisons of the measured 2D plane dose during PSQA (a) with VPMC calculation result (c), and Eclipse™ calculation result (g) at a depth of 4.0 cm for a lung cancer patient without range shifter. The corresponding 2D-3D Gamma analysis pass/fail maps are shown in (d) and (h) with 3%/2mm/10%. Subpanel (b) displays the dose profiles from VPMC and Eclipse™ in the beam direction. The black line is from the Eclipse™ calculation result, while the green line is from the VPMC calculation result. Red points are the measured results with an error bar of 2%/2mm. Subpanel (e) (f) displays the dose profile comparison between the VPMC calculated dose and the measured dose in the X direction at the Y position indicated by the horizonal line in (c) and in the Y direction at the X position indicated by the vertical line in (c), respectively. Subpanel (i) (j) displays the dose profile comparison between the Eclipse™ calculated dose and measured dose in the X direction at the Y position indicated by the horizonal line in (g) and in the Y direction at the X position indicated by the vertical line in (g), respectively. RBE = 1.1.

Figure 9.

Comparisons of the measured 2D plane dose during PSQA (a) with VPMC calculation results (c) and Eclipse™ calculation results (g) at a depth of 2.0 cm for a chest wall patient with RS. The corresponding 2D-3D Gamma analysis pass/fail maps are shown in (d) and (h) with 3%/2mm/10%. Subpanel (b) displays the dose profiles from VPMC and Eclipse™ in the beam direction. The black line is from the Eclipse™ calculation result, while the green line is from the VPMC calculation result. Red points are the measured results with an error bar of 2%/2mm. Subpanel (e) (f) displays the dose profile comparison between the VPMC calculated dose and the measured dose in the X direction at the Y position indicated by the horizonal line in (c) and in the Y direction at the X position indicated by the vertical line in (c), respectively. Subpanel (i) (j) displays the dose profile comparison between the Eclipse™ calculated dose and measured dose in the X direction at the Y position indicated by the horizonal line in (g) and in the Y direction at the X position indicated by the vertical line in (g), respectively. RBE = 1.1.

Figure 10.

Comparisons of the measured 2D plane dose during PSQA (a) with VPMC calculation results (c) and Eclipse™ calculation results (g) at a depth of 15.0 cm for a H&N patient with ERS. The corresponding 2D-3D Gamma analysis results were shown in (d) and (h) with 3%/2mm/10%. Subpanel (b) displays the dose profiles from VPMC and Eclipse™ in the beam direction. The black line is from the Eclipse™ calculation result, while the green line is from the VPMC calculation result. Red points are the measured results with an error bar of 2%/2mm. Subpanel (e) (f) displays the dose profile comparison between the VPMC calculated dose and the measured dose in the X direction at the Y position indicated by the horizonal line in (c) and in the Y direction at the X position indicated by the vertical line in (c), respectively. Subpanel (i) (j) displays the dose profile comparison between the Eclipse™ calculated dose and measured dose in the X direction at the Y position indicated by the horizonal line in (g) and in the Y direction at the X position indicated by the vertical line in (g), respectively. RBE = 1.1.