Reproducibility study of Monte Carlo simulations for nanoparticle dose enhancement and biological modeling of cell survival curves

Nanoparticle-derived radiosensitization has been investigated by several groups using Monte Carlo simulations and biological modeling. In this work we replicated the physical simulation and biological modeling of previously published research for 50 nm gold nanoparticles irradiated with monoenergetic photons, various 250 kVp photon spectra, and spread-out Bragg peak (SOBP) protons. Monte Carlo simulations were performed using TOPAS and used condensed history Penelope low energy physics models for macroscopic dose deposition and interaction with the nanoparticle; simulation of the microscopic dose deposition from nanoparticle secondaries was performed using Geant4-DNA track structure physics. Biological modeling of survival fractions was performed using a local effect model-type approach for MDA-MB-231 breast cancer cells. Physical simulation results agreed extraordinarily well at all distances (1 nm to 10 μm from nanoparticle) for monoenergetic photons and SOBP protons in terms of dose per interaction, dose kernel ratio (often labeled dose enhancement factor), and secondary electron spectra. For 250 kVp photons the influence of the gold K-edge was investigated and found to appreciably affect the results. Calculated survival fractions similarly agreed well within one order of magnitude at macroscopic doses (i.e. without nanoparticle contribution) from 1 Gy to 10 Gy. Several 250 kVp spectra were tested to find one yielding closest agreement with previous results. This highlights the importance of a detailed description of the low energy (< 150 keV) component of photon spectra used for in-silico, as well as in-vitro, and in-vivo studies to ensure reproducibility of the experiments by the scientific community. Both, Monte Carlo simulations of physical interactions of the nanoparticle with photons and protons, as well as the biological modelling of cell survival curves agreed extraordinarily well with previously published data. Further investigation of the stochastic nature of nanoparticle radiosenstiziation is ongoing.


Introduction
Radiation therapy is widely used in the treatment of cancer with approximately half of all patients having indications for its use [1,2]. The effectiveness of radiation, or tumor control probability, however, is limited by its similarly toxic effects on normal tissue. The window of radiation doses that are sufficient to achieve some tumor control effect while minimizing normal tissue complications is called the therapeutic window. Widening of this window can be achieved by using modern radiation delivery techniques or by modifying radiation effects on a cellular level using radioprotectors for normal tissue [3] or enhancing radiation damage using radiosensitizers or combination therapy [4][5][6]. The latter has been investigated extensively, primarily using metallic nanoparticles (NP) with most investigators studying gold nanoparticles (GNP) due to their high biocompatibility, high mass energy absorption coefficient and atomic number. In-vitro and in-vivo studies have shown a radiation sensitization or enhancement effect with the use of GNPs [7][8][9][10]. Theoretical studies evaluating the dose enhancement on macroscopic scales, however, suggest that GNP concentrations of 5-20 mg g −1 , i.e. 25-100 mM of gold atoms, are necessary to achieve macroscopic dose enhancement factors not exceeding a factor of four [11,12]. These concentrations, however, are orders of magnitude larger than those used in the studies by e.g. Jain et al, who used 1.9 nm diameter GNPs with up to 20 μM concentration and found a radiosensitizing effect in MDA-MB-231 cancer cells [8].
Local dose enhancement on cellular scales and smaller has been studied by many investigators using Monte Carlo techniques [13][14][15][16][17][18][19], most of which used condensed history particle transport algorithms for computational efficiency. Results from these studies were used to derive models to describe and predict the radiosensitization in the presence of metallic NP [20][21][22][23][24]. Two groups [21,23] applied a simplified local effect model (LEM) [25] to the highly heterogeneous dose distribution around NPs on cellular scales to derive an average number of lethal lesions that result in a given survival fraction. Lin et al [17] used a threestage approach to simulate the microscopic contribution to absorbed dose distribution around GNPs with the last simulation stage performed using more precise and computationally intensive track structure simulations. Results from that study were used by that group to perform biological modeling for kilovoltage photons and spread-out Bragg peak (SOBP) protons [24].
Despite the number of publications on simulation and modeling of radiosensitization, the reproducibility of those studies has not been sufficiently addressed. It is the authors' opinion however that reproducibility is a key quality feature of published research and that both, confirmatory and dissenting studies are of scientific value. During the benchmarking of simulations for related research the authors found key information to be missing from existing studies which made an evaluation of the authors' simulation configuration and setup difficult. Hence, the goal of this study was to replicate and reproduce in parts the simulations and modeling carried out by Lin et al [17,24], as a sample of published studies in the field of nanoparticle radiosensitization. In this process we evaluated important information that should be included in reports of nanoparticle radiosensitization studies to allow for feasible validation and reproducibility of simulation results.

