A novel microdosimetry-based formalism for cell survival modelling applicable to hypofractionated radiotherapy

TL;DR

This paper introduces a new model for predicting cell survival in high-dose radiotherapy by combining microdosimetry with a non-linear cell damage approach.

Contribution

A novel microdosimetry-based formalism that determines sensitive volume size through optimization and improves cell survival modeling for hypofractionated radiotherapy.

Findings

The model successfully matched 24 proton survival curves across various LET ranges and irradiation conditions.

The sensitive volume size varies with linear energy transfer (LET) and fluence spectrum, unlike standard microdosimetry.

The approach simplifies multi-particle microdosimetric spectra calculations and eliminates arbitrary assumptions about sensitive volume size.

Abstract

Hypofractionated radiotherapy requires reliable cell survival models for doses much higher than the standard 2 Gy, for which the linear-quadratic (LQ) model is not applicable. We developed an alternative approach applicable to both low doses and high doses used in hypofractionated treatments and radiobiological experiments.

We combined a standard microdosimetric technique with a recently introduced non-LQ cell survival model. Our formulation accounts for cell damage by multi-track events involving any number of particles. This is necessary for modelling cell survival at therapeutic doses. We characterise each cell type by the size R of the sensitive volume (SV) and biological response function B(q), where q is the total energy deposited in the SV after a given dose is delivered. q is a random quantity characterised by a probability density (microdosimetric spectrum) calculate with Geant4. We determine R and B(q) through an optimisation procedure that minimises differences between model predictions and cell survival measurements that cover an appropriate linear energy transfer (LET) range.

Our method eliminated a serious flaw of the standard microdosimetric approach—arbitrary SV size. We determine SV size by solving the above optimisation problem. Furthermore, our method drastically simplifies calculations of multi-particle microdosimetric spectra. We applied our approach to 24 proton survival curves for three cell lines with various irradiation conditions and LET range of 0.589-19.6keVμm-1 with good agreement between all these measurements and the model. R for a given cell type depended on fluence spectrum and increased with increasing LET owing to variations in the development and spatial spread of damage triggered by initial physical impact. This differ from the standard microdosimetry where SV size is constant.

Our model is relatively simple and suitable for implementation in a treatment planning system potentially improving treatment plan optimisation, calculation of RBEs and biologically equivalent doses.

Figures

Full text

Introduction

Hypofractionation and the use of accelerated charged particles are growing trends in radiotherapy (Durante and Formenti 2018). One of the challenges of using high doses per fraction is the well-documented limitations of the linear-quadratic (LQ) model (Elkind and Sutton 1959, Kirkpatrick et al 2008, Iwata et al 2009, Qiu et al 2020). For megavoltage x-rays the LQ model is ‘inappropriate’ when fractional dose is higher than 8–10 Gy (Qiu et al 2020). Above these doses, the survival curve (ln $[eqn]$ vs dose) becomes a straight line. This threshold is approximate and subject to biological variability. In addition, there are no data to support these values for protons. The shape of proton survival curves depends on linear energy transfer (LET). The transition point where a curved shape that can be approximated by the LQ model becomes a straight line, which is inconsistent with the LQ dependence, shifts to lower doses when LET increases. At sufficiently high LET values this point can be at zero dose. LET distributions in proton treatment plans are rather complex. Adding to this problem, large spatial dose variations in hypofractionated treatments make the LQ model a poor choice for planning proton treatments. It needs to be replaced by a better founded model applicable to a broad range of doses, for any realistic LETs.

Several methods have been developed for improving predictions of cell survival at high doses (for example, Curtis (1986), Scholz and Kraft (1994), Guerrero and Li (2004), Carlone et al (2005), Beuve et al (2008), Park et al (2008), McKenna and Ahmad (2009), Guerrero and Carlone (2010), Hanin and Zaider (2010), Wang et al (2010), Vassiliev (2012), Zhao et al (2020), Li et al (2021)). However, these studies were largely focused on x-rays and provided only limited validation data. One notable exception is the study by Scholz and Kraft (1994). Their method became a part of the local effect model (Kramer and Scholz 2000) that has been extensively tested. The downside of that method is that it is not based on any real mechanism. Survival curves are modified essentially by hand to give them a more realistic shape.

