# Validation of a Monte Carlo framework for out-of-field dose calculations in proton therapy

Marijke de Saint-Hubert, Nico Verbeek, Christian Bäumer, Johannes Esser, Jörg Wulff, Racell Nabha, Olivier van Hoey, Jérémie Dabin, Florian Stuckmann, Fabiano Vasi, et al.

► To cite this version:

Marijke de Saint-Hubert, Nico Verbeek, Christian Bäumer, Johannes Esser, Jörg Wulff, et al. Validation of a Monte Carlo framework for out-of-field dose calculations in proton therapy. Frontiers in Oncology, 2022, 12, pp.882489. 10.3389/fonc.2022.882489. cea-04125591

HAL Id: cea-04125591

https://cea.hal.science/cea-04125591

Submitted on 12 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# Validation of a Monte Carlo Framework for Out-of-Field Dose Calculations in Proton Therapy

Marijke De Saint-Hubert<sup>1†</sup>, Nico Verbeek<sup>2,3,4†</sup>, Christian Bäumer<sup>2,3,5,6</sup>, Johannes Esser<sup>2,3,7</sup>, Jörg Wulff<sup>2,3</sup>, Racell Nabha<sup>1</sup>, Olivier Van Hoey<sup>1</sup>, Jérémie Dabin<sup>1</sup>, Florian Stuckmann<sup>2,7,8</sup>, Fabiano Vasi<sup>9</sup>, Stephan Radonic<sup>9</sup>, Guillaume Boissonnat<sup>10</sup>, Uwe Schneider<sup>9</sup>, Miguel Rodriguez<sup>11,12</sup>, Beate Timmermann<sup>2,3,4,5,13</sup>, Isabelle Thierry-Chef<sup>14,15,16</sup> and Lorenzo Brualla<sup>2,3,4*</sup>

## OPEN ACCESS

**Edited by:**
Beatriz Sanchez-Nieto,
Pontifical Catholic University of Chile,
Chile

**Reviewed by:**
Silva Bortolussi,
University of Pavia, Italy
Francesco Tommasino,
University of Trento, Italy

**\*Correspondence:**
Lorenzo Brualla
Lorenzo.Brualla@uni-due.de

†These authors have contributed equally to this work and share first authorship

**Specialty section:**
This article was submitted to Radiation Oncology, a section of the journal Frontiers in Oncology

**Received:** 23 February 2022

**Accepted:** 04 May 2022

**Published:** 08 June 2022

**Citation:**
De Saint-Hubert M, Verbeek N, Bäumer C, Esser J, Wulff J, Nabha R, Van Hoey O, Dabin J, Stuckmann F, Vasi F, Radonic S, Boissonnat G, Schneider U, Rodriguez M, Timmermann B, Thierry-Chef I and Brualla L (2022) Validation of a Monte Carlo Framework for Out-of-Field Dose Calculations in Proton Therapy. *Front. Oncol.* 12:882489.
doi: 10.3389/fonc.2022.882489

<sup>1</sup>Research in Dosimetric Applications, Belgian Nuclear Research Center (SCK CEN), Mol, Belgium, <sup>2</sup>West German Proton Therapy Centre Essen WPE, Essen, Germany, <sup>3</sup>West German Cancer Center (WTZ), Essen, Germany, <sup>4</sup>Faculty of Medicine, University of Duisburg-Essen, Essen, Germany, <sup>5</sup>Radiation Oncology and Imaging, German Cancer Consortium DKTK, Heidelberg, Germany, <sup>6</sup>Department of Physics, TU Dortmund University, Dortmund, Germany, <sup>7</sup>Faculty of Mathematics and Science Institute of Physics and Medical Physics. Heinrich-Heine University, Düsseldorf, Germany, <sup>8</sup>Klinikum Fulda GAG, Universitätsmedizin Marburg, Fulda, Zurich, Germany, <sup>9</sup>Physik Institut, Universität Zürich, Zürich, Switzerland, <sup>10</sup>CEA, Université Paris-Saclay, Palaiseau, France, <sup>11</sup>Hospital Paltilla, Panama City, Panama, <sup>12</sup>Instituto de Investigaciones Científicas y de Alta Tecnología INDICASAT-AIP, Panama City, Panama, <sup>13</sup>Department of Particle Therapy, University Hospital Essen, Essen, Germany, <sup>14</sup>Radiation Programme, Barcelona Institute for Global Health (ISGlobal), Barcelona, Spain, <sup>15</sup>University Pompeu Fabra, Barcelona, Spain, <sup>16</sup>CIBER Epidemiología y Salud Pública, Madrid, Spain

Proton therapy enables to deliver highly conformed dose distributions owing to the characteristic Bragg peak and the finite range of protons. However, during proton therapy, secondary neutrons are created, which can travel long distances and deposit dose in out-of-field volumes. This out-of-field absorbed dose needs to be considered for radiation-induced secondary cancers, which are particularly relevant in the case of pediatric treatments. Unfortunately, no method exists in clinics for the computation of the out-of-field dose distributions in proton therapy. To help overcome this limitation, a computational tool has been developed based on the Monte Carlo code TOPAS. The purpose of this work is to evaluate the accuracy of this tool in comparison to experimental data obtained from an anthropomorphic phantom irradiation. An anthropomorphic phantom of a 5-year-old child (ATOM, CIRS) was irradiated for a brain tumor treatment in an IBA Proteus Plus facility using a pencil beam dedicated nozzle. The treatment consisted of three pencil beam scanning fields employing a lucite range shifter. Proton energies ranged from 100 to 165 MeV. A median dose of 50.4 Gy(RBE) with 1.8 Gy(RBE) per fraction was prescribed to the initial planning target volume (PTV), which was located in the cerebellum. Thermoluminescent detectors (TLDs), namely, Li-7-enriched LiF: Mg, Ti (MTS-7) type, were used to detect gamma radiation, which is produced by nuclear reactions, and secondary as well as recoil protons created out-of-field by secondary neutrons. Li-6-enriched LiF: Mg,Cu,P (MCP-6) was combined with Li-7-enriched MCP-7 to measure thermal neutrons. TLDs were calibrated in Co-60 and reported on absorbed dose in water per target dose (µGy/Gy) as well as thermal neutron dose equivalent per
target dose (µSv/Gy). Additionally, bubble detectors for personal neutron dosimetry (BD-PND) were used for measuring neutrons (>50 keV), which were calibrated in a Cf-252 neutron beam to report on neutron dose equivalent dose data. The Monte Carlo code TOPAS (version 3.6) was run using a phase-space file containing 10¹⁰ histories reaching an average standard statistical uncertainty of less than 0.2% (coverage factor $k = 1$) on all voxels scoring more than 50% of the maximum dose. The primary beam was modeled following a Fermi–Eyges description of the spot envelope fitted to measurements. For the Monte Carlo simulation, the chemical composition of the tissues represented in ATOM was employed. The dose was tallied as dose-to-water, and data were normalized to the target dose (physical dose) to report on absorbed doses per target dose (mSv/Gy) or neutron dose equivalent per target dose (µSv/Gy), while also an estimate of the total organ dose was provided for a target dose of 50.4 Gy(RBE). Out-of-field doses showed absorbed doses that were 5 to 6 orders of magnitude lower than the target dose. The discrepancy between TLD data and the corresponding scored values in the Monte Carlo calculations involving proton and gamma contributions was on average 18%. The comparison between the neutron equivalent doses between the Monte Carlo simulation and the measured neutron doses was on average 8%. Organ dose calculations revealed the highest dose for the thyroid, which was 120 mSv, while other organ doses ranged from 18 mSv in the lungs to 0.6 mSv in the testes. The proposed computational method for routine calculation of the out-of-the-field dose in proton therapy produces results that are compatible with the experimental data and allow to calculate out-of-field organ doses during proton therapy.

**Keywords:** proton therapy, anthropomorphic pediatric phantom, Monte Carlo simulation, out-of-field dosimetry, TLD, bubble detector

# 1 INTRODUCTION

Proton therapy (PT) enables to deliver highly conformed dose distributions owing to the characteristic Bragg peak and the finite range of protons. Nevertheless, PT is unavoidably accompanied by the production of secondary high-energy neutrons in the patient and structural materials of the beamline (1). Neutrons are of particular concern, as they are capable of traveling large distances to deposit out-of-field doses in organs located far from the primary treatment field and with a relatively high biological effectiveness (2). Furthermore, non-elastic nuclear reactions will also produce secondary protons, heavier ions, and gammas. As a result, the out-of-field radiation field in PT comprises a mixed field of radiation (including photons, neutrons, protons, and other charged particles) all with different potentials to induce biological damage. Moreover, the out-of-field radiation field, and hence the secondary dose delivered to healthy tissues, largely varies with position (close to field versus far) and depends on specific treatment parameters such as patient size and positioning, beam angles, proton energies, field size, modulation width, presence of range shifters (RSs), and the use of apertures.

