## Abstract

The biologic response of tissue exposed to radiation emitted by internal radioactivity is often correlated with the mean absorbed dose to a tissue element. However, studies show that even when the macroscopic absorbed dose to the tissue element is constant, the response of the cell population within the tissue element can vary significantly, depending on the distribution of radioactivity at the cellular and multicellular levels. These variations are believed to be the consequence of nonuniform distributions of activity among the cells or subcellular compartments that comprise the tissue element. Furthermore, the self-dose received by a cell containing radioactivity can be more radiotoxic than the cross-dose from neighboring cells. To study how the nonuniformity of activity among cells can affect the dose response, a 3-dimensional model of cells in a heterogeneous carbon scaffold was used to assess response. **Methods:** A theoretic model of a 3-dimensional cell culture was constructed, and Monte Carlo radiation transport was performed to assess self- and cross-doses for each cell nucleus in a population of 10^{6} cells. On the basis of these individual doses and on empiric models of radiation-induced cell death (i.e., reproductive failure), survival curves were simulated with different electron energies and activity distributions among the cells. **Results:** Nonuniformity of cell activities are responsible for nonuniformity of the dose at the cellular level, which in turn causes a change in the surviving fraction of the cell population from that expected on the basis of uniform activity and dose. **Conclusion:** The macroscopic dose received by a tissue cannot be used to anticipate its biologic response. The dose distribution among individual cells, because of both their nonuniform activity and geometric environment, is an important factor in determining biologic response of the tissue at the macroscopic level.

There are many variables that dictate the overall biologic response of tissues that contain radioactivity. Among the variables are tissue radiosensitivity; distribution of radioactivity at the macroscopic (1,2) and multicellular (3–7), cellular, and subcellular levels (8–10); relative biological effectiveness of radiations emitted (e.g., α, β, Auger electrons) (11); pharmacokinetics; dose rate (12); repair; oxygen tension (13); and bystander effects (14–16). Considerable progress has been made in correlating biologic response with many of these variables (17). However, we still have only a limited understanding of the correlation of biologic effects with nonuniform dose distributions that result from nonuniform distribution of radioactivity at the cellular and subcellular levels (7,18–23). Accordingly, when the distribution of activity (or sources) is uniform at the macroscopic level but nonuniform at the multicellular level, the mean absorbed dose to a tissue element may not be a suitable quantity for use in predicting biologic effects (20).

In previous work, it has been pointed out that improvement in our capacity to predict the biologic effects of incorporated radionuclides will require not only sophisticated theoretic models but also robust experimental models that can generate data to validate the theoretic approach (17). Considerable effort has been devoted to developing paired theoretic and experimental multicellular dose–response models, and these efforts have led to improvements in our understanding and modeling capabilities (7,18,24). Reports by Humm and Chin (18) and Howell and Neti (24) describe 3-dimensional (3D) experimental multicellular cluster models that afford precise control over the cellular dosimetry variables (7,24). Recently, our laboratory has devoted considerable efforts toward developing a new 3D tissue culture model that enables us to more closely simulate human tissue in vivo. This model uses a Cytomatrix carbon scaffold (Cell Sciences PTE), and it has been used to study radiation-induced adaptive responses and bystander effects elicited by nonuniform distributions of radioactivity and concomitant quantification of targeted drug delivery and biologic response in individual cells (25–27).

The purpose of our study is to develop a computer model that represents this new 3D tissue culture model. Coupling a Monte Carlo transport code to this computer model should further our understanding of the role played by each independent variable on the biologic response of the tissue by recording both the cellular self-dose (from radiation emitted by radioactivity within the cell) and the cross-dose (from radiation emitted by radioactivity in other cells). To model the complex geometry of the Cytomatrix scaffold, a representation of the scaffold with simple geometric shapes is not feasible. Instead, a 3D CT image is used to provide a faithful representation of the scaffold. Similar developments were previously implemented for bone marrow dosimetry (28,29). In these previous implementations, a sample of human trabecular bone was harvested from a cadaver and imaged at high resolution with either nuclear magnetic resonance or micro-CT. Image processing was then applied to produce a smooth 3D model of the trabecular bone microstructure. This model was then coupled with a Monte Carlo radiation transport code (30,31). In more advanced models, several image datasets were combined to provide different levels of complexity for the geometry (32). The present work, for the first time to our knowledge, marries these bone dosimetry methods to a multicellular approach to modeling the surviving fraction of cells within a Cytomatrix scaffold after exposure to nonuniform distributions of radioactivity. The difference in dose response between the self- and cross-doses is addressed (24).