The present study builds on a recently proposed non-LQ model (Vassiliev 2023). The model is based on two very general principles: (1) during a dose fraction delivery cell radiosensitivity increase owing to the accumulation of latent damage (i.e. all initial damage that is subject to transformation by relatively slow biological processes, and ultimately is either repaired or manifests as more permanent damage including cell death); and (2) radiosensitivity is bounded from above and stops increasing with dose when its maximum is reached. At this maximum point, the curved part of the survival curve (ln $[eqn]$ vs dose) ends and the linear part begins. The straight line reflects inability to repair any additional damage, which means a predominantly single-hit cell kill mechanism. For this reason, the slope of the line is determined by a single parameter: the frequency average specific energy $[eqn]$ that strongly depends on sensitive volume (SV) size. Vassiliev (2023) reviewed properties of the non-LQ model and evaluated its performance using hundreds of survival curves, for protons, x-rays and heavy ions. A particularly useful aspect of the model is that its mathematical formulation is perfectly compatible with the microdosimetric formalism, as we demonstrate in the next section.

We applied our approach to a set of 24 proton survival curves for three cell lines with different radiosensitivities: AG01522, U87 and V79 from four sources: Blomquist et al (1993), Robertson (1994), Chaudhary et al (2014), Wouters et al (2015). The data cover an LET range of $[eqn]$ . Out of abundant data for V79 cells we chose 12 experiments that provide almost uniform LET sampling over the range: $[eqn]$ , which is most relevant for treatment planning. Using data from multiple sources and representing a variety of beam properties, helps us reduce the risk of systematic errors due to selection bias and allows us to see the magnitude of real uncertainties. We present step by step our optimisation procedure, which is substantially different from the standard microdosimetric approach to survival modelling. Then we review the properties of model parameters and outline how our method can be implemented in a treatment planning system.

Methods

Combining microdosimetry with the non-LQ cell survival model

2.1.

The purpose of the study by Vassiliev (2023) was to explain at a fundamental level the general shape of cell survival curves ln $[eqn]$ versus dose, in particular, the transition to linear dependence at high doses. That study also provided a formula for cell survival based on the single-target multi-hit concept. This formula produces the correct survival shape at high doses. The formula is also accurate in the low-dose limit. The first two nonzero terms in the power series expansion in dose coincide with the LQ formula. An innovative aspect of that study is that for a given cell type the SV size is not fixed. It depends on beam quality, because in this approach, SV size characterises the range of interaction between damaged DNA sites. Thereby it accounts for the spread of damage around the particle path by physical, chemical and biological processes triggered by the initial radiation impact. SV size increases with increasing LET. This conclusion was based on a review of extensive experimental data. The same trend was also observed for inactivation cross sections.

A note on the terminology: in microdosimetry, the term ‘site’ is often used instead of ‘target volume’; we will use the term ‘sensitive volume’ (SV) for consistency with our previous work.

To model the LET, or more generally beam quality, dependence of cell survival we employ the formalism of microdosimetry. We start with a standard technique (ICRU 1983, Vassiliev 2017). Then, we combine it with the non-LQ formula (Vassiliev 2023). This results in substantial improvements compared with the standard microdosimetry. In microdosimetry, radiation impact after a dose $[eqn]$ is delivered is characterised by the probability density $[eqn]$ of the total energy $[eqn]$ deposited in the SV by all particles. It is normalised to 1 excluding the singular point $[eqn]$ :

$$[eqn]$$

Then we introduce the biological response function $[eqn]$ that maps distribution $[eqn]$ onto cell survival

$$[eqn]$$

$[eqn]$ is the probability of a cell surviving energy $[eqn]$ deposited in SV. $[eqn]$ is undefined for $[eqn]$ . Elsewhere, it is bounded $[eqn]$ , nonincreasing $[eqn]$ , and satisfies $[eqn]$ .

Equation (2) is the basis of the standard microdosimetric approach to modelling cell survival. To model cell survival for a given cell type, an optimisation procedure is executed that finds such $[eqn]$ that minimises differences between a set of survival curves (left side, equation (3)) and model predictions (right side):