Materials and methods
2.1. Simulation of physical dose distributions Monte Carlo particle transport simulations were performed using TOPAS version 3.8.p1 [26], which is built against Geant4 version 10.07.p3 [27,28]. To efficiently model the processes occurring on macroscopic (millimeters and larger) to nanoscopic (nanometer) scales these simulations were performed in separate stages [17,18,20,29]: using condensedhistory physics for macroscopic dose deposition and interaction with the nanoparticle, and using trackstructure physics for microscopic dose distribution [30].
The simulation workflow described by Lin et alet al [17] was followed and applied to photons with energies from 50 keV to 250 keV, including different spectra mimicking a 250 kVp source, as well as a proton spread-out Bragg peak (SOBP) with 133 MeV maximum proton energy and 5 cm modulation width, which was created using time features [31]. Spectra mimicking a 250 kVp source were created using discrete photon energies ranging from 50 keV to 250 keV with different weighting to create one hard and three soft spectra, which were modified to include photon energies just below and above the gold K-edge (Table S1, figure S7).
For macroscopic simulations of average dose contribution per track, the scoring volume was a 50 mm diameter and 1 mm thick water cylinder, positioned at 2 mm depth in a water phantom for photons and at three depths (8 cm, 10 cm, 12 cm) within the proton SOBP. Additionally, for proton SOBP simulations the phase space of primary and secondary particles, including neutrons and heavy ions, was scored for subsequent nanoscopic simulations.
Nanoscopic simulations were performed with the 50 nm diameter spherical GNP in vacuum and irradiated by monoenergetic or a spectrum of particles. In addition to the energy absorbed by the nanoparticle, the phase space of secondary particles leaving the nanoparticle surface together with the number of tracks entering and leaving the nanoparticle without interaction were scored. These simulations were also repeated for a water sphere of the same dimensions to calculate the dose kernel ratio at distance r from the nanoparticle center.
For microscopic simulations, the phase space of secondary particles leaving the nanoparticle or water sphere was placed at the center of a 20 μm water volume (material: G4_WATER) and dose to water was scored in radial bins with 1 nm, 10 nm, and 100 nm bin widths for radii of 0 to 200 nm, 200 nm to 2 μm, and 2 μm to 10 μm. From this, the contribution to absorbed dose per primary track interaction was calculated using the interaction probability described in 2.1.1.
Both, macro-and nanoscopic simulations used condensed history Penelope low energy electromagnetic models, g4em-penelope, which in general performs similar to the other often used model, g4emlivermore [32]. Select simulations using monoenergetic photons were repeated for comparison using the Livermore low energy models. For protons, we further included g4h-phy_QGSP_BIC_HP, g4-elas-tic_HP, g4decay, g4stopping, and g4ion-binarycascade. Microscopic simulations used g4em-dna with all additional deexcitation processes activated. The production cuts were set to 1 μm in macroscopic and 1 nm in nanoscopic and microscopic simulations.

