A Geant4 Based Simulation Platform of the HollandPTC R&D Proton Beamline for Radiobiological Studies

A Geant4 based simulation platform of the Holland Proton Therapy Centre (HollandPTC, Netherlands) R&D beamline (G4HPTC-R&D) was developed to enable the planning, optimisation and advanced dosimetry for radiobiological studies. It implemented a six parameter non-symmetrical Gaussian pencil beam surrogate model to simulate the R&D beamline in both a pencil beam and passively scattered field configuration. Three different experimental proton datasets (70 MeV, 150 MeV, and 240 MeV) of the pencil beam envelope evolution in free air and depth-dose profiles in water were used to develop a set of individual parameter surrogate functions to enable the modelling of the non-symmetrical Gaussian pencil beam properties with only the ProBeam isochronous cyclotron mean extraction proton energy as input. This refined beam model was then benchmarked with respect to three independent experimental datasets of the R&D beamline operating in both a pencil beam configuration at 120 and 200 MeV, and passively scattered field configuration at 150 MeV. It was shown that the G4HPTC-R&D simulation platform can reproduce the pencil beam envelope evolution in free air and depth-dose profiles to within an accuracy on the order of $\pm$5% for all tested energies, and that it was able to reproduce the 150 MeV passively scattered field to the specifications need for clinical and radiobiological applications.


Introduction
Proton radiotherapy is a cancer treatment modality that has seen increased use over the last two decades due to its ability to improve local conformity of radiation dose to target tumours whilst sparing the surrounding healthy tissue [1][2][3][4][5].The Holland Proton Therapy Centre (HollandPTC) is one of three proton radiotherapy centres in the Netherlands and has treated cancer patients that qualify for proton radiotherapy since 2018 [6].HollandPTC is a ProBeam (Varian, a Siemens Healthineers Company) isochronous cyclotron-based facility that features pencil beam scanning into two gantry rooms, one fixed horizontal eye treatment beamline, and one fixed horizontal Research & Development (R&D) beamline capable of producing proton energies ranging from 70 MeV up to 250 MeV [7].HollandPTC is the only proton radiotherapy facility in the Netherlands that possesses a dedicated clinical R&D beamline, and therefore maximising the scientific output from each beamtime session is of great importance.
The HollandPTC R&D beamline, designed and developed specifically for HollandPTC, consists of equipment to build a passive scattering Australia.
system in order to produce large fields of varying sizes to facilitate radiobiological studies with clinically relevant proton energies [7].To ensure that every allocated experimental beamtime is used efficiently significant planning, workflow optimisation, and pre-irradiation preparation is required.One useful way to ensure that each experimental configuration will yield the desired proton field shape, intensity, and incident energy spectra is to undertake in silico trials through the use of Monte Carlo radiation transport modelling toolkits such as Geant4 [8][9][10], FLUKA [11,12] and MCNP [13,14].This approach is standard for lower energy proton beamlines [15][16][17][18][19], and has been shown to be crucial for clinical energy proton passive scattering beamlines to enable accurate determination of delivered dose/LET in radiobiological studies [20,21].This work presents the development of a Geant4 based simulation platform of the HollandPTC R&D beamline (G4HPTC-R&D) to enable the planning, optimisation, and advanced dosimetry for radiobiological studies in both a pencil beam and passively scattered field configuration.In contrast to past studies, this work implements a six parameter https://doi.org/10.1016/j.ejmp.2023.102643Received 22 February 2023; Received in revised form 1 June 2023; Accepted 5 July 2023 non-symmetrical Gaussian pencil beam surrogate model to simulate the pencil beam properties with only the ProBeam isochronous cyclotron mean extraction proton energy as input.Section 2 describes the development of the G4HPTC-R&D simulation platform and its individual parameter surrogate functions for non-symmetrical Gaussian pencil beam model through optimisation with respect to three proton energy experimental datasets, and G4HPTC-R&D's benchmarking with respect to an additional three independent experimental datasets.The results of this process and an accompanying discussion can be found in Section 3, with an overall conclusion following in Section 4.