$$[eqn]$$

In this equation, $[eqn]$ is $[eqn]$ th survival curve in a cell survival data set. The total number of survival curves in the data set is $[eqn]$ . Each survival curve is a set of measurements at doses $[eqn]$ , where $[eqn]$ is the curve number and $[eqn]$ is the number of a dose point. The total number of dose points in $[eqn]$ th curve is $[eqn]$ .

The above approach has two serious problems. First, distribution $[eqn]$ strongly depends on the size of SV. However, ‘there is no clear indication of the relevant site size for specific biological endpoints’ NCRP 2018, 181. Second, $[eqn]$ needs to be calculated for each dose in the set $[eqn]$ . These calculations must account for the fact that the number of particles that contribute to the total deposited energy $[eqn]$ is random in interval 0 to $[eqn]$ . Such calculations are time consuming. We eliminate both these problems by applying equations (2) and (3) to the non-LQ cell survival formula proposed by Vassiliev (2023). The formula is simple and most general

$$[eqn]$$

where

$$[eqn]$$

is the average number of hits at a given dose, $[eqn]$ is the frequency average specific energy, and $[eqn]$ is the probability of a cell surviving $[eqn]$ hits (i.e. energy deposition events). $[eqn]$ characterise the ability of a cell to repair damage. Obviously, $[eqn]$ and we can expect that $[eqn]$ .

We can also represent $[eqn]$ by a similar sum (Vassiliev 2017):

$$[eqn]$$

The above is simply the total probability formula, in which we assumed that the number of hits is Poisson-distributed. Here $[eqn]$ is a conditional distribution of $[eqn]$ . The condition is that the number of hits is $[eqn]$ . All $[eqn]$ have the same normalisation as $[eqn]$ , equation (1). Next, we use equation (4) in the left side of equation (2) and use equation (6) in the right side:

$$[eqn]$$

This equation must be satisfied for any $[eqn]$ . Therefore

$$[eqn]$$

Equations (4) and (8) drastically simplify the optimisation problem. Equation (3) is now replaced with

$$[eqn]$$

where $[eqn]$ is given by equation (8). The first simplification is that $[eqn]$ does not depend on dose, it is a function of the number of hits, $[eqn]$ . Second, for all proton cell survival curves analysed in the present study, we can set $[eqn]$ for $[eqn]$ or $[eqn]$ depending on the cell type. This means that instead of the dose-dependent distribution $[eqn]$ we need to calculate no more than four dose-independent distributions $[eqn]$ . The most important improvement is that, in contrast to standard microdosimetry, in our method SV size is not chosen arbitrarily, it is now a fitting parameter. This is because $[eqn]$ strongly depends on SV size, $[eqn]$ :

$$[eqn]$$

where $[eqn]$ is the frequency average lineal energy, LET_F_ is the frequency average LET, and $[eqn]$ is the mass density of material within SV. $[eqn]$ can be approximated by LET_F_ when $[eqn]$ (Vassiliev et al 2019). Variable SV size increases the amount of calculations, because distributions $[eqn]$ has to be calculated for each size. These variations, however, are a real effect that must be accounted for. This effect was overlooked in the standard microdosimetry.

Microdosimetric Spectra

2.2.

This study is limited to protons. We also consider photon data as a reference radiation. For photons we use a method that does not require microdosimetric spectra. For protons, first, we calculate the single-event spectrum, $[eqn]$ , using Geant4 and a Monte Carlo algorithm previously described and validated (Vassiliev et al 2019). The algorithm accounts for the nanoscale dose buildup and indirect events, in which a proton misses SV but secondary electrons reach it and deposit energy. Next, using $[eqn]$ as input, we calculate multi-event spectra $[eqn]$ by numerical integration, according to the recurrence relation (Vassiliev 2017):

$$[eqn]$$

