Introduction

The availability of proton therapy has dramatically increased in recent years. According to the 2021 statistics of the Particle Therapy Co-Operative Group(PTCOG)^[1], 98 proton therapy facilities are operational worldwide, 19 facilities are under construction, and 4 facilities are being planned. Owing to the sharp falloff in the dose distribution (Bragg Peak), a proton beam can deliver most of the dose to tumors, thus protecting the surrounding healthy tissues. In addition, with the development of new types of beamline systems, such as Huazhong University of Science and Technology Proton Therapy Facility (HUST-PTF)^[2] and treatment techniques, such as flash therapy^[3], proton dose calculation requires much higher accuracy than that in conventional photon radiotherapy^[4]. Most current commercial proton-dose calculation engines used in treatment planning systems are based on analytical algorithms which can rapidly yield accurate dose results for homogeneous tissues. However, the predicted results are inaccurate when highly heterogeneous tissues are considered　^[5]. In addition, proton therapy produces many secondary particles, the information on which is of great importance for analyzing their physical and biological effects　^[6]. The analytical algorithms fail to provide detailed information on secondary fragments.

In contrast, Monte Carlo (MC) algorithms fully simulate particle reactions and model transport geometries, achieving accurate total and fragment dose results. Thus, MC simulations are considered the gold standard in dose calculations, although they are time-consuming. Currently, widely used MC tools, such as Geant4, Gate, Topas, and FLUKA, are based on Central Processing Units (CPUs). However, Graphics Processing Units (GPUs) can integrate more computing units, which can simulate transport processes in parallel　^[7-9] and significantly reduce calculation time. Jia et al. developed the first GPU-based MC proton-dose calculation engine, called gPMC ^[10]. The physics model was adapted from Fippel and Soukop ^[11], where proton transport was simulated in a class II condensed history scheme with a continuous slowing down of approximations. Ionization and elastic and inelastic interactions between the protons and materials were included. Chan et al. ^[12] applied a more precise nonelastic nuclear reaction model that incorporated the Bertini cascade and evaporation kernels by sacrificing speed. Both models are CUDA-based, and their applications are limited to Nvidia GPUs. New versions of gPMC, gPMC v2.0, and goMC were realized in the OpenCL framework ^[13,14], facilitating cross-platform applications. Although OpenCL supports multiple GPU hardware, it is mainly developed using an Application Programming Interfaces (API) with a specific hardware development interface.

In 2019, Hu et al. proposed a new parallel programming language called Taichi for high-performance numerical computations ^[15-17]. It supports most mainstream GPUs, and thus works well in cross-platform environments. Taichi is embedded in Python and uses just-in-time (JIT) compiler frameworks to offload computing-intensive Python code to native GPU or CPU instructions. Using the @ti.kernel decorator, Taichi's JIT compiler automatically compiles Python functions into an efficient GPU or CPU machine code for parallel execution. These features enable easy programming that significantly reduces the total number of lines of code compared to OpenCL and CUDA. Moreover, Taichi can be well integrated into mainstream Python deep learning ecosystems such as PyTorch, which integrate GPU high performance computing with deep learning neural networks　^[18-20]. Consequently, Taichi has been widely used in high-performance physics studies ^[16] and machine learning fields ^[21]. However, Taichi has not yet been applied to dose calculations in radiotherapy. In this work, we developed a new cross-platform GPU-based MC proton dose calculation engine for the Taichi platform, named Taichi Proton Monte Carlo (TPMC), with a model of proton transport similar to those in ^{[11, 14]}, except for a different inelastic nuclear reaction model.

