Multi-scale Monte Carlo simulations of gold nanoparticle-induced DNA damages for kilovoltage X-ray irradiation in a xenograft mouse model using TOPAS-nBio

Gold nanoparticles (AuNPs) are considered as promising agents to increase the radiosensitivity of tumor cells. However, the biological mechanisms of radiation enhancement effects of AuNPs are still not well understood. We present a multi-scale Monte Carlo simulation framework within TOPAS-nBio to investigate the increase of DNA damage due to the presence of AuNPs in mouse tumor models. A tumor was placed inside a voxel mouse model and irradiated with either 100-kVp or 200-kVp X-ray beams. Phase spaces were employed to transfer particles from the macroscopic (voxel) scale to the microscopic scale, which consists of a cell geometry including a detailed mouse DNA model. Radiosensitizing effects were calculated in the presence and absence of hybrid nanoparticles with a Fe2O3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{Fe}_2\text{O}_3$$\end{document} core surrounded by a gold layer (AuFeNPs). To simulate DNA damage even for very small energy tracks, Geant4-DNA physics and chemistry models were used on microscopic scale. An AuFeNP-induced enhancement of both dose and DNA strand breaks has been established for different scenarios. Produced chemical radicals including hydroxyl molecules, which were assumed to be responsible for DNA damage through chemical reactions, were found to be significantly increased. We further observed a dependency of the results on the location of the cells within the tumor for 200-kVp X-ray beams. Our multi-scale approach allows to study irradiation-induced physical and chemical effects on cells. We showed a potential increase in cell radiosensitization caused by relatively small concentrations of AuFeNPs. Our new methodology allows the individual adjustment of parameters in each simulation step and therefore can be used for other studies investigating the radiosensitizing effects of AuFeNPs or AuNPs in living cells.


Background
The optimization of radiotherapy is a wide and thoroughly studied field. A major limitation of external beam radiation therapy is the fact that it cannot be applied without harming healthy normal tissue around the tumor. One promising approach to reduce healthy tissue toxicity is the local enhancement of the radiation dose within the tumor using radiosensitizers, such as gold nanoparticles (AuNPs). AuNPs are among the most investigated nanomaterials for medical application, due to their favorable biocompatibility and their potential use as a contrast agent (Dai et al. 2017;Chen et al. 2016;Xin et al. 2017). The radiosensitizing effect of AuNPs is partially based on their potential to release high numbers of secondary electrons upon irradiation, which has been shown to be especially potent for kilovoltage photon beams around 100 kVp. In this energy range, the dose enhancement is mostly caused by secondary electrons with low energy such as Auger electrons with a limited range (Lechtman and Pignol 2017). Because of the highly localized effects of the radiosensitization and the possibility of AuNPs to accumulate also in healthy tissues, reliable targeting agents such as antibodies are important to make AuNP-based therapies applicable in clinical practice (Schuemann et al. 2016).
Sensitization effects of AuNPs can be studied via Monte Carlo simulations by computing the dose enhancement effect caused by AuNPs in appropriate scenarios. However, it is still not possible to precisely calculate the impact of a single AuNP on its direct microenvironment, although multiple studies have been dedicated to this topic (Li et al. , 2020Lin et al. 2014;Sung et al. 2016;Rabus et al. 2019;Haume et al. 2018). When it comes to scenarios with multiple AuNPs, a frequently used approach was the employment a local effect model (LEM)-based estimation (Kraft et al. 1999). These models usually calculate the enhancement of the effectiveness of irradiation through the addition of multiple AuNPs by extrapolating radial dose enhancements obtained from MC simulations with a single AuNP. Although they add additional assumptions and uncertainties, LEMs were commonly used to describe the dose enhancement effect of multiple AuNPs, since MC simulations on their own have been shown to underestimate the effects of radiation when compared to experimental results (McMahon et al. 2011). One reason for that is probably the lack of chemical interactions in MC simulations, which usually only compute the effects of physical interactions with a focus on secondary electrons (Kuncic and Lacombe 2018). Chemical radicals produced by water radiolysis have, however, been shown to play a significant role in the radiosensitizing effects of AuNPs (Xie et al. 2015;Rosa et al. 2017).
In this study, we introduce a multi-scale method to directly compute the DNA damage induction of ionized radiation with and without nanoparticles with a gold surface using TOPAS-nBio (Schuemann et al. 2019a). TOPAS-nBio is a radiobiology focused extension to TOPAS (Perl et al. 2012), which is a user friendly Monte Carlo toolkit based on Geant4 (Agostinelli et al. 2003). TOPAS-nBio includes a chemistry interface that employs an updated version of Geant4-DNA water radiolysis models, which were applied in our nano-scale cell simulations (Karamitros et al. 2014;Ramos-Mendez et al. 2018). Rudek et al. (2019) already observed AuNP-induced dose enhancement assuming these models. The effect on DNA damage in a detailed human nucleus model was investigated by Zhu et al. (2020a). We aim to build on both results to gain a better understanding of the radiosensitization effect of AuNPs in tumor cells by introducing a new methodology of simulating according preclinical in vivo experiments.
Our calculations were performed in three steps: mouse model, tumor and cell. On the macro-scale, a voxel mouse model including a mammary gland tumor was irradiated by an X-ray beam. Phase spaces were employed to move between the scales until detailed DNA damage could be calculated in a nano-scale cell model.