Interaction rate
The number of tracks traversing the macroscopic scoring volume of cross-section S to deposit 1 Gy of dose, N , 1Gy is given by the inverse of the average dose contribution per track. With this, the number of tracks traversing the nanoparticle of cross-section s is then ⋅ s S N , 1Gy which we call P , hit or hit probability [24]. The interaction or ionization probability P ion is defined as the ratio of tracks with at least one interaction in the volume to the total number of tracks traversing the volume, which is scored at the nanoscopic stage. Finally, the interaction rate per Gray deposited macroscopically is then given by = ⋅ P P P .
tot hit ion With this, changes in P hit are due to number of tracks necessary to deposit 1 Gy in water, whereas changes in P ion will be due to nanoparticle size and beam energy dependent interaction crosssections.

Survival curve modelling
Dose response curves were modeled following the LEM approach used by McMahon et al and others [20,24,33], specifically following the procedure outlined by Lin et al for photons and SOBP protons. Briefly, the number of lethal lesions is calculated locally for small volumes throughout the cell nucleus, is the linear quadratic model adjusted for overestimation at high doses exceeding the threshold dose D t [24]. The linear quadratic model parameters a and b are those obtained using low-LET radiation without nanoparticles. The average number of lethal lesions is then given by ( ) ; , d and the survival fraction by = -S e . N The microscopic dose distribution ( ) D x y z , , depends on the macroscopic dose D , M the number and distribution of nanoparticles inside the cell N , NP and interaction rate. For comparison and computational reasons this modeling was performed in 2D, assuming a finite but infinitessimally thin cell, which we will refer to as a LEM approach in 2D. The average number of lethal lesions is then calculated as an average over the area,

Simulation of physical dose distributions
The number of photon tracks needed to macroscopically deposit 1 Gy in a 50 mm diameter cylinder was highest for 50 keV with53 10 12 continuously decreasing to15 10 12 at 250 keV. For the generated 250 kVp (hard) spectrum it was29 10 , 12 while for the 250 kVp (soft) spectrum and those with 80 keV or 83 keV energies it was approximately33 10 . 12 For SOBP protons5 10 9 to11 10 9 tracks deposited 1 Gy depending on the position in the SOBP. The interaction probability P ion generally decreased with increasing photon energy from´-0.46 10 3 at 50 keV to´-0.04 10 3 at 250 keV. For 83 keV photons it was calculated to be´-0.53 10 . 3 For the distal and central point in the proton SOBP the value of P ion was 0.067 to 0.069, for the proximal point it was 0.052. Values for all beams as well as the total interaction rate per Gray are given in table 1.
The average dose per interaction from secondary particles leaving the nanoparticles is shown in figure 1. It is highest at the nanoparticle surface with 193 Gy to 310 Gy for photons and 51.2 Gy to 52.4 Gy for protons. At 100 nm from the nanoparticle center this reduces to 7.3 Gy to 13.6 Gy for photons and about 2.7 Gy for protons. At 1 μm distance this again reduces to 31 mGy to 47 mGy for photons and 6 mGy for protons. Values for all beams are given in table S2. At distances greater than 200 nm the dose per interaction from 100 keV photons is higher than the other beams, whereas the dose per interaction for 83 keV photons is noticeably higher closer to the nanoparticle surface  (table S2).
The dose kernel ratios (DKR) for photons (figure 3, left) ranged from 300-to 1600-times at the surface of the nanoparticle. For energies above 50 keV the DKR decreases with increasing distance from the nanoparticle with some intermittent increases at approximately 60 nm and 900 nm. For 50 keV the DKR increases after the first local maximum at 60 nm to reach a value of 1600 at distances beyond 2 μm. For protons the DKR increases proportional to the logarithm of the distance from the nanoparticle center from just above 2 at the surface to 13 at 1 μm from where it slowly increases towards 14 at 10 μm (figure 3 right, figure S2).
Differential secondary electron spectra per incident photon or proton fluence are shown in figure 4 for 200 keV photons (left) and SOBP protons (right). For both, the spectra for interactions with water are included. For photons, we further plot the spectrum of electrons from photoelectric proccesses only, which includes Auger electrons ('Au: PE'), and the spectrum of all secondary electrons when disabling or enabling Auger emission processes ('Au: Auger Off/On'). In the secondary electron spectrum from 200 keV photon interactions, one can identify the energy of the photoelectrons associated with various gold photo-absorption edges at approximately 120 keV (K-edge), 186 keV (L-edges) 197 keV (M-edges), as well as the Compton plateau below the Compton edge at 88 keV. The contribution of Auger electrons is included in energy bins 12 keV. These features can also be observed in the other photon spectra ( figure S4). For protons, the high energy edge is highest at over 200 keV for emission from GNP in the proximal portion and lowest in the distal portion of the SOBP. This high energy edge is consistent for both the WNP and GNP located in the center of the SOBP. The differential secondary electron spectra are similar between the different SOBP positions for GNP, but notably lower at all energies for WNP.

