Biological modeling of gold nanoparticle enhanced radiotherapy for proton therapy

Monte Carlo simulations were carried out using TOPAS (tool for particle simulation; version beta 12) (Perl et al 2012), which is layered on top of Geant 4, version 9.6p02 (Agostinelli et al 2003). The simulations were split into three steps as follows.

Step 1. A 20   ×   20   ×   40 cm³ water phantom was irradiated with three types of particle sources and the particle spectrum at a given depth (indicated in figure 1) was recorded as a phase space profile passing a 50 mm diameter surface perpendicular to the beam axis. The phase spaces contain all primary and secondary particles. The simulation details for each particle source are:

Zoom In Zoom Out Reset image size

Spread out Bragg peak (SOBP) proton beam.

Pristine Bragg peak proton beam.

6 MV photon beam.

Kilovoltage photon beam.

A tungsten target with a 2 mm aluminum filter was simulated which produces a 150 and 250 kVp photon spectrum impinging on the water phantom. The depth at 1 mm below water surface was selected to obtain the particle phase space, as indicated in figure 1(d).

Step 2. Four sizes of spherical GNPs were investigated with diameters of 50, 20, 10 and 2 nm. To address the low chance of hitting a GNP in the water phantom, the phase space file obtained from Step 1 was modified by adjusting both the length scale and angular distribution (see (Lin et al 2014) for more details) to have the same beam diameter as a single GNP as shown in figure 2. The electrons produced through ionizing interactions in the GNP were recorded in a phase space on the outer surface of the GNP. The Geant4-Penelope physics list was used for this simulation step and electrons were tracked down to an energy of 250 eV in gold (Bernal and Liendo 2009). The secondary electron production threshold was set to 800 eV. The number of particles interacting within the GNP was calculated as the difference between the number of particle tracks that entered the GNP and the number of tracks that emerged from the GNP without interacting. This value was then used to calculate the probability of having an ionization event inside the GNP as described in section 2.2.

Zoom In Zoom Out Reset image size

Step 3. For this step, we used the phase-space recorded in Step 2 as a source, which was placed at the center of a 40   ×   40   ×   40 µm³ water phantom. The Geant4-DNA physics package (Incerti et al 2010) was used, which allows tracking of electrons down to an energy of 7 eV in water. The dose was scored in a voxelized grid with a voxel size of 1   ×   1  ×  1 nm³ for the first 200 nm from the center, 10   ×   10   ×   10 nm³ from 200 nm to 2 μm and 100   ×   100   ×   100 nm³ from 2 to 20 μm. For a more detailed description of the setup and simulation parameters see (Lin et al 2014).

We can calculate the interaction probability per Gray with the geometrical concept illustrated in figure 2.

For 1 Gy prescribed dose, the chance of hitting a single 50 nm GNP is given by:

$$\begin{eqnarray}&&{{P}_{\text{hit}}}={{N}_{1Gy}}\times {{\left(\frac{50\text{nm}}{50\text{mm}}\right)}^{2}}\end{eqnarray}$$

where N_1Gy represents the number of particles passing a 50 mm diameter circle that deposit 1 Gy in a water phantom. The dose used in proton dose calculation is 1 Gy physical instead of 1 Gy(RBE) biological equivalent dose. The same number of particles is used to irradiate a single GNP as described in Step 2. The probability of having an ionization event by a single incoming particle is calculated as:

$$\begin{eqnarray}&&{{P}_{\text{ionization}}}=\frac{{{N}_{\text{tracks}}}}{{{N}_{1Gy}}}\end{eqnarray}$$

where N_tracks is the number of particle tracks for N_1Gy incoming particles that interacted in the GNP volume and resulted in ionized atoms. N_tracks is also defined as the number of ionization events in GNP. The interaction probability is then defined as the product of the probability of hitting a single GNP and the probability of causing an excitation event for one particle:

$$\begin{eqnarray}&&{{P}_{\text{total}}}={{P}_{\text{hit}}}\times {{P}_{\text{ionization}}}\end{eqnarray}$$

