Using Geant4 Application for Tomographic Emission (GATE), we first simulated three different single cell models with activity in the cytoplasm (to represent internalized activity) or membrane (to represent surface bound activity) of each model to determine the impact of cell geometry and source location on nucleus absorbed doses as an indicator of deoxyribonucleic acid (DNA) damage. We expanded two of the single cell models to a cell layer because of the extra-cellular range of some particulate emissions, and created nucleus absorbed dose kernels that can be combined in different ratios of cytoplasm to membrane uptake to model radiopharmaceuticals with different uptake patterns. The kernels were applied to 2D multicellular prostate cancer tumour maps (with varying amounts of normoxic, hypoxic, and necrotic regions) to obtain absorbed dose rate maps within the tumours. The impact of cell geometry, internalization ratio, and tumour morphology on the absorbed doses and absorbed dose distributions within the tumours was assessed for both 177Lu and 225Ac.

Simulation validation

To validate our simulations, we used MIRDcell V3.10, a multicellular dosimetry tool [19]. MIRDcell analytically calculates the absorbed doses from beta and alpha particles using the continuous slowing down approximation in simple spherical cellular geometries. In MIRDcell, we created two spherical cells with 14 µm diameters directly adjacent to each other, each with a 4 µm radius nucleus within a 7 µm radius cytoplasm. Homogeneous activity of 225Ac or 177Lu was placed in the cytoplasm of one cell and the self-dose per decay to the nucleus and the cross-dose per decay to the neighbor nucleus was calculated. In GATE, two identical water density cells were defined, and the same absorbed dose values were scored in the nuclei to compare against MIRDcell. Details of the GATE simulations are described below.

Monte Carlo simulations

The Monte Carlo simulation software GATE version 9.0, based on Geant4 version 10.6.2, was used for all simulations. For simulations of 225Ac, the Geant4-DNA emDNAphysics physics list was used which allows the simulation of physical interactions of charged particles (e.g. beta and alpha particles) down to very low energies (7.4 eV) in liquid water [20]. The daughter emissions of 225Ac were included in the simulations. For 177Lu, the emlivermore physics list was used which allows the simulation of internal conversion electron and beta particle interactions with matter down to about 250 eV. While the emlivermore list enabled us to obtain accurate results with 177Lu in less time than the emDNAphysics list, it was not satisfactory for 225Ac due to not simulating alpha particles [21]. For both 225Ac and 177Lu, the GATE “ion source” was used, which simulates the full photon spectrum and all charged particle emissions of the simulated radionuclide and any daughter radionuclides.

All simulations were split into 100 jobs which were run in semi-parallel. No cuts (cut-off values below which particles are no longer tracked, which are converted to energies by GATE) were used for 225Ac as the emDNAphysics list already simulates sufficiently low energies, while cuts of 0.1 µm were used for simulations of 177Lu. The source and geometry definitions for the cell geometry comparisons and nucleus absorbed dose kernels are found below in the respective sections.

Cell geometry comparison

To determine the effect of cell geometry on absorbed doses to the nucleus, we modelled three cell geometries: a spherical, ovoid, and cubic cell. All models had a nucleus, nuclear membrane, cytoplasm, and cell membrane which were water density, and the cells were surrounded by water. The total cell and nucleus volumes were approximately identical in all models (1287 µm3 for the total cell and 377 µm3 for the nucleus) and were based on cross-sectional areas of human lymph node carcinoma of the prostate (LNCaP) cells derived from human prostate cancer cells [22].

The three different geometries were modelled as such:

The spherical model consisted of concentric spheres with radii of 4.483 µm and 6.74 µm representing the nucleus and cytoplasm respectively.

The cubic model had a nucleus and nuclear membrane identical to the spherical model while the cytoplasm was a 10.87 µm length cube.

The ovoid model had an ellipsoidal nucleus with dimensions of 7.7 µm × 5.4 µm × 17.4 µm and a cytoplasm with dimensions of 12.0 µm × 8.0 µm × 25.6 µm.

In all models, the nuclear and cellular membranes were 5 nm thick and enveloped the nucleus or cytoplasm respectively. Figure 1 shows each cell model.

Fig. 1

The three cell geometries analyzed in this work, the spherical (left), cubic (middle), and ovoid (right) cell. Images were generated using visualization options in GATE