Energy distribution of protons entering an SV corresponding to an experimental survival curve $[eqn]$ is given by the fluence spectrum $[eqn]$ . Fluence spectra are rarely reported in experimental papers. Usually, only experimental setup, basic beam parameters, and depth-doses are reported. We used this information to setup Monte Carlo calculations of all $[eqn]$ . We used TOPAS (Perl et al 2012, Faddegon et al 2020) for these calculations. We adjusted our setup parameters, including those of the range modulation wheel, until agreement between Monte Carlo and reported depth-doses was achieved. These adjustments were performed using optimisation software. Given the large number of adjustable parameters, we were able to reproduces all the experimental depth-doses with good accuracy. Using best-fit values of the parameters, we generated fluence spectra at each depth at which survival curves were measured.

Optimisation methods

2.3.

For implementation of our formalism into a treatment planning system, in addition to the distribution of physical dose $[eqn]$ , beam quality needs to be characterised at every point $[eqn]$ of the irradiated volume. For this purpose we chose fluence spectrum, $[eqn]$ . This fundamental quantity can be found by solving the Boltzmann transport equation by Monte Carlo or a deterministic method. If $[eqn]$ is known, many other beam quality characteristics can be found: average LET, LET spectra, microdosimetric spectra, and others.

Specific energy zF and SV radius, R

2.3.1.

To determine $[eqn]$ we develop a regression model and find its parameters using an optimisation procedure. $[eqn]$ is represented as an integral of weighted normalised fluence spectrum

$$[eqn]$$

The weighting function $[eqn]$ is approximated by a third order polynomial of $[eqn]$ .

$$[eqn]$$

We chose this particular form based on a review of the relevant data.

To account for variable SV size we employ the following five-step optimisation procedure.

All survival curves for a given cell type are fit separately using equation (9). The fitting parameters are: $[eqn]$ and $[eqn]$ . Microdosimetric spectra are not used at this step. The number of parameters $[eqn]$ is infinite. However, we found that it suffices to fit only a subset of 2–3 consecutive $[eqn]$ . All other $[eqn]$ are set to 1 (small $[eqn]$ ) or to 0 (large $[eqn]$ ). If $[eqn]$ or $[eqn]$ , then we do not need to calculate the corresponding $[eqn]$ at any further steps. We fit each survival curve several times using different $[eqn]$ subsets, until we identify one that provides best fit. In half of all cases in our proton data set $[eqn]$ is the optimal choice. All $[eqn]$ are identified, and the first iteration of $[eqn]$ set is generated.We use the set of best-fit $[eqn]$ to determine parameters $[eqn]$ of the regression model for $[eqn]$ , equations (12) and (13). We need this model, because it allows us to calculate $[eqn]$ for an arbitrary proton beam, $[eqn]$ . Then, for each experimental survival curve, we generate a new $[eqn]$ , now using the regression model. The first iteration of the regression model and the second iteration of $[eqn]$ set are generated.For each updated $[eqn]$ we determine SV size $[eqn]$ . We cannot use for that equation (10), because $[eqn]$ is unknown. Instead, we use the try-and-error method. We pick a reasonable $[eqn]$ , then calculate with Geant4 $[eqn]$ and the corresponding $[eqn]$ . Then we adjust $[eqn]$ using the $[eqn]$ dependence given by equation (10) until Monte Carlo $[eqn]$ agrees within ±1% with $[eqn]$ predicted by the regression model. The $[eqn]$ dependence is quite accurate, so it usually takes only 2–3 iterations to determine $[eqn]$ . After completing this step we will also have calculated all $[eqn]$ for all correct sizes $[eqn]$ . All $[eqn]$ for the current iteration of $[eqn]$ set are determine and corresponding $[eqn]$ are calculate. We now have the initial state for global optimisation.First global optimisation. We fit simultaneously all survival curves for a given cell type. We optimise parameters $[eqn]$ and biological response function $[eqn]$ . We discuss the latter part in more detail in the next section. Re-optimisation of $[eqn]$ produce new set of $[eqn]$ . Then, as described above, we determine new $[eqn]$ and a new set of distributions $[eqn]$ . Final iterations for the regression model and $[eqn]$ set are derived; $[eqn]$ and $[eqn]$ are updated accordingly; the first iteration of $[eqn]$ is derived.Final global optimisation. This step is similar to step 4, except now only $[eqn]$ is optimised. All other parameters remain fixed. The final iteration of $[eqn]$ is produced.