The development of a validated Monte Carlo (MC) framework forms an important aspect in the assessment and characterization of out-of-field doses in PT. Nevertheless, the use of general-purpose MC simulations in out-of-field dosimetry is

often restricted to detector calibration, and it has highlighted important differences between MC codes and models (3). Moreover, the coupling of MC to advanced measurement and proper benchmarking of the MC codes and models are still unknown today. Once validated, MC simulations will allow to fully describe the out-of-field radiation field and foster accurate calculations of appropriate dosimetric quantities needed for the assessment of radiation damage and risks.

Out-of-field dosimetry is especially important for the radiation protection of children who might develop radiation-induced second primary tumors during their lifetime. Nowadays, the challenge for clinicians is to increase the survival rate while treating with fewer secondary effects. There is a critical need to understand the long-term health and quality of life (QoL) challenges in these populations and to assess the potential health effects of the treatment modalities to improve the survival and health of the patients.

Some medical physicists are cautious that the existing knowledge and understanding of the out-of-field doses and associated risk of inducing secondary malignant neoplasms (SMNs) is not sufficiently mature to justify the use of modern techniques, such as PT, for treating children or pregnant women (4). Therefore, a full characterization of the out-of-field doses, particularly at the PT field edge, requires special attention for the radiation protection and prevention of SMNs (5).
This study aims to set up and optimize a computational MC framework for out-of-field dosimetry in PT, through validation measurements with advanced dosimetry techniques. The study herein presented has been conducted in the framework of a European Horizon 2020 project, HARMONIC, which is addressed at improving the knowledge of the health effects of medical exposure during childhood. A central task within HARMONIC is to set up a cohort of pediatric patients treated with modern radiotherapy, including the computation of whole-body doses. Ultimately, the goal is to estimate the risk of late health effects (including the risk of second primary cancers) after pediatric radiotherapy exposures, which relies on organ dose estimation obtained from validated tools.

# 2 MATERIAL AND METHODS

## 2.1 Experimental Setup

For this study, an anthropomorphic phantom (ATOM, Computerized Imaging Reference Systems (CIRS), Inc., Norfolk, VA, USA) representing a 5-year-old child (type 705D) was used. The phantom consists of tissue equivalent (TE) materials, and 180 dosimeters can be inserted in different organ positions. For the insertion of bubble detectors, six tissue slabs were replaced by polymethyl methacrylate (PMMA) slabs manufactured at SCK CEN with dedicated inserts for this type of detector (see Figure 1).

Aiming to simulate a realistic proton treatment of a brain tumor, a clinically applied treatment plan was transferred to the conditions of the experiment. The corresponding patient should feature a cranial size and shape with a reasonable resemblance to the corresponding features of the anthropomorphic phantom. A 7-year-old female patient was selected with a diffuse midline glioma (WHO grade IV). The patient received combined radiotherapy and chemotherapy after R3 resection. A median dose of 50.4 Gy(RBE) with 1.8 Gy(RBE) per fraction was prescribed to the initial planning target volume (PTV), which

was located in the cerebellum and had a volume of 195.2 cm³. The treatment plan consisted of two ipsilateral oblique fields and a contralateral oblique field. The proton fields were delivered in a gantry room in pencil beam scanning (PBS) delivery mode employing a lucite RS with a physical thickness of 4.44 cm and a water-equivalent thickness of 5.14 cm. The concerned patient was enrolled in the prospective registry study "KiProReg" (German Clinical Trials Register: DRKS-ID: DRKS00005363) after consent from her legal guardians was obtained. This study was approved by the local ethics committee.

The treatment planning of the phantom case was conducted in the treatment planning system (TPS) RayStation (version 7, RaySearch Laboratories, Stockholm, Sweden). The X-ray CT image set was contoured in the following way: the boundaries between dissimilar TE material types were delineated; then the mass density of the volumes of the individual TE materials was overwritten with the values from the datasheet of the phantom. Furthermore, the volumes containing the thermoluminescent detector (TLD) inserts were contoured. A deformable image registration (6) was established between the clinical CT and the anthropomorphic phantom. Then the contours of the target volume and organs at risk (OARs) were mapped to the CT of the phantom. After that, the contour of the PTV was fine-tuned by alignment to the boundaries between brain substitute tissue and bone substitute tissue. The proton kinetic energies of the fields ranged from 100 to 165 MeV. The air gap between the RS and phantom surface was on average 11.0 cm regarding the minimum distance and an average of 14.6 cm regarding the distance on the central axis. The dose distribution was calculated with the MC-based dose engine (version 4.1) of RayStation. The objectives of the spot fluence optimization were similar to those of the clinical plan concerning, i.e., the right cochlea, the brainstem, and the PTV.

The delivered dose in the experiments was adjusted to the sensitivity and location of the corresponding detectors. The phantom irradiation plans for TLD measurements applied 100.8 Gy(RBE) to the PTV in a single fraction. For bubble

<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Pictures of the experimental setup. On the left are the slabs of the anthropomorphic phantom for insertion of TLDs including one of the designed PMMA slabs for BD-PNDs. On the right is the mounted 5-year-old anthropomorphic phantom (loaded with BD-PNDs) positioned in the gantry room for PT pencil beam scanning. TLDs, thermoluminescent detectors; PMMA, polymethyl methacrylate; BD-PNDs, bubble detectors for personal neutron dosimetry; PT, proton therapy.</figcaption></figure>
detectors, the delivered dose to the PTV ranged between 0.5 Gy (RBE) for close-to-field measurements and 6 Gy(RBE) for far out-of-field positions. Large changes in the dose level by modification of the overall number of monitor units (MU) are impossible in PBS because the allowed MU per beamlet ("spot") is subject to machine limitations. Thus, dedicated optimizations of the spot fluence were conducted per targeted dose level. This also included an adaptation of the spot spacing. As a result, the shape of the corresponding dose distributions was not exactly identical. This concerns, however, only the high dose region, which was not the subject of the current study.

The experiments were conducted in the West German Proton Therapy Centre Essen (WPE), which is based on the ProteusPlus proton machine (IBA PT, Louvain-La-Neuve, Belgium). Protons were accelerated in an isochronous cyclotron and subsequently slowed down in the energy selection system to adapt to the required range in the phantom. The protons were guided to a gantry-mounted, evacuated nozzle, which operated in a spot-by-spot type PBS delivery mode. The RS was mounted in a snout holder, which could be moved along the central beam axis with a linear translation stage. The cranial part of the phantom was put on a BoS Headframe (Qfix, Avondale, PA, USA), which in turn was attached to a short Patlog table (IBA PT, Schwarzenbruck, Germany). The phantom was aligned prior to the mock treatment with the lasers of the positioning system. The uncertainty was about ±2 mm. Although it could have been reduced with the X-ray-based verification system, this was not done to avoid a contribution of X-rays to the detector signal.

## 2.2 Dosimetry Systems

### 2.2.1 Thermoluminescent Detectors

TLDs, produced by IFJ-PAN (Krakow, Poland), are small cylindrical chips with a diameter of 4.5 mm and a height of 0.9 mm. Detectors of Li-7-enriched LiF : Mg, Ti (MTS-7) type were used. During PBS therapy, MTS-7 mainly detects gamma radiation, which is produced by nuclear reactions, and primary or secondary as well as recoil protons created out-of-field by neutrons. MTS-7 sensitivity to neutrons is very limited. However, Li-6-enriched LiF : Mg,Cu,P (MCP-6) detectors are very sensitive to thermal neutrons due to their high $^6$Li(n,$\alpha$)$^3$H cross section for thermal neutrons. MCP-6 was used in combination with Li-7-enriched MCP-7 to quantify thermal neutrons.

TLD detectors were read in a Thermo Scientific Harshaw 5500 reader following a preheat for 30 min at 120°C to avoid the effects of signal fading and low-temperature anomalies in the glow curves (Parisi 2018). A heating rate of 10°C/s was used to heat up TLDs up to 340°C for MTS type and up to 255°C for MCP type.

