A novel approach to seamless simulations of compact hadron therapy systems for self-consistent evaluation of dosimetric and radiation protection quantities

Hadron therapy installations are evolving towards more compact systems that require higher-quality beams for advanced treatment modalities such as proton flash and arc therapy. Therefore the accurate modelling of present and next-generation systems poses new challenges where the simulations require both magnetic beam transport and particle-matter interactions. We present a novel approach to building simulations of beam delivery systems at a level suitable for clinical applications while seamlessly providing the computation of quantities relevant for beam dose deposition, radiation protection assessment, and shielding activation determination. A realistic model of the Ion Beam Applications (IBA) Proteus® One system is developed using Beam Delivery Simulation (BDSIM), based on Geant4, that uniquely allows simulation using a single model. Its validation against measured data is discussed in detail. The first results of self-consistent simulations for beam delivery and equivalent ambient dose are presented. The results show that our approach successfully models the complex interactions between the beam transport and its interactions with the system for relevant clinical scenarios at an acceptable computational cost.

Introduction. -The growing acceptance of proton therapy and the increasing number of centers in clinical operation worldwide are driving the development of novel treatment modalities [1], such as spot-scanning proton arc therapy [2] and flash therapy [3,4] to further boost the clinical efficiency and reduce the treatment time. Those new modalities pose new challenges for the accelerator and the beam delivery system as they require a tighter control of the beam spot sizes at the isocenter, an improved control for the beam loss levels and higher beam currents. Furthermore, single-room compact centers with (a) E-mail: cedric.hernalsteens@ulb.be a small building footprint now represent a significant fraction of the newly built systems [5]. In particular, the IBA Proteus R One is a single-room compact system that was first used to treat patients in 2016. There are now 12 clinically operating installations worldwide. These advancements pose challenges for the development of numerical models able to simulate proton therapy systems globally both in terms of beam properties and in terms of secondary radiation fluences. The prediction of beam properties with a clinically relevant accuracy enables the evaluation of new beam delivery schemes from numerical simulations, in conjunction with Treatment Planning Systems (TPS). At the same time, the accurate modeling 50004-p1 of secondary radiation fluence allows computing radiation protection and shielding activation quantities. Simultaneously considering these two aspects is of particular importance when evaluating new beam delivery methods in view of upgrading existing centers or when preparing specialpurpose beams for experimental irradiation setups. It is also crucial for the conception of next-generation compact systems, where the beam characteristics, beam losses, and radiation shields are tightly coupled due to the reduced size of the installation [5]. Numerical modeling of proton therapy installations has so far used tools that can be split into three categories: Monte Carlo simulations for the beam delivery system "nozzle" (in particular for scattering techniques [6,7]), beam transport codes for simulations of the beam transfer line [8], and Monte Carlo simulations for the evaluation of the shielding requirements [9]. As active scanning beam delivery methods have become prevalent, Monte Carlo models, including the magnetic beam steering in the nozzle, have been developed and validated with experimental data [10].
In this paper, we propose a novel approach to the seamless simulation of both the beam transport system and the beam-matter interactions in the beamline and treatment nozzle, as regards the evaluation of the clinical beam properties, dosimetry, and radiation protection quantities, such as residual ambient dose and concrete shielding activation. We highlight the need and advantages of building a complete model allowing the self-consistent computation of these quantities for next-generation systems. It has been recognized [11,12] that beam transport codes that do not include detailed beam-matter interactions models (e.g., Transport [8] or MAD-X [13]) are not able to satisfactorily reproduce key beam properties at a precision compatible with the clinical requirements, especially where the typical beam size is highly restricted by aperture and collimators and can be non-Gaussian. Codes have been developed with limited particle-matter interaction capability in discrete places to model beam lines such as those for the proton therapy facility at the Paul Scherrer Institute (PSI) [14]. However, these use simplified physics models in select places. Monte Carlo codes dedicated to nozzle simulations such as TOPAS [15] are unable to simulate the full accelerator as required.
Our model of the IBA Proteus R One is built using Beam Delivery Simulation (BDSIM) [16], a Monte Carlo program that simulates the passage of particles through detailed 3D models of a beamline. Using the Geant4 [17,18], ROOT [19], and CLHEP [20] libraries, BDSIM programmatically builds geometries for the beamline components complete with electromagnetic fields, allowing transport through the entire beamline whilst benefiting from the physics models of a general-purpose Monte Carlo toolkit. The particles are propagated using a combination of accelerator tracking routines and common numerical integrators provided with Geant4. The magnetic tracking of BDSIM has been extensively benchmarked against MAD-X and other particle tracking codes [16]. Highly flexible, BDSIM models can include customised field maps, magnet pole-face rotations, detailed fringe field models and geometry allowing the inclusion of bespoke magnetic beamline elements crucial for accurate tracking, as well as permitting the inclusion of concrete shielding and non-standard machine components (e.g., energy degrader and dosimetric phantom). Together with the physics processes available within Geant4, it allows us to account for all the particlematter interactions required for accurate beam transport and simulation of clinical dosimetric properties but it also paves the way for seamless and self-consistent simulation of radiation protection quantities such as the equivalent ambient dose for the whole treatment facility.
These are distinguishing and novel features of our approach that contrast with the typical methodology where several assumptions are used to describe the beam losses in separate codes that are then used as source terms in Monte Carlo codes that do not include the magnetic tracking of the primary beam (see, e.g., [21,22]). Such assumptions break down with the recent advances of tightly integrated proton therapy systems where discrete point-like source terms and a decoupled approach lack validity. The results discussed in this paper show that a self-consistent approach is possible and within the reach of a single multipurpose simulation code.
The flexibility of BDSIM allows for further applications not demonstrated in this paper. For example, these computationally expensive simulations can be integrated within a machine learning approach such as the one described in [23] to provide optimization capabilities useful at the design stage where codes such as MAD-X or Transport are traditionally used. Also, Monte Carlo codes dedicated to radiotherapy nozzle simulations such as TOPAS are able to provide inputs for the TPS; this is a work-inprogress for the Proteus R One using BDSIM.
Model construction. -We apply our approach to the IBA Proteus R One system (P1). The P1 is a singleroom therapy center equipped with a superconducting 230 MeV synchro-cyclotron (S2C2), followed by an extraction beamline comprising an energy degrader mounted on a rotating wheel and a rotating compact gantry (CGTR), as depicted in fig. 1. The energy degrader system produces a large number of losses, especially when the beam is degraded toward low clinical ranges. This also increases the beam emittance and consequently the beam envelope in the beamline, in turn increasing distributed losses. As P1 is a compact system (the total footprint including the shielding walls is 27.4 × 12.8 × 9.6 meters), the different parts of the system cannot be treated as truly independent. It is also a flexible system where new treatment modalities such as proton arc-therapy are in development [2], making use of treatment plans requiring beams with multiple characteristics and intensities. Figure 1 schematically represents the accelerator, the energy degrader beamline, and the rotating gantry. A short extraction beamline focusses the extracted beam 50004-p2  fig. 2). The superconducting synchro-cyclotron (S2C2), the extraction line (containing a pair of quadrupoles interspersed with aperture limiting slits), the energy degrader system and the rotating gantry are shown. The concrete shielding wall (dark grey) separating the cyclotron vault from the treatment area is visible with the beamline elements fit through a cylindrical cut. A cubic water volume is placed at the isocenter for energy deposition scoring. The location of the measurement devices -beam profile monitors (P1-4G), the ionization chambers (ICs) and the isocenter-are also indicated (green). The last magnet (B3G) uses converted CAD geometry for detailed aperture and magnetic fields.
to the energy degrader, equipped with a downstream circular collimator. The input beam of the model originates from a detailed simulation of the S2C2 beam dynamics [24]. The main beamline consists of three dipole magnets and seven quadrupole magnets. Horizontal and vertical beam limiting slits are installed upstream and downstream of the energy degrader. Further slits are also present upstream of the second dipole (B2G), at a location of large dispersion, and act as the energy selection system. BDSIM generates a 3D model of the beamline from a simple textual description of the sequence of magnetic elements. Generic geometries provided by BD-SIM are used for the quadrupoles in the beamline, for the beam limiting slits and for the steering magnets. To make the model more realistic, their shapes are adjusted to be close to the actual magnets, e.g., the accurate cylindrical shape of the Q1G and Q2G magnets located in the rotating assembly within the wall separating the cyclotron vault from the treatment room. The dipole magnets have beam pipes and yokes of more complicated shapes. In order to create a realistic model, the geometry obtained directly from the CAD software is used for the third bending element (B3G). Its detailed field model and geometry are of crucial importance due to the upstream scanning method employed for beam delivery. The conversion uses pyg4ometry [25], a Python library allowing the conversion to and from different geometry description formats. In particular, the conversion from STEP files to GDML (the Geant4 geometry persistency format [26]) is processed automatically. To navigate the geometrical structure and ensure correct tracking the physical volumes in Geant4 must not overlap; however STEP geometries do not have this requirement. The geometries have been adjusted manually and pyg4ometry was subsequently used to automatically check for overlaps. In addition, field maps obtained with Opera3D are attached to the dipole magnets so that the fringe fields are taken into account for the primary beam, while the field in the return yoke is accurately described for the tracking of the charged secondaries.
The model is completed with the details of the concrete shielding surrounding the cyclotron vault and the treatment room. This complete model is shown in fig. 2. The maze at the entrance of the room is realistically modelled with tapered blocks. In addition, the model makes use of the importance sampling variance reduction technique for some crucial shielding zones such as the wall separating the vault from the treatment room and the maze [11].

