PaperlayerView exact revision on arXiv

Source-linked AI summary

A fast GPU-based Monte Carlo simulation of proton transport with detailed modeling of non-elastic interactions

H. Wan Chan Tseung, J. Ma, C. Beltran

arXiv:1409.8336v1physics.med-phphysics.comp-ph

TL;DR

Fast, accurate proton dose calculations matter in radiation therapy, and some applications require detailed nuclear-process modeling. This paper presents a GPU-based proton transport Monte Carlo simulation; its validation showed agreement with Geant4/TOPAS, and it is used in clinical treatment-plan workflows.

  • Problem

    Fast, accurate dose calculations are important in proton therapy, while some applications require detailed modeling of nuclear processes.

    Source: §I · pp. 3–4

  • Method

    A GPU-based proton transport Monte Carlo simulation models non-elastic collisions with any nucleus event by event and also models elastic collisions with nuclei other than hydrogen and oxygen.

    Source: §I · p. 4

  • Results

    98.2% 3-D gamma pass rate at 2%-2 mm was achieved for one head-and-neck plan versus TOPAS; secondary-production cross-sections showed reasonable agreement with Geant4 models and published data.

    Source: A. Non-elastic interactions · p. 15, C. Treatment plan calculations · p. 19

  • Takeaways & Limitations

    The simulation is used clinically for treatment-plan verification and as the dose engine for Monte Carlo-based IMPT optimization.

    Source: §V · p. 23

    Neutrons and de-excitation gammas from nuclear interactions in the nozzle and patient are not propagated, contributing primarily to a small dose deficit.

    Source: C. Treatment plan calculations · p. 20

Abstract

from arXiv · show

Purpose: Very fast Monte Carlo (MC) simulations of proton transport have been implemented recently on GPUs. However, these usually use simplified models for non-elastic (NE) proton-nucleus interactions. Our primary goal is to build a GPU-based proton transport MC with detailed modeling of elastic and NE collisions. Methods: Using CUDA, we implemented GPU kernels for these tasks: (1) Simulation of spots from our scanning nozzle configurations, (2) Proton propagation through CT geometry, considering nuclear elastic scattering, multiple scattering, and energy loss straggling, (3) Modeling of the intranuclear cascade stage of NE interactions, (4) Nuclear evaporation simulation, and (5) Statistical error estimates on the dose. To validate our MC, we performed: (1) Secondary particle yield calculations in NE collisions, (2) Dose calculations in homogeneous phantoms, (3) Re-calculations of head and neck plans from a commercial treatment planning system (TPS), and compared with Geant4.9.6p2/TOPAS. Results: Yields, energy and angular distributions of secondaries from NE collisions on various nuclei agree well with the Geant4 Bertini and Binary cascade models. The 3D-gamma pass rate at 2\%-2 mm for treatment plan simulations is typically 98\%. The net calculation time on a NVIDIA GTX680 card, including all data transfers, is $\sim$20 s for $1\times10^7$ proton histories. Conclusions: Our GPU-based MC is the first of its kind to include a detailed nuclear model to handle NE interactions of protons with any nucleus. Dosimetric calculations are in very good agreement with Geant4/TOPAS. Our MC is being used to perform fast routine clinical QA of pencil-beam based treatment plans, and has also been adopted as the dose engine in a clinically-applicable MC-based IMPT TPS. The detailed nuclear modeling will allow us to perform very fast linear energy transfer and neutron dose estimates on the GPU.

I. INTRODUCTION

GPU proton-transport Monte Carlo can be fast, but earlier implementations simplified or omitted nuclear interactions. This work adds detailed elastic and non-elastic modeling while retaining comparable calculation times.