TLDs were calibrated with a Co-60 source in terms of kerma "free in air" ($K_{air}$), which was then converted to absorbed dose to water ($D_W$) using conversion factor $D_W/K_{air} = 1.12$ determined by the ratio of the mass energy absorption coefficient for water to air for the energy of Co-60 (7). MTS-7 data were expressed in absorbed dose in water per target dose (physical) [µGy/Gy]. To quantify the thermal neutron dose, the data from MCP-7 were subtracted from MCP-6. Next, we applied $K_{air}$ to neutron dose equivalent conversion coefficients for thermal neutrons (1.24 *

$10^{-02}$ mSv/mGy), as described in (8). However, it should be noted that this conversion coefficient has been reported to have an uncertainty of up to a factor 2, which is related to the uncertainties on energy and angular distribution of neutrons as well as on the light collection of the TLD reader (8). Finally, data were expressed as thermal neutron dose equivalent data (µSv), which were normalized to the physical target dose (Gy) and expressed as [µSv/Gy]. The estimated uncertainty was 100% (8).

Uncertainties on the out-of-field absorbed dose were assessed considering dosimeter reproducibility (1.8%), batch reproducibility (1.9%), Co-60 calibration uncertainty (2.4%) (9), and background uncertainties, which were dependent on the measured dose and reached up to 11% (coverage factor $k = 1$) for the farthest positions. For the energy response of MTS-7 detectors, the energy dependence for both photons (10) and protons (11) was considered. For the energy dependence of photons, a characteristic spectrum was used in the study of (12). Photon energies ranged between 30 keV and 10 MeV. Assuming a flat energy dependence above 1 MeV (10, 13), the calculated uncertainty was below 1% ($k = 1$). For the proton energy dependence, a uniform distribution over the proton energies was assumed for energies up to the maximum proton energy used in this study (165 MeV). This resulted in an uncertainty on the proton energy response of 5% ($k = 1$).

### 2.2.2 Bubble Detectors

Bubble detectors for personal neutron dosimetry (BD-PNDs) (Bubble Technology Industries, BTI, Chalk River, ON, Canada) were used to measure neutrons of energies above 50 keV. These cylindrical detectors are 15 cm in length and 2 cm in diameter, but the sensitive part, where bubbles are created, is only 7 cm in length and 1.6 cm in diameter (see **Figure 1**). BD-PNDs were calibrated with a Cf-252 source to obtain neutron dose equivalent by applying fluence-to-dose equivalent conversion factors derived from kerma factors $k(E)$ and a quality factor as a function of neutron energy ($Q(E)$) for ICRU tissue, as described previously (14). The final data were expressed in neutron dose equivalent per target dose (physical) [µSv/Gy]. Uncertainties of BD-PNDs are estimated to be on average 20% ($k = 1$).

## 2.3 Monte Carlo Framework

The well-established Geant4 (15–17) wrap-up MC code TOPAS v3.6 (Geant4) (18), in conjunction with the Matlab (The Mathworks, Inc., Natick, MA, USA)-based matRad v2.10.1 (19) project to create a DICOM-based dose verification system, was used to simulate the out-of-field absorbed dose distribution. For this purpose, matRad was extended by including the possibility to process DICOM RTIon files. With this feature, it was then possible to create the TOPAS input files with the treatment room-specific radiation parameters employing matRad as the TOPAS syntax parser. The simulations for the determination of the neutron equivalent dose at a point and the proton and gamma out-of-field dose could then be conducted.

### 2.3.1 Beam Model

To simulate the anthropomorphic phantom irradiation, it was necessary to run TOPAS simulations that reproduce the
commissioned beam. Mean energy and spread have been adjusted to reproduce the measured depth dose curves following the methods of (20, 21), and (22) in 5 MeV steps from the lowest energy available (that is, 100.0 MeV) up to the highest energy available (that is, 226.7 MeV). The commissioned beam data for the corresponding reference values were yielded by measurements with the plane parallel Bragg peak chamber (PTW, Freiburg, Germany) (23). Simulated and measured depth doses agreed within ±0.01 cm at R₈₀. In addition, the Fermi-Eyges parameters from the beam model implemented in RayStation were used to fully characterize the proton pencil beam.

Furthermore, an MU/ion calibration was performed to determine the number of protons in TOPAS corresponding to the respective MUs. For this purpose, reference fields, consisting of 1,681 spots with 0.25-cm spacing arranged in a symmetrical square around the isocenter, were simulated with 5 × 10⁵ protons per spot for the 27 different energies that make up the beam parameter database. The protons started 50 cm upstream of the isocenter at the nozzle exit of the treatment head. The Fermi-Eyges parameters were back-projected in vacuum to the nozzle exit (24–26). The method of (27) was applied to obtain the spot positions of the protons, taking into account the deflection of the protons at the two foci from the scanning magnets of the pencil beam dedicated nozzle. Downstream, a water tank with a volume of 50 × 50 × 50 cm³ was created. The isocenter was at 3-cm depth of the water tank. Simulation conditions thereby correspond to the conditions for beam-monitor calibration at the WPE, generally following reference dosimetry according to TRS 398 and DIN 6801-1 (28, 29). In the simulations, a cylindrical tally with a diameter of 1 cm and a thickness of 0.5 cm was then placed at the isocenter. The following physics models were employed (30): *g4em-standard_opt4*, *g4h-phy_QGSP_BIC HP*, *g4decay*, *g4ion-binarycascade*, *g4h-elastic HP*, and *g4stopping*. In accordance with (31), the mean excitation energy of G4 water was set to 78 eV. The density was set to 1 g/cm³. The overall allowed maximum step size for the condensed history algorithm in TOPAS was set to 0.1 cm. The production cut for all secondary particles was set to 0.05 cm. No variance-reduction techniques were employed. These simulations reached a standard statistical uncertainty of less than 0.6%.

The MU/ion calibration was then performed employing the simulation result of the tally in dose [Gy/MU] and using the dose meter set in [Gy(RBE)/MU]. The calibration yielded was then added to the beam parameter database so that the treatment room-specific machine parameter file included 27 mean energies, energy spreads, Fermi-Eyges parameters, and the number of protons per single MU.

## 2.3.2 RT Integration in TOPAS

The matRad code was extended to read DICOM RTIon files including both RT Plan and RT Struct. The RS used was considered in the matRad configuration files with the 4.44-cm-thick material lucite comprising a mean excitation energy of 74.0 eV, a density of 1.19 g/cm³, and material composition of 8.05% H, 59.98% C, and 31.96% O, corresponding to the material definition as included in the PSTAR database given by the

National Institute of Standards and Technology (NIST). Since RayStation specifies spot positions at the isocenter, these were back-calculated analogously to the procedure for MU/ion calibration. Based on the machine parameter file, the required beam data were linearly interpolated starting from the energies given by the RayStation RT Plan file for the individual energy layers. The DICOM CT images were incorporated into TOPAS using the *TSImageCube* function. The device and scan protocol-specific density corrections were applied, as well as the full Schneider model comprising 25 different stoichiometric tissues (32). The grid size of the inserted CT was in accordance with the one employed in RayStation, that is, 0.2 × 0.2 × 0.2 cm³.

The matRad and TOPAS build was validated with a simulation in which the absorbed dose to water was tallied and compared with the RayStation simulation by means of a 3D gamma test. The 10⁹ simulated histories were distributed based on the respective MU weights per spot. The simulations were run on four Intel® Xeon® CPU E5-2670 v3 @ 2.30GHz (48 cores) with 64.0 GB RAM each. The three applied fields were divided into four runs each so that a total of 12 individual simulations were necessary to obtain the full dose distribution. The number of histories relative to the distribution of MUs for each spot was sufficient to avoid undersampling by more than 0.01%. The results from the simulations were then obtained by summing the individual runs of the respective fields. Since the number of histories corresponding to the median dose of 50.4 Gy(RBE) is known, it can be multiplied by the appropriate factor to obtain a comparable dose within the matRad-based analysis of the results.

## 2.3.3 Out-of-Field Dose Calculations