The biological response function, B(q)

2.3.2.

For each cell line we determine a unique function $[eqn]$ . Our procedure is based on equation (9). Simultaneous fitting of multiple cell survival curves poses two challenges. First, low-dose points contribute much more to standard loss functions, such as the sum of squared residuals, than do high-dose points. Second, there are large variations in uncertainties between individual data points. The common practice of weighting each point by inverse standard error is wrong. Uncertainties of estimated errors are higher than those of the survival data. Also, the estimated errors are random numbers, and dividing by random numbers usually results in bad statistical properties. In short, with this weighting, the optimiser is largely driven by noise.

To achieve a better balance of contributions to the loss function $[eqn]$ from high- and low-dose points, we fit ln $[eqn]$ instead of $[eqn]$ . If survival was exponential, this method would produce the standard linear regression. To mitigate the influence of large variations in uncertainties of survival data, we calculate for each survival curve the Huber loss $[eqn]$ (Huber, 1964). In this method, based on statistics of residuals, points are identified that are likely to have a high uncertainty, and contributions from such points to $[eqn]$ are reduced. For fitting simultaneously multiple survival curves, we minimise the following loss function

$$[eqn]$$

The biological response function $[eqn]$ must have certain properties that we listed previously. A properly scaled hyperbolic tangent satisfies this requirement. Hence, we chose the following form for $[eqn]$ :

$$[eqn]$$

This model has three independent parameters: $[eqn]$ , and $[eqn]$ . Parameter $[eqn]$ can be found from the condition $[eqn]$ :

$$[eqn]$$

Implementation into a treatment planning system

2.4.

It is impractical to use microdosimetric spectra in treatment planning software. Fortunately, there is a method for replacing spectra $[eqn]$ with fluence spectra $[eqn]$ . The method requires a few pre-tabulated 2D functions that map the $[eqn]$ space onto $[eqn]$ . For an explanation of the algebra, we refer to Vassiliev et al (2019). First, we need to introduce microdosimetric spectra normalised per one source particle, $[eqn]$ . These functions have a singularity at $[eqn]$ , because the probability of a source particle missing SV entirely is a finite number. We also introduce the probability of a non-zero energy deposit $[eqn]$ . From the formula for conditional probability, we can write:

$$[eqn]$$

Further, for polyenergetic beams, we have

$$[eqn]$$

and a similar formula is valid for $[eqn]$ . Next, using the above expressions, we transform equation (8):

$$[eqn]$$

where $[eqn]$ is a biological response function that depends on proton energy

$$[eqn]$$

Then, to implement our method in treatment planning software, a small number of these functions, for example, $[eqn]$ and $[eqn]$ need to be tabulated. These will be two-dimensional tables, because $[eqn]$ also depend on $[eqn]$ . With these tables probabilities $[eqn]$ are calculated as a simple one-dimensional integral. Generating tables of $[eqn]$ involves extensive Monte Carlo simulations. However, for any given cell type represented by a $[eqn]$ , these calculations are performed only once.

We will also need to determine $[eqn]$ for any given fluence spectrum. First, we determine $[eqn]$ using the regression model. Then, for a given $[eqn]$ we solve this equation for $[eqn]$

$$[eqn]$$

This requires pretabulated $[eqn]$ derived using the method we used above for $[eqn]$ .

Results and discussion

Cell survival data-set

3.1.

All cell survival data analysed in the present study are listed in table 1. Chaudhary et al (2014) results were chosen because they extend to high LET values. We also chose V79 cells, because abundant data are available. We were able to select 12 experiments that almost uniformly sample an important for treatment planning LET range: $[eqn]$ . Throughout the present study ‘LET’ means frequency average LET. This is for the following reasons: (1) consistency with the microdosimetric terminology; (2) ‘track average LET’, when it was first introduced, meant initial kinetic energy of a heavy charged particle divided by its total range, and this is not what we and many recent studies calculate; (3) dose average quantities do not appear anywhere in our formalism; (4) frequency average LET has a good predictive power whereas dose averaging, if done correctly, can produce value strongly contradicting relative biological effectiveness (RBE) data (Vassiliev 2021).