I. INTRODUCTION

  • The paper’s stated motivation is to make Monte Carlo dose calculations practical for routine clinical QA, which existing CPU-based packages can make computationally prohibitive.A TOPAS calculation with ∼1×10^7 histories in a 150 cm3 tumor volume typically takes around one hour on a 100-node CPU cluster.Source: §I · p. 3
  • The GPU proton-transport methods cited here report 2–30 s per 1 × 10^7 histories, but omit or simplify proton–nucleus interactions, often considering only 16O.The paper identifies this as a major approximation in earlier GPU work.Source: §I · pp. 3–4
  • Detailed nuclear modeling matters for metallic implants, secondary-particle production, neutron-dose estimates, and dose-averaged LET calculations.The paper describes non-elastic interactions as an intranuclear cascade followed by nuclear evaporation.Source: §I · p. 4
  • The proposed class II GPU Monte Carlo models event-by-event non-elastic collisions with any nucleus and elastic collisions beyond hydrogen and oxygen, with very good agreement against Geant4/TOPAS.It reports a net time of ∼20 s for 1 × 10^7 histories on a NVIDIA GTX680, including CPU operations and transfers.Source: §I · p. 4

1. Geometry

The transport model maps CT voxels to material compositions, estimates proton energy loss and scattering, and advances histories through voxelized geometry until a collision or termination condition.

1. Geometry

  • CT voxels from −1000 to 3060 HU are represented by 26 materials containing combinations of up to 12 elements; values above 3060 HU are treated as metallic implant, titanium by default.The HU-to-density conversion follows Schneider et al..Source: §1 · p. 5

2. Stopping powers, multiple scattering, energy straggling

  • Stopping-power and range tables support energy-loss calculations, while the model uses scaled Highland scattering and Vavilov-Gaussian energy straggling.The Highland formula is scaled by 1.07 for comparable angular spread with Geant4.9.6p2; the full Molière alternative is about five times slower.Source: §2 · pp. 5–6

3. Transport mechanics

  • Each step is limited by the next nuclear interaction, voxel boundary, or proton range; histories stop when the proton ranges out, falls below 0.1 MeV, or exits the geometry.A random-hinge algorithm accounts for transverse displacement after multiple-scattering steps.Source: §3 · p. 6

4. Non-elastic nuclear interactions

The model combines intranuclear-cascade and evaporation simulations with coherent nuclear elastic scattering and TOPAS-derived beam-spot phase space. Its scope includes explicit approximations in post-cascade transport and nozzle modeling.

4. Non-elastic nuclear interactions

  • Non-elastic interactions use tabulated cross-sections to select a target, then track nucleons through a two-component Fermi-gas nucleus with Pauli blocking until an allowed collision produces outgoing nucleons and a residual nucleus.The nucleus uses a Woods-Saxon density approximated by 15 constant-density shells; residual nuclei de-excite through particle evaporation.Source: §4 · p. 7
  • Evaporation uses GEM with isotropic emission in the nuclear rest frame, but photon evaporation, pre-equilibrium emission, fission, and Fermi break-up are omitted; neutrons are treated as dosimetrically irrelevant.Only protons from the cascade and evaporation stages are propagated, while heavier charged particles are locally deposited.Source: §4 · p. 7

5. Coherent nuclear elastic interactions

  • Coherent elastic collisions use relativistic kinematics and Geant4 total cross-sections, with scattering angles sampled from Ranft’s differential-cross-section parameterization.The parameterization depends on target mass number and invariant momentum transfer.Source: §5 · pp. 7–8

6. Simulation of proton beam spots from scanning nozzle

  • TOPAS nozzle simulations generate phase-space lists near the nozzle exit for all beam energies and configurations, storing x, y, u, v, and kinetic energy.Central-axis lists can represent off-axis spots because their kinematic distributions are essentially identical; range-shifter spots receive an energy correction.Source: §6 · pp. 8–9

B. GPU implementation

The implementation uses CUDA and weighs two GPU parallelization paradigms: one step per thread and one particle history per thread. The paper chooses the latter while batching non-elastic events to limit divergence.

1. Parallelization approach

  • GPU transport Monte Carlo can parallelize either basic track components across threads or assign each thread a particle history.The first approach can reduce branching, while the second keeps a history’s kinematic variables in registers between steps.Source: §1 · pp. 9–10
  • One history per thread can reduce global-memory transactions, although differing histories cause thread divergence, especially when warps handle multiple particle types or secondaries.The alternative step-per-thread approach requires global reads and writes for kinematic variables at each kernel launch.Source: §1 · p. 10
  • The authors select one history per thread and batch non-elastic events for parallel processing to limit divergence.They cite a neutron-transport comparison in which one step per thread was ∼10 times slower, with memory latency outweighing divergence-free execution.Source: §1 · p. 10

