A preliminary Monte Carlo study for the treatment head of a carbon-ion radiotherapy facility using TOPAS

In medical physics it is desirable to have a Monte Carlo code that is less complex, reliable yet flexible for dose verification, optimization, and component design. TOPAS is a newly developed Monte Carlo simulation tool which combines extensive radiation physics libraries available in Geant4 code, easyto-use geometry and support for visualization. Although TOPAS has been widely tested and verified in simulations of proton therapy, there has been no reported application for carbon ion therapy. To evaluate the feasibility and accuracy of TOPAS simulations for carbon ion therapy, a licensed TOPAS code (version 3_0_p1) was used to carry out a dosimetric study of therapeutic carbon ions. Results of depth dose profile based on different physics models have been obtained and compared with the measurements. It is found that the G4QMD model is at least as accurate as the TOPAS default BIC physics model for carbon ions, but when the energy is increased to relatively high levels such as 400 MeV/u, the G4QMD model shows preferable performance. Also, simulations of special components used in the treatment head at the Institute of Modern Physics facility was conducted to investigate the Spread-Out dose distribution in water. The physical dose in water of SOBP was found to be consistent with the aim of the 6 cm ridge filter.


INTRODUCTION
Over the past 30 years, ion beam radiotherapy using particles like proton, helium, carbon ions has been gradually and slowly increasingly in popularity in radiation oncology [1][2][3].The inverted depth dose profile ⎯the increase of dose along the penetration depth (the Bragg Peak) and a sharp drop after the Bragg Peak ⎯ make ion beams an ideal tool for treatment of the deepseated tumours [1][2][3][4][5][6][7].For carbon ions, the excellent dose profile is further boosted by an additional increase in the relative biological effectiveness (RBE) along the depth till the end of the particle range [2-4, 6, 8].Because of the distinguished physical and biological properties, carbon ion radiotherapy (CIRT) has become a unique research focus around the world.In China, after almost 10 years of research about heavy ion therapy at the Institute of Modern Physics (IMP), Chinese Academy of Sciences (CAS), sophisticated and advanced techniques of beam delivery and treatment planning system have been established [5,6].Since 2006, a total of 213 patients have been treated using carbon ions and most of them are observed with impressive therapeutic effect [5,6].Encouraged by the superior clinical results, two dedicated heavy-ion therapy facilities are now under construction in Lanzhou and Wuwei, China, respectively.Both of them are expected to start treatment in the next few years [6].
Monte Carlo technique has become ubiquitous in medical physics [8].Monte Carlo particle transport codes are practical and useful in calculating the radiation fields and dose distributions.They have also been widely used in hadron therapy as an alternative to fast but less accurate analytical algorithms or to provide data when it is impossible to conduct a physical experiment [9][10][11].A number of studies have shown that Monte Carlo simulation is the most accurate method to calculate doses in radiotherapy [9][10][11][12].Indeed, general purposed Monte Carlo codes such as Geant4, FLUKA, PHITS, and MCNPX are reportedly used for radiotherapy to validate the clinical treatment plans that are mostly calculated by using analytical algorithms [12,13].Clearly, the accuracy of the Monte Carlo codes is of the foremost importance as it will ultimately impact the quality of the treatment plan.In reality, however, the special components of the treatment head used in ion beam radiotherapy are too complicated to be faithfully modelled using those general purpose Monte Carlo codes [12].TOPAS ("TOol for PArticle Simulation") was developed as a new Monte Carlo research platform that is based on the toolkit Geant4 ("GEometry ANd Tracking 4"), and was specially designed for particle radiotherapy simulations, calculations and visualizations [12,14,15].The aim of TOPAS code is to provide an easier-to-use, efficient but still accurate applications for medical physicists [12].Although this tool has been widely tested and verified in simulation of proton therapy, TOPAS utilization in carbon ion therapy has not been reported [12,14,15].
To evaluate the feasibility and accuracy of TOPAS simulations in carbon ion therapy, a licensed TOPAS code (version 3_0_p1) was used to perform dosimetric study for carbon ion therapy.