Reference radiation

3.2.

The reference radiation for AG01522 and U87 cells was 225 kVp x-rays. For all V79 cells measurements, the reference radiation was ^60^Co. A simple formula is available (Vassiliev et al 2020) for the conversion of RBE values from one reference photon radiation to another. Here, we did not do this conversion. We used equation (9) to fit these data. The fitting parameters were $[eqn]$ and survival at zero dose, $[eqn]$ . $[eqn]$ was not constrained to be 1 in order to account for uncertainties in plating efficiency. However, it was bounded within ±0.05. We do not have sufficient data to determine $[eqn]$ for photons. Hence, we did not use microdosimetric spectra in this case. Figure 1 shows results of the fit. The significant scatter of V79 data points reflects discrepancies between data from different sources. The best fit values of all the parameters are given in table 2. Two observations can be made. First, most cells survive single hit, even when irradiated by 225 kVp x-rays with an $[eqn]$ (Vassiliev et al 2020). Second, radiosensitivity is a result of the interplay between $[eqn]$ , which is related to the range of interaction between damaged sites, and $[eqn]$ that characterise the ability of a cell to repair damage.

To apply our cell survival model to photons, $[eqn]$ are not necessary. It is easier for each beam energy, for example 6 MV, 6 MV FFF, 10 MV, etc determine three parameters, $[eqn]$ and two $[eqn]$ values, as table 2 shows.

Microdosimetric spectra

3.3.

Calculation of $[eqn]$ spectra with Geant4 is the most time consuming part of our approach. For each of the 24 survival curves, these calculations are repeated about five times using different SV sizes. Transforming an $[eqn]$ spectrum into a multi-particle spectra $[eqn]$ , takes a few seconds. This transformation includes integration that reduces statistical uncertainties. All microdosimetric spectra in the present study are for polyenergetic proton beams. Figure 2 shows an example of calculated spectra from one of the cases used in the present study.

zF and SV size R

3.4.

Table 3 presents best-fit values of parameters $[eqn]$ of the regression model for $[eqn]$ . The corresponding weighting functions $[eqn]$ are plotted in figure 3. The logarithmic scale was chosen for a better view of the most variable parts of the curves.

Figure 4 shows $[eqn]$ as a function of $[eqn]$ for the three cell lines. $[eqn]$ determines cell survival at high doses. The asymptotic expression for all cell survival curves in the high-dose limit was derived previously (Vassiliev 2023). The result is that ln $[eqn]$ is a linear function of dose with the slope $[eqn]$ :

$$[eqn]$$

This means that a larger $[eqn]$ value (and therefore a smaller $[eqn]$ value) indicate higher radioresistance at high doses. The corresponding SV sizes are shown in figure 5. For all three cell lines, $[eqn]$ tends to increase with increasing LET_F_. This trend and the numerical values are consistent with our previous study and with the behaviour of inactivation cross sections (Vassiliev 2023). $[eqn]$ is a single biological parameter that determines the slope of cell survival curves at high doses. The slope depends strongly on $[eqn]$ . From equations (21) and (22), $[eqn]$ . For $[eqn]$ depends only weakly on $[eqn]$ and, in fact, can be approximated by LET_F_ Vassiliev et al (2019). Hence, $[eqn]$ . These observations emphasise the high importance of variable SV size that our model has introduced.

Biological response function, B(q)

3.5.

Best-fit values of all the parameters that define $[eqn]$ are given in table 4. These functions for the three cell types are plotted in figure 6. They all originate at $[eqn]$ . Then, they have an almost flat part indicating that most cells will survive $[eqn]$ (i.e. three ionisations). Then, a steep, approximately exponential, decrease begins. The interpretation of figure 6. is complicated by the fact that $[eqn]$ alone does not determine $[eqn]$ . Instead, we need to consider the product $[eqn]$ , according to equation (8), where $[eqn]$ depends strongly on SV size, $[eqn]$ . For illustration purposes only, let us assume that $[eqn]$ . Then, figure 6 suggests that U87 cells are more radioresistant at low doses than AG01522. Further, we assume that $[eqn]$ . Assuming that all cells receive the same dose, this means that the average energy deposited in SV(V79) is 1/8 of that deposited in SV(U87). This difference is reflected in figure 6 by a substantial shift of the V79 curve towards lower energies. The reason for the shift is that $[eqn]$ for V79 does not reach energies $[eqn]$ as high as those reached by $[eqn]$ for VU87.