## MATERIALS AND METHODS

### Overview

In brief, a 3D multicellular dosimetry model was created comprising a carbon scaffold, cell nuclei randomly distributed in the pores of the scaffold, and medium that fills the remaining space (representing a combination of the cytoplasm, cell membrane, and extracellular matrix). Cellular radioactivity was restricted to the nucleus. The relative activity in each cell was assigned according to a lognormal probability distribution. The Monte Carlo transport code EGSnrc (33) was then used to assess both the self-dose and the cross-dose for each individual cell nucleus. These absorbed doses provided the input data that were needed to calculate the probability of survival for each individual cell and the overall surviving fraction of the population of cells within the Cytomatrix scaffold.

The present study focused on monoenergetic electrons emitted upon decay of a hypothetical radionuclide. It has been shown that inclusion of photons for a radionuclide that emits both electrons and photons does not alter significantly the results, provided that the volume of the cell population is small enough that the photon-absorbed fraction is negligible in comparison to the electron-absorbed fraction (34). Furthermore, the initial use of monoenergetic electrons facilitates benchmarking of the computer model. Our objective was to further our understanding of the interrelationship between electron energy and activity distribution in terms of their effect on the surviving fraction of cells in a 3D multicellular model of heterogeneous composition.

### Micro-CT Model of Carbon Scaffold

The Cytomatrix scaffold was imaged under dry conditions, to provide adequate contrast between ligaments and pores, using a desktop cone-beam micro-CT scanner (Scanco Medical AG). Cubic voxels with edge lengths of 6 μm were characterized by spatial coordinates and a pixel intensity proportional to its average radiographic absorption coefficient. To reduce noise, the image dataset was median-filtered using a kernel size of 3 × 3 × 3 voxels. After filtration, the threshold intensity to distinguish carbon from air was determined using a gradient-magnitude technique (30). Such binary segmentation provides a faithful representation of volume fractions between the 2 media. However, because of the jagged representation of the boundary between voxels, this segmentation does not preserve the surface area of the interface between the 2 media. This voxel effect has consequences for dosimetric assessment using Monte Carlo transport and is especially important for low-energy electrons (35). At a 6-μm voxel size, the effect is pronounced for electron energies below 20 keV (e.g., for electrons emitted by ^{3}H or ^{125}I). To overcome this problem, a smoothing technique suitable for Monte Carlo transport code was used (31). This technique facilitates representation of the boundary between carbon and air with a smooth and continuous surface that eliminates most of the voxel effects. Figure 1 is a rendering of a micro-CT image of the Cytomatrix disk (diameter, 8.5 mm; thickness, 2 mm). Voxel analysis yields an average pore diameter of approximately 0.7 mm, ligament diameter of approximately 50–100 μm, and ligaments occupying 4.1% of the total Cytomatrix scaffold. These parameters depend on the manufacturer's average pore size, available in 10–100 pores per inch.

### Construction of Cell Field within Cytomatrix Scaffold

For modeling the biologic environment, the air in the pores of the scaffold image was computationally replaced with cells. Cell nuclei were represented by spheres (diameter, 8 μm). The assumption of spheric geometry yields self-doses within 15% of those calculated with extreme elliptic geometry over the range of electron energies considered here (8). The cytoplasm was not explicitly modeled; rather, nuclei were randomly positioned within the pores but were not allowed within 1 μm of the ligaments, and nuclei were separated by at least 2 μm. With these restrictions, 1.5 × 10^{8} nuclei fit in the scaffold, with a mean distance of 2.48 ± 0.45 μm between neighboring nuclei and a packing density of 16%. A 31% cell packing density and mean distance of 0.48 ± 0.45 μm between neighboring cells is obtained if a 1-μm-thick cytoplasmic shell surrounds each nucleus. Because our personal computer cluster cannot support computations with 10^{8} cells, this study was limited to 10^{6} cells in the central portion of the Cytomatrix device and specifically a cylinder with a diameter of 1.322 mm and height of 1.252 mm.