The strategy for GPU acceleration in Taichi is similar to those in OpenCL and CUDA. Primary and secondary stacks were built to store primary protons and secondary particles, respectively. The GPU calculation threads first processed the primary protons to start transport and then saved all the secondary particles produced into the secondary stacks. After all the primary protons stopped or escaped from the region of interest, the secondary protons were transported on the GPU threads again and treated as primary protons, except that the information of the produced secondaries was not documented. Using the Taichi platform reduced the number of lines of code to only a few hundred. The computation time was reduced by allocating the deposited dose to different local dose counters to avoid the atomic addition effect. In our study, we prepared 110-, 160-, and 200 MeV proton pencil beams and four water and human tissue phantoms for dose calculation. The gamma passing rates, calculation speeds, and total number of lines of code were evaluated and compared with those obtained from Geant4-TOPAS and other GPU-based simulations.

The remainder of this study is organized as follows. In Section 2, we describe the physical process of the TPMC, parallel simulation strategy, and Taichi GPU memory arrangement. The results of the dose calculations for Taichi and the other platforms are compared in Section 3. Finally, in Section 4, we conclude our study and discuss its current and future applications.

Methods

In this study, we employed the class II condensed history simulation scheme ^[22] and continuous slowing down approximations for proton transport. A simplified physics model was established to simulate the nuclear reactions of the primary and secondary protons. The physical process of proton transport was almost the same as those in Refs ^[10,11], which will be briefly introduced in the following section. On the other hand, the proton transport was accelerated in the Taichi framework on GPUs, where parallel threads were applied and a new dose counter model was developed to solve thread-block problems.

2.2

GPU acceleration

As mentioned in Sect. 1, the implementation of CUDA is limited to the GPU of NVIDIA. Although OpenCL can solve this problem ^[13], it requires an additional application programming interface and lacks functional libraries, thereby creating higher programming barriers for medical physicists and developers. However, Taichi can overcome the drawbacks of these two languages ^[15]. It is based on Python and can therefore conveniently use available packages such as random, numpy, and pandas. When applied to GPU acceleration, an easy outermost scope for-loop command automatically parallelizes work items on different threads in the Taichi kernel, whereas in CUDA and OpenCL, appropriate GPU machine codes must be designed ^[15]. Moreover, switching from Nvidia cards to other GPU cards can be achieved simply by changing only the head files in Taichi, instead of (almost completely) rewriting the code in OpenCL. Finally, the compilation time and memory space are saved in Taichi because it is performed on the LLVM-based compiler Clang ^[15], and Taichi also tunes the parameters to best explore the GPU architecture. Table 1 briefly compares these three languages.

Table 1Full-screen View

Development environment and bases of different platforms.

Platform	Base	GPU	Library
CUDA	C++	Nvidia	Nvidia
OpenCL	API	Multiple	OpenCL
Taichi	Python	Multiple	Python

Show more

The cross-sectional data of proton-nucleus reactions and the density and stopping power distribution of the target were first loaded by the CPU and stored in the GPU global memory, which was called constantly during the simulations. The transport of primary protons was randomly simulated stepwise on different GPU threads until they stopped or moved out of the target. The automatic parallelization of particles on a GPU can be achieved using @ti.for-loops command in TPMC. Secondary particles were generated after inelastic nuclear interactions, and their information, including the location, energy, and direction, was recorded in the GPU global memory and pushed into the secondary particle stack. The transport of secondary particles was blocked until all the primary protons were simulated. Subsequently, we randomly allocated secondary protons to different threads for transport.

Although the particles were transported in parallel on the GPUs, the dose deposition process in each thread was not independent. Individual dose deposition on the global dose counter causes a thread block, where the dose information of other threads must wait until the current process of proton information transfer ends, known as atomic addition. This treatment was the primary reason for the time-consuming nature of the GPU-based MC dose engine. In this study, we applied the same strategy as in gPMC V2.0^[14] to mitigate the memory conflict problem. Eight local dose counters were prepared to balance the calculation efficiency and memory size. The deposited energy was randomly assigned to a dose counter after each step, and the total doses were summed for the eight counters after all the particle transports were completed.