Survival curve modelling
Survival fractions (SF) were calculated using all sources with GNP distribution in the cytoplasm, nucleus, and cell for macroscopic doses in absence of GNP, D , M from 0 Gy to 10 Gy (figures 5, 6, S6) and tabulated for 2 Gy and 8 Gy (table S3). Briefly, without nanoparticles the survival fractions were 0.78 and 0.03 at 2 Gy and 8 Gy, respectively. For nanoparticles distributed only within the cytoplasm, SF(2 Gy) ranged from 0.53 to 0.78 for both photons and protons, while SF(8 Gy) ranged from 0.0004 to 0.028 for photons and was approximately 0.027 for all protons.
For nanoparticles distributed only within the nucleus, SF(2 Gy) was  -10 10 for 50 keV and 83 keV, -10 11 for 100 keV, approximately -10 6 for 80 keV and 150 keV, and approximately 0.1 for 200 keV and 250 keV. Using 250 kVp spectra SF(2 Gy) was lowest for the spectrum including 83 keV photons at´-1.7 10 , 7 followed by the soft spectrum, the spectrum including 80 keV photons, and finally the hard spectrum, each with approximately one order of magnitude difference. For protons the nucleus SF(2 Gy) was between 0.20 and 0.35 while SF(8 Gy) ranged betweeń -3.0 10 6 for at the proximal and´-8.5 10 5 at the distal end of the SOBP. Table 2 lists the survival fractions at 2 Gy and 8 Gy for SOBP protons and 250 kVp (w/80 keV) from this work as well as those from Lin et al [24]. Finally, survival curves for all beams are shown in figure 5 for 250 kVp spectra, figure 6 for SOBP protons, for all others the reader is referred to figure S6.

Discussion
In this work we investigated the physical interactions of GNPs with mono-and polyenergetic photon beams as well as SOBP protons, followed by modeling of cell survival curves assuming different cellular distributions of the nanoparticles. We followed the procedure outlined by Lin et al [17,24] and employed a multi- Table 1. Number of tracks traversing a 50 mm diameter cylinder that result in 1 Gy absorbed dose (N 1Gy ), interaction probability for a 50 nm GNP (P ion ), and total interaction rate (P tot ) for photons and SOBP protons. stage approach consisting of macroscopic, nanoscopic, and microscopic Monte Carlo particle transport simulations followed by biological modeling using a LEM-style approach in 2D.