### Distribution of Activity

The distribution of activity among the cells was assumed to be uniform or lognormal. For the latter, each of the 10^{6} cells was assigned a relative activity using a lognormal probability density function (PDF) normalized to unity (22). Lognormal shape parameters (σ = 0.6, 1.0, and 2.0) were considered. The PDF and cumulative probability function are shown in Figure 2. The relative probability of a decay occurring in a given cell is also given by the PDF. Thus, during a simulation, cell nuclei were randomly selected for a decay based on probabilities obtained from the lognormal PDF. Source particles were generated arbitrarily in random nuclei, with activity uniformly distributed in individual nuclei.

### Monte Carlo Radiation Transport

Radiation-absorbed doses to cell nuclei within the Cytomatrix scaffold were assessed using the EGSnrc Monte Carlo radiation-transport code (33). The transport geometry was limited to a rectangular solid region (1.380 × 1.380 × 1.308 mm) that contained the cylinder defining the field of nuclei. Every time a particle was transported outside the region, it was discarded and not followed further. This reduced the central processing unit (CPU) load for the simulations. Monoenergetic electrons (10, 30, 100, 300, and 1,000 keV) were used as source radiations. Photons produced during the electron transport were transported by EGSnrc. A history is defined as the process of radiation transport of a single particle along its ionization track. For each history, energy deposition events in each cell nucleus were recorded and stored according to self- and cross-deposition of energy in each nucleus.

The compositions of the components of the model were assumed to be carbon graphite with a 2.0 g/cm^{3} density for the ligaments, cell nuclei as defined by the International Commission on Radiation Units and Measurements (36), and liquid water for the remaining volume. The cytoplasm was not considered a source or target region to reduce CPU time. Other transport parameters were Rayleigh scattering (important at low photon energy), isotropic electron emission, and 1-keV cutoff energy for both photons and electrons.

### Cell Killing and Calculation of Surviving Fraction

The surviving fraction of radiolabeled cells (SF_{labeled}) within a multicellular cluster can be described by Equation 1,_{self} and D_{cross} are the self- and cross-doses to the cell nucleus (24). Similarly, for a given radionuclide, D_{37,self} and D_{37,cross} are the self- and cross-doses required to achieve 37% survival, respectively. The need to account for differing radiotoxicities for self- and cross-doses arises because energy deposited in a cell by a radionuclide incorporated into its DNA has been shown to be more radiotoxic than energy deposited by emissions arising from decays in neighboring cells (24).

Equation 1 can be used to calculate the probability that an individual cell survives a given radiation insult. For the current study, the parameters D_{37,cross} and D_{37,self} were assigned values of 4.0 and 1.2 Gy, respectively. These values were experimentally determined for ^{131}I-iododeoxyuridine incorporated into the DNA of V79 cells maintained in 3D clusters (21). The surviving fraction of the cell population was determined as follows. For each cell, a survival probability was calculated by substituting its self- and cross-dose into Equation 1. A random number 0 and 1 was computer-generated and compared with the survival probability. If the random number was smaller than the generated probability, the cell was scored as a survivor. Otherwise, it was scored as reproductively dead. The fraction of survivors among the 10^{6} cells that compose a given simulation represents the surviving fraction of the cell population.

### Theoretic Survival Curves

The calculation of the surviving fraction provides only the fraction for a specific simulation with a specific number of electron histories. To build a complete survival curve, it was necessary to execute many simulations, each with an increasing number of histories. To achieve several logs of kill, simulations consisting of up to 5 × 10^{10} histories were required to deliver 100 Gy. However, because of the large number of cells and complex geometric configuration of the nuclei within the scaffold, a single history can be so CPU-intensive that the computer cluster cannot run enough histories to deliver 100 Gy within a reasonable time frame. Furthermore, even when millions of histories were run, it was often observed that an increase in the number of histories resulted in an appreciable change in the mean absorbed dose to a given cell nucleus per decay in the scaffold. However, such large fluctuations did not necessarily affect the surviving fraction because of the large number of cells involved. Accordingly, an empiric technique was used to determine the optimum number of histories. Preliminary simulations with 10^{4}, 3 × 10^{4}, 10^{5}, 3 × 10^{5}, and 10^{6} histories were first performed. The self- and cross-energy depositions in each cell nucleus were scored for each simulation and used to calculate the self- and cross-doses to each cell nucleus, and surviving fractions were determined for each simulation. These 5 surviving fractions were plotted as a function of their corresponding mean absorbed doses to the population of nuclei. Additional points on the graph were obtained by multiplying the self- and cross-doses (at 10^{6} histories) for each nucleus by a scaling factor and resimulating the surviving fraction. Thus, a complete survival curve was constructed up to 100 Gy. A new survival curve was then created by adding a sixth point with 3 × 10^{6} histories and then adding points using factors to scale up from doses based on 3 × 10^{6} histories. This process was repeated several times with an increasing maximum number of histories (up to 1,000 × 10^{6}) until the new curve was not appreciably different from the preceding curve.