According to (33) and (34), for the simulation of the neutron equivalent dose at a point and the gamma dose, respectively, the geometry of the treatment room and the gantry pit must be implemented in the MC code in order to obtain a comprehensive neutron spectrum including thermal neutrons and a gamma dose to water as accurately as possible. A fully rotating gantry around the isocenter of the treatment room was modeled in TOPAS. The geometrical models of the scanning and bending magnets, as well as of the counterweight, were simplified. The geometry also included a rotating table, a maze, a rolling floor, and a gantry pit including all walls, ceilings, and floor with the corresponding materials of the treatment room and gantry pit. The simulation environment thus had a volume of 20.0 × 7.0 × 20.0 m³. The world material in these simulations was altered from vacuum to air to account for ionization occurring in the air relevant for thermal neutrons and gammas. The physical properties of air were modeled according to (31). However, to ensure that the initial protons hit the upstream side of the RS as accurately as possible with respect to the Fermi-Eyges parameters, the initial proton transport started in a vacuum box with the lateral dimensions of the RS, and a longitudinal extension of 50 cm upstream of the isocenter to the upstream side of the RS. The additional energy loss due to air scattering downstream the RS, up to the surface of the anthropomorphic phantom, was accounted for *via* a slight increase in the mean energy of each spot. This increase was calculated by an interpolation employing the data provided in the PSTAR database for the continuous
slowing down approximation in air based on a logarithmic cubic spline fit and the predicted track length of the protons. The influence of the additional air scattering on the energy spread is negligible, as well as the number of initial protons absorbed in air subsequent to the interaction in the RS.

The out-of-field dose simulations used 10¹⁰ primary protons. The simulation parameters were set analogously to those described above for the MU/ion calibration. The gamma and proton dose to water, as well as the neutron dose equivalent at a point, were tallied. To determine the neutron dose equivalent at a point, fluence-to-dose conversion factors within the TOPAS *AmbientDoseEquivalent* scorer were employed based on a logarithmic energy binning. The built-in conversion factors were adopted to the appropriate fluence-to-neutron dose equivalent conversion factors as described in the study from (14) by using tissue kerma factors *k(E)* and quality factors as a function of neutron energy (*Q(E)*). In addition, simulations under identical conditions were performed in which only the dose to water of secondary protons was scored to elaborate on the effect of high-energy protons close to the irradiation field.

The spatial grid size in the TOPAS simulations was chosen analogously to the grid size in the RayStation simulations, namely, 0.2 × 0.2 × 0.2 cm³. The doses for the neutron dose equivalent at a point, as well as the gamma and proton dose to water, were obtained from reading out the RTStructs of the contoured TLD positions via averaging of the voxels comprising the respective areas.

## 2.4 Calculation of Total Dose Equivalent and Organ Doses

Following validation of the MC framework, TOPAS simulations were used to calculate the total dose equivalent. The contribution of the following particles was considered for the out-of-field

dose, and their contributions were added to obtain the total dose equivalent:

1. Primary protons and assuming a generic RBE of 1.1.

2. Secondary gammas considering an RBE=1.

3. Neutrons as calculated according to the method described in Section 2.3.3.

Finally, organ doses were calculated relying on the specific locations of TLD inserts in the CIRS phantom, which correspond to certain organs to allow for organ dosimetry. In a total of 28 organ doses per target, the dose was calculated as well as the total dose to the child's organs considering a total target dose of 45.8 Gy (physical dose), i.e., 50.4 Gy(RBE).

# 3 RESULTS

## 3.1 Measured Out-of-Field Doses

### 3.1.1 MTS-7 Thermoluminescent Detectors

Absorbed dose in water per target dose, as measured with MTS-7 detectors, revealed doses ranging from 2,842 ± 181 to 7.9 ± 0.5 µGy/Gy at respectively 7.8 and 50 cm from the isocenter (see **Figure 2**). Compared to the target dose, the out-of-field dose was more than 2 orders of magnitude lower close to the field and decreased to 5 orders of magnitude lower doses beyond 35 cm out-of-field. MTS-7 results report on the absorbed dose in water from non-neutron contributions, which are dominated by protons and also gammas that contribute to their signal. Close to the field, one can expect a larger contribution of primary protons, while further away from the field, these primary protons will not be measured, as recoil protons will dominate the field.

<figure><img src="image_4.png" /></figure>

**FIGURE 2** | Overview of experimental data. The absorbed dose in water per target dose [µGy/Gy] is plotted as a function of distance for MTS-7 (left figure). The uncertainty bars in this figure are not displayed since they are smaller than the symbol size. In the right figure, the neutron dose equivalent data are plotted as a function of distance for BD-PND and MCP-6 in combination with MCP-7. Uncertainty bars ($k = 1$) are plotted for all BD-PNDs and for MCP-6/MCP-7 for every 5th data point to maintain readability of the plot. Notice that the abscissas and ordinates axes of both figures are the same for comparison purposes.
### 3.1.2 MCP-6 and MCP-7 Thermoluminescent Detectors

Thermal neutron dose was obtained by subtracting MCP-7 from MCP-6 detectors. The neutron dose equivalent per target dose ranged between 20.3 and 0.08 µSv/Gy for 7.7 to 50 cm from the isocenter, respectively (see **Figure 2**). It should be noted that the uncertainty ($k = 1$) on the calibration factor is 100%, which is due to a number of different contributions as described in (8). These data can also be expressed in gamma equivalent neutron doses (not plotted), as it is sometimes referred to by other groups to quantify the thermal neutron dose in terms of gamma dose equivalent. This easy approach is obtained by MCP-6 minus MPC-7 doses as calibrated in Co-60 dose in water ($D_W$). Hence, gamma equivalent neutron doses in this study ranged between 1846.0 and 6.9 µGy/Gy.

### 3.1.3 Bubble Detector for Personal Neutron Dosimetry

Finally, it was possible to obtain the neutron dose equivalent (neutron energies greater than 50 keV) in 6 different positions. Bubble detector data doses were between 915 ± 183 and 16 ± 3 µSv/Gy for respectively 6.7 cm and 43 cm from the isocenter (see **Figure 2**). When comparing these results with the data obtained

by MCP-6 in combination with MCP-7, as described in Section 3.1.2, it was noted that the thermal neutron dose equivalent data were much lower compared to the neutron dose equivalent measured by BD-PNDs (neutron energies greater than 50 keV). More specifically, the thermal neutron dose contribution to the total neutron dose was only 2.5% close to the field and decreased to 0.9% at 50 cm.

## 3.2 TOPAS Simulations

### 3.2.1 Target Dose Distribution

**Figure 3** shows the result of the simulated absolute absorbed dose distribution from TOPAS. The dose to water was scored in order to compare to the RayStation results employing a 3D gamma test. On the lower left side in **Figure 3**, it can be seen that for a global 3D gamma test with passing criteria of 1%, 2 mm produces a gamma pass rate of 99.338%. This makes it evident that the implementation of the verification system was successful. To ensure that the differences between the overwritten materials of the anthropomorphic phantom in RayStation and the Schneider model did not cause discrepancies, the $R_{80}$ range of each simulated irradiation field from TOPAS was analyzed with the corresponding $R_{80}$ ranges of the fields from RayStation, where discrepancies were smaller than ±0.02 cm.

<figure><img src="image_5.png" /></figure>

**FIGURE 3** | Comparison between the absorbed dose distribution calculated by TOPAS (top left) and the absorbed dose distribution calculated by RayStation (top right), both shown in axial plane. A 3D gamma test was performed once for a global gamma of 1%, 2 mm (bottom left), and for a local gamma of 1%, 1 mm (bottom right), both excluding the dose outside the anthropomorphic phantom so as not to bias the test result.
On the lower right side in **Figure 3**, it can be seen that in a local 3D gamma test with 1%, 1 mm, the obtained gamma pass rate was 42.565%. It is recognized that most of the deviations occur outside the prescription range in this sensitive gamma test.

### 3.2.2 Out-of-Field Dose Simulations

In the top part of **Figure 4**, results demonstrate absorbed doses from different contributors including protons, all protons as well as secondary protons and gammas. Results show that the proton dose was approximately 2 mGy/Gy close to the field and drops to 1 µGy/Gy, at the very end of the phantom. Looking into the secondary protons, they match the total proton dose from 30 cm and higher, which points toward the fact that at only beyond 30 cm, primary protons will not contribute to the out-of-field doses. Secondary protons can be generated from the primary beam (which will be absorbed close to the field), but most likely they are created by neutrons as recoil protons. The gamma-induced dose was always lower than the proton dose and ranged on average between 1% and 54% of the proton dose close to the field and far out-of-field (50 cm), respectively. On the contrary, when neglecting the primary protons and only considering secondary protons, the gamma dose was on average 30% and 70% close to the field and far out-of-field, respectively. It should be noted that the uncertainty bars become wider the further the dose has been simulated out-of-field. This is related to the decreased particles in these regions hitting the small TLD volumes. To avoid biasing the data, no variance reduction techniques have been applied in the simulations.

Neutron dose equivalent considering all neutrons and only thermal neutrons are shown in the right plot in **Figure 4**. Data demonstrate a neutron dose equivalent of 528 ± 41 µSv/Gy at a 7.8-cm distance, which decreases to 11.2 µSv/Gy at 50 cm from the isocenter. Clearly, the contribution from thermal neutrons to the total is very small and remained below 1%. Assuming the