Structure of TOPAS code
TOPAS is a particle transport simulation code developed with the specific aim of making Monte Carlo simulations user-friendly for users in the particle therapy community.Currently in the 3_0_p1 version, TOPAS is designed as a "user code" layered on top of Geant4, which includes the standard Geant4 toolkit, plus additional codes to access Geant4 functionality [12].TOPAS can model quantities such as passive scattering or scanning beam treatment heads, model a patient geometry based on CT (computed tomography) images, score dose and fluence, save and restart phase space [12].It provides advanced graphics and is fully four-dimensional (4D) views to handle variations in the beam delivery and patient geometry during treatment [12].Initial TOPAS validation results show agreement with measured data at the University of California San Francisco Comprehensive Cancer Center involving the eye treatment head, MGH stereotactic alignment in radiosurgery treatment head and MGH treatment heads, which are simulated in scattering and scanning modalities [12].
The heart of TOPAS is the parameter control system which allows the user to control nearly everything by constructing txt-files without the need of detailed knowledge of Geant4 [12,16].In TOPAS, it is the user input "parameter files" that specify everything: geometry, particle source, fields, motion, scoring, graphical output and physics settings [12,16].Every command line in a txt-file has the same structure: "Parameter-Type: Parameter-Name = Parameter-Value" and the order of command lines does not matter [16].With the correct format of "parameter files", one can run TOPAS as a command-line program with the name of the top level parameter file.That file includes whatever other parameter files are needed.Each parameter file is a simple text file consisting of one or more lines, specifying either an include file or a parameter definition.Each parameter definition line has the same easily mastered format that specifies a parameter type, parameter name, parameter value [16].

Physics model for ion simulations in TOPAS
TOPAS provides a list of default physics models for proton simulation and it has been carefully validated by clinical proton beam measurements at Massachusetts General Hospital [12,14,17].Now for the new version of TOPAS (3_0_p1), the physics models in the default list includs: "g4em-standard_opt4" "g4h-phy_QGSP_BIC_HP" "g4decay" "g4ion-binarycascade" "g4h-elastic_HP" and "g4stopping" [16].This list includes models that handle not only protons but also all secondary particles (neutrons, helium ions, deuterons, tritons, photons, electrons) produced by proton interaction with targets.The default physics list provided by TOPAS is the kind 'modular physics list' in Geant4, which allows that multiple physics 'processes' and 'models' can be assigned as alternatives to the same physical process in the user's application, the actual physics setting being decided upon at run time [17].Unlike Geant4, the users do not have to specify their choice of physics models and settings in a piece of C++ code.There is no single default physics in Geant4 and the variety of physics models and many adjustable settings within each model may lead to the dilemma that various groups in different fields are using their own physics lists which may or may not be properly validated.Of course, users can also further adjust settings from their parameter file, turning off various processes, adjusting the step size or the range cut, etc. Adjustments may be set independently in different geometry components, allowing, for example, more detailed simulation in an ionization chamber [12].Parameters also allow the user to replace the TOPAS physics list with any of the reference physics lists included in Geant4.As a result, TOPAS provides user an accurate but still flexible default physics models for proton therapy simulation in case some users are not expert enough.But can the default physics models be directly used for carbon ion simulation?Maybe corresponding changes and adjustments are needed when it comes to carbon ions with higher energy, since the default physics models have only been verified for therapeutic proton simulation.
The default modular physics list can be mainly characterized by addressing the following types of interactions: (1) electromagnetic process, (2) elastic scattering, (3) inelastic scattering of protons and neutrons and (4) inelastic scattering of heavier ions [17].For example, the standard electromagnetic package is an analytical model that derives directly from QED calculations and describes the interactions of photons and all charged particles down to 1 keV [17,18].The energy loss of hadrons is calculated by the Bethe-Bloch formula down to 2 MeV, below which a parameterization based on ICRU 49 is implemented [17].Scattering is considered by a condensed history algorithm with functions to calculate angular and spatial distributions that match the Lewis theory, validated with an electron scatter benchmark [12,19].After development by years, some of these models have recently undergone major improvements and corrections.For example, according to a benchmark work about hadron therapy both in measurements and Monte Carlo simulation, it is revealed that G4QMD model is preferable over the BIC model for the correct prediction of fragment yields [10].G4QMD and simultaneous inclusion of de-excitation models will result in a better agreement with the experimental data compared with the BIC model [10].While in TOPAS, the intranuclear cascade propagating interactions is addressed by the Geant4 BIC model.As a result, for carbon ion radiotherapy simulation in TOPAS, it still needs to investigate if the BIC model in default physics models should be replaced by the G4QMD model.