## RESULTS

### Number of Track Histories Required for Simulations

Twenty survival curves were built with an increasing number of histories until the shape of the curve did not change appreciably. The convergence process for 3 conditions—uniform distribution of a hypothetical radionuclide that emits electrons with initial energies of 10, 100, and 1,000 keV—is shown in Figure 3. On the basis of Figure 3A, 1,000 × 10^{6} histories are required for convergence of the survival curve for 10-keV electrons. About 300 × 10^{6} and 100 × 10^{6} histories are needed for the 100- and 1,000-keV electrons (Figs. 3B and 3C), respectively.

Table 1 provides the number of histories required to achieve convergence for all initial electron energies considered for a uniform activity distribution. The corresponding mean doses are also given. Because only a small dose is required to achieve convergence, only the low-dose region of the converged curves in Figure 3 are built with simulations that have a sufficient number of histories to produce the actual dose. In contrast, the high-dose portion of the curves (>1–5 Gy depending on the initial particle energy) are obtained by the extrapolated manner detailed in the previous section. It was also determined that the distribution did not substantially affect the number of histories required to reach convergence.

### Dependence of Survival Curves on Energy and Activity Distribution

The dependence of the shape of the survival curves on electron energy is shown in Figures 4A–4D for uniform and lognormal distributions of activity among the entire cell population. Similarly, the dependence of curve shape on activity distribution is shown in Figures 5A–5D for initial electron energies of 10, 30, 100, and 1,000 keV, respectively. Exponential survival curves corresponding to D_{37} = D_{37,self} = 1.2 Gy and D_{37} = D_{37,cross} = 4.0 Gy are also shown. They represent survival curve limits that result when every cell receives the same dose in the form of only self-dose or only cross-dose, respectively. They will be used as references to interpret further results. In keeping with our experimental observations for ^{131}I with clonogenic survival as the biologic endpoint (21), these reference curves have no shoulder that is typical of low–linear energy transfer radiations. More complex dose–response models may be needed for other radiations, endpoints, and tissue environments in which repair processes may be more robust (37).

## DISCUSSION

The model presented herein is a Monte Carlo approach that uses multicellular dosimetry to obtain radiation-absorbed dose and biologic response of individual cells located in a geometrically complex and heterogeneous environment. The complex array of variables that were studied produced a series of survival curves with a number of implications.

### Track Histories

Figure 3 shows that the number of electron track histories used to create a survival curve can have a sizeable impact on curve shape. Simulations of 10^{6} 10-keV electron histories result in a survival curve that saturates at 37% regardless of the mean absorbed dose to the nucleus. This result arises as a consequence of statistical fluctuations that are inherent in the simulation. The range of a 10-keV electron in liquid water is about 2.5 μm, and few of the electrons emitted from within a cell nucleus can escape to deliver a cross-dose to neighboring cells. Therefore, cells containing radioactivity that emits monoenergetic 10-keV electrons receive only self-dose. Given that there are only 10^{6} cells in our population, and only 10^{6} histories are simulated, each cell should emit an average of 1 electron. When the activity is uniformly distributed among the cells, selection of the source cell for a given decay (and corresponding history) follows a binomial distribution, with each cell having a selection probability of 10^{−6}. In this scenario consisting of few histories, 36.8% of our cells are expected to emit no electron at all and, as a consequence, receive no self-dose. When the doses to the cell nuclei based on 10^{6} histories are extrapolated (multiply each cell nucleus dose by a dose-scaling factor) to build the high-dose portion of the curve, the same cells remain unirradiated and survive regardless of the mean dose to cell nucleus. With 3 × 10^{6} trials, the binomial distribution gives 4.98% of the cells with no source, which agrees with our second curve in Figure 3A. These percentages give confidence in our approach to determining the optimal number of histories because in the limit of 10-keV electrons, the model behaves exactly as one would expect from simple probabilistic arguments.