G4HPTC-R&D simulation platform development and proton pencil beam model optimisation
Geant4 version 10.06.p01 was utilised to develop the G4HPTC-R&D simulation platform based on the experimental geometry of the HollandPTC R&D beamline in its passively scattered field configuration seen in Fig. 1.A total of eight key geometric elements outlined in Table 1 were implemented in G4HPTC-R&D that could be enabled or disabled depending on the beamline operational mode (i.e.pencil beam or passively scattered field configuration).The scattering foil and dual ring seen in Fig. 1 facilitate the expansion and shaping of the initial proton pencil beam to generate a uniform intensity proton field.The expanded and shaped beam is collimated at two stages along its path through the use of a first and second stage variable open cross-section brass collimator to produce square fields of up to 200 mm × 200 mm.At the irradiation/measurement stage where different radiobiological endstations are placed in Fig. 1, a Lynx ® detector (IBA Dosimetry, Schwarzenbruck, Germany) can be seen which is used to assess the shape and quality of the proton field.Finally, an additional geometrical element that is not show in Fig. 1 was implemented: a water box that mimics the properties of a QUBEnext (DE.TEC.TOR, Turin, Italy) detector. 1Table 1 outlines the dimensions and materials of the different elements, with additional information relating to their design and orientation along the beamline's path outlined in Rovituso et al. [7].
The non-symmetrical Gaussian pencil proton beam of the Hol-landPTC R&D beamline was implemented in G4HPTC-R&D using a six parameter surrogate model emerging after the Kapton vacuum pipe exit window. 2 These six parameters model three important characteristics of the proton beam at the Kapton vacuum pipe exit window: the 2D Gaussian lateral beam spot size (  ,   ), the 2D Gaussian beam spot angular deviation (  ,   ), and the 1D Gaussian energy spread of the cyclotron generated proton beam with initial mean energy  0 and energy spread .Through the optimisation workflow outlined in Fig. 2, a set of individual surrogate functions were developed for each of these parameters to enable the modelling of the Gaussian pencil beam properties using the ProBeam isochronous cyclotron mean extraction proton energy with respect to experimental measurements at 70, 150, and 240 MeV.For each combination of these six parameters and detection media (Lynx ® detector or surrogate QUBEnext detector water box) investigated with G4HPTC-R&D, a total of 10 6 primary protons were run and the transport of all particles was simulated using a combined Geant4 ''Standard EM Option 4'' and ''QGSP_BIC_HP'' physics list [10,22] with atomic deexcitation enabled, a particle production range cut of 200 μm, and a low-energy cut off of 250 eV.

Table 1
Name, dimensions, and materials of the geometric elements that were implemented in G4HPTC-R&D to mimic the experimental configuration seen in Fig. 1.In the first stage, the lateral beam spot size (  ,   ) and angular deviation (  ,   ) were optimised with respect to the proton pencil beam envelope cross-section measured with a Lynx ® detector in free air at 230, 530, 911, 1230, 1530, 1830 and 2045 mm down-stream from the Kapton vacuum pipe exit window.Both sets of experimental and simulated Lynx ® detector data 2D beam profiles were fitted with a 2D Gaussian function to obtain the Full Width at Half Maximum (FWHM) in  and .The Full Width at Tenth Maximum (FWTM) was extracted from the central -and -axis planes to investigate the tails of the distributions as a second figure of merit.The agreement between experimental and simulated FWHM values as a function of the pencil beam energy was assessed through the us of the Sum of Squared Errors (SSE) metric:

Name
where FWHM sim, is the simulated and FWHM exp, is the experimental FWHM summed over  = 7 distances.Starting with (predefined) initial values for   ,   ,   and   from Rovituso et al. [7], a two step optimisation with the SSE metric was undertaken.In the first step a ±10% offset sweep around initial values was explored, and the combination that resulted in the smallest SSE was selected as an initial estimates for each parameter.A second parameter sweep of all four of the initial estimate for each parameter with a ±5% offset was then undertaken to fine tune and ensure that each parameter value was not a local minimum in the optimisation.
In the second stage, the two remaining beam parameters ( 0 and ) that modelled the initial mean energy and energy spread of the cyclotron generated proton pencil beam were optimised with respect to the experimental proton pencil beam depth-dose profiles at 70, 150, and 240 MeV.Experimental measurements of the proton pencil beam depth-dose profiles at 70, 150, and 240 MeV were obtained with a QUBEnext detector placed at the beam isocentre 911 mm from the Kapton vacuum pipe exit window. 3Using the optimised lateral beam spot  size (  ,   ) and angular deviation (  ,   ) values at each energy, the mean energy ( 0 ) was varied in steps of 0.1 MeV around the ProBeam isochronous cyclotron mean extraction proton energy and the energy spread () in ±5 steps of 0.05 MeV with respect to initial estimates taken from Rovituso et al. [7] to compare to the experimental proton pencil beam depth-dose profiles.Each proton pencil beam depth-dose profile data was fitted with a Bortfeld function [23], and the position of the 80% dose in the distal falloff (R80), the distance between the distal position of the 80% and 20% dose values (R80-R20 distal falloff), and peak-to-entrance ratio were extracted as figures of merit [24].Comparison of the simulation and experimental results for these three figures of merit were used to optimise  0 , , and to provided a general figure of merit.
The set of six parameter values obtained in the first and second stages of the optimisation workflow outlined in Fig. 2 form the basis of the developed individual parameter surrogate functions (third stage).They were solved through the mapping of each parameter value at 70, 150 and 240 MeV to second-order polynomial functions.With these surrogate functions G4HPTC-R&D can model the non-symmetrical Gaussian pencil beam properties at the Kapton vacuum pipe exit window of the R&D beamline with only the ProBeam isochronous cyclotron mean extraction proton energy as input.

G4HPTC-R&D independent experimental benchmarking
The refined non-symmetrical Gaussian pencil beam model and G4HPTC-R&D simulation platform was benchmarked with respect to three independent experimental datasets.Two of these experimental datasets were of the HollandPTC R&D beamline operating in its pencil beam configuration at 120 and 200 MeV, and the other was the HollandPTC R&D beamline operating in its passively scattered field configuration at 150 MeV to generate a 100 mm × 100 mm field at the irradiation/measurement stage.Experimental measurements and G4HPTC-R&D simulations were undertaken at 120 and 200 MeV in an identical manner to that outlined above to obtain beam envelope evolution in free air and depth-dose profiles in water datasets.These experimental and simulated FWHM, FWTM, R80, R80-R20 distal falloff, and peak-to-entrance ratio results at each energy were compared to assess the validity of the refined non-symmetrical Gaussian pencil beam model and G4HPTC-R&D simulation platform.
Experimental measurement of the 100 mm × 100 mm field generated at the irradiation/measurement stage for the HollandPTC R&D beamline operating in its passively scattered field configuration was undertaken using the Lynx ® detector at 150 MeV.Large field simulations were also performed at 150 MeV with all beam elements implemented as shown in Fig. 1 and for an inner open cross-section of the final stage beam defining brass collimator set to 100 mm × 100 mm.A total of 4 × 10 7 primary protons were run, and the transport of all particles was simulated using the same physics configuration outlined above for the pencil beam configuration simulations.Three figures of merit were utilised to assess the validity of the refined G4HPTC-R&D simulation platform operating in a passively scattered field configuration: (1) the field uniformity  , (2) the -index mean value, and (3) the -index global pass rate.The field uniformity ( ) figure of merit across the field for both the experimental and simulation data was calculated by: where   is the maximum intensity across the field, and   is the minimum intensity across the field.It should be noted that a uniformity of above 97% is required for radiobiological experiments [7,20].Whereas for the -index mean value and global pass rate figures of merit, the -index can be written: with: where   is the spatial location in the experimental field,   is the spatial location in the simulated G4HPTC-R&D field, (  ,   ) is the distance between the two locations, (  ,   ) difference in dose of between the two locations,  is the distance difference criterion, and  is the dose difference criterion [25,26].A total of three different  and  criterion combinations were assessed respectively: (1) 3 mm and 3%, (2) 4 mm and 4%, and (3) 5 mm and 5%.