50004-p3
Model validation. -The model is validated with experimental data for three key observable quantities: beamline transmission efficiency, energy deposition at isocenter (Bragg peak) and beam profile measurements. This allows confident use of the model for two types of simulations: on one hand results regarding the beam losses, their interactions with the surrounding beamline elements and the concrete shielding, on the other hand the precise evaluation of beam delivery and dose deposition at the isocenter.
First, the beamline transmission is compared with measured values. This provides a global assessment of local features that cannot be measured directly, e.g., aperture restrictions and emittance increases. The emittance increase induced by the energy degrader and the effect of the circular collimator on the beam transverse and energy spectra have been reported in [27]. The total efficiency is defined as the product of two terms. The first term represents the extraction beamline efficiency extr . The second one is the beamline efficiency beam that represents the transmission from the degrader to the isocenter.
where I P1E specific is the beam current measured at the Beam Profile Monitor (BPM) P1E, placed on the degrader wheel, when the slits of the extraction beamline are set to their specific opening corresponding to the requested range; I P1Eopen is the beam current measured at P1E when the extraction slits are completely open and I IC is the beam current measured at the ionization chambers (IC) located between B3G and the isocenter. The results are shown in fig. 3 with comparison to measurements at two independent installation sites with their own calibration. One can observe the excellent agreement at all ranges even in the transition zones between the degrader materials. This result validates the global BDSIM model of the P1 and its degrader. It should be noted that at high energies the horizontal slits in the extraction beamline are purposely closed to artificially reduce the beamline transmission, in order to avoid too large variation with energy. The depth-dose profile in a water phantom located at isocenter serves as an observable to assess the energy spectrum of the degraded beam. The energy deposition is scored using a 6 cm × 6 cm × 1 mm mesh. The depthdose is thus an integral measurement over the transverse beam size, even as it scatters with depth, as provided by the StingRay probe [28] used for the measurements. The results exhibit an outstanding agreement between the experimental data and the simulated Bragg peaks over the full span of ranges at isocenter, as shown in fig. 4 for an example requested ranges of 8.82 cm ( fig. 4(a)) and 28.65 cm ( fig. 4(b)) where experimental data were available. As an experimental absolute degrader angle to range calibration is unavailable, a virtual calibration of the degrader wheel using BDSIM was performed for 123 degrader positions each with 2 × 10 7 events. With linear interpolation, Fig. 3: Absolute error on the total efficiency t (eq. (1)) as a function of the beam range at isocenter. The error is given with respect to two sets of experimental data independently measured on two IBA P1 sites. The three materials composing the degrader wheel are indicated. the model can be queried for a given range at isocenter. Therefore, the measured and simulated R 90 proton ranges (90% dose in the distal region of the peak) are equal by construction. The possibly different overall shape of the peaks and the dose at skin are predicted within a 2-3% error margin. The distal fall-off, defined as the distance 50004-p4   between the distal 80% and 20% normalised dose, is also predicted with an error of less than 100 mm. A systematic error analysis was performed by offsetting the simulated data in range and comparing the χ 2 , where the minimal χ 2 was found for a range offset of 80 micron. This is far below the measurement resolution and we therefore conclude that the virtual calibration does not contribute to the residual error on the dose points. The transverse beam profiles along the beamline have also been validated against measured data. Whilst the beamline transmission validation already provides a global insight, the validation of the beam profiles provides a local validation required to provide clinically relevant information about the beam spot sizes. As shown in fig. 1, four BPMs are located along the beamline (P1-4G). The intercepting BPMs have unequal strip sensor widths and a Gaussian profile is fitted to the histogram by the machine control system. Due to aperture restrictions and collimation, the beam profiles are not expected to be Gaussian, however, only this fitted σ from the control system is available for comparison. The post-processing of the BDSIM simulations reproduces these steps using data collected in the simulation from passive "sampling" planes immediately after the position of each BPM. An example comparison between experimental and simulated results (1M initial protons) is shown in fig. 5 for P2G and a range of 15.5 cm at the isocenter. These comparisons include the non-Gaussian input beam distribution, interaction in the degrader and aperture restrictions, so provide a combined assessment of the transverse emittances and beam optics in a realistic method. Measured data have been acquired for three different P1 centers (labelled I, II and III). The systems of each center differ by their respective imperfections and are individually calibrated to meet the specifications at the isocenter. The simulated profiles are in good agreement with the experimental data and the results lie within the natural variation between centers of 5 to 10%.
The BDSIM model does not include magnetic or alignment imperfections and we use the optics configuration of center III as the nominal configuration. Similar discrepancies have been already observed with the OPAL model of the Proscan facility [14]. These conclusions apply similarly to the other monitors and all ranges.
We conclude that the detailed Geant4 physics models together with the realistic beamline geometry and field descriptions implemented in BDSIM are key factors to reproduce the measured data with a high accuracy. It is suitable to produce validated simulation results for dose deposition and radiation fluences to be used for residual dose computations and material activation either directly within BDSIM or for external codes, e.g., FISPACT-II [29].
Results and H * (10) ambient dose maps. -The interactions of the beam with the degrader and the different elements present in the beamline induce the emission of ionizing particles all around the machine. This radiation has an impact on the activation of the concrete shielding, on the beamline elements and electronics (ageing and activation) and can contribute extra dose to the patient. Preliminary studies of the shielding activation using BDSIM can be found in [11] following a benchmark of concrete shielding activation computation with GEANT4 against other Monte Carlo codes (MCNPX and PHITS) [30]. The benchmark shows that the use of the Lige Intranuclear Cascade model (G4 INCL) or the Binary Cascade Model (G4 BIC) physics neutron transport model provide a reasonable agreement with the MCNPX results taken as reference. 3D per-particle species fluence and dose maps form key contributions of the present work and can be used to evaluate the above-mentioned phenomena. The compound quantity H * (10) has been chosen to evaluate these crucial aspects of our model for its relevance in the context of radiological protection for personnel. Other fluence-related quantities can be obtained similarly.  radiological protection purposes, protection quantities are defined by the International Commission on Radiological Protection (ICRP) and are used to impose limits on radiation exposure in medical facilities [31]. One of the most used protection quantities is the equivalent dose, H T , in a tissue or organ T given by H T = R w R × D T,R where D T,R is the average absorbed dose from radiation R, in tissue T, and w R the radiation weighting factor of radiation R [31]. The protection quantities cannot be directly measured but can be reasonably estimated by operational quantities, defined by the International Commission on Radiation Units (ICRU), which are related to the radiation field characteristics [31]. The operational quantity suggested for area monitoring of strongly penetrating radiation is the ambient dose equivalent, H * (10), defined as the dose produced in the ICRU sphere at a depth of 10 mm on the radius opposing the direction of the irradiation field [32], where the ICRU sphere is defined as a 30 cm diameter, tissue-equivalent sphere of unit density. The ambient dose equivalent can be computed by making the product between the particle fluence φ (cm −2 ) and the fluence-to-ambient dose equivalent conversion coefficients h 10 (Sv × cm 2 ). The resulting expression for the ambient dose is H * (10) = φ(E, particle) × h 10 (E, particle). The h 10 conversion coefficients of ref. [32] (computed using the Monte Carlo transport code FLUKA) are parametrised using a double logarithmic cubic spline fit to increase the precision in energy [33] and can be queried in BDSIM. A specific ambient dose scoring capability has been developed through the use of fluence scorers and the conversion factors and is now fully integrated within BDSIM.
The projected dose maps are shown in fig. 6. The H * (10) ambient equivalent dose is integrated over a volume |y| < 50 cm on both sides of the cyclotron mid-plane. A model of the cyclotron losses, obtained from prior detailed simulations of the S2C2 extraction has been implemented statically using distributed point-like sources. It accounts for the losses during acceleration and at extraction. The resulting dose map is represented in fig. 6(a) which shows, as expected, that the H * (10) dose in the treatment area is five orders of magnitude lower than the peak value in the cyclotron vault. The large contribution of the cyclotron extraction is also visible. The dose maps obtained from the fluence induced by the beamline losses are shown in fig. 6(b). The major contribution to the losses and resulting fluence, at such low energy, comes from the degrader and the circular collimator. Distributed losses are present along the beamline, while few losses are present in the final part of the gantry. The interactions of the primary beam with the air column at the exit of the gantry are also visible. The total maps are presented in 50004-p6 fig. 6(c) and fig. 6(d) for a range of 8.82 cm and 28.65 cm, respectively. A direct conclusion is reached if one refers to the large difference in beamline transmission ( fig. 3) between these two ranges. For low ranges, the global transmission is low, but the losses occur mainly in the energy degradation section, while, for higher ranges, the beam transmission is higher, but the losses are more distributed along the gantry. This also comes from the fact that at higher energies, the emittance increase induced by the degrader is not immediately cut by the circulator collimator but downstream by the transverse slits and by the beam pipe.
Conclusion. -The model and results presented in this paper show that a self-consistent approach, considering the full coupling between all components of a hadron therapy treatment system is possible and within the reach of a single multi-purpose simulation code. The seamless simulation from the cyclotron extraction to the treatment room allows computing all the major quantities of interest. The results show excellent agreement with the design performances of the system. In particular, the simulations provide a truly realistic beam distribution, including the energy spectrum, and secondaries at the end of the beamline that can be used for detailed modeling of the beam delivery systems to assess the dose deposition at the patient. The equivalent ambient dose is obtained from the beam losses arising from the detailed modeling and simulation of the beam transport system. It paves the way for further use of simulated radiation fluence data in the cyclotron vault and treatment room. The availability of a complete model reproducing key experimental figures, as provided by BDSIM in the form of a start-to-end selfconsistent model provides a complete numerical toolbox for machine studies and preparation of beam configurations compatible with novel treatment modalities. * * *