secondary proton production is mainly from recoil protons generated from neutrons, we also calculated the average quality factor from these data. We divided the calculated neutron dose equivalent data by the absorbed dose quantity from secondary protons and derived an average Q-factor of approximately 10, which did not vary significantly between locations out-of-field and is in line with literature data. Finally, the left plot in **Figure 4** also shows the total dose equivalent per target dose, which considers the contribution of the primary protons close to the field edge. At larger distances, it is clear that the total dose equivalent is dominated by the neutron dose, as the gamma dose does not contribute significantly to the dose equivalent (1%-6%).

### 3.3 Comparison of Out-of-Field Doses

In **Figure 5**, a comparison between the dose obtained from the MTS-7 detectors and TOPAS, from both protons and gamma contributions, revealed a good agreement. On average, the TOPAS doses, including proton and photon contributions, were 18% lower compared to those from the MTS-7 detectors with a slightly better agreement at larger distances. Nevertheless, the last 3 data points reveal a lower experimental dose compared to the TOPAS dose, which was within uncertainties, due to the larger uncertainties of the calculations at these positions as well as the higher uncertainties of the measurement points, as measured doses are closer to the background, and background uncertainties are 11%. In the right plot of **Figure 5**, the comparison between the experimental data obtained with BD-PNDs and simulated neutron dose equivalent data shows a good agreement close to the field and far from the field. Nevertheless, at 23 and 28 cm, the measured data were 50% lower compared to the simulated ones. One should keep in mind the use of PMMA slabs for the insertion of bubble detectors, which may have an impact on the out-of-field doses, particularly in regions where

<figure><img src="image_6.png" /></figure>

**FIGURE 4** | Overview of TOPAS results. Absorbed dose in water per target dose [µGy/Gy] is calculated for all protons, only secondary protons and gammas (left figure). Neutron dose equivalent data are plotted as a function of distance for all neutron energies considered, and when only considering the thermal neutrons (right figure). Uncertainty bars ($k = 1$) are plotted only once every 10th data point to maintain readability of the plot.
<figure><img src="image_7.png" /><figcaption>**FIGURE 5** | Comparison between experimental data and TOPAS. On the left, the absorbed dose in water per target dose [µGy/Gy] is plotted for MTS-7 measurements and TOPAS simulations summing proton and gamma doses. On the right, neutron dose equivalent data are plotted for BD-PND and TOPAS simulations of neutron doses. BD-PND, bubble detector for personal neutron dosimetry.</figcaption></figure>

the density of the phantom material should be lower such as the lung region. Moreover, bubble detectors have an uncertainty of 20% (k = 1), which does not include the uncertainty related to their energy response, as this is not known for energies above 20 MeV.

The measured thermal neutron doses are plotted in **Figure 6** together with the MC calculated thermal neutron doses. The experimental doses are higher than the simulated data. The difference is the largest close to the field where a 15-fold higher dose was measured, while out-of-field, the experimental data were a factor of 2 higher but within uncertainties.

## 3.4 Organ Dose Calculations

Dose calculations were grouped per organ in **Figure 7**, using a total target dose of 45.8 Gy (50.4 Gy(RBE)). The average thyroid dose per target dose was found to be 2,673 µSv/Gy corresponding to a total dose of 120 mSv. For organs in the chest region such as the lungs and thymus, average organ doses of 18 and 32 mSv were calculated for the total target dose, respectively. The breast dose was 17 mSv, while the heart dose was 8.3 mSv. For the liver and stomach, the obtained average doses were 4.1 and 3.4 mSv, respectively. Gonad doses were 1.1 and 0.6 mSv for the ovaries and testes, respectively.

<figure><img src="image_8.png" /><figcaption>**FIGURE 6** | Comparison of thermal neutron doses per target dose [µSv/Gy] for experimental data from MCP-6/MCP-7 data and TOPAS simulations. Uncertainty bars (k = 1) are plotted once every five points for clarity reasons.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**FIGURE 7** | Organ dose calculations produced by TOPAS at various TLD positions of the phantom. Data are grouped per organ, and the dose is reported as total dose equivalent per target dose [µSv/Gy] in lower abscissas and total dose equivalent for a target dose of 45.8 Gy in the upper abscissas. Uncertainty bars ($k = 1$) are given for each position.</figcaption></figure>

## 4 DISCUSSION

TOPAS was chosen in this study, as it is capable of calculating the dose distribution, including the aura consisting of neutrons and gammas (25, 35). Moreover, it was shown that TOPAS is able to simulate multiple Coulomb scattering in the lucite RS used in the present work with sufficient accuracy (26). As a first verification step of the MC framework, the in-field dose distribution was compared to RayStation, showing an excellent agreement in the 3D global gamma test as shown in **Figure 3**. However, the more sensitive local 3D gamma test revealed discrepancies outside the volume, which, as described in the work of (25), are due to the fact that RayStation simplifies the transport of secondary particles except for protons. For this reason, TOPAS was used to determine the out-of-field doses. Detailed modeling of the beam was performed, as well as detailed modeling of the geometry of the treatment room, including the gantry pit, to account for all possible sources of secondary

particles resulting from interactions between particles scattered and generated in the phantom, as well as the environment (33, 34).

A good agreement was observed between TLD measurements, MTS-7 type measuring the non-neutron component of the mixed radiation field, and TOPAS simulations tallying the proton and gamma component. A slight overestimation by MTS-7 detectors was however observed, which could be due to the slightly increased response of MTS detectors for protons (11). We assigned an uncertainty to this potential error of 5% ($k = 1$) in TLD data. Nevertheless, we did not use the spectra as input and assumed a uniform distribution of proton energies, which could be off for certain locations in the field, particularly when protons reach the end of their range. Nevertheless, we did not correct MTS-7 data for non-linearity in the energy response for protons and photons. Even though this can be considered a limitation, we do believe that the impact will be very small, and we considered it in the uncertainty of the detector data.
Moreover, the mentioned 12% overresponse of LiF-based TLDs for out-of-field measurements in photon beams does not apply to PT, as out-of-field spectra are much softer than in proton beams (36, 37). Furthermore, it should be noted that in the three outmost TLD positions, there is an underestimation of the dose with TLDs, which is most likely due to the very low doses measured in these positions and uncertainties associated with the background signal (11%).

Another reason for the slight overestimation of TLDs could be the contribution from heavy charged particles or fragments created by secondary neutrons in the phantom itself. This was not considered in the simulation due to the very low statistics of such particles but could have resulted in a signal in MTS-7 detectors. Still, heavy charged particles generated outside the TLD volume and created in the CIRS phantom could contribute, but due to their short range, the likelihood of reaching the TLD detector is small. Moreover, the sensitivity of TLD detectors for heavy charged particles is low, as the TL efficiency depends largely on the ionization density, as studied experimentally and using microdosimetric models (11, 38).

The use of Li-6-enriched LiF : Mg,Cu,P (MCP-6) detectors allowed us to assess the thermal neutron's contribution due to the high ⁶Li(n, α)³H cross section for thermal neutrons. However, it should be noted that the quantification of thermal neutrons and the use of appropriate conversion factors are subject to large uncertainties as previously described (8). The method assumes that the dose is deposited only by thermal neutrons, and neutrons are isotropic. Moreover, corrections are applied for superficial dose deposition in TLDs and, therefore, a decreased light attenuation. No information about the energy spectrum and angular distribution was obtained, and therefore, the calculated conversion factor has a large uncertainty (100%). In the current study, the experimental results are always higher as compared to the simulated data, which were the highest close to the field where a 15-fold higher dose was measured. This large discrepancy can be explained not only by the large measurement uncertainty but also by uncertainties in the TOPAS simulations related to the detailed modeling of beam, gantry, bed, walls, and other room components. As thermal neutrons are created by neutrons slowing down during collisions, it is very challenging to model these accurately. However, when moving out-of-field, the experimental data were a factor of 2 higher than simulated data, which could be explained by the more isotropic nature of the thermal neutrons when moving further away from the isocenter. Still, it was observed that the contribution of thermal neutrons to the total neutron dose equivalent is very small (within 1%). Similar findings, with thermal neutron doses contributing <1% to the total neutron dose equivalent, have also been observed in (39).

In general, the agreement between the measured and simulated neutron doses was good, despite the large uncertainties of the measurements performed with bubble detectors (20%, k = 1). Moreover, this uncertainty did not consider potential dependence on their response as a function of neutron energies. In fact, their response is well characterized until 20 MeV, above which response is not fully described. Nevertheless, previous data using BD-PNDs in PT have shown good agreement with $H_p$ (10) reference