G4HPTC-R&D proton pencil beam model optimisation
Table 2 presents the optimal values for the non-symmetrical Gaussian pencil beam surrogate model parameters   ,   ,   and   obtained through the optimisation workflow outlined in Fig. 2 at 70, 150 and 240 MeV.In the case of   ,   ,   ,   and E, all of these parameters decrease with increasing initial mean proton energy as expected [27][28][29].Whereas for the E 0 , the obtained values are within less than 2% of the requested ProBeam isochronous cyclotron mean extraction proton energy.This small difference between experimental and G4HPTC-R&D value for the E 0 parameter can be attributed to: (1) generation of the non-symmetrical Gaussian pencil beam spot after the Kapton vacuum pipe exit window, (2) uncertainties in the alignment and the resolution of the QUBEnext detector, and (3) uncertainties in relevant Geant4 proton cross-sectional data and physics models which are on the order of ±10% [22,29].
Fig. 3 presents a comparison between the fitted FWHM and FWTM extracted values from the experimental and G4HPTC-R&D proton pencil beam envelope cross-section measurements in free air at 70, 150, and 240 MeV using the optimal values for   ,   ,   , and   outlined in Table 2.The dotted lines represent a ±5% deviation with respect to experimental results to aid in assessing the accuracy of the G4HPTC-R&D results.For all energies the G4HPTC-R&D FWHM values in both  and  were reproduced to within less than 5% indicating a high level of correlation.With the FWTM values, the G4HPTC-R&D 70 and 150 MeV data is also within 5% of experimental results indicating a high level of correlation.However, at 240 MeV the G4HPTC-R&D 2D proton pencil beam envelope cross-section FWTM values increase as a function of distance in air with the furthest two distances having a difference of up to 9% relative to their respective experimental values.It should be noted that this level difference is still within the established uncertainties in relevant proton Geant4 proton cross-sectional data and physics models which are on the order of ±10% [22,29], and the fact that this difference is still on the order of 5% indicates a high level of correlation between the experimental and G4HPTC-R&D data.
The depth-dose distributions of experimental and G4HPTC-R&D results obtained with these optimised values of Table 2 are shown in Fig. 4, and their corresponding proton range (R80), distal fall-off (R80-R20) and peak-to-entrance ratio values are displayed in Table 3.The depth-dose profiles are normalised on the maximum dose, and show agreement in R80 and R80-R20 to within less than 0.5 mm (or 2%) for all energies.As for the peak-to-entrance ratio, a difference on the order of 5% can be observed for all three energies.Again given that the difference for all figures of merit are on the order or less than 5%, a high level of correlation is present between the experimental and G4HPTC-R&D depth-dose profiles at 70, 150, and 240 MeV.

Refined G4HPTC-R&D proton pencil beam model independent experimental benchmarking
Table 4 presents the second-order polynomial function constant values of the six refined non-symmetrical Gaussian pencil beam surrogate model parameters obtained via the optimisation workflow outlined in Section 2. A comparison between the fitted FWHM and FWTM extracted values from the experimental and refined G4HPTC-R&D proton pencil beam envelope cross-section measurements in free air for a mean ProBeam isochronous cyclotron extraction proton energy of 120 and 200 MeV using these values can be seen in Fig. 5.As in Fig. 3, the dotted lines represent a ±5% deviation with respect to experimental results to aid in assessing the accuracy of the G4HPTC-R&D results.Inspection of Fig. 5 shows that the G4HPTC-R&D FWHM values in both  and  were reproduced to within 5% of the experimental data at 120 and 200 MeV.With the exception of a single outlier in the 120 MeV data at a distance of 1830 mm (7.6% difference in the ), it can be seen that the G4HPTC-R&D FWTM values in both  and  were reproduced again to within 5% of the experimental data at the 120 and 200 MeV.
The depth-dose distributions of experimental and refined G4HPTC-R&D results for a mean ProBeam isochronous cyclotron extraction proton energy of 120 and 200 MeV are shown in Fig. 6, and their corresponding proton range (R80), distal fall-off (R80-R20) and peakto-entrance ratio values are displayed in Table 5.For the 120 MeV depth-dose profile data the G4HPTC-R&D R80 value agrees to within 0.5 mm (1%) of the experimental data, and R80-R20 value is the same as the experimental data.Whereas for the 200 MeV depth-dose profiles, the G4HPTC-R&D R80 value is within 1.5 mm (1%) to the experimental data and R80-R20 is on the order of 0.5 mm (11%) less then the experimental data.Furthermore, Table 5 also shows that a difference on the order or less than 5% is present in the experimental and G4HPTC-R&D peak-to-entrance ratio values at both the 120 and 200 MeV.