2. Software structure

The simulation uses staged GPU kernels to generate proton spots, transport histories, process nuclear interactions and score dose. It keeps particle data on the GPU across kernel calls and calculates dose statistics after processing the proton groups.

2. Software structure

  • Five kernels handle beam-spot initialization, proton transport, intranuclear cascade, nuclear evaporation and voxel dose statistics.Transport continues until a non-elastic collision or an end condition; secondary protons are processed through the transport and interaction kernels.Source: §2 · p. 10
  • Ten groups of N/10 proton histories are used to estimate dose statistical errors, with each group split into subgroups to fit GPU memory constraints.Subgroup size n is an integer multiple of the number of GPU threads in a block.Source: §2 · p. 11
  • The kernel sequence batches nuclear interactions for cascade and evaporation, repeats the kernels for secondary protons, then computes total and RMS dose per voxel.Simulation inputs are initialized before the primary-group loop; dose statistics are calculated after all groups are processed.Source: §2 · pp. 11–12

3. GPU memory allocation and management

  • Particle structure-of-arrays remain on the GPU between kernel calls, while host-GPU transfers occur mainly during initialization and at simulation end.GPU global memory stores particle stacks and dose tallies; structure-of-arrays support coalesced transactions, and particle data are read and updated once per kernel call.Source: §3 · pp. 12–13

C. Treatment plan calculations

The plan workflow preprocesses TPS exports into inputs for the GPU Monte Carlo, including a reduced CT volume and structure masks. Validation compares nuclear predictions, phantom dose distributions and treatment-plan calculations with Geant4/TOPAS.

C. Treatment plan calculations

  • TPS DICOM plan, structure and image files are converted into Monte Carlo inputs, including a trimmed CT volume, beam and nozzle parameters, and a structure bit mask.The bit mask supports user-directed HU overrides and restricting dose output to selected structures.Source: C. Treatment plan calculations · p. 13

D. Validation

  • Nuclear-model validation compares non-elastic collisions on therapeutically relevant nuclei at 70–230 MeV with Geant4 Bertini and Binary predictions and published measurements when available.Source: D. Validation · p. 14
  • Phantom validation uses 1-mm voxelized water, tissue, bone, calcium and titanium, while treatment-plan validation compares three head-and-neck cases against TOPAS using native CT voxel size.The treatment-plan comparison uses the 3-D gamma index and percentage dose-difference distributions; the authors note gamma may miss systematic biases smaller than its pass criterion.Source: D. Validation · p. 14
  • The percentage dose-difference distributions are assessed alongside statistical-error distributions, with peak shifts and RMS differences indicating biases or other model discrepancies.The comparison includes GPU-versus-TOPAS dose differences and the statistical errors of each simulation.Source: D. Validation · pp. 14–15

III. RESULTS

The reported secondary-production cross-sections generally agree reasonably with Geant4 cascade models and published measurements. The GPU Monte Carlo is reported to run about 50 times faster than Geant4 Bertini for the cited interaction calculation.

A. Non-elastic interactions

  • Secondary proton and neutron differential cross-sections show reasonable agreement with Geant4 Bertini and Binary models and published measurements across the displayed collision cases.The cases include 90 and 113 MeV p–12C, 200 MeV p–16O and 200 MeV p–40Ca collisions at various angles; values are normalized to Cugnon et al. total cross-sections.Source: A. Non-elastic interactions · p. 15, Figure 4 · p. 16
  • The GPU-based Monte Carlo is approximately 50 times faster than the Geant4 Bertini model for the reported interaction calculations.The comparison concerns 2.6 × 10^6 proton–16O and proton–40Ca interactions at 200 MeV.Source: A. Non-elastic interactions · pp. 15–16, Table I · p. 15
  • The illustrated pencil-beam dose comparisons include depth-dose, isodose contours and dose profiles for comparison with Geant4.9.6p2.The top-right figure compares iso-dose contours; the Geant4 contours are dotted and overlap the simulation contours.Source: A. Non-elastic interactions · p. 17, Figure 5 · p. 17