Simulation of physical dose distributions
For photons, the decrease in number of tracks needed to macroscopically deposit 1 Gy in a 50 mm diameter cylinder of water with increasing energy is consistent with the maximum amount of energy that can be deposited by a photon (table 1). With the Compton effect becoming dominant in water at energies above 50 keV, a 250 keV photon can transfer as much as 124 keV to a Compton electron, while a 50 keV photon can only transfer 8.2 keV. With a higher probability for photoelectric effect with oxygen at lower energies, however, this translates into only a 3.5-fold difference between N 1Gy for 50 keV and 250 keV. For interactions with GNP, changes in photon energies in the kiloelectronvolt range profoundly affect the interaction probability, due to the higher absorption edges for photoelectric effect in gold. This results in rapidly decreasing values for P ion from 0.464 at 50 keV to 0.135 at 80 keV; with the gold K-edge at 80.7 keV, P ion increases sharply for 83 keV photons and continuously decreases thereafter to 0.035 for 250 keV photons. The impact of the K-edge can also be seen in the values for the soft 250 kVp spectra (soft, +80 keV, +83 keV): while the spectrum with 83 keV photons has the highest average energy of the three, its  Figure 2. Contribution to absorbed dose per interaction at up to 50 nm from the nanoparticle center for different generated 250 kVp spectra (left) and SOBP protons (right). The photon spectra were generated to include contributions from photons at 80 keV (blue/ circle), 83 keV (orange/triangle), and to be harder (green/star) or softer (square/red).
interaction probability exceeds that of the soft 250 kVp spectrum.
For protons, N 1Gy and P ion depend on the stopping power, which increases with decreasing proton energy and thus, with depth in the SOBP; higher stopping power corresponds to more energy transferred to delta electrons per proton. However, with the maximum energy transfer to delta electrons in a single collision limited to approximately 2.2% of the proton's kinetic energy, the maximum energy of delta electrons decreases with increasing depth in the SOBP (figure 4, right). This in turn leads to a higher stopping power for the delta electrons, which corresponds to an increased interaction probability with the GNP, explaining the moderate increase in P ion with depth from 0.52 to 0.69 between the proximal and distal end of the SOBP. The total interaction rate per 1 Gy macroscopic dose, P , tot follows the combined behavior of its factors, P ion and P , hit where the latter is the product of N 1Gy and the cross-section ratio, which is ⎛ ⎝ ⎞ ⎠ s S = = -50 nm 50 mm 10 .