The GPU implementation of the aforementioned process is schematically illustrated in Fig. 1. Protons were simulated in batches with a batch size of $N_{\text{batch}}=65536$ protons in all cases, except for the last batch in our simulation. After $N_{\text{p}}$ primary protons were loaded into a batch, the transport simulation was initiated. It is noteworthy that $N_{\text{p}}$ is equal to $N_{\text{batch}}$. In each batch, the dose information of $N_{\text{p}}$ primary protons, $N_{\text{s}}$ secondary protons, other secondary fragments, and all tertiary particles was calculated following the physics model and saved in the local dose counter $D_{\text{r}} \left(r=1-8\right)$. After all the particle transports were performed in a batch, the dose information in $D_{\text{r}}$ was sent to the CPU memory and added to the total dose, $D_{\text{total}}=D_{\text{total}}+{\displaystyle \sum }_r{D}_{\text{r}}$. The simulation process ends when the number of simulated primary protons $N_{\text{pri}}$ is equal to the number of total primary protons $N_{\text{total}}$.

Fig. 1Full-screen View

Schematic of the GPU acceleration of TPMC

In this study, all TPMC GPU dose calculation processes were performed on an NVIDIA Tesla V100 PCle 32GB card system with 5120 processor cores. The operating frequency was 1230 MHz with 14 TFLOPS in single-precision mode. All the CPU procedures were performed on an Intel Xeon 8260 CPU with 24 processor cores system, processor base frequency of 2.4 GHz, and 1 TB of maximum memory.

2.3

Test model

We used a Gaussian pencil beam with a spot size of 5 mm and scattering angle of $0^{\circ }$. The number of primary protons simulated was 10⁷ and the beam energies were 110, 160, and 200 MeV. Four phantoms of different levels of heterogeneity, containing water, lung tissue, and bone tissue were designed for testing, as shown in Figure 2. The size of the phantom was 10 cm × 10 cm × 30 cm, with a 1 mm × 1 mm×1.5 mm scoring grid. The densities and elemental compositions of the materials are listed in Table 2.

Table 2Full-screen View

Densities and elemental compositions of the materials in the testing phantoms. The densities and composition values are in units of g/cm³ and mass weighted percentage, respectively.

Material	Density	H	C	N	O	Na	P	S	Cl	K	Mg	Ca
Water	1.0	11.2			88.8
Lung	0.3	10.3	10.5	3.1	74.9	0.2	0.2	0.3	0.3	0.2
Bone	1.9	13.5	16.0	4.2	44.5	0.3	9.5	0.3			0.2	21.5

Show more

Fig. 2Full-screen View

Size, shape, and materials of the four testing phantoms.

Results

3.1

Dose validation

The physics setting of TOPAS included the following Geant4 modules: g4em-standard_opt3, g4hphy_QGSP_BIC_HP, g4decay, g4ion-binarycascade, g4helastic_HP, and g4q-stopping, following Testa et al.[27].

3.1.1

Depth dose curve

Figure 3 shows the depth-dose curves obtained for TPMC and TOPAS in the pure water phantom (phantom in Figure 2). The dose values for beams with different energies are normalized by the maximum dose (Bragg peak) of 110 MeV. The dose differences between TPMC and TOPAS mainly occur around the Bragg peak regions.

Fig. 3Full-screen View

(Color online) Depth dose distribution (top) and relative dose difference (bottom) for the pure water phantom for monoenergetic proton pencil beams of 110 MeV (black peak at 91 mm), 160 MeV (red peak at 175 mm), and 200 MeV (blue peak at 258 mm), calculated by TPMC and TOPAS. The dose value for different energy beams is normalized with respect to the maximum dose (Bragg peak) of 110 MeV.

3.1.2

Lateral dose verification

