The linear-quadratic (LQ) parameterization of survival fraction [SF(D)] inherently assumes that all cells in a population receive the same dose (D), albeit the distribution of specific energy z over the individual cells f(z,D) can be very wide. From these microdosimetric distributions, which are target size dependent, we estimate the size of the cellular sensitive volume by analyzing its influence on the LQ parameterization of cell survival. A Monte Carlo track structure code was used to simulate detailed tracks from a 60Co source as well as proton and carbon ions of various energies. From these tracks, f(z,D) distributions were calculated for spherical targets with diameters ranging from 10 nm to 12 μm. A cell survival function based on f(z,D) was fitted to experimental LQ α values, revealing an intrinsic limitation that target size imposes on the usage of f(z,D) to describe the linear term of the LQ parameterization. The results indicate that such threshold volume arises naturally from the relationship between the particle's probability of no-hit and the probability of cell survival. Further analysis led to the proposal of a radiobiological property , defined as the mean lineal energy corresponding to the target size that allows equivalence between the mean inactivation dose (MID) and the mean specific energy . The fact that is an increasing continuous function of target size within the range of biological targets of interest in radiobiology, ensures the uniqueness of for any radiation quality, thus, its potential usefulness in modeling. In conclusion, an accurate estimation of such threshold volumes may be useful for improving modeling of cell survival curves.
The clonogenic cell survival assay is a straightforward and inexpensive in vitro assay that quantifies the proportion of cells that retain reproductive integrity, and thus proliferate to create viable colonies, after irradiation with incremental amounts of dose. The resulting cell survival ratio has been and still is fundamental to the study of biological effects from ionizing radiation. The method most used for describing cell survival data is the so-called linear-quadratic (LQ) parameterization (1), where the probability for survival at a dose (D) is given by:
In addition to the trivial Poisson interpretation, that the parameters α and β are proportionality factors for lethal events, a satisfactory bio-mechanistic explanation is still missing. Equation (1) provides a quantification of the integrated response of the cell where many complex intra- and extracellular processes are involved. Indeed, to decouple each process is a colossal task. However, with the advent of particle radiotherapy, there has been a renewed interest in the development of semi-mechanistic mathematical models for prediction of the LQ parameters. Implementation of these models in radiation treatment planning systems would allow optimization of biological effectiveness throughout the treated volume. This is essential, particularly for treatment with heavier ions (e.g., carbon), since the large and rapid change in the relative biological effectiveness (RBE) at the end of the Bragg Peak may pose a serious threat to any healthy tissue located behind (2, 3). In proton therapy, it is common practice to apply a generic constant RBE value of 1.1 (4), however, there are also strong indications of the benefits of implementing variable RBE in clinical planning (5–7).
Micro- and nanodosimetric-based semi-mechanistic models assume the existence of a so-called “sensitive volume,” which is generally conceptualized as the volume within a cell that will most likely induce cell death if affected by the ionizing radiation (8). It is natural to assume that such volume is found within the cell nucleus, as it harbors the DNA molecule, yet there is increasing evidence suggesting induced cell death by compromising other cellular structures (9). The underlying truth is that without knowledge of the subcellular composition of the sensitive volume, it remains unmeasurable with today's biotechniques, and therefore an issue open for debate. For these types of models, knowledge of the size of such a target is crucial, since it represents a way towards characterizing radiation according to either the distribution of energy deposition, or the spatial distribution of the points of interaction (e.g., ionization clusters) within such targets.
In this work we explore the influence of target size variation on microdosimetric energy deposition distributions and the effect it may have on cell survival LQ parameterization. The dose-dependent frequency distribution of specific energy f(z,D) provides a full microscopic description of the distribution of energy deposition in a population of targets as a function of the average dose (10). This makes f(z,D) a candidate for use in the estimation of biological response at a subcellular level.
MATERIALS AND METHODS
Track Structure Simulations
A set of 500 proton and carbon-ion tracks (where a track consists of the primary particle together with its entire secondary particle cascade) was simulated in a water phantom with the updated version (11) of the Monte Carlo (MC) track structure code LIonTrack (12). To cover the thickness variation among healthy and cancerous mammalian cell nuclei, each particle track, starting from a unidirectional point source, was scored along a 12-μm-long segment without lateral restriction. Longitudinal equilibrium of the released secondary electrons along this segment was provided by starting the particle at the top of a padding layer upstream and stopping the particle at the end of a second padding layer downstream of the 12-μm scoring layer (see Fig. 1A). At the center of the scoring layer the nominal energies were 0.91, 1.40, 1.72, 3.18, 3.59 and 4.97 MeVu–1 for the protons and 5.27, 10.95 and 76.9 MeVu–1 for the carbon ions. Such energies correspond to those at which the referenced cell survival experiments (see Table 1) were performed. The transport cut-off energy of all secondary electrons was set at 50 eV, whereas the electron production cut-off was 1 eV, which is the lowest energy in the tabulated cross sections used by LIonTrack. Thus, any electron produced with or reaching a kinetic energy below the transport cut-off was not transported further and its energy was deposited on the spot. The energy deposition (ED) profile consisting of the coordinates where the EDs (ionizations and/or excitations) took place was tallied in separate files, one per simulated ion track. For the 60Co source, a monoenergetic 1.25 MeV photon beam was used for the release of electron tracks in liquid water using an isotropic point source. No scattered photons were considered because in in vitro irradiation conditions their production is minimal, thus making their contribution to biological damage negligible.
Published LQ Alpha (α) and Beta (β) Parameters from Experimentally Obtained Survival Curves for V79 Cell Line
The distribution of a sum of uncorrelated stochastic variables equals the convolution of their distributions. Therefore, the accumulation of dose in a microscopic subregion affected by a certain number of tracks can be derived by repeatedly convolving the frequency distribution of specific energy for one track f1(z) with itself until the number of tracks v is reached. The mean of f1(z) represents the mean specific energy that one track would deposit in the target volume and therefore the mean number of tracks n that would yield a dose is estimated by . Assuming that the number of tracks is Poisson distributed, the dose-dependent frequency distribution of specific energy is given by Kellerer and Chmelevsky (10) as:Equation (2) has two distinct components: 1. The first is a Poisson distribution e–nnv/v!, giving the probability that v tracks will contribute to the total dose D when the mean number of tracks is equal to n; 2. The frequency distribution of specific energy for v tracks , which gives the probability for a certain amount of energy to be imparted by exactly v tracks in a specified target. By definition, the mean of f(z,D) is equal to D.
Importantly, f1(z) and all quantities derived from it are dependent on target size.
Calculation of f1(z) and f(z,D)
The f1(z) distributions for protons and carbon ions were calculated by dividing the scoring layer described above, into layers of thickness equal to the diameter of the target volume. Each layer was filled with a grid of spheres with matching diameters while making sure that the track's core would traverse the center of the “middle” sphere at each layer, as shown in Fig. 1A. As a mean of variance reduction technique, the generated tracks were translated laterally to random positions across the projected area of the middle spheres (as shown in Fig. 1B) to reproduce the correct chord length distribution of tracks traversing a convex target at different positions. Each of the 500 simulated tracks was translated to 100 random positions, achieving a statistical uncertainty less than 1.5%. According to the reciprocity theorem [e.g., Attix (13)], this setup allows scoring of energy deposition not only by the track's core in the middle spheres but also from secondary electrons that traverse any of the surrounding spheres. Repeating the scoring procedure after each translation of the track yields a histogram which, after conversion of the deposited energy to specific energy, results in the f1(z) distribution. This process was repeated for spheres with diameters ranging from 0.01 to 12 μm.
A different approach was needed for the 60Co source because of the tortuous trajectories of the produced electrons. We applied the same procedure as used by Villegas et al. (14) for calculation of f1(z) distributions where the volume of a bounding box enclosing each track (starting from a point source immersed in water) was divided into a grid of cubic sub-volumes of the desired size and for which the energy deposition was scored separately. However, two changes were applied. The first consisted of exchanging the scoring volumes from cubes to spheres. The second was to randomize the starting position of the first electron interaction within the spherical sub-volume in which it was found, thus avoiding bias in energy deposition within such sub-volume. Translation of the entire track was done accordingly. A total of 20,000 electron tracks were used to reach an acceptable statistical uncertainty.
For each of the spherical volumes, f(z,D) was determined according to Eq. (2) for doses between 1–6 Gy at 1 Gy intervals, thus covering the range of doses at which the published experimental survival assays were performed and whose LQ parameters will be used for later analysis (see Table 1). The convolutions of f1(z) to yield fv(z) were implemented using Mathematica™ version 11.1 (Wolfram Research Inc., Champaign, IL) cf. Villegas et al. (14). For some energies, the Poisson probability for v = 0 had the dominating weight factor, and to maintain normalization of f(z,D), Eq. (2) was implemented with f0(z) set to a Dirac delta distribution at z = 0.
Cell Survival Fraction and Microdosimetry
The clonogenic cell survival assay provides a simple method for quantifying survival fractions, SF(D). Essentially, the fraction of cells that manages to proliferate from independently exposed populations to various dose levels is quantified. Recalling from the microdosimetry theory that the energy deposited in each individual cell is governed by f(z,D), where D is the dose given to the cell population, and assuming the existence of a survival fraction SF(z), where cells receiving equal z have equal probability of mortality, we can hypothesize that the SF(D) can be described by:Eq. (3) is indeed the most basic in terms of modeling, namely that biological response is governed by the specific energy received and that the biological effect exhibited by the irradiated cells is independent of the effect on the others (i.e., no bystander effect). Yet, for the purpose of this work, it is ideal to express the survival fraction as in Eq. (3) because it opens the possibility to directly study the influence of target size on biological response, as f(z,D) hinges on target size. Following Poisson statistics, SF(z) can be described by e–az, where parameter a represents the mean number of lethal events per specific energy, yielding:
To study the behavior of the parameter a as a function of target size, we make use of experimentally derived SF(D) approximated by LQ parameterization. For the following analysis we use only the linear term, as it describes SF(D) at fraction doses that are clinically relevant in ion therapy. Values for the parameter a can thus be fitted by minimizing the following objective function:
In practice, the f(z,D) frequencies present intrinsic noise, which makes brute-force integration challenging. To facilitate computations, we instead make use of the cumulative dose-dependent specific energy function , which enables application of the integration by parts yielding:Eq. (6) includes full integration the Dirac pulse at z = 0 in f(z,D). Mathematically, the lower limit should be marked as 0–, however, since negative values of z have no physical significance, we kept the notation with 0 as the limit. The right-hand side integral of Eq. (6) can be further divided into two separate ranges; the first is [0,zmax], where zmax is the maximum tabulated value of F(z,D) and the second being [zmax,∞] for which the evaluation of F(z,D) is by definition equal to unity. Therefore, Eq. (6) becomes: Eq. (7) into Eq. (5) yields the final objective function
A data set of 10 in vitro experimental survival curves for the Chinese hamster lung cell line V79 was selected from the literature. Each survival curve corresponds to exposure to different radiation types varying in energy: 6 protons, 3 carbon ions and 60Co photons. Both of the LQ parameters are shown in Table 1.
The Distributions f1(z) and f(z, D)
Examples of the calculated f1(z) distributions are shown in Fig. 2. At a diameter of 0.01 μm the f1(z) rapidly decreases with z in more or less the same pace for all radiation species and energies; however, for the lowest carbon ion (5.27 MeVu–1) a distinct peak is observed at approximately 3 × 105 Gy. This peak becomes distinctively triangular as the target volume increases. Likewise, f1(z) for all ions (protons and carbons) presents such a characteristic triangular shape at some given target size. This behavior is a hallmark of the chord length distributions resulting from straight track segments traversing a sphere. Thus, the triangular part of f1(z) is due to the track cores, whereas the structures at low specific energies are due to the delta electron tracks. On the other hand, the absence of the triangular pattern at any target size in the f1(z) of the 60Co comes as a result of the tortuous path of the electron tracks. The triangular part of f1(z) shifts to lower z values with increasing particle energy.
The resulting f(z,D) for D equal to 1 and 5 Gy are shown in Fig. 3 for spheres with diameters of 0.01, 0.5 and 10 μm. The distinctive triangular shape seen in some f(z,D) curves are reminiscent of the triangle observed in f1(z) for protons and carbon ions. For these cases, the mean number of tracks n is approximately one, causing f(z,D) to be essentially equal to f1(z) renormalized by the Poisson factor of Eq. (2). It is important to point out that when n≪1, the dominant term of the f(z,D) distribution is the no-hit term i.e., f(z = 0,D). At this point it is much more likely that the particle track will miss the target entirely than to hit it. But if a hit occurs, the deposited specific energy can exceed up to 5 orders of magnitude of the mean dose (see the cases of the 0.01- and 0.5-μm diameter spheres shown in Fig. 3).
If either the target size or the mean dose is increased, the number of convolutions in Eq. (2) also increases, eventually causing f(z,D) to become normally distributed around the mean regardless of radiation quality. The actual target size or mean dose at which such transition occurs can vary from one quality to another. In Fig. 3 it can be observed that 1 Gy given to a 10-μm diameter sphere, the f(z,D) is normally distributed (note that the log-log plot causes a visual skewness of symmetric distributions) when irradiated by 60Co or 0.91 MeVu–1 protons but not when irradiated by 5.27 MeVu–1 carbon ions. If the dose increases above 5 Gy then f(z,D) becomes normally distributed for these carbon ions.
Fitting of the LQ Parameter
The behavior of the parameter afit in Eq. (4) as a function of target diameter for different radiation species was explored by fitting to available experimental values of α. The results expressed as the ratio afit/αexp for 60Co, 0.91 MeVu–1 protons, and 76.90 and 5.27 MeVu–1 carbon ions are shown in Fig. 4. As expected, this ratio tends towards unity with increasing diameter regardless of radiation species, but the target size at which the ratio reaches unity varies with radiation quality. For 60Co, a 2-μm diameter sphere is sufficient to provide equality (less than 1% difference) between parameter afit and α, whereas for the lowest carbon energy the ratio is approximately 1.3 for a 12-μm diameter sphere.
More interestingly, our results demonstrate the existence of a threshold in target size below which there is no satisfactory fit of Eq. (5) as the residuals become too large (>>10-fold). This is indicated in Fig. 4 by the dashed lines. It can also be observed that such threshold size varies among radiation qualities. Our results show that the threshold size for 60Co [lowest linear energy transfer (LET) radiation in this work] is between 0.1 and 0.5 μm. As LET increases, so does the threshold size, reaching between 6 and 7 μm for 5.27 MeVu–1 carbon ions (highest LET radiation in this work).
Microdosimetric Restriction on the Sensitive Volume
The existence of a threshold size discussed in the previous section is closely related to an intrinsic geometrical limit imposed by the microdosimetric formalism [e.g. (19)]. Under the assumption that cells are autonomous, the survival fraction is the conditional probability of survival P(S) given the probability of no-hit P(NH) (the probability that a cell is missed by a particle and/or its secondaries) and the probability of hit and repair P(HR) (the probability that a cell hit by the particle and/or its secondaries survives as a result of a successful damage repair), i.e.,
Recalling that n is the mean number of tracks passing through a target exposed to a mean dose D, then P(NH) is numerically equal to the Poisson probability with v = 0 (no tracks), i.e., P(NH) = e–n. It is worth noting that P(NH) can also be extracted by integrating the surviving term of f(z,D) when setting v = 0, i.e., . Since f0(z) is a Dirac delta distribution at z = 0, the integral will yield the same result.
Comparison between the no-hit probability as a function of dose for different target sizes with the SF(D) clearly illustrates the existence of a lower limit in size of the sensitive volume within a cell. Namely, target sizes whose P(NH) curves are greater than the experimentally measured SF(D) violate Eq. (9) and cannot be used to describe the linear term of SF(D). Figure 5 shows such comparison for V79 cells. As expected, the slope of the P(NH) curves increases with decreasing volume indicating that smaller targets are easier to miss, regardless of radiation quality. Yet, the limit in sensitive size varies with the microdosimetric characteristics of the ionizing radiation tracks, e.g., increasing limit size with decreasing initial kinetic energy of the particle.
The range of limits described in this section is equivalent to the range of threshold sizes found by the fitting process in the previous section. This similarity is clear evidence that the single track microdosimetric distribution of energy deposition in relationship to the size of the available sensitive size within a cell affects, to some extent, the probability of survival of the cell.
Mean Inactivation Dose
The mean inactivation dose (MID), defined as the area under the survival curve,19). Dividing MID by (mean specific energy deposited by the passage of one track) will yield the mean number of inactivation tracks, m. In other words, m is the fraction of n (defined in the microdosimetric background section), that will cause (on average) a lethal event. It follows that the volume at which one track is expected to become lethal can be estimated by looking for the target size at which equals the MID (i.e., m = 1). Following the results of the previous sections, it is anticipated that this “single-track lethal-volume” also varies with radiation quality. Knowledge of this volume enables the calculation of a special frequency-mean lineal energy for which the sizes are tuned to yield m = 1 and can thus be seen as a radiobiological property of the radiation quality. From here on, we will use the sub-index MID to emphasize that those values of have been calculated with such target sizes.
In this work, calculation of values started by extracting the functions and (d being the diameter of spherical targets) from the set of f1(z) frequencies calculated via MC simulations (see Materials and Methods) for the various radiation qualities. The MID values were calculated with Eq. (10) employing the experimental LQ parameters given in Table 1. Then, the function was searched for the condition m = 1 to obtain the diameter of the “single-track lethal-volume” per radiation quality. Finally, the corresponding functions were evaluated at such diameters to acquire the desired values.
Figure 6 shows the diameter and the deposited energy of the “single-track lethal-volume” as a function of . The smooth increase of both curves observed across radiation species is noteworthy. Both curves can be easily parameterized by a relationship of the type , which suggests that could be a potential candidate for radiation quality characterization for modeling. It is worth pointing out that such characterization is limited to the biological conditions under which the clonogenic assay has been performed, thus, the results in Fig. 6 apply only for V79 cells.
All microdosimetric quantities are target size dependent. This represents a hindrance for radiobiological models based on micro- and nanodosimetry as the actual sensitive volume sizes within a cell cannot be experimentally determined. However, radiobiological models like the fourth version of the local effect model (LEM IV) (20), Katz's track structure model (21) and the modified microdosimetric kinetic model (mMKM) (22) assume that the critical target is found within the cell nucleus (micrometer size) and that this target can be further divided into small independent nanometer-sized sub-volumes. In these models, cell survival depends on the accumulated severity of the lesions within the sub-volumes, often referred to as domains. Accordingly, the key to modeling survival curves is to find a relationship between the lesions and distribution of energy deposition within the domain. Therefore, most modeling work would benefit from correct estimates of the domain's volume. Meanwhile, the parameter for “cell nucleus target” volume plays a minor role, appearing as a fixed value in the above mentioned, as well as other more recent, modeling approaches (23–25). The results of this work suggest that a more accurate estimation of the size of this target can be calculated and that perhaps such volume plays role that is equally as important as that of the domain because of its relationship with the ionizing radiation's probability of interaction. Indeed, the relevance of such probability (in the form of the action cross section) is hinted at in Katz's track structure model (21), as it directly determines the “ion-kill” component of the cell survival fraction. The “gamma-kill” component is, however, dependent on a fitted value of the geometrical cross section of the cell nucleus target.
The hypothesis behind this study is that the expected biological effect (e.g., measured by SF(D)) is a result of the averaging of the individual effects caused by the radiation's intrinsic variation of energy deposition (e.g., specific energy) at the level of subcellular volumes. Thus, the natural candidate for expressing the SF(D) under these terms is the f(z,D) distribution because it explicitly describes the probability of energy deposition for a given target size and given mean dose [see Eq. (4)]. The biological effect in Eq. (4) is carried by the parameter a and it is expected to exert its effect solely on the sensitive volume of the cell. To find the parameter value of a, a fitting process to published experimental survival curves was performed. For each radiation quality tested here, a range of target sizes yielded suitable afit values through the application of Eq. (8). Figure 4 shows that as target size increases, the value of afit approximates the experimental LQ value of α for all radiation qualities, albeit at different rates. Indeed, at larger volumes, f(z,D) becomes normally distributed around the mean dose D with a small enough spread that translates, within our hypothesis frame, into an almost negligible difference between the cell's individual biological effect and the mean overall effect.
Although no further mechanistic conclusions can be drawn from this behavior, our study revealed the existence of a threshold in target size below which the fit fails. More interestingly is that these threshold sizes may significantly vary between radiation species (spheres with diameters between 0.1 and 0.5 μm for 60Co to diameters between 6 to 7 μm for carbon ions). The mathematical reason for the lack in the goodness-of-fit is that at the threshold size the first term of the f(z,D) summation [i.e., f0(z)e–n] becomes the dominating factor, especially at low doses (generally below 3 Gy for the ion energies tested here). Consequently, there is no valid afit value that can satisfy Eq. (8) across the dose levels commonly used in cell survival experiments.
At first, this finding seems to challenge the assumption that cells have a unique sensitive volume. However, the notion of a limit in sensitive volume that varies with radiation quality comes out naturally from the mathematical constraint that the probability of no-hit at each dose has to be less than or equal to the probability of cell survival (see Results). Essentially, the relationship between the microdosimetric distribution of energy deposition (which is a property of radiation quality that changes with target size) and the measured cell survival sets a minimum (threshold) volume within a cell that the ionizing track should affect, to cause a lethal effect. This concept directly influences the modeling of the initial slope of the survival curve (e.g., α parameter) because it is closely related to the lethal potential of the single-track hit. Published models for heavy ions such as LEM IV [e.g., (26)], the mMKM model [e.g., (27)] and multi-hit single-target model (28) agree well with experimental survival data when using cell nucleus target volumes corresponding to spheres with diameters between 6 and 8 μm. These volumes are in agreement with our estimated threshold volume for the lowest-energy carbon ion. This suggests that regardless of the model's theoretical framework, to model data for lower-energy particles, it is necessary to introduce a target that is at least the size of its threshold volume. Data for particles of higher energy will be naturally covered as the threshold volume decreases with increasing particle energy. It is also important to clarify that the cell nucleus target does not necessarily equal the visible size of the cell's nucleus, as shown by Hawkins (29) for V79 cells in different stages of the cell cycle.
Additionally, we calculated the magnitude of such threshold volume by finding the target size at which the mean specific energy for one track equals the MID. Thus, the threshold volume can be interpreted as the volume that, if hit by one track, will provoke a lethal event in the cell. The diameters of these spherical volumes were then used to calculate a special kind of mean lineal energy , which can be regarded as a radiobiological property characteristic of radiation quality. The surprisingly smooth curves that result from characterizing radiation quality by its values (see Fig. 6) hints at the property's potential in modeling. Microdosimetric measurements of mean specific energy as a function of target size can be readily obtained for hadron-therapy beams by using state-of-the-art detectors. These measurements, coupled to a consistent cell survival database, would yield values that can characterize mixed field beams at various depths. Moreover, such experimental output can provide a reliable benchmarking of Monte Carlo codes.
The spatial distribution of energy deposition as described by f(z,D) is notably influenced by changes in target size. The no-hit peak becomes the dominant factor as volume decreases for any radiation quality. This behavior imposes a natural limit on the possible sensitive size (at the cell nucleus level) of the cell when its relationship with cell survival data is studied. Moreover, such a threshold volume is radiation quality-dependent, increasing with decreasing particle energy. This finding led to the derivation of a potentially useful microdosimetric property , which characterizes ionizing radiation across radiation modalities in relationship to their cell inactivation capacities.
We gratefully acknowledge the Swedish Radiation Safety Authority (SSM) for their financial support. The computations were performed on resources provided by SNIC through the Uppsala Multidisciplinary Centre for Advanced Computational Science (UPPMAX) under project nos. p2011144 and snic2016-7-92.