The local effect model (LEM) originally developed for carbon ion therapy seeks to relate the dose response curve following high linear energy transfer (LET) radiation to a reference x-ray radiation (Scholz and Kraft 1994, Scholz et al 1997, Elsässer and Scholz 2007). In order to consider the secondary delta electrons, the radial dose distribution produced by secondary electrons around the primary particle track is used as input. The dose deposited as a function of the distance to the ion path (r) has a 1/r² dependence, indicating a high dose near the ion track followed by a rapid dose decrease as r increases.

The dose in close vicinity of a GNP is particularly high due to the large number of secondary electrons ejected from the GNP. The dose decreases rapidly as the distance from the GNP surface increases (McMahon et al 2011a, Lin et al 2014). The LEM concept thus provides a framework to predict the dose response curve for GNP radiosensitization by replacing the radial dose deposition around an ion track with the radial dose deposition around a GNP. In our approach, 'GNP-LEM', three inputs are needed, 1) a spatial dose distribution around a GNP for a given particle source, 2) the Linear-Quadratic (LQ) model to parameterize the x-ray dose response curve into α and β components and 3) a geometrical model of the cell to characterize the volume of the nucleus and the distribution of the GNPs.

Within the LEM framework, the number of lethal events N as a function of dose (D) follows a Poisson distribution, thus the surviving fraction as a function of photon dose (S_x(D)) is determined by

$$\begin{eqnarray}&&{{S}_{x}}(D)={{\text{e}}^{-N(D)}}\end{eqnarray}$$

Because the LQ model overestimates the response at high doses (Astrahan 2008) above a threshold dose D_t, the dose response curve reaches a maximum slope given by ${{S}_{\text{max}}}=\alpha +2\beta {{D}_{t}}.$ As a result, the dose response curve can be represented as a linear quadratic linear (LQL) model:

$$\begin{eqnarray}s(D)=\left\{\begin{array}{c} {{\text{e}}^{-\alpha D-\beta {{D}^{2}}}} & :D\leq {{D}_{t}} \\ {{\text{e}}^{-\alpha {{D}_{t}}-\beta {{D}_{t}}^{2}}}{{\text{e}}^{{{S}_{\text{max}}}\left(D-{{D}_{t}}\right)}} & :D>{{D}_{t}} \\ \end{array}\right.\end{eqnarray}$$

For a homogeneously irradiated cell nucleus volume V, the lethal event density ν(D) is

$$\begin{eqnarray}&&\upsilon (D)=\frac{N(D)}{V}=\frac{-\ln \left({{S}_{x}}(D)\right)}{V}\end{eqnarray}$$

For a subvolume ΔV, the number of lethal events is

$$\begin{eqnarray}&&{{\bar{N}}_{\text{lethal}}}=\upsilon (D)\times \Delta V=N(D)\times \frac{\Delta V}{V}=-\ln {{S}_{x}}(D)\times \frac{\Delta V}{V}\end{eqnarray}$$

For an inhomogeneous radiation field, the total number of lethal events N_total can be calculated by integrating the lethal events over the whole target volume based on the distribution of local doses within the nucleus, which gives

$$\begin{eqnarray}&&{{\bar{N}}_{\text{total}}}={\int}_{V}\frac{-{{S}_{x}}\left(\text{d}\left(x,y,z\right)\right)}{V}\text{d}V\end{eqnarray}$$

where d(x,y,z) represents the distribution of local doses. The dose response curve in the presence of GNPs can then be calculated by

$$\begin{eqnarray}&&{{S}_{\text{GNP}}}={{\text{e}}^{-{{{\bar{N}}}_{\text{total}}}}}\end{eqnarray}$$

In the presence of GNPs, the local dose distribution is represented by a homogeneous dose from the delivered treatment field in addition to the dose deposited in the vicinity of GNPs (dose_GNP). The dose_GNP is determined by multiplying the dose from a single ionizing event by the number of GNPs, the interaction probability per Gray and the prescribed dose.

The GNP-LEM developed in this study was implemented in 2D, where the volume integration is reduced to an area integration over the cell nucleus. Cells appear to be flattened out when they attach to a petri dish as compared to in cell suspension and the thickness of the cell is set to be 2 µm.

The dose response relationship of MDA-MB-231 breast cancer cells was considered in this work. The parameters for the linear-quadratic equation were α = 0.019 and β = 0.052 (Jain et al 2011). In addition, two other α and β values representing a high and low α/β ratio cell line were considered, which are DU145 prostate cancer cells with α = 0.13 Gy⁻¹, β = 0.022 Gy⁻² and α/β = 6 Gy, and human brain tissue with α = 0.1 Gy⁻¹, β = 0.05 Gy⁻² and α/β = 2.0 Gy (Van der Kogel 1986).

The size of the cell and nucleus varies depending on cell lines and cell cycle. For instance, MDA-MB-231 cells have an average size of 13.5 µm (Coulter et al 2012) while PC-3 cells have an average cell size of 13.1 µm and nucleus size of 8.2 µm (Lechtman et al 2013). For this study, we used 13.5 µm as cell diameter with an 8 µm diameter radiation sensitive nucleus at the center, which is in good agreement with the range of many cancer cell lines (Zhao et al 2008, Jiang et al 2010). An additional 2.5 µm thickness was added outside of the cell to represent the extracellular media region.

The following five GNP geometries were investigated:

• (a)  
  Nucleus: GNPs are randomly distributed in the cell nucleus only (figure 3(a)).
• (b)  
  CellHomo: GNPs are randomly distributed in the whole cell including the nucleus and cytoplasm space (figure 3(b)).
• (c)  
  Cytoplasm: GNPs are randomly distributed in the cytoplasm only (figure 3(c)).
• (d)  
  Media: GNPs are randomly distributed in the extracellular media only (figure 3(d)).
• (e)  
  Complex: GNPs are randomly distributed inside of the cell and extracellular media with the same concentration. Within the cell, 90% of GNPs are distributed in the cytoplasm and 10% of GNPs are distributed throughout the whole cell volume including the nucleus (figure 3(e)).

Zoom In Zoom Out Reset image size

The threshold dose D_t is an empirical and cell line dependent parameter, and may be approximated for some cell lines by experimental or clinical data (Astrahan 2008). A D_t of 20 Gy is used for this study but the effect of several D_t values ranging from 10 to 30 Gy as described in the literature (Krämer and Scholz 2000, Paganetti and Goitein 2001, Astrahan 2008) was also considered.

Dose response curves were calculated from 0–8 Gy for scenarios covering a range of GNP distributions, sizes and concentrations for the three irradiation modalities. The sensitizer enhancement ratio (SER) was used to quantify the enhancement. SER was calculated from the mean inactivation dose (MID), defined as the area under the dose response curve. The SER was obtained by dividing the MID for the control sample by the MID in the presence of GNPs. The MID for the control sample is defined as the MID of the survival curve without GNPs added. The MID for gold-treated samples was obtained from the dose response curve calculated by GNP-LEM.

• (a)  
  For each particle source, we investigated the GNP SER for the five GNP geometries as described in section 2.4. A concentration of 1 µM using 50 nm diameter GNPs results in 1.4   ×   10⁵ GNPs for the Nucleus, CellHomo and Cytoplasm geometries (based on a cylindrical volume of 13.5 μm diameter and 2 μm thickness). For the Complex case, there were 1.26   ×   10⁵ GNPs in the cytoplasm and 1.4   ×   10⁴ GNPs homogenously distributed outside the cell.
• (b)  
  We then repeated a) using 20, 10 and 2 nm GNPs for the CellHomo geometry, once assuming the same number of GNPs and once assuming the same weight of GNPs.
• (c)  
  The SER as a function of GNP concentration was investigated using 50 nm GNPs. We investigated two GNP geometries (the CellHomo and Cytoplasm) for each modality. For the kV photons, the considered range of GNP concentration was between 10 nM and 1 µM, while for the proton source and 6 MV Linac photon source, the considered concentration ranged from 100 nM to 10 µM.