The GATE DoseActor tool was attached to each nucleus to score the deposited energy and ultimately calculate the absorbed radiation dose within each voxel. The DoseActor outputs 3D images of the deposited energy (“edep”) in units of MeV and the associated uncertainty. The absorbed dose was calculated from the deposited energy outputs as described below in Sect. "S-value calculation". For each single cell geometry, the DoseActor had a voxel size of 0.2 µm × 0.2 µm × 0.2 µm. 1 × 109 and 1 × 107 primaries were simulated for 177Lu and 225Ac respectively. Simulations with the emDNAphysics list are more time intensive (taking approximately 6 days to simulate 1 × 107 primaries of 225Ac, compared to approximately 30 min to simulate 1 × 109 primaries of 177Lu with the emlivermore list) and 1 × 107 primaries was the maximum amount that was feasible with 225Ac.

Nucleus absorbed dose kernels

While human tissue consists of layers or clusters of cells, simulating thousands of cells with GATE is computationally intensive and time consuming. To best replicate this biological environment with our simulation constraints, we created a 2D multi-purpose nucleus absorbed dose kernel which describes the absorbed dose deposited in each cell nucleus in a layer of cells when activity is present in the cell at the center of the grid. This kernel can be convolved with simulated or real tissue activity images to move from activity values to absorbed dose rate values, enabling dosimetry and analysis of tissue heterogeneity. We created kernels for both 177Lu and 225Ac, both source locations, and with both spherical and cubic cell geometries to investigate the effect of cell clustering (i.e. the presence of gaps in between the cells) on kernel values.

To create the kernels, cells were replicated in the positive X and Y directions every cell length (in this case 13.498 µm or 10.875 µm for the spherical and cubic cells respectively) using a Bash script. This created a layer of cells with the cell containing activity placed at the left most bottom quadrant (see Fig. 2). To reduce the simulation memory, each of the replicated cells only contained a nucleus and a cytoplasm; however 10 µm was added to the replicated cytoplasm to account for the missing membranes and maintain consistent cellular spacing. The original source cell had all four cell regions (nucleus, nuclear membrane, cytoplasm, and cellular membrane). The areas in between the cells were water density.

Fig. 2

The 225Ac cell distribution for the cubic cell kernel modelled in GATE. The cell with the green nucleus is the source cell. The 177Lu cubic kernel was identical, however was 56 × 56 cells instead of 10 × 10

For each radionuclide and geometry (spherical and cubic), two kernels were created (see Fig. 3): i) when activity was placed homogenously in the cytoplasm of the central cell to simulate radiopharmaceutical internalization, and ii) when activity was placed on the cell membrane to simulate a radiopharmaceutical that is cell surface bound. The GATE DoseActor was attached to each nucleus in the cellular grid and was 1D, meaning the “edep” output file ascribed a singular value corresponding to the total energy deposited in each cell nucleus (compared to the single cell simulations where the DoseActor was 3D, resulting in a 3D matrix as output). The pixel size was one cell size (see Table 1) to enable convolution with cellular activity maps. Despite this, because the DoseActor is attached to the nucleus it only describes the energy deposited in the nucleus. 1 × 108 and 1 × 106 primaries were simulated for 177Lu and 225Ac respectively.

Fig. 3

A 2D, 5 cell × 5 cell subset of the input for the simulation of the spherical absorbed dose kernels (not to scale) with one full cell geometry in the bottom left quadrant. All other cells contain only the cell nucleus and extranuclear region. The white gaps are water density and the black lines represent the pixel boundaries

Table 1 The final dimensions and lengths of all kernels post-processing

Using Python, the GATE DoseActor “edep” outputs for each cell nucleus were combined to create a 2D one-quarter kernel of the deposited energy with the source activity in the bottom left corner. This quarter kernel was transformed appropriately to create a full kernel, with one of the center rows and columns removed to obtain a kernel describing the deposited energy in each cell nucleus from activity in the source region of the center cell. The kernel was converted to units of S-value (absorbed dose per decay) as described in Sect. "S-value calculation". The final sizes of each kernel are presented in Table 1. While simulating only a one-quarter kernel and extrapolating to a full kernel omits the contribution from scattered radiation off of cells outside the quarter region, we believe this effect would be minimal.

While the maximum range of beta particles from 177Lu is around 2.0 mm in tissue, the uncertainty in S-value of the kernel at 2.0 mm from the center cell was very high, given the very small pixel sizes of these kernels. Because the X90 (the radius of a sphere in which the beta particles have released 90% of their energy) of 177Lu in water is approximately 600 µm [23], we did not encompass the full range of beta particles from 177Lu (from the kernel center to the end in each direction) as that was not computationally feasible.

To assess the impact of different ratios of internalized versus extracellular bound activity, we combined the kernels in three different ratios of cytoplasm to membrane activity: 3:1, 1:1, and 1:3. This was done using the following equation:

$$_=\frac_+ m\cdot _}$$

(1)

where c and m are the ratios of activity in the cytoplasm to membrane respectively, and Kcytoplasm and Kmembrane are the cytoplasm and membrane kernels respectively. The ratios of 3:1, 1:1, and 1:3 were chosen to represent the difference between radiopharmaceuticals with cytoplasm dominant uptake, membrane dominant uptake, or equal uptake. However, these kernels can be summed in any ratio to match the internalization pattern of a specific radiopharmaceutical that can be measured with internalization assays.

S-value calculation

S-values were calculated according to the MIRD schema, which defines the S-value as the mean absorbed dose to a target region per radioactive decay in a source region [24]. The “edep” matrices obtained from the GATE DoseActor were used to calculate all S-values. For the single cell simulations, the deposited energy in each voxel of the matrix was summed, converted from MeV to Joules, and divided by the nucleus mass and number of simulated primaries (also known as decays or histories) to obtain the total S-value in the cell nucleus from the source region in units of Gy/Bq/s as shown in Eq. 2. The same calculation was performed for each pixel of the extended nucleus absorbed dose kernels, however no summation was necessary as each pixel already represented the total energy deposited in each nucleus.

$$_= \frac^_}_\cdot _}$$

(2)

where Ei is the energy deposited in the ith voxel, n is the total number of voxels in the cell nucleus, Mnucleus is the mass of the nucleus, and Nprimaries is the number of primaries simulated with GATE.

2D tumour maps

2D multicellular tumour maps containing normal (healthy) cells, blood vessel cells, normoxic cancer cells, hypoxic cancer cells, and necrotic cancer cells were used to study the impact of hypoxia and necrosis on the absorbed dose distributions. These maps were created by Ahn et al., and a description of their creation and validation is described in the cited work [25]. Three distinct tumour morphologies were used in this study, composing mostly normoxic cancer cells, hypoxic cancer cells, or necrotic cancer cells. Two sizes of each morphology were studied, corresponding to approximately 5.5 mm and 10 mm in diameter.

To convert the tumour grids to activity maps, we assumed that the activity uptake within each cell depended on its oxygenation level, as the amount of oxygen was assumed to correlate to blood vessel availability and the ability for the radiopharmaceutical to reach the cell. According to Ahn et al., cells in the tumour model became hypoxic or necrotic when they had oxygen levels of 0.08 to 0.5% and < 0.08% respectively [25] and we assumed a value of 6% O2 for normoxic cells based on typical tissue oxygen levels [18, 26]. Based on these numbers, we assumed 91.2%, 7.6%, or 1.2% of the activity of each radionuclide went to each normoxic, hypoxic, and necrotic cell respectively. This was implemented by assigning 91.2 Bq, 7.6 Bq, or 1.2 Bq to each pixel of the tumour map depending on the cancer cell in that pixel, and then normalizing each pixel by the tumour’s total activity (obtained by summing the activity in each pixel) to ensure a total activity of 1 Bq per tumour. This enabled absorbed dose rate maps to be normalized in units of Gy/Bq/s. No activity was assumed to accumulate in the normal cells or blood vessel cells, under the assumption that there would be no uptake in off-target cells.

The summed 3:1, 1:1, and 1:3 cytoplasm to membrane ratio nucleus S-value kernels were convolved with the tumour activity maps to produce a 2D image of the absorbed dose rate per unit activity (DR/A) within the tumours in units of Gy/Bq/s. To perform the convolution, we used the ndimage.convolve function available in the Python SciPy package [27].

Data analysis

For each tumour morphology and size, the mean DR/A to each cell type (normal cells, blood vessel cells, and the three cancer cell types) was calculated. Dose-volume histograms (DVHs) for each cancer cell type and tumour morphology were created as a measure of absorbed dose heterogeneity within the tumours.

To assess the differences between the three different internalization ratios and the two different cell geometries used in the kernels on the tumour DR/A maps, a statistical analysis of variance (ANOVA) test was done on the mean DR/A to each cell type for each morphology using the scipy.stats function by SciPy in Python [27]. Two sets of ANOVA tests were done: the first compared the kernel geometries (the mean DR/A across all cells after the tumours were convolved with the kernels of different geometries, while the internalization ratio was kept the same) and the second compared the internalization ratios (the mean DR/A across all cells after the tumours were convolved with the kernels of different internalization ratios, while the geometry was kept the same). These tests were evaluated for both 177Lu and 225Ac.