measurements (39). Still, an underestimation was measured in some positions, which could be assigned to the use of PMMA slabs instead of the anthropomorphic materials containing tissue materials. The discrepancy is the largest, and up to 50% lower for measured data compared to simulated data, in the area where there is lung material. As the TOPAS model was based on the phantom's CT and did not model the PMMA slabs and the exact material densities during measurements, this could be the cause of the discrepancy. Although this could be considered a limitation of the study, the dose data, namely, tissue dose, as determined in the TOPAS simulations are highly relevant for translation to the clinic and epidemiology. Another reason for the discrepancy could come from the uncertainty in neutron dose simulations due to missing cross-section data above 20 MeV for which many codes need to rely on the use of nuclear models. Several models are available; however, it is still open which models are more suitable. The choice of the binary intra-nuclear cascade (BIC) model in GEANT4 within this study was based on previous data demonstrating a good agreement with experiments (25, 34, 40).

Once validated, the MC simulation framework allowed to calculate the total dose, which was assessed as the dose quantity *dose equivalent*, considering the biological effectiveness of the radiation, which is dependent on the radiation type and energy. In previous studies, usually, the total dose in out-of-field positions is calculated considering only the neutron and gamma dose, for example, in (12, 14). Our study, however, demonstrates that closer than 30 cm from the isocenter, the contribution from primary protons is significant. We therefore added this contribution by applying a generic RBE of 1.1 for protons. We did not calculate the proton's energy, and we do recognize that this calculation may be a simplification, as the RBE will depend on the proton's linear energy transfer (LET). Still, the RBE-LET relationship is under investigation and would require more extensive calculations of (micro)dosimetric quantities, which we considered out of the scope of this paper.

Unfortunately, the calculation of the total dose equivalent was not possible from the experimental detectors, as TLDs will measure, apart from the gamma dose, both the primary and secondary protons. When summing these to the neutron doses from BD-PNDs, this will lead to an overestimation due to the double counting of recoil protons. Research is ongoing on how to combine different detector systems, with various response functions, to overcome the challenges of mixed radiation fields and to allow for an accurate experimental measurement of the total dose equivalent in the future. Similar issues are encountered for dosimetry in space where also complex mixed radiation fields consisting of neutrons, photons, protons, and heavier ions exist (41). In space, often silicon telescopes and other spectroscopic devices are used in order to be able to separate the different radiation field components and to obtain an estimation of the total dose equivalent. However, such detectors are too bulky to be used in phantom measurements in PT.

Previous works mainly evaluated out-of-field doses in PT through experimental measurements in water phantoms (1) or anthropomorphic phantoms (42–44). Only few studies have modeled the PT beam in detail to allow MC calculations of out-
of-field doses in PT (45–47). Studies connecting MC and experimental data in PT are mostly limited to measurements with ambient monitors or Bonner sphere systems in the room (33, 34, 48). Such studies are lacking in phantom measurements, which are of utmost importance for patient care. Our study provides a first step in the development of an MC framework that allows us to fully characterize the out-of-field radiation field and eventually could lead to tools for dose and risk optimization in children.

It is important to make the framework less computationally demanding, as now a very detailed beam model is used. In the same way as it was done in the study from (34), it will be important to identify the origin of the secondary radiation in the beamline component and to allow simplification of the beam model. Therefore, in the next developments, it must be considered which components of the simulation the complexity can be reduced without changing the simulation results to such an extent that it no longer coincides with the present result within the standard statistical uncertainty. For this purpose, a traceable approach is to determine which part of the treatment room has the largest share in the secondary particle generation or scattering. In this way, the dimensions of the room can presumably be reduced, which is synonymous with a reduction in the simulation time. Likewise, the origin of the secondary radiation depending on the direction of flight and momentum of the protons escaping the nozzle, as well as the scattering in the body and phantom, needs to be analyzed with the scope of simplifying the applied irradiation field, which in turn reduces the effort required to generate the phase space files and ultimately eliminates the time-consuming step function feature with the equally time-consuming phase space sampling of each spot.

The closest comparison of our experimental data to literature could be made to two studies performed within EURADOS WG9 (42, 43), measuring out-of-field doses during PT in the same anthropomorphic 5-year-old phantom treated for a brain tumor (6-cm diameter). One study described the response of passive detector systems in PT out-of-field dosimetry (42), where no RS (no RS) was used, while another study verified the impact of using an RS or 3D-printed beam compensator (BC) on the out-of-field doses (43). At 12 cm, the neutron dose equivalent data were 120 µSv/Gy (no RS) in the study from (42) versus 130 µSv/Gy (BC) and 180 µSv/Gy (RS) in the study from (43). Our data reported a dose of 260 µSv/Gy, which could be due to the larger volume in the current study (195.2 cm³) compared to the previous studies (65 cm³) (42, 43). Interestingly, these studies also compared data to photon plans for the same phantom, tumor size, and location, revealing that intensity-modulated radiation therapy (IMRT) and 3D conformal radiation therapy (3D-CRT) (49) are at least one order of magnitude higher than PT at 30 cm. As part of another paper within this special issue, the same case was treated with IMRT and volumetric modulated arc therapy (VMAT), and we noted a reduction in out-of-field dose of a factor of approximately 5, close to the field for both IMRT and VMAT, while at 30 cm, the difference was a factor of 35 and 20 lower for PT as compared to IMRT and VMAT, respectively.

We reported on organ doses for a target dose of 50.4 Gy(RBE) by multiplying the normalized doses with the physical target

dose (45.8 Gy). Nevertheless, in the clinical treatment course, a second beam set delivered 3.6 Gy(RBE) in two fractions to a PTV, which was cut at the inferior side to protect the spinal cord. The corresponding reduction in the volume of the modified PTV amounted to 2.6%. Thus, the out-of-field contributions of the second beam set can be approximated by the corresponding values of the initial beam set.

Looking into the organ doses, the thyroid dose was the highest, yielding 120 mSv, while other organ doses ranged between 18 mSv for the lungs and 0.6 mSv for the testes. According to BEIR VII (50), the lifetime attributable risk (LAR) for cancer incidence and for an exposure at 5 years old is most elevated for the breast, lungs, and thyroid with LAR values of 914, 608, and 419/10⁵/0.1 Gy for women, respectively, while for men, the values were 261 and 76/10⁵/0.1 Gy for the lung and thyroid, respectively. Applying these risk factors to our data, we estimated a risk for secondary thyroid cancer of 0.6% for women and 0.1% for men. The risk for breast cancer was 0.2%. However, we should be aware that these risk models are mainly for low doses and low dose rates, and they cannot easily be extrapolated to radiotherapy where the dose is fractionated and organ dose is heterogeneous. These estimations should be considered with even more caution in the context of PT, as the effects of high-LET particles (i.e., protons and heavier ions) are outside the scope of the BEIR VII report.

Knowledge about potential long-term sequelae of treatment modalities needs precise data on the oncologic treatment, related to not only radiotherapies, such as the dose-volume histogram for every OAR (in the field and out of the field), but also the cumulative dose of every drug of chemotherapy (including new molecules and corticosteroids) and precise information on surgeries. This information needs to be complemented with a long period of follow-up.

The HARMONIC project was set up to provide direct evidence of the late health effects of low, moderate, and high radiation doses from modern radiotherapy techniques using protons or photons. Following up pediatric patients treated with PT will strengthen the epidemiological basis for assessing radiation risk in pediatric patients and will provide complementary information to the contribution from the large historical childhood survivor cohorts treated prior to 2000, which did not include new treatment modalities (51). The HARMONIC project therefore builds the infrastructure and instruments to evaluate the potential health, QoL, and social impacts of medical exposures to ionizing radiation in children, with potential for advanced patient-specific dose reconstruction, as presented here, and mechanistic investigations. It aims at providing evidence on the magnitude of possible cancer and non-cancer effects (including neurovascular, cardiovascular, and endocrine system effects), which may arise following cancer treatment with modern techniques, including PT in pediatrics.

# 5 CONCLUSION

As the role of proton beam therapy is continuously increasing, particularly when very young children are concerned, the
understanding of out-of-field doses and their impact on secondary
cancer induction is essential. Since the HARMONIC project aims to
investigate the incidence of secondary cancer, a reliable calculation
of the out-of-field dose is of crucial importance. In this framework,
the development of a validated MC system forms an important
aspect in the assessment and characterization of out-of-field doses
in PT. Once validated, MC simulations allow to fully describe the
out-of-field radiation, permitting calculations of appropriate
dosimetric quantities needed for the assessment of radiation
damage and risks. In this study, the coupling of MC to advanced
measurements with different detector types enabled the
performance of a proper benchmarking of a widely used MC
code, for use in out-of-field dosimetry.