The interaction probability per prescribed Gray as defined by P_total in section 2.2 for the three particle sources is shown in table 1. Figure 4(a) shows the dose per ionization event as defined in section 2.2 by the average dose deposited per incoming particle track that interacts in a GNP. It has been shown that the GNP dose using 6 MV photons has a considerable contribution from electron ionization as well as photon ionization (McMahon et al 2011a, Lin et al 2014). For the 6 MV photon beam, the dose per ionization is plotted for the electron and photon contribution individually. Within 10 nm from the GNP surface, the dose per ionization for the 6 MV particle spectrum electron component is the highest, but drops off below the photon component after 100 nm. The dose deposited at the GNP per ionization agrees within 20% for the proton beam, the 250 kVp photon beam and the photon component of the 6 MV photon beam. As the distance from the GNP surface increases, the dose from proton irradiation drops off fastest. At a distance of 1 μm from the GNP surface the dose from 250 kVp photons and the photon component of the 6 MV photon beam is ten times higher than the dose from proton irradiation. The difference increases further and at a distance of 10 μm from the GNP surface and beyond, the dose from 250 kVp photons is 30 times higher than the dose from proton beam irradiation at the center of the SOBP (SOBP2).

Zoom In Zoom Out Reset image size

Beside the differences in spatial dose deposition per ionization interaction for different treatment modalities, the interaction probability also differs by up to two orders of magnitude as shown in table 1.

The kilovoltage photon source has the highest interaction probability with GNPs per Gray, and it is energy dependent. 150 kVp photons have a 30% higher interaction probability than 250 kVp photons.

The MV photon source has a much lower interaction probability compared to kV photon and the probability increases with increasing depth. The interaction probability at PDD50 is nearly twice as high as at PDD100 depth. This increase is mainly due to the beam softening due to scattering, as shown in our previous study (Lin et al 2014).

The interaction probability for the considered proton source is on the same order of magnitude as for the MV photon source. For a SOBP proton beam, the interaction probability is independent of depth, and about 20% higher than in a pristine (mono-energetic) proton beam.

The interaction probability per prescribed Gray for four GNP sizes at the center of the proton SOBP (SOBP2) is shown in table 2. The dose per ionization event for a proton beam impinging on GNPs of various sizes is shown in figure 4(b). For one ionization event, 2 nm GNPs produce nearly a 215 times higher dose near the surface of the nanoparticles compared to 50 nm GNPs. The difference is reduced to less than 10% as the distance from the surface increases to 1 µm. For the same energy absorbed by a single GNP, the secondary electrons generated in a large GNP are more likely to lose their energy before reaching the surface. Such self-absorption contributes to the lower dose deposited around the GNP by one ionization event for larger GNPs.

The results are shown in figure 5 and table 3. The relative biological effect (RBE) of protons and different photon radiation types is not considered in the current study, and the enhancement effect is calculated solely based on the addition of GNPs. The survival curve for 'noGNP' for all cases is produced from a homogeneous dose in the nucleus area using α and β values from the photon survival curve.

Zoom In Zoom Out Reset image size

For all three irradiation modalities, GNPs distributed according to the Nucleus geometry resulted in the highest GNP enhancement. In general, GNPs located inside of the cell nucleus contribute most significantly to the cell killing.

3.3.1. For protons.

3.3.2. For kV photons.

Kilovoltage photons cause significant GNP enhancement at as low as 100 nM GNP concentration. The SER is 5.97 and 4.88 for Nucleus and CellHomo geometry (figure 6(a)). For kV photon, GNP sensitization is more effective for a 150 kVp photon source (SER = 4.88) than a 250 kVp source (SER = 4.04) for the CellHomo geometry (see figure 6(c)). Even when no GNPs are internalized in the cell nucleus, the SER is 1.77 and 1.56 for 150 and 250 kVp x-rays, respectively. In addition, at a GNP concentration of 100 nM in the media, the SER is only 1.08 and 1.06 for 150 and 250 kVp x-rays, respectively, but increased to 1.82 and 1.62 when the concentration increased to 1 µM.

Zoom In Zoom Out Reset image size

For the proton SOBP at the depth of SOBP2, the SER only increases from 1.001 to 1.008 when the GNP concentration in the media increases from 1 to 10 µM. As shown in figure 4(a), at a distance of 1 μm from the GNP surface the dose induced by the GNP from 250 kVp photons is 10 times more than when using proton irradiation. For instance, the GNP dose from kV photon irradiation is approximately 0.05 Gy per ionization events at 1 µm from GNP surface, while the GNP dose is about 0.005 Gy per ionization for proton irradiation. Therefore, the GNPs present in the media contribute to the cell damage more effectively for kV photons than protons. When the cells are treated with a high concentration of GNPs, sensitization can be achieved for kV photons even without cellular uptake of GNPs.