Interpretation of the convergence of the survival curves for the 100- and 1,000-keV electrons (Figs. 3B and 3C) is more difficult because of the added contribution of cross-dose. However, even in these cases, at low history numbers, there are cells that do not receive a dose because their nucleus is never traversed by a particle. The corresponding survival curves will saturate regardless of the magnitude of the mean dose. Thus, determination of the optimal number of histories is essential.

### Dependence of Survival Curves on Energy

There are 2 interesting findings in Figure 4A. First, as the energy increases, the curves shift from an exponential response with D_{37} = D_{37,self} = 1.2 Gy toward the exponential response with D_{37} = D_{37,cross} = 4.0 Gy. This shift occurs because when 10-keV electrons are emitted, the cells do not receive cross-dose and the simulated curve exactly follows the self-dose exponential curve with the expected D_{37} = 1.2 Gy. At 30 keV, with an electron range of 17.5 μm in soft tissue, the cells start to receive a significant cross-dose and the curve starts to shift toward that for response due only to cross-dose. At 1,000 keV, the self-dose is insignificant, compared with the cross-dose, and the survival curve behaves nearly as the cross-dose exponential curve with D_{37} = 4.0 Gy.

The second interesting point in Figure 4A relates to the shape of the curves. At low energy they follow the straight line of an exponential model with a D_{37} corresponding to D_{37,self}. Because of the short range of the electrons, the geometric position of the cells has little to do with the amount of dose they receive. As the energy increases slightly, the local surrounding cells begin to contribute a cross-dose and all cells continue to receive a dose close to the mean dose. Therefore, calculation of the surviving fraction with Equation 1 yields an exponential curve. As the energy is increased further, the cross-dose becomes ever more important. However, the higher electron stopping power of the carbon ligaments relative to water provides a shielding effect for those cells that are close to ligaments. Furthermore, cells near ligaments have fewer neighbors because of geometric constraints. On average, these cells receive about half the cross-dose that the more predominant and centrally located cells receive. That is why at 100 keV and 300 keV, the surviving fraction is higher than if all cells received the same dose as the mean dose. This effect is apparent in Figure 4A because the surviving fraction starts to diverge from the straight line at high dose. When the electron energy increases to 1,000 keV, the electrons easily traverse the ligaments and no longer provide significant shielding as they do at intermediate energy. Therefore, at 1,000 keV, the dose distribution among the cells is again more uniform—except for those cells that are on the edge of the simulation model. Thus, at high energy, the survival curve almost reverts to a straight cross-dose exponential curve (D_{37} = 4.0 Gy).

The survival curves for 1,000-keV electrons show a small degree of upward curvature (Fig. 4), likely due to an edge effect wherein the cells at the edge of the modeled cylinder receive, on average, a lower dose than the cells at the center of the cylinder. In fact, for a uniform distribution of activity in a sphere, the dose rate at the edge of a sphere is half the dose rate at the center (2). Edge effects also include the 2-fold-lower doses received by cells adjacent to ligaments.

### Dependence of Survival Curve Shape on Activity Distribution

Many of the energy-dependent features of the survival curves found in Figure 4A are also observed in Figures 4B–4D; however, interpretation of these curves first requires an understanding of their dependence on the activity distribution. Figure 5A shows that, for 10-keV electrons, the shape of the survival curve is highly dependent on activity distribution. The range of 10-keV electrons is shorter than the dimensions of the cells, and therefore each cell receives only self-dose (no cross-dose). With a uniform distribution of activity, each cell receives the same dose, and the surviving fraction depends solely on the mean dose. Thus, the survival curve for the short-range 10-keV electrons follows an exponential response commensurate with D_{37,self}. With a lognormal distribution, however, a significant fraction of the cells have low activity, compared with the mean (Fig. 2B). For example, when σ equals 2.0, 74% of the cells contain less than the mean activity and, as a consequence of being irradiated only by themselves, also receive correspondingly less dose than the mean. Thus, they are more likely to survive than their neighbors with higher activity. The cell nuclei that receive a larger dose than the mean are more likely to die but they represent only 26% of the population and thus do not compensate for the excess survival. Even when σ equals 0.6, only 40% of the cells receive a higher dose than the mean (Fig. 2B), and they do not balance the 60% that are more likely to survive. As depicted in Figure 5A, this effect is substantial.