12
The contribution to absorbed dose per interaction in the GNP is overall similar between different photon energies and between different positions within the SOBP (figures 1, 3). An exception are beams that include 83 keV photons where the dose at distances up to approximately 175 nm from the NP surface is increased ( figure 1 left, 3 left). While higher energy photons also undergo K-edge absorption, the energy  of the photoelectrons is lowest for K-edge absorption of 83 keV at about 2.3 keV, which has a higher electron stopping power and will thus, deposit more energy closer to the NP surface compared to higher energy photoelectrons [34]. This in turn explains both, the lower contribution to absorbed dose for 83 keV at distances closer to 10 4 nm due to the lack of higher energy photoelectrons and the increased contribution to absorbed dose for 100 keV at these distances, whose photoelectrons have a larger range of approximately 2.5 μm. Higher energy photon beams on the other hand produce photoelectrons with energies >20 keV with corresponding ranges in water >10 μm and lower electron stopping power, which thus, deposit less dose within the first 10 μm from the NP center. For protons the impact of reduced delta electron energy with increased depth in the SOBP can be observed, too (figure 2, right). However, this effect is less pronounced than for photons due to the comparatively high energy and thus range of delta electrons ( figure 4, right), making the dose per interaction effectively the same for all positions within the SOBP.
The dose kernel ratio for photons is large, ranging from 200 to 2000 at distances less than 1 μm from the NP center, which is due to the presence of Auger electrons at short distances and photoelectrons with intermediate to high energies ( figure 4, left). For photons the DKRs shows a strong energy dependence with lower DKRs at higher energies, except for 100 keV photons at short distances, which is due to the addition of K-edge photoelectrons compared to 50 keV. The DKR for 50 keV photons increases at distances past 1 μm due to the presence of gold L-shell photoelectrons with energies of approximately 35 keV [35]. A similar effect is expected for the higher energies, but only at distances past 10 μm due to the higher energy of their L-shell photoelectrons. These results are in line with those reported for 125 I, 169 Yb, and 192 Ir [15], whose gamma ray spectra have average energies of approximately 35 keV, 93 keV, and 380 keV.
For protons, the DKR is similar for all positions in the SOBP due to the similarity of NP secondary electron spectra: it increases linearly with the logarithm of distance from the NP up to 1 μm from the NP center where it plateaus at a value of 13 to 14 (figure 3 right, S3 right). The proton stopping power of gold is approximately 5 to 10 times higher than that of water in the range of 10 MeV to 50 MeV, meaning that more delta electrons are produced per track in gold compared to water, assuming similar delta electron energy distributions. The increase in DKR from the nanoparticle surface to about 1 μm (figure S3 right) may be explained by differences in the secondary electron spectra between the gold and water particles (figure 4, right). Differences between the spectra are within one order of magnitude up to approximately 125 keV, resulting in the limited increase in DKR. The higher energy and thus, longer range electrons from proton interactions with gold then result in the sharper DKR increase to about fourteen-fold thereafter.
Comparing Penelope and Livermore low energy physics models, the spectra of secondary electrons ( figure S4) are overall similar with only minor differences in the distributions between different photoelectron peaks and the Compton edge. While the contribution to absorbed dose per interaction follows the same overall behavior with beam energy and distance, the Livermore physics models lead to an increased contribution to absorbed dose per interaction of approximately 16% at short distances (up to 50 nm) from the NP ( figure S2).
Overall, the physical stage simulation results show good agreement with those by Lin et al [17,24] for monoenergetic photons and SOBP protons. The additional dose per ionization (figures 1, 2) was converted to additional dose per incident photon for monoenergetic photons and to additional dose per MeV absorbed in the NP for SOBP protons to allow for direct qualitative comparison (figure S1) with figures 4 and 5 from their work [17]. In addition to the dose per incident photon at the nanoparticle surface which was approximately 10 −1 Gy for 50 keV to just less than 10 −2 Gy for 250 keV, our results also include the 'cross-over' of the 50 keV and 100 keV dose curves at large distances from the NP surface (figure S1).
This agreement indicates that for GNP the method of obtaining the proton SOBP may be less important: Whereas we generated a SOBP using spot scanning of proton pencil beams with different energies, Lin et al utilized a commissioned Monte Carlo simulation of a clinically used passive scattering system, using an explicitly modeled modulator wheel to change the the proton energy [36]. Due to the induced scattering, the proton energy distribution at depth is usually different compared to of pristine Bragg peaks. This in turn would have the potential to affect P ion and P hit , though we expect it to be minimal given the good agreement of the physical results, and similarly expect the impact on the biological modeling to be limited. Furthermore, the secondary neutron fluence is substantially larger in a double scattering system compared to pencil beam scanning. Thus, if neutron interactions with the GNP contributed significantly to the changes in local energy deposit, we would have expected a substantial deviation of our work from that of Lin et al, indicating that for GNP charged particle interactions are dominant. This, however, may not be generalizable to all metals, especially those with high inelastic neutron cross-sections like gadolinium [37].
The DKRs for monoenergetic photons (figure 3) are similarly consistent with those reported by Lin et al: specifically the increase in the DKR for 50 keV above 10 3 at approximately 1 μm distance. This increase is due to a substantial decrease in radial dose around the WNP which occurs as electrons with maximal kinetic energy from Compton scattering reach the end of their range. For 50 keV the energy transferred to electrons is at most 8.2 keV, which corresponds to a range of less than 2.5 μm [34]. While this range limit also applies to Compton electrons from GNPs, more dose is deposited by gold L-shell photoelectrons (for 50 keV photons), which have a kinetic energy of approximately 35 keV [35], corresponding to a range in excess of 10 μm [34]. The spectra of secondary electrons (figure 4) are similarly consistent with those reported by Lin et al, noting that we are reporting the differential energy spectra per incident photon or proton fluence whereas Lin et al reported the secondary electron fraction which scales each individual curve such that the sum of all bins equals unity. For example, in their work this leads to an offset between the spectra from GNP interactions with and without Auger emissions included, most easily identified in an offset in the Compton plateau below 100 keV, which is unphysical and makes it difficult to compare the different spectra. Regardless, the impact of Auger emissions is clearly identified in both figures in the peaks at electron energies below 12 keV and at approximately 50 keV in the Compton plateau. It is important to point out that any Auger emissions secondary to photoelectric effects are created by the same process in Geant4 that also generates photoelectrons, and hence, are included in the 'Au: PE' spectrum. Additional contributions to the complete 'Au: AugerOn' spectrum include atomic deexcitation processes like particle induced x-ray emission (PIXE) [38].
In contrast to the agreement described above, the reproduction of the DKR for polyenergetic photon spectra was complicated by a missing description of those spectra. Attempting to generate these spectra yielded to DKR that were approximately one order of magnitude higher (figure S3, left) than those reported by Lin et al in their figure 6 [17]. This is due to the large impact of gold's K-edge at 80.7 keV, which we have discussed above using monoenergetic photon beams as well as different 250 kVp spectra. These mock spectra can be considered as to describe a lower and upper bound with respect to the radial dose around the nanoparticle from a 250 kVp source with 2 mmAl filtration, given that the radial dose from a SpekPy-generated spectrum of these parameters falls between the curves obtained using spectra with and without 83 keV photon contributions (figure S8) [39]. These differences in dose per interaction and DKR also translate into differences in the biological modeling of survival curves which will be discussed next. This points to the need for a detailed description of energy spectra to allow for validation and reproducibility of the experiments, not just for in-silico but also in-vitro and in-vivo studies using kilovoltage energies.