Materials and methods
All MC simulations were carried out in TOPAS v3.2 for Linux machines and Geant4 v10.5p1. For reproducibility, all developed extensions and parameter files were made generally accessible via Github (Link: https:// github. com/ AKlap proth/ Multi Scale_ AuNP_ TOPAS). Simulations were performed in three major steps: mouse model, tumor and cell. The particles were transferred from one step to the next via phase space files. These files take a snapshot of all tracks passing a given surface and include all the necessary information for each scored track to be recreated in the next stage of our simulations. The information included were particle type, energy, location and momentum at the time of scoring. The initial beam spectra were calculated with SpecCalc employing the parameters from the small animal radiation research platform (SARRP) device to simulate a realistic experimental setup (Wong et al. 2008).
For mouse model and tumor simulations, the standard Geant4 electromagnetic physics list option 3 was used. It is focused on electron, hadron and ion tracking outside a magnetic field and the most accurate standard model for medical applications in the micrometer range or above (Physics Lists EM constructors in Geant4 10.4 2021; Allison et al. 2016). The use of more accurate physics lists on these scales would have caused an immense increase in computation time with no significant benefit. This is because most of the low energy tracks, which are the most time-consuming to simulate, occur in locations, where they have no effect on the outcome. For the cell simulations, we applied a combination of two more accurate physics lists to make the outcome on the nanometer scale as accurate as possible. The exact physics setup in the cells is explained in the subsection "Cells".
The nanoparticles, which were used in the cell simulations, were hybrid nanoparticles. While their surface consists of a 1-nm-thick gold layer, they feature an Fe 2 O 3 core with a diameter of 2nm, adding up to an overall diameter of 4nm. We selected these nanoparticles to describe an experimental setup that was conducted for theranostic purposes (Kang et al. 2020). The iron oxide core offers advantages for both therapy and diagnostics, since it enables the gold-Fe 2 O 3 hybrid nanoparticles (AuFeNPs) to be detectable via MRI scans.

Mouse model
The voxel model of a 21 g mouse developed by Xie and Zaidi (2013) was used and embedded in a 120 × 120 × 160 mm 3 box filled with air (cf. Fig. 1A). The model consists of 200 × 200 × 512 cubical voxels with a border length of 200 µm. A new TOPAS extension was developed that allows the insertion of an ellipsoid shaped tumor at any location inside the model. Normal tissue, skin and air in the defined area is then replaced by tumor tissue. In case the tumor extends outside the body, exposed tumor tissue is covered by a layer of skin with a thickness of one voxel. For this work, a tumor with diameters 5 mm× 4 mm× 5 mm was planted near the left hind leg to simulate a mammary gland carcinoma. Different materials were defined for each component of the model, either a predefined material available in Geant4 or water with an adjusted density (cf. Table 1). The respective densities were chosen according to previous studies with this mouse model (Xie and Zaidi 2013). Since the tumor consisted of water in the following simulation steps due to limitations in the availability of physics parameters for subcell scale simulations, standard G4_WATER material was used for the tumor in this step as well to ensure consistency. Two different types of kVp photon beams were used for the simulated treatment, both based on the radiation source used in the SARRP, when set to either 100 kVp or 200 kVp. The respective photon spectra were calculated with the SpekCalc toolkit and are shown in Fig. 2 (Poludniowski et al. 2009). Particle tracks originated from a flat, circular particle source placed 20 cm below (i.e., shifted along the y-axis) the tumor center. Its diameter was set to 5 mm, according to the tumor's X-and Z-diameters (cf. Fig. 1B). All particles entering the tumor were scored in a phase space file and multiplied 750 times, which means for each particle scored during the simulation 750 entries were made in the phase space file, each with the same particle type and energy. To increase variance each duplicated particle was positioned on a randomly generated point on the tumor surface, while its momentum was rotated accordingly around the center. Particles that still possessed the initial momentum along the y-axis were assigned a random position on the tumor surface at its bottom or Y-Minus side. All other particles were rotated randomly around the y-axis going through the tumor center.

Tumor
For simulations inside the tumor the entire geometry setup consisted of water to ensure consistency throughout the simulation steps. An ellipsoidal water phantom with diameters of 5 mm× 4 mm× 5 mm was defined as the tumor. Spheres with a diameter of 100 µm were placed at 3 different positions inside the tumor to represent cellular regions of interest (cf. Fig. 3). For each of the spheres, all particles entering the volume were scored in a second phase space file. Parameters were chosen to emulate a mammary carcinoma in a xenograft mouse model. A popular cell line for this purpose is 4T1 mouse tumor cells. Since they tend to metastasize in the mammary gland, when transplanted on another location, and show some similarities to human mammary carcinomas, they are an ideal choice for Fig. 1 Sketch of the voxel model simulations. A Displays a visualization of the implemented mouse model, where each color represents a different organ. B Shows the approach of the simulations. The source consists of parallel photon beams (displayed in green) and is aimed directly at the mammary gland tumor (displayed in orange). The coordinate system has been added for orientation, not to mark the actual origin of the simulation geometry Fig. 2 X-ray photon spectra. The 100-kVp spectrum ranges from 10 to 100 keV with a bin size of 0.05 keV and the 200-kVp spectrum from 20 to 200 keV with bin size 0.1 keV. The energy for each created particle is chosen randomly based on the percentages defined by the respective spectrum. Both spectra possess the same four intensity peaks mammary xenograft tumors (Pulaski and Ostrand-Rosenberg 2001). 4T1 cells are triple-negative breast cancer cells, which means they test negative for both the expression of estrogen receptors, progesterone receptors and an overexpression of the HER2 protein. Patients with triple-negative breast cancer usually have a relatively poor prognosis, since the tumors do not respond to many commonly used treatment methods (Foulkes et al. 2010). These tumor cells are therefore of high interest regarding new and optimized treatment methods (Takai et al. 2016).

Cells
Cells were built as spheres with a diameter of 20 µm and a nucleus diameter of 13.8 µm according to parameters of 4T1 tumor cells (cf. Fig. 4) (Oelze et al. 2004). The nucleus was placed at the center of the cell and filled with DNA using a model based on the human DNA model published by Zhu et al. (2020a, b). The model was modified to conform to the male mouse genome by adjusting the voxel count, voxel size and the number of chromatin fibers per voxel. The resulting nucleus has a diameter of 13.8 µm and each chromosome contains a realistic number of base pairs. The voxels and base pairs per chromosome can be found in Table 2, and a sketch of the resulting nucleus structure Fig. 3 Sketch of the tumor simulations. The particle source, which is shown in green, is the phase space file recorded during the mouse model simulation. Most particles enter the tumor with their initial momentum without previous interaction, but particles can enter from any direction. Cell-sized spheres are defined at three different positions. All particles entering a cell are stored in their respective phase space    Fig. 6A. Figure 6B provides a more detailed depiction of each fiber and how the DNA double helix is wrapped around the histones. The nucleus consists of 24,464 cubical voxels with a border length of 3.833 µm, which are divided into 40 chromosomes. All in all, the nucleus includes around 5.19 Giga base pairs of DNA. The phase spaces from the tumor simulation were used as source. Each particle's distance to the cell center was reduced to 10 µm, while its momentum, particle type and energy stayed unchanged. Particles were then multiplied 100 times, which means for each particle scored during the tumor simulation, 100 particles originated in the respective cell simulations. To increase variance, position and momentum were rotated randomly around the y-axis going through the cell, while the other parameters remained unchanged. All simulations were performed twice, once without AuFeNPs and once with 1 Mio AuFeNPs placed randomly around the nucleus with a maximum distance of 100 nm (cf. Fig. 4). This results in an AuFeNP concentration of around 0.225% by weight in relation to the cell. The cytoplasm and nucleus consist of water in the simulations, so the Geant4-DNA-based physics and chemistry lists , which simulate the full particle track structure and are part of TOPAS-nBio, could be used in these regions. The standard G4_WATER material was used in all regions of the nucleus except the DNA backbone where the density of water was adjusted to 1.407 g/cm 3 (Smialek et al. 2013). Due to limitations in available cross sections in Geant4-DNA, the Geant4 Livermore physics model was applied in gold and Fe 2 O 3 . Production cut and electromagnetic range were both set to 10 eV and 1 MeV and the cut for electrons was defined as 0.1 nm. To achieve the combination of two different models two regions were defined, one containing all AuFeNPs and one all other parts within the cell. The simulation switched between the two electromagnetic physics models depending on the region a particle traversed. Accordingly, only physics interactions (condensed histories) were tracked and generated inside AuFeNPs, while in the rest of the cell (water) we also simulated the propagation of chemical species. For the computation of DNA damage based on our track structure simulations, we chose standard assumptions and parameters used in previous Monte Carlo studies Rudek et al. 2019;Lampe et al. 2018;Meylan et al. 2017;Sakata et al. 2019). Direct damage was only produced by interactions in the DNA backbone including the adjacent hydration shells. Indirect damage was only produced in the DNA backbone. If at least 17.5 eV (Rudek et al. 2019;Lampe et al. 2018;Sakata et al. 2019) of deposited energy caused by physical interactions in a single back bone including its hydration shell was scored during one history, it was considered as a direct strand break (SB). The chemistry model was applied according to the work of Zhu et al. (2020a) with a chemical stage time of 1.0 ns. Six different particle types that are included in the standard Geant4-DNA chemistry extension were quantified, namely hydrogen radicals (H · ), molecular hydrogen ( H 2 ), hydrogen peroxide ( H 2 O 2 ), hydronium ( H 3 O + ), solvated electrons ( e aq ) and the hydroxyl group (OH), which includes hydroxide (OH-) and hydroxyl radicals ( · OH). Chemical species produced inside DNA regions were immediately killed to emulate that no water radiolysis occurs there. H · , e aq and · OH diffusing into DNA regions were also immediately terminated. Only · OH radicals that had just entered the DNA backbone were able to induce an indirect SB with a probability of 40% (Meylan et al. 2017;Sakata et al. 2019).
Damage was estimated in form of direct and indirect single-strand breaks (SSBs) and double-strand breaks (DSBs). Direct damage is caused by physical interactions and indirect damage by chemical radicals. Two SBs that occur on opposing DNA strands with a distance of 10 or fewer base pairs were defined as a DSB. When scored, DSBs were classified depending on the types of their underlying SBs into direct DSBs (two direct SBs) indirect DSBs (two indirect SBs) and hybrid DSBs (one of each). The results were displayed in the standard DNA damage (SDD) data format (Schuemann et al. 2019b). In addition, the number of produced chemical species was counted and classified by type.

Results
Our simulation results show that the cell depth as well as the presence of AuFeNPs affects the simulation results. The exact values of the results and the respective standard deviations can be found in Additional file 1: Tables S1-S12.
As shown in Fig. 7, the dose applied to the nucleus is increased in scenarios, where AuFeNPs are present for each cell location inside the tumor. The position itself also plays a role for the results, as the dose decreases the further back a cell lies in relation to the particle source. This location effect was only significant for simulations with 200-kVp photons.
The SB count shows the same behavior and is increased in simulations including AuFeNPs around the nucleus and decreases with increasing depth only for 200-kVp photons (cf. Fig. 8). Indirect SBs generally account for around twice the number of direct SBs, which conforms to previous findings (Zhu et al. 2020a). Notably, the dose and damage increase was observed in all scenarios, although the concentration of AuFeNPs was relatively low in comparison to previous AuNP studies (Sung et al. 2017;Sung and Schuemann 2018).
In accordance with the proportion between direct and indirect SBs, indirect and hybrid DSBs outweigh direct DSBs in all simulations (cf. Fig. 9). A clear shift in DSBtype ratios depending on AuFeNP presence or cell location could not be detected.
The effect of AuFeNPs is especially prominent when examining the count of produced chemical species for different scenarios (cf. Fig. 10). A clear enhancement in the presence of AuFeNPs can be observed for each investigated location and photon spectrum. The increase of total numbers is mainly attributed to the very high enhancement ratios of OH and H 2 . However, not all chemical radical counts are increased, when adding AuFeNPs. The count of H 3 O particles and e aq decrease significantly (cf. Table 3). We exclusively considered · OH for the scoring of indirect DNA damages, as hydroxyl radicals are generally considered to be the mediator of much of the induced DNA damage (Balasubramanian et al. 1998). Therefore, the increase of produced OH radicals is especially relevant and directly connected to the observed increase in indirect DNA damage. Analogically to the already mentioned results, we also observe an overall decrease of produced chemical species as a function of depth for 200-kVp photons.

Discussion
We successfully developed a method for multi-scale simulations of mouse tumor irradiations in radiobiological experiments, from large scale (mouse model) to small scale (DNA). This is a major step away from LEM-based models, allowing for a better understanding and more accurate modeling of in vivo experiments. This was in part achieved by accounting for realistic effects on the radiation field of tumors within the mouse geometry and by including chemical reactions that are induced by radiolysis.
In our simulations, we detected an enhancement of deposited energy and chemical radicals produced by water radiolysis in simulations including AuFeNPs, resulting  Zhu et al. 2020a, b). When irradiated, the depth of a cell inside the tumor affects the dose that is applied to its nucleus. This naturally also affects the amount of DNA damage, which is illustrated by the number of strand breaks. The addition of gold nanoparticles causes a small but significant increase in both direct and indirect strand breaks, which is caused by secondary electrons and their adjunctive chemical species. However, some limitations have to be kept in mind with regard to the results.
In the MC simulations, many physical, chemical and biological parameters were used, and these parameters are mostly evaluated from experimental data, and they are subject to large uncertainties. The physical parameters such as cross sections for electrons have an uncertainty of 5% at the lower energy of 1 keV, 17-20% uncertainty at very low energy of 100 eV in water (Thomson and Kawrakow 2011). In contrast, the chemical and biological parameters are mostly derived from in vitro cellular and molecular experimental reactions and radiobiological effects, and they compose even larger uncertainties (Zhu et al. 2020b).
In the present study, we did not perform a detailed uncertainty analysis of the results of SSBs, DSBs and chemical species shown in Figs. 7,8,9,10. Such an analysis would require significant further work on parameters analysis as performed in the work of Zhu et al. However, as Zhu et al. found out, the threshold energy of 17.5 eV used in this study and the probability of radicals for generating an indirect damage, account for most of the uncertainty of SSBs and DBSs. Zhu et al. (2020b) estimated differences of up to 34% and 16% for SSB and DSB yields, respectively, caused by all Fig. 9 Total number of computed DSBs with and without AuFeNPs. Results are subdivided into direct DSBs caused by physical interactions, indirect DSBs caused by chemical reactions and hybrid DSBs caused by one physical and one chemical reaction the parameters used in TOPAS-nBio. These uncertainties are applicable to our work as well.
The major limitation for the simulations is the lack of models for chemical interactions inside and on the surface of AuFeNPs. This means that each chemical track  encountering an AuFeNP is instantly eliminated (i.e., equivalent to reacting with the AuFeNP without consequence) and can no longer produce DNA damage. Rudek et al. found that AuNPs could thus lead to a reduction of availability of chemical species when modeled with track structure MC codes (Rudek et al. 2019). This is a considerable discrepancy between simulations and reality, where interactions between chemical radicals and gold atoms exist and can cause downstream reactions. Cheng et al. (2012) detected an enhancement of chemical activity around AuNPs that could not only be explained by water radiolysis caused by X-rays. They attributed this effect to the activation of gold atoms on the surface of AuNPs by superoxides. In addition, there are still no reliable models for the surface chemistry of gold layered nanoparticles themselves, which have been shown to enhance reactive oxygen species (ROS) production (Pan et al. 2009;Liu et al. 2014;Sicard-Roselli et al. 2014;Seo et al. 2017). Possible charge accumulations in the AuFeNPs and their surface could further affect the reactions or influence which species are attracted (or repelled) by the AuFeNPs. Another factor to consider regarding in vivo studies is the nanoparticle coating necessary for targeting and biocompatibility. Xiao et al (2011) showed that such coatings can decrease the radiosensitization effect of AuNPs significantly. Such coatings can further impact the chemical reactions around the GNPs. Another limitation is the currently available physics models of Geant4-DNA, which do not yet include cross sections for gold. Thus, standard Livermore models were employed, which show good accuracy for keV electrons, but show limitations below the 100 eV range . The detected effect of AuFeNPs was narrow and no increase in the fraction of SBs involved in DSBs was detectable, which might be expected due to the agglomeration of damaging events around the nanoparticles. A possible reason might be the distance of AuFeNPs to regions that contain DNA and thus are relevant for the results, due to the relatively low concentration of DNA inside the nucleus. There was still a definite AuFeNP-related increase in SBs, dose and chemical species detected for all depths, which is in agreement with previous results Rudek et al. 2019). The release of additional physics models will be a big help in improving and validating simulation results for the physical effects of AuFeNPs.
The next important step will be the validation of simulation results in biological experiments. As the study of Zhu et al. (2020b) showed, simulations are sensitive to the chosen parameters and physics models. Experimental data will therefore play a crucial role in deciding which simulation parameters can accurately describe observed phenomena. The output in form of SSBs and DSBs can be used to compare cell survival outcomes from in vitro experiments and in vivo studies with mice. The described multi-scale approach can further be used to investigate the effect of AuFeNPs at the involved scales and on the efficiency of AuFeNP-enhanced radiotherapy.