Description of the treatment head at IMP
After verification and adjustment in physics models for carbon ion simulation, a study was conducted to model the treatment head which are practically used at IMP treatment terminal.At the IMP carbon ion oncology center, the ion beam is firstly deflected by scanning magnet to form a uniform irradiation field.Then carbon ions will go through the passive scattering modality treatment head in order to produce the desired dose distribution in PTV (Planning Target Volume). Figure 1 shows the simplified layout of the passive scattering treatment head at IMP.The treatment head is a complicated system which usually consists of various kinds of special components such as collimator, range shifter, ridge filter, ion chambers, aperture, etc. [11,20].For the dose calculation and verification, four important components (As seen in figure 1) are practically constructed in this simulation.The primary collimator is a component made of 4 identical copper bars with each thickness of 7.5cm, whose function is blocking the unnecessary particle and controlling the original size of the treatment field.Range shifter, which is constituted by various dimensioned Lucite slices (mass density ȡ= 1.19kg/m3), is used to modulate the depth of the pristine Bragg Peaks.The ridge filter is a stationary plate shaped device made of aluminum typically (mass density ȡ=2.7 kg/m3), and consists of about a total of 40 ridge bars abreast placed every 5 mm in this modelling.The beam energy will be degraded by bring in the ridge filter [21,22].After the ridge filter, the Spread-Out Bragg peak will be produced to cover the target in depth with a relatively homogeneous RBE-weighted dose [2][3][4].The SOBP dose distribution span depends on the period of the ridge filter structure and the angular straggling of the ridge filter placement.Multi-Leaf Collimator (MLC), which is made of the tungsten copper alloy (W 90%, Cu 10% and mass density ȡ= 16.75kg/m3), is a component specially chosen and customized for the individual case to block the unwanted particles.In this study, the MLC arrangement of each slice is consistent with the previous treatment case at IMP.All these special components have been carefully optimized and refined over years at IMP to achieve a high degree of dose conformation to the target volume while protect the OAR (Organs At Risk).

Verification of TOPAS default physics model for proton and carbon ions simulation
Figure 2 shows the pristine Bragg Peaks for monoenergetic proton and carbon ions calculated using TOPAS.Figure 2a is the calculation results of proton beam with energy from 50 to 250 MeV. Figure 2b is the analogous results for carbon ions.Both of the simulation were performed based on the TOPAS default physics models.As can be seen from the Figure 2a, when the primary energy of proton beam is lower, for example less than 150 MeV, the proton Bragg Peak in water is sharp and the dose level is higher at the peak region.But for proton with higher energy of 200 MeV or greater, the Bragg Peak region will be broadened due to severe energy straggling for protons with such high energies.But for carbon ions of 400 MeV/u, the pristine Bragg Peak is still sharp confirming that carbon ions are much heavier than protons and therefore the energy straggling effects are less important.Depth in water(cm) Although the TOPAS default physics models for proton therapy simulation have been testified and validated by many studies [14,15].As discussed earlier in section 2.2, TOPAS will evolve following the Geant4.In TOPAS v3_0_p1, the default physics models have been slightly changed.For example, the electromagnetic model has been modified from opt3 to opt4.For comparison, Figure 3 provides the verification for proton ranges using the PSTAR range data from National Institute of Standard and Technology [23].The mean projected range (or range) of a proton beam, ܴ0, is defined as the penetrating depth in the material where half of the protons have been stopped by undergoing electromagnetic interactions only.It also corresponds to 80% of the maximum dose after the Bragg peak [24].It can be seen from the Figure 3 that the default physics models are fairly accurate for proton calculations with negligible differences (less than 2%).For carbon ions simulations, TOPAS does not provide corresponding physics models, so the precision of result using the TOPAS default physics models still needs further investigation.Furthermore, several studies reported that the QMD model are preferable than the BIC model, while the TOPAS default physics models utilized the BIC model and the QMD model is abandoned [10,16].As a result, it is necessary to compare and determine the more accurate physics model.Figure 4 compares the results obtained by QMD and BIC models.Both calculations are verified by the measurements data from IMP experiment.It was found that both the QMD and BIC models can give highly accurate results, especially for carbon ions with energy less than 300 MeV/u.But when carbon ion energy reaches 400 MeV/u, slight differences are observed before and after the Bragg Peak.As seen from Figure 4 (d), the QMD result curve agrees better with the measurement, suggesting that the QMD model is more accurate than BIC model when carbon ions are at a high energy.Noted that the vertical axis of Figure 4 is the relative dose rather than the absolute dose.