Survival curve modelling
The differences in dose per interaction and interaction rate per macroscopic dose discussed above translated into differences in survival fraction for the different beams (figures 5, 6, S6). As expected, higher additional dose per interaction at short distances from the NP resulted in greatly reduced survival fractions even at low macroscopic doses when NP were simulated inside the nucleus (table S3 & figure S5, 80 keV & 83 keV). For NP distribution within the cytoplasm, survival fractions at low macroscopic doses were comparable (figure 5, center, figure S6). Any enhancement with NP distributed in the cytoplasm only would be due to dose at several hundred nanometers or more from the NP center. At these distances the additional dose per interaction is not only approximately three orders of magnitude lower than at the surface but also more consistent between different sources (figure 1) reducing differences in dose to the nucleus, which is the active target. With more interaction events at 8 Gy macroscopic dose, the probability of an event closer to the nucleus is increased, explaining the decrease in SF(8 Gy) for 50 keV, 83 keV, and 100 keV compared to 80 keV or 250 keV. The survival fraction for uniform nanoparticle distribution within the cell averages the behavior of nucleus and cytoplasm geometry as one would reasonably expect. Of note, survival fractions for SOBP protons are highest at the distal end of the SOBP ( figure 6). This is due to the lower number of tracks needed to deposit a given macroscopic dose, Differences between the distal end of the SOBP and the other depths become large, i.e. more than a factor of three, only for 8 Gy macroscopic dose and only in the cell and nucleus geometry. This is due to the lower dose per interaction (figures 1(a)nd (3)), inducing fewer lethal lesions locally, as well as the overall lower interaction rate per 1 Gy (table 1), meaning that only at higher doses the number of interactions and hence, additional local dose become significant compared to the homogeneous dose throughout the nucleus.
The modeled survival fractions in this work are in line with those previously reported (table 2) [24]. Differences for SF(2 Gy) were less than a factor of two for both center SOBP protons and 250 kVp photons when using our 250 kVp (w/80 keV) source for comparison. At higher macroscopic doses, differences for SOBP protons were larger with survival fractions in this work 3.5 (cell) to 80 (nucleus) times lower. Given the order of magnitude of the effect, 10 −3 for the cell and 10 −6 for the nucleus geometry these differences may still be reasonable. Data from Lin et al were requested but were not available so survival fractions were estimated based on the published figures, adding uncertainty to these numbers. For example, SF(8 Gy) for photons reached 10 −10 for the cell geometry and several orders of magnitude smaller for the nucleus geometry, making it not comparable as their figures did not extend further. However, one may also consider survival fractions  -10 10 to be negligible and effectively zero. A limitation of this survival curve modeling is the simplification of the problem into 2D, ignoring the volumetric distribution of GNPs, which would reduce the spatially varying dose contribution from secondary electrons stemming from GNP interactions and hence, increase the survival fractions. A 3D dose accumulation and calculation of average number of lesions, however, is computationally difficult to perform due to memory requirements and time needed to perform matrix operations.
In a multicentric review of Monte Carlo calculated dose enhancement ratios, Rabus et al highlight that what we reported as dose enhancement or contribution to absorbed dose is rather enhancement or changes in the average microdosimetric specific energy [40]. We attempted to include the stochastic nature of GNP interactions and local 'dose enhancement' by sampling collections of GNPs that undergo an interaction given a uniform dose deposited by the