Cell survival and RBE

3.6.

In figures 7–10 best-fit survival curves are compared with experimental data. Good agreement was achieved for AG01522 and U87 cells. These were single source data (Chaudhary et al 2014). In contrast, we selected V79 data from several indendent sources that used a variety of beam energies and SOBP parameters. As a result, we observed discrepancies, sometimes substantial, between different studies. We, however, consider this as a positive property of our V79 set. These discrepancies reflect uncertainties in experimental data of this type. We divided our V79 results into two subsets based on $[eqn]$ values (given in the captions). $[eqn]$ is the dose at which 10% of cells survive. These results are shown in figures 9 and 10. Our model correctly outlines variations of the data with dose and beam quality. However, accurate fit was not achieved for each individual curve.

In figure 11, we show RBE as function of LET_F_ calculated from the best-fit survival curves. The RBEs were calculated at proton doses of 2 Gy and 10 Gy. To help understand this figure, we note that at 10 Gy, RBE is largely determined by $[eqn]$ values. At 2 Gy the damage repair ability (i.e. $[eqn]$ ) becomes a strong factor. A particularly interesting result is for V79 cells. In this case, the RBE at 2 Gy suggests a high probability of damage repair. However, at 10 Gy owing to a rather high $[eqn]$ , the RBE decreases to values close to and below 1 (relative to ^60^Co). Please note that all V79 data are associated with high uncertainties. Table 2, figures 1, 9 and 10 give a good idea of their magnitude.

Finally, figure 12 summarises mechanisms of cell-kill by protons as a function of dose. This figure shows that the main mechanism of killing AG01522 cells at 2 Gy and 10 Gy involves four or more protons. For V79 cells the main mechanism at 2 Gy involves two or three protons, and more than three protons mechanism dominates at 10 Gy. The smaller numbers of protons for V79 cells can be attributed to a smaller SV size. The single-proton mechanism dominates below 0.38 Gy (AG01522) and 1.22 Gy (V79). These results have serious implications for now popular studies of DNA damage that use microscopic Monte Carlo simulations. In a typical setup, a single proton is launched and DNA damage is recorded. Then, before the next proton is launched, the DNA is reset to its original undamaged state. Such techniques can predict only a single-proton mechanism, which is a minor factor at therapeutic doses.

Conclusions

The present study addresses the need for a cell survival model that is accurate over a broad range of doses, including high doses used in hypofractionated radiotherapy and in radiobiological experiments. We combined two previously existing approaches: a non-LQ cell survival model that predicts the dose dependence and satisfies the above requirement, and the most general version of the microdosimetric approach for predicting beam quality dependence. This combination has eliminated two serious problems of the standard microdosimetry: arbitrary choice of SV size and the need to calculate dose-dependent multi-particle microdosimetric spectra $[eqn]$ . In our model, SV size is essentially a fitting parameter. We found that SV size increases with increasing LET. In the standard microdosimetry SV size is constant. Variable SV size reflects real physical, chemical and biological processes that cause spatial spread of radiation-induced damage. Damage spreads by up to several hundred nanometres from the particle path.

We have tested our model by applying it to a set of 24 proton survival curves for three cell lines. The set includes a variety of irradiation conditions and covers a broad LET range. We achieved good agreement between all these measurements and the model. Our analysis shows that the shape of survival curves is a result of an interplay between cell ability to repair damage and the extent of spatial spread of damage. At sufficiently high doses, the cell’s ability to repair damage becomes a negligible factor. Our data also show that at therapeutic doses cell kill results predominantly from a combined damage from several protons. The single-proton mechanism is a minor factor and becomes negligible above approximately 10 Gy.

Our model is relatively simple and suitable for implementation in a treatment planning system. This model has the potential to improve treatment plan optimisation, calculation of RBEs and biologically equivalent doses.