Lateral dose verification was performed by comparing the results for the four phantoms calculated by the TPMC and TOPAS at the same incident proton energy, E = 160 MeV. The relative error maps of the dose distributions are shown in Fig. 4. In most of the regions with energy deposition, the results of TPMC are consistent with those of TOPAS. These differences mainly appear in the Bragg peak areas for phantoms (a), (b), and (c). Interestingly, for the pure water phantom (a), our predicted dose values are larger than those of TOPAS at energies before the Bragg peak, and TPMC decreases faster than TOPAS after the Bragg peak. These differences increase further when lung tissue is inserted into the proton beam pathway. However, the opposite results are obtained when the inserts were replaced with bone tissue, as shown in panels (b) and (c) of Figure 4. For the heterogeneous phantom, large differences appear along the extension line of the boundary between the two inserted materials, apart from the two Bragg peaks. Despite this difference, the Bragg peak regions of TPMC and TOPAS are consistent in all the test cases.

Fig. 4Full-screen View

(Color online) Lateral dose difference of TPMC and TOPAS in the four test phantoms for a 160 MeV monoenergetic proton pencil beam: (a) difference in water phantom, (b) difference in lung phantom, (c) difference in bone phantom, and (d) difference in heterogeneous phantom.

The influence of the heterogeneous medium was further studied by simulating two additional incident proton energies, 110 MeV and 200 MeV, for phantom (d). The lateral dose contour maps of the TPMC and TOPAS in the three cases, the corresponding lateral profiles at a depth of 6 cm (material boundary), and the two Bragg peaks are shown in Figures 5 (a-c). The results show that the lateral dose distribution in TPMC is consistent with that in TOPAS, regardless of the proton beam energy, phantom depth, or material difference.

Fig. 5Full-screen View

(Color online) Lateral dose difference of TPMC and TOPAS in the heterogeneous phantoms for: (a) 110, (b) 160, and (c) 200 MeV mono-energetic proton pencil beams.

3.1.3

Quantitative analysis

To evaluate the accuracy of the TPMC quantitatively, we calculated the gamma passing rates for the four testing phantoms with three incident proton beam energies. The criteria were set to 1 mm/1% and 2 mm/2%. A mask was added to each phantom whose height and width were set to 5 cm and 30 cm, respectively. The depth of the mask was 1.1 times the depth of the Bragg peak. A deeper Bragg peak was selected for the heterogeneous phantom. The gamma passing rates in various situations are listed in Table 3. We note that all the values are greater than 0.99 and 0.95 under the criteria of 2 mm/2% and 1 mm/1%, respectively. These quantitative results indicate that the TPMC is sufficiently precise for proton treatment planning. In addition, the gamma passing rates of the TPMC gradually decrease with increasing proton beam energy and complexity of the phantom.

Table 3Full-screen View

Masked Gamma passing rates of TPMC for four testing phantoms with 110, 160, and 200 MeV proton beams.

Energy	Phantom	Criteria
		1 mm/1%	2 mm/2%
200 MeV	Water	99.2	99.9
	Water + Lung	99.0	99.9
	Water + Bone	98.0	99.9
	Heterogeneous	96.6	99.3
160 MeV	Water	99.1	99.9
	Water + Lung	98.7	99.7
	Water + Bone	98.1	99.9
	Heterogeneous	97.5	99.5
110 MeV	Water	99.4	99.9
	Water + Lung	99.6	99.9
	Water + Bone	98.3	99.6
	Heterogeneous	98.1	99.7

Show more

3.1.4

CT Image Validation

The accuracy benchmark of TPMC against TOPAS was then performed on a CT-acquired patient phantom. A square beam with an energy of 110 MeV was incident. The dose distributions are shown in Fig. 6. The TPMC dose calculation results for the horizontal and inclined injection beams are shown in Figs. 6a and 6b, and the TOPAS results are plotted in Figs. 6c and 6d. Figures 6e and 6f depict the dose differences between the two MC results. For regions with doses above 50% of the maximum dose, the average dose difference is less than 1.5%. We also performed a gamma-index test (2 mm/2%), and the passing rate exceeds 99%. Focusing on the dose points within the 10% isodose lines, the passing rate is greater than 96% for both cases.

Fig. 6Full-screen View

(Color online) Dose calculation result in the case of the patient. The first two rows represent the dose results from TPMC and TOPAS respectively. The third row shows the relative dose difference between the two methods