Reporting recommendations
We have shown that small changes in the photon spectrum, especially in the kilovoltage range, may have a non-negligible impact on the physical simulation and biological modeling of radiosensitization due to the presence of metallic nanoparticle within a cell. In line with the AAPM task group report 61 on 40-300 kV x-ray beam dosimetry [41], reporting only the maximum photon energy (kVp) is insufficient. However, given the presence of metallic NPs also the combination of half value layer and maximum photon energy as suggested by that report remains insufficient given the differences observed in the results between the 250 kVP spectra with 80 keV or 83 keV added, as they would have overall similar half value layer if measured in thickness of aluminum or copper. If a published tool like SpekCalc or SpekPy [39,42] were used to generate the spectra, the parameters passed to these tools should be reported. Otherwise, we recommend to include a graphical representation of the spectrum at a minimum, e.g. as shown in figure 8 in the comparison study by Rabus et al [40], or ideally a tabulated histogram. Such a histogram should feature sufficiently small bin widths of no larger than 1 keV close to cross-section resonances of materials of interest for radiosensitization. Since this may include elements with relatively low atomic numbers like iron [43] or germanium [44] one should start to consider this for photon energies starting at 55 keV to include the K-edge of iron. The upper limit for a detailed description of the spectrum may be chosen as energies exceeding 88 keV as K-edge of lead; however, one may extend this to 107 keV to include the K-edge of actinium, given its potential use in radionuclide therapy and hence, potential for interactions with external beam radiotherapy. Outside of these limits a larger bin width may be used. For proton SOBPs, where differences within the SOBP were found to be small compared to those seen with photons, we nontheless recommend that the procedure to obtain the SOBP be described in more detail to ensure reproducibility.
Modeling of the cell survival curves can be performed differently, given the relatively sparse description of this process in most publications. In this work we described the accumulation of dose from ionizations around a nanoparticle. However, one may also calculate the dose due to nanoparticles per macroscopic dose and distribute that dose kernel at all nanoparticles present. This will potentially underestimate the enhancement by averaging out the local dose increase, if close to or within the nucleus. Alternatively, one may also calculate the number of lethal lesions induced by ionizations at the nanoparticle location and aggregate those. While this ignores synergistic effects from multiple nanoparticle ionizations and thus, yields larger survival fractions, one may find it equally reasonable as it excludes damage from multiple tracks.

Conclusions
We replicated the simulations performed by Lin et al [17,24] using a newer version of TOPAS and Geant4 and without access to some in our opinion critical information. Results agreed extraordinarily well for physical stage simulations and reasonably for biological modeling making us confident in suggesting the use of up-to-date simulation software and data files. We also gave recommendations on how to include detailed descriptions of photon spectra in reports of radiosensitization via metal nanomaterials to increase reproducibility.
The stochastic nature of nanoparticle interactions was partially considered in this work. Its further investigation and implications on biological outcomes remains warranted to understand the mechanisms and sensitivities of nanoparticle radiosensitization [40]. Finally, we believe that reproducibility studies remain warranted.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.