Table 5
Proton range R80 (mm), distal fall-off R80-R20 (mm) and peak-to-entrance ratios of the experimental and G4HPTC-R&D of the 120 and 200 MeV depth-dose distributions seen in Fig. 6.Fig. 7 presents the experimental and refined G4HPTC-R&D 100 mm × 100 mm fields and their respective central -axis lines profiles generated at the irradiation/measurement stage of the HollandPTC R&D beamline operating in its passively scattered field configuration.
Both fields appear similar in shape and intensity, with their central axis lines profiles illustrating that the maximum difference in between the two is less than 5%.Assessment of the field uniformity in the axis yields a value of 97.4% and 96.7%, and in the -axis yields a  value of 97.7% and 96.9% for the experimental and G4HPTC-R&D data respectively.This less than 1% field uniformity difference in along both axis indicates a good level of agreement between the experimental and G4HPTC-R&D data, and that the field quality is sufficient for radiobiological experiments [7,20].Table 6 presents the G4HPTC-R&D 100 mm × 100 mm field index mean and -index global pass rate values with respect to the experimental 100 mm × 100 mm field for three different  and  criterion combinations.It can be seen that as the  and  combinations increase in value, the -index mean value decreases and the -index global pass rate increases.Under the most strict  = 3 mm and  = 3% criterion combination, the -index global pass rate is 96.6% and exceeds the clinically accepted pass rate threshold of 90% for this criterion combination [30,31].For  = 4 mm and  = 4% the -index global pass rate increases to 99.9%, and then to 100% for  = 5 mm and  = 5% criterion combination.
These three independent experimental benchmarking trials of the refined G4HPTC-R&D simulation platform illustrated that it is able to reproduce the physical characteristics of the HollandPTC R&D beamline operating in both its pencil beam and passively scattered field configurations to within an acceptable level of agreement for clinical and radiobiological applications.The slight differences observed between the refined G4HPTC-R&D and experimental data for these three independent experimental benchmarking trials can be attributed to two primary factors: (1) uncertainties in the alignment each beamline element and the resolution of the Lynx ® /QUBEnext detector, and (2) uncertainties in relevant Geant4 proton cross-sectional data and physics models which are on the order of ±10% [22,29].The impact of beamline geometric element alignment uncertainties is particular relevant in the case of the 100 mm × 100 mm passively scattered field data where misalignment on the level of a millimetre or two can cause excessive ''haloing'' around the edges of the field.Even under optimal alignment conditions, this ''haloing'' effect is still present and can be observed in the both 2D field maps and central -axis line profiles seen in Fig. 7 via the increased in relative intensity at the field edges and corners.Further work is already underway to explore the impact of these element alignment uncertainties on the 2D proton Linear Energy Transfer (LET) spectra distribution at the surface of cell/tissue culture for the different radiobiological endstations under development at the HollandPTC R&D beamline [7].