Conclusions
In this work, a new Monte Carlo-based method was established to simulate radiationinduced DNA damages in AuNP or AuFeNP-assisted radiotherapy. This multi-scale approach covers the whole range of preclinical in vivo studies and can therefore be valuable for parameter optimization and analyzing results in clinical cancer radiotherapy settings in the future. Including the AuFeNPs in the simulation geometry, as opposed to using a LEM-based approach, utilizes the ongoing increase in both processing power of computers and advancement of Monte Carlo models to produce more accurate and traceable results. The detailed nucleus model allows direct counting and classification of SSBs and DSBs and the output in the SDD format makes the results comparable to similar studies and experimental data (Schuemann et al. 2019b). The results show that even low concentrations of gold can cause a noticeable increase in DNA damage after kV irradiations and highlights the importance of taking chemical interactions into account.
The inclusion of the tumor as a separate step allows the consideration of tumor heterogeneity in future studies. It is well known that the tumor landscape can differ vastly for different locations of the same tumor, with changing AuFeNP uptake and penetration, cell type, radiosensitivity and microenvironment (Alfonso and Berk 2019;Dagogo-Jack and Shaw 2018;Marusyk and Polyak 2010;Zhivotovsky et al. 1999). These effects can be taken into account, for example, by including different SB probabilities for different cell locations or adjusting repair parameters, when calculating cell survivability.
This methodology is not restricted to X-ray or metal nanoparticle studies. It can be easily adjusted to cover all investigations of radiation-induced DNA damage. The inclusion of a complete mouse model enables dosimetry applications in organs of interest. Human cell studies can be performed by replacing the mouse model with a phantom displaying the respective area of interest, adjusting cell size, and replacing the mouse DNA model with the already established human equivalent (Zhu et al. 2020a). The kVp photon beam may be replaced by any other radiation source supported by TOPAS, e.g., MV photon beams (LINAC) or proton beams for studies with or without nanoparticles. This work can be used as a template for future multi-scale radiation studies in different settings.