Search

Collections

Journals

About

Contact us

My IOPscience

Dependence of Monte Carlo microdosimetric computations on the simulation geometry of gold nanoparticles

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2013 Phys. Med. Biol. 58 7961 (http://iopscience.iop.org/0031-9155/58/22/7961) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 128.114.34.22 This content was downloaded on 28/11/2014 at 19:15

Please note that terms and conditions apply.

IOP PUBLISHING

PHYSICS IN MEDICINE AND BIOLOGY

Phys. Med. Biol. 58 (2013) 7961–7977

doi:10.1088/0031-9155/58/22/7961

Dependence of Monte Carlo microdosimetric computations on the simulation geometry of gold nanoparticles Piotr Zygmanski 1 , Bo Liu 1,2 , Panagiotis Tsiamas 1 , Fulya Cifter 3 , ¨ Hesser 2 and Erno Sajo 3 Markus Petersheim 2 , Jurgen 1 Department of Radiation Oncology, Brigham and Women’s Hospital, Dana-Farber Cancer Institute and Harvard Medical School, Boston, MA 02115, USA 2 Department of Radiation Oncology, University Medical Center Mannheim, University of Heidelberg, Germany 3 Department of Physics and Applied Physics, Medical Physics Program, University of Massachusetts Lowell, Lowell, MA 01854, USA

E-mail: [email protected]

Received 15 May 2013, in final form 28 September 2013 Published 30 October 2013 Online at stacks.iop.org/PMB/58/7961 Abstract Recently, interactions of x-rays with gold nanoparticles (GNPs) and the resulting dose enhancement have been studied using several Monte Carlo (MC) codes (Jones et al 2010 Med. Phys. 37 3809–16, Lechtman et al 2011 Phys. Med. Biol. 56 4631–47, McMahon et al 2011 Sci. Rep. 1 1–9, Leung et al 2011 Med. Phys. 38 624–31). These MC simulations were carried out in simplified geometries and provided encouraging preliminary data in support of GNP radiotherapy. As these studies showed, radiation transport computations of clinical beams to obtain dose enhancement from nanoparticles has several challenges, mostly arising from the requirement of high spatial resolution and from the approximations used at the interface between the macroscopic clinical beam transport and the nanoscopic electron transport originating in the nanoparticle or its vicinity. We investigate the impact of MC simulation geometry on the energy deposition due to the presence of GNPs, including the effects of particle clustering and morphology. Dose enhancement due to a single and multiple GNPs using various simulation geometries is computed using GEANT4 MC radiation transport code. Various approximations in the geometry and in the phase space transition from macro- to micro-beams incident on GNPs are analyzed. Simulations using GEANT4 are compared to a deterministic code CEPXS/ONEDANT for microscopic (nm–μm) geometry. Dependence on the following microscopic (μ) geometry parameters is investigated: μ-source-toGNP distance (μSAD), μ-beam size (μS), and GNP size (μC). Because a micro-beam represents clinical beam properties at the microscopic scale, the effect of using different types of micro-beams is also investigated. In particular, a micro-beam with the phase space of a clinical beam versus a plane-parallel beam with an equivalent photon spectrum is characterized. Furthermore, the 0031-9155/13/227961+17$33.00

© 2013 Institute of Physics and Engineering in Medicine

Printed in the UK & the USA

7961

7962

P Zygmanski et al

spatial anisotropy of energy deposition around a nanoparticle is analyzed. Finally, dependence of dose enhancement on the number of GNPs in a finite cluster of nanoparticles is determined. Simulations were performed for 100 nm GNPs irradiated in water phantom by various monoenergetic (11 keV–1 MeV) and spectral (50 kVp) sources. The dose enhancement ratio (DER) is very sensitive to the specific simulation geometry (μSAD, μS, μC parameters) and μ-source type. For a single GNP the spatial distribution of DER is found to be nearly isotropic with limited magnitude and relatively short range (∼100– 200 nm for DER significantly greater than 1). For a cluster of GNPs both the magnitude and range are found much greater (∼1–2 μm). The relation between DER for a cluster of GNPs and a single GNP is strongly nonlinear. Relatively strong dependence of DER on the simulation micro-geometry cautions future studies and the interpretation of existing MC results obtained in different simulations geometries. The nonlinear relation between DER for a single and multiple GNPs suggests that parameters such as the number of adjacent nanoparticles per cell and the distances between the GNPs and the cellular target may be important in assessing the biological effectiveness associated with GNP. S Online supplementary data available from stacks.iop.org/PMB/58/7961/mmedia

(Some figures may appear in colour only in the online journal) 1. Introduction Recently, Monte Carlo (MC) simulations have been used to investigate the dose enhancement from gold nanoparticles (GNPs) within microscopic (nm–μm) distances in response to clinical x-ray beams (Jones et al 2010, Lechtman et al 2011, McMahon et al 2011, Leung et al 2011). The dose enhancement was characterized via considering numerous interaction phenomena and quantities, such as photoelectric absorption, properties of secondary electrons ejected from GNPs, photoelectric energy conversion, as well as dose ratios. These studies included both external (linac, x-ray tube) and internal (brachytherapy) x-ray sources. Dose enhancement was determined for various simplified simulation geometries approximating clinical macroscopic conditions. However, the intrinsic dependence of the computations on the parameters of the microscopic simulation geometry has not been recognized and comprehensively studied. Similarly, the spatial distribution of energy deposition due to a single versus multiple nanoparticles has not been fully characterized. Due to various limitations in the MC method, simulations of energy deposition at submicroscopic (or nanoscopic) levels must be carried out separately from the simulations at macroscopic (mm–cm) levels. Typically, in the context of GNPs, macroscopic simulations of coupled photon–electron transport for clinical beams are performed using MC codes such as EGSnrc (Jones et al 2010) or MCNP5 (Lechtman et al 2011) that employ the condensed history electron transport approximation. Sub-micron (nanoscopic) simulations are carried out using codes whose energy cutoff permits the simulation of shorter electron ranges, such as Penelope (Lechtman et al 2011), GEANT4 (Leung et al 2011), or others that can consider electron track structures, such as NOREC (Jones et al 2010). In the first (macroscopic) stage of the simulation the x-rays are transported across macroscopic depths (mm to cm) to a location near the target GNP in a homogeneous medium where the phase space is determined. Phase space incorporates all particle types, their numbers, energy and direction of motion. In the second stage, this phase space (or sometimes spectrum) is used as the source for a

Dependence of MC microdosimetric computations on the simulation geometry of GNPs

7963

new MC simulation to determine the microscopic dosimetric quantities. Therefore a bridge between the two levels (macro–micro) has to be established. Modeling the transition between the macroscopic and nanoscopic levels of simulation is essential in order to ensure that the magnitude of the computed dose enhancement effects is not misrepresented. This is because dose enhancement is very sensitive to beam quality (energy content) (Tsiamas et al 2012, 2013a, 2013b) and therefore it may depend on the specific simulation geometry as well as on the angular component of the microscopic phase space used in the second stage. In the present paper, we investigate the dependence of the computed dose enhancement on the simulation geometry and on the angular distribution of the microscopic source in relation to the macroscopic phase space at a given depth in the medium. The dose enhancement ratio (DER) is used as the metric of clinical relevance, which is the ratio of energy deposited in the water medium with and without GNPs, i.e., DER(x, y, z) = DGNP/Dwater. In addition, we study the relation between the DER for a single versus multiple GNPs. In particular, we focus on the behavior of DER as a function of distance from the GNPs and of the number of GNPs in a cluster. The anisotropy of the DER is also considered. 2. Materials and methods 2.1. General We focus our analysis on the computation of the DER in the vicinity of GNPs using the GEANT4 MC model simulating various geometries. Geant version 4.9.4 with Penelope-based model of interactions (photons, electrons, positrons) was applied. Specifically, the ‘atomic de-excitation’ and ‘ActivateAuger’ modules were turned on (Chauvie et al 2004, Guatelli et al 2007). DER is defined as the local energy deposited within a microscopic voxel in the medium near the GNP (or multiple GNPs) divided by the energy deposited in the same voxel without the presence of GNP(s). In general, the macroscopic phase space at a given point in the medium can be determined for any heterogeneous composition with macroscopic concentrations of nanoparticles therein. However, we focus on a homogeneous medium with a single GNP or a cluster of GNPs. Here, we use the terms ‘microscopic’ and ‘micro-’ to denote effects occurring at nanometer-to-micrometer distances from the GNP for which nanometer spatial resolution is required. Thus the microscopic scoring voxel size in the proximity of GNPs will be in nanometers. In order to better understand the effect of the spectral characteristics of the incident beam on the induced and computed DER, and to improve the computational efficiency without using aggressive non-analogue variance reduction, the MC simulations of x-ray beams are performed in two stages: macroscopic and microscopic. This is a well-known technique in multi-scale radiation transport simulations, which has also been used by the nanoparticle computational literature (Jones et al 2010, Lechtman et al 2011). In this work, unless otherwise stated, the macroscopic source is an external monoenergetic x-ray beam of size S at the surface of a water phantom (figure 1—step 1). The entering x-ray beam is plane-parallel, normally incident on the phantom, and it is transported to a macroscopic depth d in water where the macroscopic phase space is determined using GEANT4. This includes the photon and electron spectra and angular fluences. In the second stage, the central axis portion of the macroscopic phase space within a 10 cm × 10 cm area at depth d is used as a source to obtain a microscopic phase space for a micro-beam of lateral size μS. This is a microscopic source (μ-source) whose origin is at depth d in the phantom, and it has the same angular and energy distribution as the macroscopic phase space within μS about the central axis, except its spatial resolution matches the small size of GNP (figure 1 – step 2).

7964

P Zygmanski et al

Figure 1. Schematic diagram of the macroscopic and microscopic simulation geometries for a single gold nanoparticle (GNP) of varying size μC × μC × δ (case 1) and a square array of varying number of GNPs (N × N) separated by distances C (case 2).

Next, the microscopic beam (μ-beam) is transported from depth d to a single or multiple GNPs and micrometers beyond them. The distance between the origin of the μ-source and the GNP is denoted μSAD. Simulations are preformed for spherical, square, and rectangular GNPs. Unless otherwise stated, we present DERs computed for rectangular GNPs of thickness δ = 100 nm and cross-section μC × μC perpendicular to the direction of the primary x-ray beam. The micro-beam size (μS), GNP size (μC) and the distance between the origin of micro-source and the GNP (μSAD) are varied, and their impact on the spatial distribution and the magnitude of the DER, as well as on its angular anisotropy is determined. Additional cases are also presented using spherical GNPs in spherical geometry. The reasons why monoenergetic macroscopic x-ray beam and micro-geometry are adopted are several. In spectral beams the information about contributions from specific energies is lost or at least confounded. Second, DER is sensitive to beam quality (Tsiamas et al 2012, 2013a, 2013b). It depends on the beam delivery technique or clinical source type (e.g. open beam versus IMRT, standard versus FFF linac) and field parameters (depth, field size, location within the field, in-field versus out-of-field). Therefore, investigation of the dose enhancement potential and related principal characteristics of monoenergetic beams is appealing, as the effect of a spectral beam may be modeled by the superposition of monoenergetic beams. Third, the second stage simulation with nanometer spatial resolution cannot be efficiently achieved using millimeter or centimeter sized sources. In order to determine the impact of the simulation geometry on the computed DER, the following studies were performed. Case 1. Dependence of DER micro-beam size, and the source-to-GNP distance, i.e., parameters, μS, and μSAD, respectively. Case 2. Dependence of DER on the number of GNPs in a rectangular N × N array and on their relative distances (parameters N and C). Case 3. Dependence of DER on the angular distribution of μ-source with respect to the macrosource.

Dependence of MC microdosimetric computations on the simulation geometry of GNPs

7965

Case 4. Angular anisotropy of DER for a single GNP (scoring energy in three directions: backward, forward and oblique). Case 5. Comparison of DERs computed using GEANT4 (MC) and CEPXS/ONEDANT (deterministic) transport codes. Case 6. Energy dependence of DER. Case 7. Dependence of DER on the contribution of energy fluence from adjacent planar clusters. In studies 1, 2 and 7 finite slab geometry with rectangular GNPs was adopted, as explained above. Studies 3, 4, 5 and 6 were done in spherical geometry with spherical GNPs. Experiments in vitro and in vivo found that nanoparticles tend to form clusters (Zheng et al 2008, Jeong et al 2009, Rahman et al 2009, Zheng and Sanche 2009). This is supported by colloid theory that predicts a non-negligible collision density among the nanoparticles, especially at high number densities, which occur at and near the administration site, or in storage of nanoparticles. Hence the most likely morphology is a cluster of nanoparticles rather than solitary nanoparticles. Owing mostly to self-shielding and nanoparticle packing geometry within the cluster, the dose distribution about the cluster of nanoparticles cannot be obtained by superimposing dose distributions due to many single nanoparticles. The dose to the cellular or sub cellular target is therefore a nonlinear function of the number of particles and their morphology. To study the effect of single versus multiple GNPs we adopt two approaches in cases 1 and 2. First, we investigate the dependence of DER on the size of a GNP cluster forming an N × N planar array perpendicular to the micro-beam (μC × μC). Assuming that the cluster is densely packed, we consider two limiting cases: a single rectangular GNP and an effectively infinite GNP slab in the lateral directions, both 100 nm thick (figure 1, right-hand side). The lateral infinite condition is relative to the scoring voxel size in the MC simulations and to the range of secondary electrons originating in the GNP. Second, we investigate the dependence on the number of isolated GNPs of size δ × δ × δ = 100 nm × 100 nm × 100 nm in a rectangular array of a fixed size (figure 1—beam’s eye view). In the latter case, the separation distance, C, between the GNPs varies but the array size μC × μC is fixed. Thus the two limiting cases are: a collection of single GNPs separated by large distances to regard them as isolated GNPs and densely packed GNPs forming a GNP slab. These geometries allow us to characterize the dependence of DER on the size of the GNP cluster, the distances between the nanoparticles, and the number of nanoparticles in the cluster, as well their relation to a single GNP and a large cluster of GNPs. To permit cross-comparisons, in these two cases we simulate rectangular instead of spherical GNPs, with practically no bearing on the final results compared to simulations based on spherical GNPs. Other geometries considered in this work entail spherical GNPs, as detailed below. 2.2. Simulation geometry (case 1) Simulation geometries for the macro- and microscopic stages are shown in figure 1 (steps 1 and 2). The macroscopic water phantom size is 100 × 100 × 30 cm3. Single or multiple GNPs are located at depth of d = 2 cm, where minimal changes in the beam energy and angular distribution are observed compared to larger depths. Even though any depths are allowed, the goal of using the shallow 2 cm is to keep the photon spectrum closer to the original distribution for the aforementioned reasons. The energy of the external photon beam is E. The field size of the external macroscopic photon beam incident on the phantom is S × S = 50 × 50 cm2, while the micro-beam has varying size, ranging from μS × μS = 100 nm × 100 nm to 5 μm × 5 μm. The phase

7966

P Zygmanski et al

space of the micro-beam is generated from the macroscopic beam by extracting the particles within a 10 cm × 10 cm area and creating an equivalent phase space contained within the μS × μS area of the μ-source about the central axis at location d depth in the phantom. The μ-source-to-GNP distance is also varied (μSAD = 50 nm–5 μm). The size of a single GNP is denoted μC and the size GNP cluster (array) is μC × μC, which is varied up to 2 μm × 2 μm. The thickness of each GNP (single or in a cluster) is δ = 100 nm. The thickness of scoring voxel size varies from 10 to 100 nm, depending on its distance from GNP. Specifically, within the first 150 nm from the center of the GNPs the voxel size is 10 nm, which increases to 50 nm in the range extending from 150 to 600 nm, and to 100 nm from 600 nm to 5 μm distances. For μC = 100 nm a single nanoparticle is implied, and for μC = 2 μm a densely (case 1) or loosely (case 2) packed cluster (containing up to 1002 = 104 nanoparticles) are implied. When the micro-beam size is equal or smaller than the size of the GNP (μS μC) narrowbeam geometry is formed, in which there is lateral particle disequilibrium, i.e., there is a loss of particles from the micro-beam between the micro-source and the GNP. In contrast, when the micro-beam is greater than the size of the GNP (μS μC) lateral (semi-) equilibrium exists at the central axis of the GNP and the beamlet. Hence, when μS > μC the flux reaching the GNP is nearly constant as function of μS, while it is gradually decreasing when μS μC. 2.3. Single versus multiple GNPs (case 2) Figure 1 also shows the schematic diagram of the GNP cluster to study DER as a function of the number of nanoparticles in an N × N array (bottom right). Cubical GNPs of δ × δ × δ = 100 nm × 100 nm × 100 nm size are employed. The size of the GNP cluster is fixed (μC × μC = 2 μm × 2 μm) and the center-to-center distance, C, between individual nanoparticles is varied (C = 50 nm – 950 nm). Therefore, the number of nanoparticles in the cluster is inversely proportional to C2. 2.4. Angular fluence and secondary electron build-up in the μ-source (case 3) To characterize the dependence of DER on the angular fluence of the μ-source and the secondary electron build-up in the micro-beam, which are related to the depth of GNP in the phantom, simulations were performed and compared for two scenarios: when the phase space of the μ-source is equivalent to that of the macroscopic beam at depth d, and when in the micro-beam the angular photon fluence is replaced by a plane-parallel fluence and the secondary build-up electrons are removed. The energy distribution of the parallel beam is obtained from the macroscopic phase space. The setup geometry is otherwise as given in figure 1 (step 2). 2.5. Angular anisotropy of DER (case 4) To investigate the angular anisotropy of DER about a single nanoparticle, a spherical GNP of δ = 100 nm diameter was assumed as shown in figure 5. The μ-source was placed on the central axis at 2 cm macroscopic depth at μSAD = 5 μm away from the GNP, and scoring of energy deposition was along the beam axis. The distance from the surface of GNP is denoted as r, and the angle with respect to the incident micro-beam direction is θ . Dose is scored in three intervals (θ < π /4, π /4 < θ < 3π /4, θ > 3π /4) corresponding to forward-, oblique-, and backscatter. Two monoenergetic sources, at 20 and 100 keV, were considered, targeting the K- and L-shells of gold.

Dependence of MC microdosimetric computations on the simulation geometry of GNPs

7967

2.6. GEANT4 simulation of boundary crossing (case 5) As part of the present work, we also investigated the energy deposition close to the GNP boundaries with high spatial resolution and low noise. GEANT4 was used to model a 100 nm thick spherical GNP shell exposed to photons in spherical symmetry, shown in figure 6, which effectively makes this problem one-dimensional. Photons originate in a uniform spherical surface source, concentric with the 100 nm gold shell, with initial direction normal to the surface. The surface source is 1 mm away from the gold shell (μSAD = 1 mm) whose radius is 2 cm. We compared DER(r) results generated by GEANT4 to results obtained using the CEPXS/ONEDANT coupled photon–electron deterministic transport code in spherical geometry (Lorence et al 1989, Lorence and Beutler 1997). The choice of 2 cm was guided by the need to provide relatively large radius to ‘flatten’ the curved geometry of the spherical shells (both source and GNP are spherical shells). If the source-to-GNP distance radius, then the transport from the spherical shell source to the GNP is conceptually equivalent transport in slab geometry. This was done to provide a fair comparison of GEANT5 to CEPXS doses. Note also that the two geometries (spherical and rectangular) result in different scoring and obviously slightly different results. 2.7. Energy dependence of DER (case 6) Energy dependence of DER was studied for a spherical geometry in order to decrease the stochastic errors in the MC simulations. The incident photons were varied in energy between 10 keV and 1 MeV with sufficiently high resolution to correctly simulate the K- and L-shell transitions. 2.8. Dependence of DER on dose contributions from adjacent planar clusters (case 7) In our study we did not directly simulate macroscopic concentrations of GNPs or multiple GNP clusters. Instead, we determined the order-of-magnitude effect due to multiple GNP clusters or macroscopic concentrations of GNPs by performing convolution of doses based on MC dose for a single cluster. To further investigate how DER increases for multiple clusters, we assumed that there are M densely packed GNP clusters (slabs) as in figure 8(A). Energy deposition for M slabs was obtained by summing the contribution of individual slabs obtained from MC simulations by offsetting them to the required distance along the z-axis. 2.9. GEANT4 simulation parameters In this study, the GEANT4 low-energy package was selected with cutoff energy of 250 eV for electron and photon transport in both gold and water. Atomic relaxation, including Auger production, was enabled with the photoelectric process, along with the Compton and electron ionization processes. The detailed simulation of atomic relaxation cascade is very important for micro- or nano-scale calculations, and it improves the accuracy of MC results in the low-energy domain (Lorence et al 1989, Lorence and Beutler 1997). In addition, GEANT4’s low-energy package uses the condensed history algorithm with multiple scattering approximation for charged particle transport. The step size of charged particles permitted nanometer spatial resolution. The total number of histories for each simulation ranged from 108 to 2 × 109. A total of 30 MC runs were performed in various geometries. More details about packages used in GEANT4 are provided in Lorence and Beutler 1997.

7968

P Zygmanski et al

2.10. Deterministic radiation transport simulation An alternative to the stochastic method employed by MC computer models is the deterministic solution of Boltzmann’s transport equation. The main advantages that deterministic methods have over MC are computation speed, lack of statistical uncertainties (aka noise), and that the full phase space is obtained in the entire simulation geometry in a single run rather than at designated tally voxels or surfaces. They are also less sensitive to boundary artifacts when high spatial resolution is desired, i.e., when the separation of boundaries are small compared to the mean free path or range of the radiation, as we do in this study. With precise solutions, the computation burden is significantly reduced compared to the MC method (Lorence et al 1989, Lorence and Beutler 1997). Owing to these advantages, deterministic methods are more conducive to acquiring physical quantities with high spatial resolution. On the other hand, deterministic methods require more computation memory, especially in multi-dimensional calculations (Lorence et al 1989). We used the CEPXS/ONEDANT coupled electron–photon transport code package (Lorence et al 1989). The CEPXS module is a multigroup-Legendre cross-section generating code that produces consistent electron cross-section sets for one-dimensional neutral particle discrete ordinates (SN) codes in Gaussian directional sets. The multigroup-Legendre crosssections enable coupled electron–photon transport calculations using ONEDANT. ONEDANT is an advanced discrete ordinates code that solves the one-dimensional multigroup transport equation in plane, cylindrical, spherical, and two-angle plane geometries. Regular and adjoint, inhomogeneous and homogeneous problems subject to vacuum, reflective, periodic, white, albedo, or inhomogeneous boundary flux conditions are solved. General anisotropic scattering is allowed and anisotropic inhomogeneous sources are permitted. CEPXS/ONEDANT has been extensively validated (Lorence et al 1989) and benchmarked for medical physics applications (Williams and Sajo 2002). For cross-comparison with GEANT4 results in case 5 (figure 6), CEPXS/ONEDANT simulations were performed in one-dimensional spherical geometry, which is also equivalent to slab geometry via geometric transformation (Jaeger et al 1968). Therefore, this code is suitable to simulate individual and tightly packed clustered GNPs assuming a nearly uniform reaction rate (photoelectric, Compton, etc) across the volume of the nanoparticle. At practical photon energies this requirement is satisfied.

3. Results 3.1. The effects of the simulation geometry on the DER (case 1) The DER as a function of distance from the surface of nanoparticles depends on the simulation parameters μSAD, μS, and μC, as defined in figure 1. The results of the computations are collected in figures 2(A)–(D). Figures 2(A) and (B) show the DER as function of distance from the GNP and for different values of μSAD for the case when the incident energy is E = 100 keV, μS = 1 μm, and μC = 100 nm. Owing to the degradation of photon energy and the increasing availability of secondary electrons to impinge on the GNP, the DER increases with μSAD. The dependence of DER on μS is displayed in figures 2(C) and (D) for the case when μSAD = 50 nm, μC = 100 nm, and E = 100 keV. When μS > 100 nm = μC, the DER does not appear to vary significantly with the size of the micro-beam. This is expected, as in a beam that is wide compared to the size of the GNP (μS > μC) the photon fluence at the central axis of GNP is in lateral equilibrium. Therefore, in this case the fluence does not tend

Dependence of MC microdosimetric computations on the simulation geometry of GNPs

7969

(A) (B)

(C) (D)

(E)

Figure 2. (A) Dose enhancement ratio (DER) as function of distance, z, from GNP for different μSAD cases (μS = 1 μm). (B) DER dependence on μSAD (μS = 1 μm). DER at points 5 and 50 nm away from gold surface were selected. (C) DER as function of z, for different μS field sizes (μSAD = 50 nm). (D) DER dependence on μS at two selected distances, z = 5 and 50 nm, from the GNP (μSAD = 50 nm). Simulations done in slab geometry of figure 1 for E = 100 keV monoenergetic photons and μC = 100 nm, x = y = 100 nm. (E) Large artifacts in DER can be seen by selecting extreme values for the micro-source size and distance to GNP (μS = 100 nm, μSAD = 1 mm) instead of reasonable values (μS = 10 μm, μSAD = 10 μm) for 50 kVp.

to change with the micro-beam size and the DER remains nearly constant with μS. In contrast and remarkably, when μS = 100 nm = μC, the DER has significantly higher values than for any of the larger field sizes. This behavior is due to the narrow-beam geometry, which is not obvious on prima facie. In the case when the beam size equals or it is smaller than the GNP (μS μC, in the present case when μS = μC = 100 nm) the narrow beam experiences lateral scattering and a net loss of photons and secondary electrons between the μ-source and the location of the

7970

P Zygmanski et al

(A)

(B)

Figure 3. (A) DER as function of distance from GNP surface, z, for various distances, C, between the nanoparticles in the N × N cluster (E = 100 keV, μSAD = 50 nm, μS = 2 μm, x = y = 100 nm). (B) DER as function of distance from GNP center, r, for two limiting cases: single GNP versus densely packed cluster of GNPs (forming a slab of μC = 2 μm). Figure 1 shows the simulation geometry.

GNP. Hence, the photon and electron fluxes at the location in the phantom medium where the DER is computed are smaller in the narrow beam than in a wide beam where lateral scatter is in equilibrium. The net electron loss from the beam results in an underestimate of the energy deposition without GNP, while the loss of photons via lateral scattering has a smaller effect on the energy deposition with or without GNP. Therefore, computational simulation of dose deposition about the GNP using an incident beam whose size is equivalent or smaller than the size of the GNP strongly overestimates the dose enhancement. For instance, figure 2(E) shows that extremely large values of DER can be obtained if unrealistic simulation parameters are used (e.g. μS = 100 nm and μSAD = 1 mm). We stress that the under- and overestimates are not due to the inability of the MC simulation to correctly predict the dose, but they are due to inadequate simulation of irradiation conditions. Furthermore, in our study we did not directly simulate macroscopic concentrations of GNPs, but rather a single GNP or a finite cluster of GNPs. 3.2. Single versus multiple nanoparticles (case 2) Figure 3(A) displays the DER as function of distance from the GNP surface, z, for different values of C in an array of GNPs whose geometry is shown in figure 1(B). Dose is scored in a series of voxels along the central axis, with 100 nm × 100 nm size in the transversal plane. Two limiting cases are also shown in figure 3(B): When C = Cmin = 50 nm, the cluster is densely packed (contiguous planar cluster or slab) and the DER is maximum. At large separation, when C approaches 950 nm, the DER is nearly equivalent to the DER for a single nanoparticle. In this case, the DER rapidly falls off with distance, and it becomes 1 (little or no enhancement) at and beyond about 250 nm. 3.3. Angular distribution and secondary electrons in the μ-source (case 3) The μ-source used in the second stage simulation is a microscopic representation of the macroscopic beam of the first stage simulation, and it has the same energy and angular distribution. The phase space contains the angular photon and electron fluence due to interactions in the phantom medium between the external source and GNP. If we replace this μ-source with a plane-parallel photon-only beam with identical energy spectrum the

Dependence of MC microdosimetric computations on the simulation geometry of GNPs

7971

Figure 4. DER as function of distance from the GNP surface, z, with and without the effects of secondary electron build-up and angular scattering in the micro-source due to the presence of water in front of the GNP (E = 100 keV, μS = 1 μm, μC = 100 nm, μSAD = 50 nm, x = y = 100 nm). Simulations were done in slab geometry as in figure 1. Spectral input means that the micro-beam consists of plane-parallel photons whose energy spectrum is identical to the phase space photons. Phase space photons include the angular components.

Figure 5. DER as function of distance from GNP center r, and angle with respect to the incident direction (E = 100 and 20 keV, μSAD = 5 μm, μS = 5 μm). Simulation geometry 4 is shown on the left.

resulting DER is altered as seen in figure 4. In the latter case, the dose in the water medium without GNP can be underestimated due to the lack of contributions from secondary build-up electrons, which gives rise to a large increase in DER for the parallel beam. Therefore, the effect of the depth in the medium at which the GNP is located cannot be neglected. Computations of the angular anisotropy of DER were performed in spherical geometry using two incident monoenergetic beams of 20 and 100 keV. The geometry and the results are shown in figure 5. The distance to the surface of gold is denoted as r, and the angle with respect to the incident direction is θ . Dose is scored in three angular intervals, θ < π /4, π /4 < θ < 3π /4, and θ > 3π /4, corresponding to forward-, oblique- and backscatter, respectively. It can be seen that the three curves of DER, corresponding to the three angular intervals in each

7972

P Zygmanski et al

Figure 6. GEANT4 benchmarked against CEPXS/ONEDANT using the geometry shown on the left for case 5 (E = 30, 100, 300 keV, r = (r − 2 cm). Radius r = 2 cm is the central radius of the gold shell, which is 100 nm thick.

beam, are almost identical, which means that the dose enhancement is nearly but not perfectly isotropic about the spherical GNP. 3.4. GEANT4 simulation of boundary crossing and benchmark to CEPXS/ONEDANT (case 5) The comparison between Geant4 and CEPXS/ONEDANT for selected monoenergetic photon sources (30 keV, 100 keV, and 300 keV) is shown in figure 6. DER(r) is expressed in terms of distance r from the center of the GNP shell at r = 2 cm. The DER calculated with GEANT4 agrees within a few per cent with that predicted by the deterministic computations. Considering the higher spatial resolution of CEPXS/ONEDANT and the noise in GEANT4, it represents a very good agreement. 3.5. Energy dependence of DER (case 6) The magnitude of enhancement depends on the energy of the x-ray beam. The energy dependence of the DER on the surface of a spherical GNP is shown in figure 7. Spherical symmetry was used to improve the scoring statistics, which is especially important for low energy x-rays below 20 keV, which are mostly attenuated within a depth of 2 cm in the macroscopic medium, resulting in low contribution to the tally voxel. As shown in case 4 (figure 5) the angular anisotropy of DER for low-energy incident photons is not significant, therefore using spherical symmetry is a good approximation. Of particular interest is that the ratio of energy absorption coefficients below 20 keV is nearly identical to the DER in this energy range. 3.6. Dependence of DER on dose contributions from adjacent planar clusters (case 7) In this analysis it is assumed that the x-ray spectra in water is not altered and that selfscreening of GNP clusters can be neglected, which is valid for low concentration of GNPs (e.g. 10 mg g−1 = 1% by mass reported in the GNP literature). For this analysis, we select M planar GNP clusters separated by x distances as in figure 8. Each planar cluster is a rectangular array of μC × μC size, perpendicular to the beam direction. The GNPs within each cluster are separated by C distances, and the size of each GNP is δ × δ × δ. The average (macroscopic) concentration of GNPs is the ratio of the total mass of GNPs divided by the mass of water:

Dependence of MC microdosimetric computations on the simulation geometry of GNPs

7973

Figure 7. DER(E) as a function of energy E of monoenergetic macroscopic isotropic x-rays originating from the center. The GNP shell of radius 2 cm surrounds the source. Dose enhancement ratio is shown for the proximal and distal surfaces as well as for the center of the GNP shell. The gold ratio of mass energy absorption coefficients μen /μwater is also plotted to show similarities with en the DER.

Figure 8. (A) Schematic representation of multiple GNP clusters separated by distances x along the incident beam direction and C across. (B) Dose enhancement DER(x = 0, M) on the distal surface of one of the first GNP clusters as a function of the total number M of the clusters from (A).

mGNP /mwater = 19.3 × (δ × δ × δ)/(C × C × x), where mGNP = M × 19.3 × (δ × δ × δ) × (μC × μC)/(C × C) and mwater = M × 1.0 × (μC × μC × x). Here, the density of gold is taken as 19.3 g cm−3 and that of the water is 1.0 g cm−3. It follows that for the average concentration of mGNP/mwater = 10 mg g−1 the distance between the GNP clusters must be x = 19.3 × 100 × (δ × δ × δ)/(C × C). For instance, when δ = 100 nm and C = 439 nm the distance between adjacent clusters is x = 10 μm. Similarly, for δ = 100 nm and C = 621 nm the distance is x = 5 μm. For smaller GNP sizes, e.g., δ = 10 nm, we have x = 1 μm for C = 44 nm and x = 500 nm for C = 62 nm. These simple examples show that for low macroscopic concentrations of GNPs the distances between GNP clusters are of the order of cell sizes. Thus, assuming that GNPs are distributed approximately uniformly among the different cells, for a specific cell the most important contribution must come from the GNPs within or near the same cell. Figure 8(B) shows an example of DER (x = 0, M) on the distal surface of the first GNP cluster for E = 30 keV, δ = 10 nm and x = 1 μm as a function of the number of clusters M. The total DER from all the clusters increases nonlinearly with M and reaches saturation at

7974

P Zygmanski et al

M ≈ 5–10 (equivalent to the distance of M · x = 5–10 μm from the first to the last slabs). In this case the relative increase of DER due to increasing M from 1 to 10 is only about 30%. 4. Discussion 4.1. The effect of simulation micro-geometry In order to offer a realistic representation of DER values across various MC simulations, a series of computations were carried out to evaluate the impact of the particular simulation geometry. Broader micro-beams were found to be more appropriate, although they required more computation time and weakened the statistics. Narrow micro-beam geometry, such that the micro-beam size is comparable to the size of the GNP, lead to large, often unrealistic, dose enhancement (figure 2). To correctly approximate the macroscopic beam conditions, when using narrow micro-beams large distances between the micro-source and the GNP should be avoided, otherwise particles that scatter laterally escape the region of interest and the DER values may be biased. In the case of plane-parallel spectral photon micro-sources, secondary electron build-up may lead to significantly underestimated dose in water without GNP, thereby artificially increasing the computed DER (figure 4). The aforementioned behavior can be explained as follows. While the magnitude of dose enhancement (DER) nonlinearly increases with μSAD and decreases with μS, simultaneously employing plane-parallel spectral instead of angular phase space micro-source leads to a strong overestimate of DER. The directionality of particles as well as micro-beam size and distance to GNP are related to each other. Their relative contribution to DER depends on specific parameters. It is helpful in this analysis to make a distinction between energy deposition at short versus long ranges and central axis versus off-axis contributions. Energy deposition in the short-range is dominated by low-energy electrons, whereas in the long range it is dominated by higher energy electrons and photons. Photons that travel in lateral direction can escape from the beam but they can also contribute to the dose at the GNP depending on where they originate in and in which direction they are scattered. For instance, figure 2(E) shows the geometry effect for an angular phase space microbeam of the same size as the nanoparticle (μS = 100 nm). Here, the distance from GNP is relatively large (μSAD = 1 mm). In this geometry there is very little contribution to the central axis dose from the off-axis locations in the surrounding water. Furthermore, there is a strong lateral disequilibrium. Therefore, DER is dominated by the central axis x-rays within the narrow beam that move forward and are absorbed by the GNP. Furthermore, without the GNP present the dose in water is smaller than in equilibrium. Hence the combined effect is that the DER is artificially high. The opposite happens for a broad micro-beam (μS = 10 μm μC = 100 nm). In the latter case, there are relatively large contributions to dose from off-axis points within the broad micro-beam and DER is not so much dominated by the central axis photons and therefore DER is smaller. When plane-parallel spectral photon, instead of phase space micro-source, is used, the electron build-up represents an additional confounding factor. The lack of primary electrons in the micro-beam (due to ignoring transport through the medium up to a short distance away from the GNP) leads to lower dose in water without GNP. In the presence of GNP, the electron production is solely due to photon interactions in the GNP. The combined effect is exaggerated DER. GNP DERs were calculated in various geometrical configurations irradiated with different monoenergetic photons varying from 11 keV to 1 MeV. The maximum DER (in the order

Dependence of MC microdosimetric computations on the simulation geometry of GNPs

7975

of 100) was observed for a planar densely packed GNP cluster when irradiated by low keV photons. This result is lower than the values of microscopic dose enhancement factors (mDEF) or the values of interaction enhancement ratios (IER), reported in the literature (Jones et al 2010, Leung et al 2011), which range up to 1000 and more. Our results show that the computation of the dose enhancement is sensitive to the parameters of the simulation geometry, including the distance between the nanoparticles C, which have to be explicitly considered. These discrepancies can be explained by the different simulation geometry and the different definition of DER compared to prior literature in which mDEF (Jones et al 2010) and IER (Leung et al 2011) were defined. Specifically, mDEF was defined in terms of the total secondary electron contributions from all the uniformly distributed gold atoms averaged over a macroscopic volume of water (first stage simulation) and assumed to be originating from a point source in water (second stage simulation). On the other hand, IER was defined for a plane-parallel micro-beam of the same size as the GNP diameter, located 1 mm away from the GNP, and in terms of the total number of interactions in a single GNP. Thus, mDEF is an energy deposition ratio derived from the macroscopic simulation of uniform concentration of gold atoms, and IER is a quantity not directly related to spatial distribution of energy deposition. Furthermore, the transition from macroscopic to microscopic simulation is achieved in a number of ways with various renormalizations of the scored quantities and therefore a direct comparison of these quantities to the microscopic DER presented herein is not possible. Finally, it has to be noted that in some models a linear dependence on GNP concentration was assumed (Ngwa et al 2010, Berbeco et al 2011). This assumption is based on macroscopic considerations under charged particle equilibrium and that the inter-particle effects, such as mutual shielding, are negligible, which can lead to further overestimation of the DER. In contrast, in our work we study the local DER due to energy deposition in nanoscopic voxels surrounding a single GNP or a cluster of GNPs. In this way, the magnitude of the DER increases with the number of GNPs per cluster (N × N) but it is not linearly proportional to them. This nonlinear relation asymptotically reaches saturation at a finite cluster size. This behavior seems intuitive because at microscopic distances among a cluster of GNPs the contribution from the nearest GNP is dominant and the dose fall-off with distance is nonlinear. If one recalls that the cell sizes range from about 2 μm (endothelial cells) to 10 μm (somatic cells) and that the DER from a planar cluster is approaching unity (DER→1) for distances larger than few μm, the nearest GNPs are essentially those that pertain to a given cell (GNPs inside and outside of the cell). This geometrical effect is readily seen in figure 8(B) and it has impact on cell survival as investigated in a recent radiobiological study examining stochastic effects due to uneven energy distribution among a population of cells with variable configurations of GNPs (Zygmanski et al 2013). Under non-uniform GNP distribution (the number of GNPs and their distance to specific cellular targets) such stochastic effects may lead to decreased radiobiological effectiveness of GNP therapy. Note that we simulated GNPs of various cross-sectional sizes (μC) and fixed thickness (δ = 100 nm). Specific numerical results depend on the size of the nanoparticles, however, certain trends remain irrespective of the size. For the case of a single nanoparticle, the larger the cross-sectional size and thickness the larger the probability for x-ray interaction with GNP and the larger the DER values. But the characteristic shape of DER with two distinct ranges (due to Auger and photoelectric) remains unaltered. Similarly, the dependence of MC simulations on the parameters of micro-beam geometry we address in the paper is still applicable. For the case of GNP cluster, the larger the cross-sectional size of the individual nanoparticle, the larger the overall DER of the cluster, but the dependence on DER on the relative distances between the nanoparticles C and their number N × N will show general behavior.

7976

P Zygmanski et al

4.2. Statistical uncertainties Other uncertainties are due to the stochastic nature of energy deposition and depend on the number of histories, simulation geometry parameters and dose magnitude within a given voxel. We have kept the number of histories in the range from 108 to 2 × 109. For a given number of histories, statistical uncertainty increases with micro-beam size μS, distance from the origin of the micro-source to GNP μSAD (due to attenuation and scatter of x-rays), and distance from GNP to the location of a given voxel (due to electrons from GNP). On the other hand, uncertainty decreases with dose and the size of GNP (or the number of GNPs in a cluster). Furthermore, because DER is a dose ratio its uncertainty is a geometric combination of the two sources of uncertainties for dose in pure water and water with GNPs. Typically the voxel to voxel deviation in DER was up to a few per cent. However, for large values of μS and μSAD, and distances away from GNP and a single 100 nm GNP voxel-wise variations in DER were reaching up to 7%–25%. But this was for regions that are microns away from the GNP, where DER was close to 1.0, thus had little significance from the point of view of the present study. Smaller variation in DER can be obtained by smoothing dose in water before calculation of DER ratio. This strategy is possible based on the fact that dose in water is a smooth function of distance. Dose in water is the main source of statistical uncertainties because it is much smaller than dose in water with GNP. However, our DER data does not include this type of smoothing except when we compared DER to CEPXS/ONEDANT, which was motivated by the fact that this code employs a deterministic method and does not bear statistical noise. 4.3. Energy dependence The shape of the energy dependence of DER resembles the one of the gold photoelectric cross-section, however for energies less than about 20 keV there is a sudden drop of the DER. For energies below 20keV a better resemblance is obtained with the mass energy absorption μwater is also plotted in figure 7. One can observe the K- and coefficient ratio μgold en en , which water L- edges in both DER and μgold μ en en . In charged particle equilibrium the dose can be well approximated as the product of the mass energy absorption coefficient and the energy fluence. Thus, this may indicate that for very low energy photon beams, below 15 keV, charged particle equilibrium is approached but not so for larger energies. More information about Auger and photoelectron contributions to DER is provided in the supplementary material (Online supplementary data available from stacks.iop.org/PMB/58/7961/mmedia). 5. Conclusion We used GEANT4 to study dose enhancement of GNPs at the nanometer-to-micrometer distances with high spatial resolution. Macroscopic monoenergetic beams were transported in macroscopic water medium to generate microscopic beams whose phase space inherited that of the macroscopic beam. The microscopic beams were transported over various microscopic distances away from single GNPs or GNP clusters. Results obtained using GEANT4 were compared to those using the deterministic discrete ordinates computer code CEPXS/ONEDANT, and showed very good agreement. Dose enhancement was found to be almost isotropic about a spherical GNP but its magnitude strongly depended on the simulation micro-geometry (μS, μSAD) and phase space properties (plane-parallel versus full angular phase space). Dose enhancement nonlinearly scaled with the number of nanoparticles in a given planar cluster or the number of adjacent planar clusters. Therefore, the magnitude of DER is dominated by the nearest GNPs in a planar cluster or nearest planar clusters, and is less

Dependence of MC microdosimetric computations on the simulation geometry of GNPs

7977

influenced by further neighbors. It is shown that a careful simulation of the irradiation geometry is required to avoid artifacts that tend to exaggerate the computed DER. The strong dependence of DER on the simulation micro-geometry cautions future studies and the interpretation of existing MC results obtained in different simulations geometries. Further studies are required to determine the impact of macroscopic concentrations of nanoparticles on beam hardening, however, these studies can be carried out using macroscopic radiation transport and will not alter the main conclusions of the present study. References Berbeco R I, Ngwa W and Makrigiorgos G M 2011 Localized dose enhancement to tumor blood vessel endothelial cells via megavoltage x-rays and targeted gold nanoparticles: new potential for external beam radiotherapy Int. J. Radiat. Oncol. Biol. Phys. 81 270–6 Chauvie S et al 2004 Geant4 low energy electromagnetic physics IEEE Nuclear Science Symp. Conf. Record vol 3 pp 1881–5 Guatelli S, Mantero A, Mascialino B, Nieminen P and Pia M G 2007 Geant4 atomic relaxation IEEE Trans. Nucl. Sci. 54 585–93 Jaeger R G, Blizard E P, Chilton A B, Grotenhuis M, Hoenig A, Jaeger T A and Eisenlohr H H 1968 Engineering Compendium on Radiation Shielding vol II (New York: Springer) Jeong S-Y et al 2009 Systemic delivery and preclinical evaluation of au nanoparticle containing beta-lapachone for radiosensitization J. Control. Dis. 139 239–45 Jones B L, Krishnan S and Cho S H 2010 Estimation of microscopic dose enhancement factor around gold nanoparticles by Monte Carlo calculations Med. Phys. 37 3809–16 Lechtman E, Chattopadhyay N, Cai Z, Mashouf S, Reilly R and Pignol J P 2011 Implications on clinical scenario of gold nanoparticle radiosensitization in regards to photon energy, nanoparticle size, concentration and location Phys. Med. Biol. 56 4631–47 Leung M K K, Chowa J C L, Chithrani B D, Lee M J G, Oms B and Jaffray D A 2011 Irradiation of gold nanoparticles by x-rays: Monte Carlo simulation of dose enhancements and the spatial properties of the secondary electrons production Med. Phys. 38 624–31 Lorence L J Jr and Beutler D E 1997 Radiation transport phenomena and modeling, part a: codes; part b: applications and examples Sandia Report SAND97-2135 (Albuquerque, NM: Sandia National Laboratories) Lorence L J Jr, Morel J E and Valdez G D 1989 Physics Guide to CEPXS: a multigroup coupled electron–photoncrosssection generating code Sandia Report SAND89-1685 (Albuquerque, NM: Sandia National Laboratories) McMahon S J et al 2011 Biological consequences of nanoscale energy deposition near irradiated heavy atom nanoparticles Sci. Rep. 1 1–9 Ngwa W, Makrigiorgos G M and Berbeco R I 2010 Applying gold nanoparticles as tumor-vascular disrupting agents during brachytherapy: estimation of endothelial dose enhancement Phys. Med. Biol. 55 6533–48 Rahman W N, Bishara N, Ackerly T, Cheng F H, Price J, Wong C, Davidson R and Geso M 2009 Enhancement of radiation effects by gold nanoparticles for superficial radiation therapy Nanomedicine 5 136–42 Tsiamas P, Sajo E, Cifter F, Theodorou K, Kappas C, Makrigiorgos G M, Marcus K and Zygmanski P 2013b Beam quality and dose perturbation of 6 MV flattening-filter-free linac Eur. J. Med. Phys. at press Tsiamas P, Sajo E, Cifter F, Theodorou K, Kappas C, Makrigiorgos G M, Marcus K and Zygmanski P 2012 SU-E-T-35: Optimal clinical megavoltage x-ray beam quality for contrast enhanced RT (CERT) Med. Phys. 39 3710 Tsiamas P et al 2013a Impact of beam quality on megavoltage radiotherapy treatment techniques utilizing gold nanoparticles for dose enhancement Phys. Med. Biol. 58 451–64 Williams M L and Sajo E 2002 Deterministic calculations of photon spectra for clinical accelerator targets Med. Phys. 29 1019–28 Zheng Y, Hunting D J, Ayotte P and Sanche L 2008 Radiosensitization of DNA by gold nanoparticles irradiated with high-energy electrons Radiat. Res. 169 19–27 Zheng Y and Sanche L 2009 Gold nanoparticles enhance DNA damage induced by anti-cancer drugs and radiation Radiat. Res. 172 114–9 Zygmanski P, Hoegele W, Tsiamas P, Cifter F, Ngwa W, Berbeco R, Makrigiorgos M and Sajo E 2013 A stochastic model of cell survival for high-Z nanoparticle radiotherapy Med. Phys. 40 024102