3.3.3. For MV Photons.

GNP sensitization is more effective at deeper depths. The SER increased from 2.59 at PDD100 to 4.74 at PDD50 for the CellHomo geometry at 1 µM GNP concentration (see figure 7(c)). Like protons, 6 MV photons have a lower interaction probability than the kV photons, and thus we list the SER for 1 µM GNP concentration in table 3.

Zoom In Zoom Out Reset image size

The effect of GNP sizes on radiosensitization in proton beams is shown in figure 8 and table 4. The dose response curve is calculated for proton irradiation at the center of the SOBP (SOBP2) according to the CellHomo geometry. The enhancement reduces with decreasing GNP size when using the same number of GNPs inside a cell because of the difference in the probability of GNP particle interaction for a given particle fluence (see table 2). On the other hand, if the same weight of GNPs is taken up by the cell, smaller GNPs cause more damage. Smaller GNPs have a much higher dose per ionization at their surface due to lower internal absorption of secondary electrons by the GNPs, as shown in figure 4(b). As a result, 10 nm GNPs cause about a factor 3 more damage than the same weight of 50 nm GNPs.

Zoom In Zoom Out Reset image size

We evaluated the impact of the choice of D_t. The dose response curve is calculated for proton irradiation at the center of the SOBP (SOBP2) with 1 µM GNPs distributed according to the CellHomo geometry. 50 and 2 nm GNPs are used for the calculation. The results are shown in figure 9 and table 5. Three values of D_t, i.e. 10, 20 and 30 Gy are compared for proton beams. The results are displayed for 2 nm GNPs and 50 nm GNPs.

Zoom In Zoom Out Reset image size

The choice of D_t does not affect the results for 50 nm GNPs significantly (see figure 9(a)). Changes in the SER are within 5% as shown in table 5. However, for 2 nm GNPs the SER increased from 3.24 to 4.32 when D_t increases from 10 to 30 Gy. As shown in figure 4(b), the dose per ionization near the GNP surface for the 2 nm GNPs is more than 100 times higher than the dose for 50 nm GNPs and is in the range of several hundreds or thousands of Gray within several nanometers of the GNP surface. Truncating the LQ model for a linear transition at the lower dose (small choice of D_t) therefore has a higher impact on 2 nm GNPs than 50 nm GNPs.

The radiosensitization for three cell lines with different α and β values were evaluated in the presence of GNPs, while the cell geometry is left unchanged. The three cell lines are: MDA-MB-231 breast cancer cells with α = 0.019 Gy⁻¹, β = 0.052 Gy⁻² and α/β = 0.37 Gy (Jain et al 2011), DU145 prostate cancer cells with α = 0.13 Gy⁻¹, β = 0.022 Gy⁻² and α/β = 6.0 Gy, and human brain tissue with α = 0.1 Gy⁻¹, β = 0.05 Gy⁻² and α/β = 2.0 Gy (Van der Kogel 1986). The dose response curve is calculated for proton irradiation at the center of the SOBP (SOBP2) with 50 nm GNPs distributed according to the CellHomo geometry. As shown in figure 10, the change in SER is about 7% for the three different cell lines, while the survival fraction at 8 Gy is two orders of magnitude smaller for MDA-MB231 and DU145 cells and one order of magnitude smaller for cells from the human brain.

Zoom In Zoom Out Reset image size

The results are shown in figure 11. For kV photons, the GNP enhancement increases quickly as the GNP concentration increases. A GNP concentration of 80 and 120 nM is needed to reduce the 2 Gy survival fraction to 0.5 for the Complex and Cytoplasm geometry, respectively.

For 6 MV photons, a GNP concentration of 4.4 µM is needed to reduce the 2 Gy survival fraction to 0.5 for the Cytoplasm geometry. A much lower concentration of 1.3 µM is needed to reduce the 2 Gy survival fraction to 0.5 for the Complex geometry.

For protons, a much higher concentration of GNPs beyond 10 µM is needed to reduce the 2 Gy survival fraction to 0.5 for both Complex and Cytoplasm geometry. However, if the GNP can be internalized into cell nucleus, a GNP concentration of 1.9 µM is needed to reduce the 2 Gy survival fraction to 0.5 for the Nucleus geometry.