Treatment head modelling and SOBP calculation
As a preliminary study, a simplified model of the IMP treatment head was constructed to verify the spread-out dose distribution in water.But even in a real treatment case, only the ionization chamber, collimator and the compensator are included, which would have a slight effect on the SOBP study.The full modelling of a more realistic treatment head is now under progress.Figure 5 is the visualization for the special components of the treatment head, and the detailed sketches of the ridge filter and multi-leaf collimator are supplied.Figure 6 is the calculation of SOBP by 300 MeV/u carbon ions using the default physics models.Noted that this result is obtained by 300 MeV/u carbon ions getting through a 6 cm ridge filter and 5 cm range shifter, which is also the reason that the SOBP end is not consistent with the pristine Bragg Peak of 300 MeV/u carbon ions.The function of the ridge filter has also been verified in figure 6, as can been seen, the SOBP region spans from nearly 11 cm to 17 cm, which exactly accords with the design aim of spreading the peak for 6 cm.SOBP calculated by TOPAS is the physical dose, the relative biological effectiveness has not been considered.According to the design of the ridge filter, the biological effective dose, which is the product of physical dose and RBE, should be uniform and with the plateau at the SOBP region.Studies related to biological effective dose of carbon ions are still the on-going effort.

Conclusion
To conclude, TOPAS is well-designed to make the simulation for Ion Beam Radiotherapy, both proton and carbon ions, easier and more efficient.The default physics models in TOPAS are elaborately chosen and have been validated with a high accuracy.However, when it comes to the high energetic carbon ions like 400 MeV/u or higher, the G4QMD model is recommended instead of the BIC model.Especially when using the TOPAS in other research field and the object is the heavier ions like neon, iron or uranium, etc., the G4QMD model is the preferable choice.Other outstanding features supported by TOPAS are including the import of XIO format phantom, DICOM files and CAD models.Time related feature and the support of electromagnetic field simulation are two other highlight features in TOPAS.But, still, the calculation capability based completely on CPU is limited even with the parallel run of multi-thread.One of the possible solution is transplanting the computing to a GPU heterogeneous environment where the Monte Carlo method can have a chance to be more efficient and practical for clinical usage [25].This project was supported by the NSFC (Natural Science Foundation of China) with grant (Grant No. 11575180).

Fig. 1 .
Fig. 1.Layout diagram for the simulation of the treatment head (the four components are primary collimator (No.1), range shifter (No.2), ridge filter (No.3) and multi-leaf collimator (No.4), respectively).The photographs are the multi-leaf collimator and range shifter routinely used at IMP.

Fig. 2 .
Fig. 2. Depth dose distribution of proton (a) and carbon ions (b) based on TOPAS default physics models.

Fig. 3 .
Fig. 3. Comparison of proton ranges calculated using TOPAS and NIST data

Fig. 4 .
Fig. 4. Comparison of QMD and BIC model with the measurement for carbon ion from 100 MeV/u to 400 MeV/u.