Discussion

To ensure the accuracy of our physics model, it is worth noting that the maximum depth-dose difference in the pure water phantom occurred around the Bragg peak region. This is because nuclear reactions play a dominant role in the production of many short-range secondary fragments. However, the energy of these fragments was assumed to be locally deposited without additional interactions in the TPMC. This difference increases with increasing incident proton energy, which can be understood in the same way that the higher the projectiles energy, the more the fragments produced. In addition, the lateral dose difference mainly originates from the approximate treatment of the delta electrons, whose influence on the proton transport directions is neglected in our physics model. The accuracy of TPMC decreased when other materials, such as lung or bone tissues, were inserted into water. The transport of delta electrons is probably non-negligible in low density materials for lung tissues ^[10], while for bone tissues, proton-calcium inelastic nuclear interactions should be included. Moreover, the phantom with the heterogeneous medium exhibited the lowest gamma passing rate. The cross section of protons and materials might dramatically change when protons pass through the heterogeneous interface, but we only consider the material information of the current voxel in that case. These results indicate that the accuracy of TPMC may decrease further when the structure of the human body is more complicated.

The lateral dose differences of TPMC and TOPAS in Fig. 4 exhibit opposite results when lung or bone tissues were inserted compared to those of the pure water phantom. This was attributed to the material conversion method used to calculate the water-equivalent steps in Eq. (2). For lung and bone tissues, although energy loss is assumed to be conserved in such conversion, the value of $f_{\text{s}}\left(\text{ρ},T_{\text{p}}\right)$ is not absolutely accurate, since it is obtained by fitting to experimental data [11]. Thus, the energy of the protons after passing through the inserts (other than water) changed compared with the real situation, leading to a shift of approximately 1 mm in the Bragg peak position.

Although both CUDA and OpenCL can achieve better computation performance through numerous well developed functions, the implementation of Taichi in proton dose calculations is still appealing because it supports multiple GPU devices and has a low programming barrier. Taichi also exhibits advantages in other aspects. In proton transport, Taichi continues to process all the particles on the GPUs, whereas OpenCL sends the information of the generated secondary protons from the GPU to the CPU memory after the simulation of each batch finishes. This information is returned to the GPU during secondary proton transport. Such cross-transportation may lead to a relatively low efficiency, as shown in Table 4, where the simulation time of OpenCL in each case is slightly longer than that of Taichi. In addition, Taichi supports most Python function libraries such as Pytorch and Numpy, which naturally integrate dose calculations within the framework of deep learning. This feature provides a new aspect for future work, namely, using convolution neural networks to improve the accuracy of the present study. Similar studies were conducted by Ryan Neph et al. ^[20], who used 3D Unet and high-noise dose maps (10⁵ incident particles) to predict low-noise dose distributions (10⁷ incident particles) in treatment planning. Another work by Wu et al. ^[28] improved the accuracy of the pencil beam algorithm using Unet. Both studies yielded promising results, indicating that deep learning has great potential for dose prediction. Finally, the great power of Taichi in gradient descent optimization can be applied to other aspects, such as image processing and dose optimization. Our goal is to develop a Taichi-based treatment planning system for proton therapy in the future.

Conclusion

We developed a graphics processing unit (GPU)-based Monte Carlo proton dose calculation engine within the Taichi framework. Protons are transported in a class II condensed history simulation scheme with various post step proton-nucleus interactions. The simulations were parallelized on GPUs for acceleration. The results of the testing models indicated that the accuracy satisfied the clinical requirements adequately. The simulation speed was approximately 7000 times faster than that of the full MC simulation using TOPAS, and the number of lines of code was approximately 10 times lesser than those of CUDA and OpenCL. Our study provides a fast, accurate, and easy-to-use proton dose calculation engine for algorithm developers, students, and medical physicists. Further studies are needed to implement Taichi in other components of treatment planning, such as image processing and dose optimization.