Though all curves for 10-keV electrons follow the self-dose exponential curve at low dose (Fig. 5A), they increasingly deviate from this curve as the activity distribution becomes increasingly nonuniform (i.e., as σ increases from 0.6 to 1 to 2). However, as the electron energy is increased, the effect of the nonuniformity is minimized because the cross-dose becomes the dominant factor in dictating the shape of the response curve (Figs. 5B–5D). At intermediate electron energies (30–100 keV), both the self- and the cross-doses play important roles in the shape of the curves. However, at high energies (1,000 keV), the cross-dose is fairly uniform and, because the self-dose becomes insignificant, the curves do not depend on the activity distribution. In fact, all 4 curves in Figure 5D (1,000 keV) follow the cross-dose exponential curve (D_{37,cross}), with the exception of a small deviation caused by the reduced cross-dose along the edges of the simulated geometry.

This understanding of the dependence of the shape of the survival curve on activity distribution facilitates explanation of Figures 4B–4D. When the activity distribution is uniform (Fig. 4A), the surviving fraction increases with increasing electron energy because of the corresponding increased contribution of the less radiotoxic cross-dose (D_{37,cross} = 4 Gy) relative to the more radiotoxic self-dose (D_{37,self} = 1.2 Gy). A small increase in nonuniformity causes a substantial change in the survival curves for low-energy electrons (10 and 30 keV) but little change in the shape of the survival curves corresponding to high-energy electrons (100, 300, and 1,000 keV). In Figure 4C, which represents results for a lognormal distribution when σ equals 1.0, the same pattern emerges, but the changes in curve shapes for low-energy electrons become even more significant. Finally, dramatic changes in the response to low-energy electrons are manifested when the shape parameter is increased to σ equals 2.0 (Fig. 4D). In this case, notable changes are also observed for 100-keV electrons. The responses to 300 and 1,000 keV remain unaffected by the activity distribution because of the long range of these electrons in water and their high penetrability through the ligaments of the Cytomatrix scaffold.

## CONCLUSION

Monte Carlo simulations within a realistic model of 3D tissue culture show that the character of the biologic response (cell survival) of a cell population, labeled with a monoenergetic electron emitter, depends on the initial energy of the electron and the nonuniformity of the activity distribution among the cells. Both the initial electron energy and the activity distribution, as characterized by the lognormal shape parameter σ, determine the relative contributions of the self- and cross-doses to the absorbed dose to the cell nucleus. Earlier studies have shown that, when the radionuclide is incorporated into the DNA in the cell nucleus, self-dose can be more lethal than cross-dose. This disparity is largest for Auger electron emitters. The interplay between these parameters can profoundly affect the shape of cell survival curves (Figs. 4 and 5). Changes in the shape of the survival curve are most significant at low electron energy, for which self-dose predominates. At high electron energy, cross-dose dominates and mitigates the dependence of the curve shape on energy. Therefore, the mean absorbed dose to the cell population is often not a reliable parameter for characterizing the biologic response of a tissue that contains radioactivity—an important conclusion considering that this Monte Carlo model can be readily extended to more complex dosimetry applications such as human trabecular bone. Last, it is also recognized that a range of cell sizes and packing densities are possible and that these can also significantly affect the biologic response.

## Acknowledgments

We thank Scanco (Basserdorf, Switzerland) for acquiring the micro-CT images. Microscope images were provided by John Akudugu. This study was supported in part by NIH grants R01CA083838-09 and R01CA116743. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Cancer Institute or the National Institutes of Health.

- © 2011 by Society of Nuclear Medicine

## REFERENCES

- Received for publication June 7, 2010.
- Accepted for publication December 29, 2010.