The proposed computational method for calculation of the
out-of-the-field dose in PT produces results that are compatible
with the experimental data. The validated framework allowed a
detailed characterization of the radiation field and the calculation
of out-of-field organ doses during PT. The development of such
an MC framework could lead to tools for dose and risk
optimization in children.

DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be
made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS

MDS-H: experimental design, experimental setup, data analysis, and writing. NV: Monte Carlo simulation design, setup and execution, experimental setup, data analysis, and writing. CB:

REFERENCES

1. Stolarczyk L, Trinkl S, Romero-Expósito M, Mojzeszek N, Ambrozova I, Domingo C, et al. Dose Distribution of Secondary Radiation in a Water Phantom for a Proton Pencil Beam—EURADOS WG9 Intercomparison Exercise. *Phys Med Biol* (2018) 63:085017. doi: 10.1088/1361-6560/aab469

2. Ottolenghi A, Baiocco G, Smyth V, Trott K, Consortium A. The Andante Project: A Multidisciplinary Approach to Neutron Rbe. *Radiat Prot Dosim* (2015) 166:311–5. doi: 10.1093/rpd/ncv158

3. De Saint-Hubert M, Farah J, Klodowska M, Romero-Expósito MT, Tyminska K, Mares V, et al. The Influence of Nuclear Models and Monte Carlo Radiation Transport Codes on Stray Neutron Dose Estimations in Proton Therapy. *Radiat Meas* (2022) 150. doi: 10.1016/j.radmeas.2021.106693

4. Newhauser WD, Durante M. Assessing the Risk of Second Malignancies After Modern Radiotherapy. *Nat Rev Cancer* (2011) 11:438–48. doi: 10.1038/nrc3069

5. Harrison RM, Ainsbury E, Alves J, Bottollier-Depois JF, Breustedt B, Caresana M, et al. EURADOS Strategic Research Agenda 2020: Vision for the Dosimetry of Ionising Radiation. *Radiat Prot Dosim* (2021) 194:42–56. doi: 10.1093/rpd/ncab063

6. Weistrand O, Svensson S. The ANACONDA Algorithm for Deformable Image Registration in Radiotherapy. *Med Phys* (2015) 42:40–53. doi: 10.1118/1.4894702

7. Hubbell JH. SMS. NIST Standard Reference Database, Vol. 126. (2004).

8. Van Hoey O, Parisi A. Development and Validation of a Model for Assessing Neutron Fluence With Lithium Fluoride Thermoluminescent Detectors. *Radiat Phys Chem* (2021) 188:109688. doi: 10.1016/j.radphyschem.2021.109688

treatment planning, experimental setup, and editing. JE: Monte
Carlo geometry coding. JW: preparation of Monte Carlo
configuration files and data analysis. RN: experimental setup.
OVH: experimental analysis. JD: Monte Carlo simulations. FS:
Monte Carlo data analysis and experimental setup. FV:
simulation validation, data analysis, and writing—review. SR:
analytical model and writing—review. GB: experimental setup
writing—review. US: simulation validation and writing—review.
MR: Monte Carlo data analysis. BT: clinical analysis, writing—
review, and methodology. IT-C: epidemiological review, writing
—review, and methodology. LB: conceptualization, supervision,
writing, and writing—review. All authors contributed to the
article and approved the submitted version.

FUNDING

The presented research has been funded by the HARMONIC
project. The HARMONIC project (Health effects of cArdiac
fluoRoscopy and MOderN radIotherapy in paediatrics) has
received funding from the Euratom research and training
program 2014-2018 under grant agreement No 847707. MR
acknowledges funding from the Sistema Nacional de
Investigación de Panamá. IT-C acknowledges support from the
Spanish Ministry of Science and Innovation and State Research
Agency through the “Centro de Excelencia Severo Ochoa 2019-
2023” Program (CEX2018-000806-S) and support from the
Generalitat de Catalunya through the CERCA Program.

ACKNOWLEDGMENTS

The authors are thankful to Niklas Wahl from the matRad project.

9. De Saint-Hubert M, De Angelis C, Knežević Z, Michalec B, Reniers B, Pyszka E, et al. Characterization of Passive Dosimeters in Proton Pencil Beam Scanning - A EURADOS Intercomparison for Mailed Dosimetry Audits in Proton Therapy Centres. *Phys Med* (2021) 82:134–43. doi: 10.1016/j.ejmp.2021.01.073

10. Parisi A, Dabin J, Schoonjans W, Van Hoey O, Mégret P, Vanhavere F. Photon Energy Response of LiF : Mg,Ti (MTS) and LiF : Mg,Cu,P (MCP) Thermoluminescent Detectors: Experimental Measurements and Microdosimetric Modeling. *Radiat Phys Chem* (2019) 163:67–73. doi: 10.1016/j.radphyschem.2019.05.021

11. Parisi A, Van Hoey O, Vanhavere F. Microdosimetric Modeling of the Relative Luminescence Efficiency of LiF : Mg,Ti (MTS) Detectors Exposed to Charged Particles. *Radiat Prot Dosim* (2018) 180:192–5. doi: 10.1093/rpd/ncx188

12. Jia SB, Hadizadeh MH, Mowlawi AA, Loushab ME. Evaluation of Energy Deposition and Secondary Particle Production in Proton Therapy of Brain Using a Slab Head Phantom. *Rep Pract Oncol Radiother* (2014) 19:376–84. doi: 10.1016/j.rpor.2014.04.008

13. Bruggmoser G, Saum R, Saum F, Gainey M, Pychlau C, Kapsch RP, et al. Correction Factors kE and kQ for LiF-TLDs for Dosimetry in Megavoltage Electron and Photon Beams. *Z Med Phys* (2015) 25:186–91. doi: 10.1016/j.zemedi.2014.03.010

14. Romero-Expósito M, Domingo C, Sánchez-Doblado F, Ortega-Gelabert O, Gallego S. Experimental Evaluation of Neutron Dose in Radiotherapy Patients: Which Dose? *Med Phys* (2016) 43:360–7. doi: 10.1118/1.4938578
15. Agostinelli S, Allison J, Amako K, Apostolakis J, Araujo H, Arce P, et al. Geant4—a Simulation Toolkit. *Nucl Instrum Meth A* (2003) 506:250–303. doi: 10.1016/S0168-9002(03)01368-8

16. Allison J, Amako K, Apostolakis J, Araujo H, Arce Dubois P, Asai M, et al. Geant4 Developments and Applications. *IEEE T Nucl Sci* (2006) 53:270–8. doi: 10.1109/TNS.2006.869826

17. Allison J, Amako K, Apostolakis J, Arce P, Asai M, Aso T, et al. Recent Developments in Geant4. *Nucl Instrum Meth A* (2016) 835:186–225. doi: 10.1016/j.nima.2016.06.125

18. Perl J, Shin J, Schümann J, Faddegon B, Paganetti H. TOPAS: An Innovative Proton Monte Carlo Platform for Research and Clinical Applications. *Med Phys* (2012) 39:6818–37. doi: 10.1118/1.4758060

19. Wieser H, Cisternas E, Wahl N, Ulrich S, Stadler A, Mescher H, et al. Development of the Open-Source Dose Calculation and Optimization Toolkit Matrad. *Med Phys* (2017) 44:2556–68. doi: 10.1002/mp.12251

20. Bortfeld T, Bray T. An Analytical Approximation of the Bragg Curve for Therapeutic Proton Beams. *Med Phys* (1997) 4:2024–33. doi: 10.1118/1.5981116

21. Clasie B, Depauw N, Fransen M, Gomà C, Panahandeh HR, Secol J, et al. Golden Beam Data for Proton Pencil-Beam Scanning. *Phys Med Biol* (2012) 57:1147–58. doi: 10.1088/0031-9155/57/5/1147

22. Grassberger C, Lomax A, Paganetti H. Characterizing a Proton Beam Scanning System for Monte Carlo Dose Calculation in Patients. *Phys Med Biol* (2015) 60:633–45. doi: 10.1088/0031-9155/60/2/633

23. Bäumer C, Koska B, Lambert J, Timmermann B, Mertens T, Takoukam Talla P. Evaluation of Detectors for Acquisition of Pristine Depth-Dose Curves in Pencil Beam Scanning. *Med Phys* (2015) 16:151–63. doi: 10.1120/jacmp.v16i6.5577