B. Calculations in homogeneous phantoms

Homogeneous-phantom tests show close agreement with Geant4 across beam energies and materials, including titanium; square-field simulations also agree well.

B. Calculations in homogeneous phantoms

C. Treatment plan calculations

Treatment-plan recalculations showed close dosimetric agreement with TOPAS, with 2%-2 mm gamma pass rates near 98%; the GPU implementation was substantially faster in the reported comparison.

C. Treatment plan calculations

  • The 3-D gamma pass rates at 2%-2 mm were 98.2%, 97.8%, and 98.9% across the three head-and-neck plans.For the illustrated plan, TOPAS was the dose-map reference; its DVHs and representative coronal dose wash agreed very well with the GPU MC.Source: C. Treatment plan calculations · pp. 19, 21
  • The GPU dose was on average 0.2% lower than Geant4, attributed primarily to unpropagated neutrons and de-excitation gammas from the nozzle and patient.The authors found no apparent spatial pattern in the dose difference; higher differences occurred in low-statistics voxels.Source: C. Treatment plan calculations · pp. 19–20
  • The estimated model-difference variance was σmodel ∼1%, with δ-ray handling hypothesized as its main contributor.TOPAS propagates δ-rays, whereas the GPU MC deposits electron kinetic energy locally; the TOPAS statistical error was also slightly higher, attributed to differing energy-straggling implementations.Source: C. Treatment plan calculations · p. 20
  • A single K20Xm card was up to ∼200 times faster than TOPAS on a 100-node cluster, with dosimetrically close results.The comparison included three head-and-neck plans; computational times are shown in Table II.Source: C. Treatment plan calculations · p. 21, Table II · p. 21

IV. DISCUSSION

The GPU MC is used for routine treatment-plan QA and as the dose engine in an MC-based IMPT optimization code. The authors report near-real-time QA feedback and estimate a half-hour calculation for a complex plan on a 24-GPU cluster.

A. Current applications of our GPU-based MC

  • Routine QA deployment on a dual K20Xm server makes near-real-time feedback on most commercial-TPS plans achievable.A user interface was being developed, and a Python script allows MC dose to be imported into the TPS for evaluation and display.Source: A. Current applications of our… · pp. 21–22
  • The GPU MC is the dose computation engine in a developed intensity-modulated proton therapy optimization code.The code is being upgraded to a robust MC-based IMPT optimization system.Source: A. Current applications of our… · p. 22
  • A complex three-field head-and-neck plan is estimated to take half an hour on a 24-GPU cluster, supporting the authors' view that MC-based IMPT planning can be clinically applicable on a routine basis.The estimate includes the initial dose map for all relevant beam spots.Source: A. Current applications of our… · p. 22

B. Future applications

The GPU Monte Carlo calculates physical dose and dose-averaged LET simultaneously, with LETd values higher where protons range out. Planned additions include δ-ray production and transport and secondary-neutron propagation; the conclusions also report clinical use of the simulation.

B. Future applications

  • Planned physics-model improvements are δ-ray production and transport and propagation of secondary neutrons in the patient.The authors state that neutron transport will enable fast and accurate neutron-dose calculations.Source: B. Future applications · p. 22
  • The LETd map shows the predicted distribution in the coronal plane corresponding to figure 7.The prediction covers areas receiving at least 10% of the maximum dose.Source: Figure 9 · p. 22, B. Future applications · pp. 22–23
  • LETd calculations in water agree well with TOPAS, and LETd is higher where protons are ranging out.Physical dose and LETd are calculated simultaneously without considerably extending net computational time.Source: B. Future applications · p. 23
  • Current LETd values omit small contributions from proton recoils due to secondary-neutron scattering and from heavy-ion recoils.The neutron-scattering contribution is planned for inclusion after adding a neutron-transport kernel.Source: B. Future applications · p. 23

V. CONCLUSIONS

  • The MC was used clinically for treatment-plan verification and as the dose engine for MC-based IMPT optimization.The conclusions report typical processing time of approximately 20 s per 1 × 10^7 proton histories.Source: §V · p. 23
Loading 1409.8336v1…