Conclusion
This work both developed and characterised the performance of a Geant4 based simulation platform of the Holland Proton Therapy Centre (Netherlands) R&D beamline (G4HPTC-R&D) to enable the planning, optimisation, and advanced dosimetry for radiobiological studies in both a pencil beam and passively scattered field configuration.It implemented a six parameter non-symmetrical Gaussian pencil beam surrogate model to simulate the R&D beamline in both a pencil beam and passively scattered field configuration.Three different experimental proton datasets (70 MeV, 150 MeV, and 240 MeV) of the pencil beam envelope evolution in free air and depth-dose profiles in water were used to develop a set of individual parameter surrogate functions to enable the modelling of the non-symmetrical Gaussian pencil beam properties with only the ProBeam isochronous cyclotron mean extraction proton energy as input.This refined beam model was then benchmarked with respect to three independent experimental datasets of the R&D beamline operating in both a pencil beam configuration at 120 and 200 MeV, and passively scattered field configuration at 150 MeV.It was shown that the G4HPTC-R&D simulation platform can reproduce the pencil beam envelope evolution in free air and depthdose profiles to within an accuracy on the order of ±5% for all tested energies, and that it was able to reproduce the 150 MeV passively scattered field to the specifications need for clinical and radiobiological applications.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: This work was partially funded by Varian Medical Systems (a Siemens Healthineers Company).

Fig. 1 .
Fig. 1.HollandPTC R&D beamline configured in its passively scattered field configuration (top).Seven of the eight implemented geometric elements in G4HPTC-R&D to mimic this experimental configuration, and their relative distances with respect to the Kapton vacuum pipe exit window in mm, can be seen in (bottom).Here the seven of the eight key geometric elements are: (1) the kapton vacuum pipe exit window, (2) scattering foil, (3) beam monitor, (4) dual ring, (5) first stage beam defining collimator, (6) second stage beam defining collimator, and (7) the front of the irradiation/measurement stage.Note that in (top) a Lynx ® detector can be seen at the irradiation/measurement stage, and the QUBEnext detector is not shown.

Fig. 2 .
Fig. 2. G4HPTC-R&D simulation platform six parameter non-symmetrical Gaussian pencil beam surrogate model optimisation and individual parameter surrogate function mapping workflow with respect to experimental measurement at 70, 150, and 240 MeV.This parameter optimisation and individual surrogate function mapping workflow is composed of three stages: (1) lateral beam spot size (  ,   ) and angular deviation (  ,   ) optimisation with respect to the proton pencil beam envelope cross-sectionals evolution measurements in free air, (2) initial proton beam mean energy (E 0 ) and energy spread (E) optimisation with respect to proton pencil beam depth-dose profiles in water, and (3) mapping of each optimised parameter value to an individual parameter second-order polynomial surrogate function.

Fig. 3 .
Fig. 3. Experimental and G4HPTC-R&D 2D beam cross-section FWHM and FWTM values as a function of distance from the Kapton vacuum pipe exit window at 70, 150 and 240 MeV.Here error bars representing the accuracy of FHWM data fitting cannot be resolved due to their scale being on the same order as FWHM data symbols.

Fig. 5 .
Fig. 5. Experimental and G4HPTC-R&D 2D beam cross-section FWHM and FWTM values as a function of distance from the Kapton vacuum pipe exit window at 120, and 200 MeV.Here error bars representing the accuracy of FHWM data fitting cannot be resolved due to their scale being on the same order as FWHM data symbols.

Fig. 7 .
Fig. 7. Experimental and G4HPTC-R&D 100 mm × 100 mm fields (top) and their respective central -axis lines profiles (bottom) generated at the irradiation/measurement stage of the HollandPTC R&D beamline operating in its passively scattered field configuration.Here the pink horizontal lines (top) represent where the central -axis lines profiles were taken, and the blue dotted lines outline the target field 100 mm × 100 mm of interest.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2
Optimised non-symmetrical Gaussian pencil beam surrogate model parameters at 70, 150 and 240 MeV: lateral spread  in  and  (mm) and angular spread  in  and  (radians), initial mean energy E 0 (MeV) and energy spread E (MeV).

Table 4
The second-order polynomial function constant values of the six refined nonsymmetrical Gaussian pencil beam surrogate model parameters obtained via the optimisation workflow outlined in Section 2.Here the second-order polynomial function is written Parameter(E) = E 2 + E +  where E is the mean extraction proton energy of the ProBeam isochronous cyclotron.

Table 6
G4HPTC-R&D 100 mm ×100 mm field -index mean and -index global pass rate values with respect to the experimental 100 mm ×100 mm field for three different  and  criterion combinations.