24. Gottschalk B. Techniques of Proton Radiotherapy: Transport Theory. *arXiv* (2012) 1204:4470.

25. Verbeek N, Wulff J, Bäumer C, Smyczek S, Timmermann B, Brualla L. Single Pencil Beam Benchmark of a Module for Monte Carlo Simulation of Proton Transport in the PENELOPE Code. *Med Phys* (2021) 48:456–76. doi: 10.1002/ mp.14598

26. Verbeek N, Wulff J, Janson M, Bäumer C, Zahid S, Timmermann B, et al. Experiments and Monte Carlo Simulations on Multiple Coulomb Scattering of Protons. *Med Phys* (2021) 48:3186–99. doi: 10.1002/ mp.14860

27. Kimstrand P, Traneus E, Ahnesjö A, Grusell E, Glimelius B, Tilly N. A Beam Source Model for Scanned Proton Beams. *Phys Med Biol* (2007) 52:3151–68. doi: 10.1088/0031-9155/52/11/015

28. Andreo P, Burns D, Hohlfeld K, Huq M, Kanai T, Laitano F, et al. Absorbed Dose Determination in External Beam Radiotherapy: An International Code of Practice for Dosimetry Based on Standards of Absorbed Dose to Water. *IAEA Tech Rep Ser* (2006) pp 140–141.

29. DIN-Normenausschuss Radiologie (NAR). *Procedures of Dosimetry With Probe-Type Detectors for Proton and Ion Radiation - Part 1: Ionization Chambers*. DIN 6801-1. (2019).

30. Jarlskog C, Paganetti H. Physics Settings for Using the Geant4Toolkit in Proton Therapy. *IEEE Trans Nucl Sci* (2008) 55:1018–25. doi: 10.1109/ TNS.2008.922816

31. ICRU Report 90. Key Data for Ionizing-Radiation Dosimetry: Measurement Standards and Applications. In: *Reports of the International Commission on Radiation Units and Measurements*.

32. Schneider U, Pedroni E, Lomax A. The Calibration of CT Hounsfield Units for Radiotherapy Treatment Planning. *Phys Med Biol* (1996) 41:111–24. doi: 10.1088/0031-9155/41/1/009

33. Trinkl S, Mares V, Englbrecht FS, Wilkens JJ, Wielunski M, Parodi K, et al. Systematic Out-of-Field Secondary Neutron Spectrometry and Dosimetry in Pencil Beam Scanning Proton Therapy. *Med Phys* (2017) 44:1912–20. doi: 10.1002/mp.12206

34. Englbrecht FS, Trinkl S, Mares V, Ruhm W, Wielunski M, Wilkens JJ, et al. A Comprehensive Monte Carlo Study of Out-of-Field Secondary Neutron Spectra in a Scanned-Beam Proton Therapy Gantry Room. *Z Med Phys* (2021) 31:215–28. doi: 10.1016/j.zemedi.2021.01.001

35. Gottschalk B, Cascio E, Daartz J, Wagner M. On the Nuclear Halo of a Proton Pencil Beam Stopping in Water. *Phys Med Biol* (2015) 60:5627–54. doi: 10.1088/0031-9155/60/14/5627

36. Scarboro SB, Followill DS, Howell RM, Kry SF. Variations in Photon Energy Spectra of a 6 MV Beam and Their Impact on TLD Response. *Med Phys* (2011) 38:2619–28. doi: 10.1118/1.3575419

37. Kry SF, Bednarz B, Howell RM, Dauer L, Followill D, Klein E, et al. AAPM TG 158: Measurement and Calculation of Doses Outside the Treated Volume From External-Beam Radiation Therapy. *Med Phys* (2017) 44:e391–429. doi: 10.1002/mp.12462

38. Berger T, Hajek M. TL-Efficiency—Overview and Experimental Results Over the Years. *Radiat Meas* (2008) 43:146–56. doi: 10.1016/j.radmeas.2007.10.029

39. De Saint-Hubert M, Saldarriaga Vargas C, Van Hoey O, Schoonjans W, De Smet V, Mathot G, et al. Secondary Neutron Doses in a Proton Therapy Centre. *Radiat Prot Dosim* (2016) 170:336–41. doi: 10.1093/rpd/ncv458

40. Arce P, Bolst D, Bordage MC, Brown JMC, Cirrone P, Cortés-Giraldo MA, et al. Report on G4-Med, a Geant4 Benchmarking System for Medical Physics Applications Developed by the Geant4 Medical Simulation Benchmarking Group. *Med Phys* (2021) 48:19–56. doi: 10.1002/mp.14226

41. Benton E, Benton E. Space Radiation Dosimetry in Low-Earth Orbit and Beyond. *Nucl Instrum Meth B* (2001) 184:255–94. doi: 10.1016/S0168-583X (01)00748-0

42. Kneževic Z, Ambrozova I, Domingo C, De Saint-Hubert M, Majer M, Martínez-Rovira I, et al. Comparison of Response of Passive Dosimetry Systems in Scanning Proton Radiotherapy-a Study Using Paediatric Anthropomorphic Phantoms. *Radiat Prot Dosim* (2018) 180:256–60. doi: 10.1093/rpd/ncx254

43. Wochnik A, Stolarczyk L, Ambrožová I, Davidková M, De Saint-Hubert M, Domański S, et al. Out-Of-Field Doses for Scanning Proton Radiotherapy of Shallowly Located Paediatric Tumours-a Comparison of Range Shifter and 3D Printed Compensator. *Phys Med Biol* (2021) 66:035012. doi: 10.1088/1361- 6560/abcblf

44. Hälg RA, Besserer J, Boschung M, Mayer S, Lomax AJ, Schneider U. Measurements of the Neutron Dose Equivalent for Various Radiation Qualities, Treatment Machines and Delivery Techniques in Radiation Therapy. *Phys Med Biol* (2014) 59:2457–68. doi: 10.1088/0031-9155/59/10/ 2457

45. Taddei PJ, Mirkovic D, Fontenot JD, Giebeler A, Zheng Y, Kornguth D, et al. Stray Radiation Dose and Second Cancer Risk for a Pediatric Patient Receiving Craniospinal Irradiation With Proton Beams. *Phys Med Biol* (2009) 54:2259–75. doi: 10.1088/0031-9155/54/8/001

46. Geng C, Moteabbed M, Seco J, Gao Y, George Xu X, Ramos-Méndez J, et al. Dose Assessment for the Fetus Considering Scattered and Secondary Radiation From Photon and Proton Therapy When Treating a Brain Tumor of the Mother. *Phys Med Biol* (2015) 61:683–95. doi: 10.1088/0031- 9155/61/2/683

47. De Saint-Hubert M, Tymińska K, Stolarczyk L, Brkić H. Fetus Dose Calculation During Proton Therapy of Pregnant Phantoms Using Mcnpx and Mcnp6.2 Codes. *Radiat Meas* (2021) 149. doi: 10.1016/ j.radmeas.2021.106665

48. Leite AMM, Ronga MG, Giorgi M, Ristic Y, Perrot Y, Trompier F, et al. Secondary Neutron Dose Contribution From Pencil Beam Scanning, Scattered and Spatially Fractionated Proton Therapy. *Phys Med Biol* (2021) 66:225010. doi: 10.1088/1361-6560/ac3209

49. Majer M, Stolarczyk L, De Saint-Hubert M, Kabat D, Kneževic Z, Miljanic S, et al. Out-Of-Field Dose Measurements for 3d Conformal and Intensity Modulated Radiotherapy of a Paediatric Brain Tumour. *Radiat Prot Dosim* (2017) 176:331–40. doi: 10.1093/rpd/ncx015

50. National Research Council. *Health Risks From Exposure to Low Levels of Ionizing Radiation: BEIR VII Phase 2*. Washington, DC: The National Academies Press (2006). doi: 10.17226/11340

51. Winther JF, Kenborg L, Byrne J, Hjorth L, Kaatsch P, Kremer LCM, et al. Childhood Cancer Survivor Cohorts in Europe. *Acta Oncol* (2015) 54:655–68. doi: 10.3109/0284186X.2015.1008648

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of
the publisher, the editors and the reviewers. Any product that may be evaluated in
this article, or claim that may be made by its manufacturer, is not guaranteed or
endorsed by the publisher.

Copyright © 2022 De Saint-Hubert, Verbeek, Bäumer, Esser, Wulff, Nabha, Van Hoey, Dabin, Stuckmann, Vasi, Radonic, Boissonnat, Schneider, Rodriguez,

Timmermann, Thierry-Chef and Brualla. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author (s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.