A model for relative biological effectiveness of therapeutic proton beams based on a global fit of cell survival data

We introduce an approach for global fitting of the recently published high-throughput and high accuracy clonogenic cell-survival data for therapeutic scanned proton beams. Our fitting procedure accounts for the correlation between the cell-survival, the absorbed (physical) dose and the proton linear energy transfer (LET). The fitting polynomials and constraints have been constructed upon generalization of the microdosimetric kinetic model (gMKM) adapted to account for the low energy and high lineal-energy spectrum of the beam where the current radiobiological models may underestimate the reported relative biological effectiveness (RBE). The parameters (α, β) of the linear-quadratic (LQ) model calculated by the presented method reveal a smooth transition from low to high LETs which is an advantage of the current method over methods previously employed to fit the same clonogenic data. Finally, the presented approach provides insight into underlying microscopic mechanisms which, with future study, may help to elucidate radiobiological responses along the Bragg curve and resolve discrepancies between experimental data and current RBE models.

As pointed out in the main text, the linear density of the atomistic excitations and ionizations undergo a transition to a highly compact distribution as the primary protons slow down. Therefore, the deviation from a Poisson distribution appears to be significant at the end of the range of proton. Figures 1A and 1B in the main text illustrate the ionization events within the cellular dimensions generated by traversing a proton with initial energy of 80 MeV and 1 MeV. The MC simulation has performed by using Geant4 DNA. As illustrated in these figures, the number of events in low energies (e.g., 1 MeV) is orders of magnitude greater than high energies (e.g., 80 MeV). Such difference in compactness of the ionization clearly justifies deviation from the Poisson distribution as the proton keeps losing kinetic energy.
Refs. [28,29,57] describes the original approach in mathematical formulation of the cell survival in the standard microdosimetric kinetic model (MKM), starting from the repair-misrepair fixation model (RMF) [49][50] in a cell nucleus domain. In this approach, Hawkins has shown that the time-integrated solution of the linearized RMF mass-action equations, averaged over the ensemble of the cell nucleus domains, leads to the linear-quadratic dependence of cell survival on the deposited dose, the first two terms in Eq.(1). Moreover this model predicts to be a linear function of LET and a constant and independent of LET. A similar approach with a distinction of including the non-linear (quadratic) term in the solutions of the RMF mass-action equations that accounts for the chromosome misrepair binary endjoining leads to higher order deposited dose and LET terms. These terms in Eqs. (1) to (4)  Due to large fluctuations of energy deposition in sub-micrometer volumes, the ionizing radiation is characterized by the probability distribution of specific energy, and single event specific and lineal energies and their expectation values, the frequency-mean ( , , and dose-averaged ( , , as well as their higher order statistical moments. The cell survival fraction is therefore given by , 3 where is the mean lethal lesions, averaged over specific energy distribution, in the ensemble of domains in all cell nuclei. Linear solutions and LQ cell survival: First we consider the linear approximation in RMF model, Eq. 3 , in which 0, corresponding to Eq.(7) in Ref.
[57]. In this limit, the analytical solution of Eqs. 1 and 2 can be easily found using the Green's function method . 4 Here is the solution of Eq. (A1) in the linear approximation. It is straightforward to justify that the homogeneous solution of Eq. (A1) is identical to zero, thus we do not consider it in Eq. (A4). It is also a straightforward calculation to show , where is the Heavyside function, i.e., 1 if and 0 otherwise. The steps in calculating retarded Green's function, , include converting the integral equation, Eq. (A4), to a differential equation for by substituting (A4) in (A1) and imposing the initial condition 0 for 0 where 0. Similarly we define where . Here the bar over denotes energy deposition averaging on the ensemble of cell nuclei domains, specific to a lineal-energy distribution. For an acute radiation dose, , the solution of Eq.(A4), , leads to ̅ . By averaging over the lineal-energy distribution and all cell nuclei and their domains we obtain a linear-quadratic model in cell-survival where and . Here we use the identity ̅ ̅ [58] that accounts for the spatial averaging of the energy deposition fluctuations, where , and is the single event distribution function of specific energy deposition, a counterpart distribution function of the lineal-energy . Furthermore ̅ ; where ; is distribution function accounting for all events and is mean number of events and/or tracks. Applying a relation between and (see e.g., Eq.(II.28) in Ref. [57] or Eq. (8)  ; . 6 Here , , and (the radius of a spherical domain) are the phenomenological fitting parameters.

Non-linear expansion of RMF solutions, going beyond LQ cell survival:
We now turn to perform a perturbative expansion to calculate the non-linear solution of the RMF model, Eqs.
(A1) and (A2). To go beyond the RMF linear solutions presented above, we assume to be a small parameter, hence we expand about perturbatively and linearize the resulting massaction kinetic equation to obtain the dynamics of the small fluctuations describing deviations from linear DSB solutions. We define and recall Eq. (A1) to obtain a linear massaction equation for

. 5
Here is the exact solution of the RMF model. As we defined, is the linear solution of the RMF model, hence describes the difference between exact and linear solutions. It is more convenient to transform Eq. (A5) into a more compact form The last two terms in Eq. (A11) are the contribution of the terms omitted in Eq. (A8) due to linearizing in the limit where is negligible. To transform (A11) to a form similar to the linearquadratic model we must calculate the statistical fluctuations in microscopic dose deposition throughout the averaging over cell nucleus domains, assuming equivalence between the ensemble averaging over the domains and the spatial averaging of the energy deposition fluctuations over the cell nuclei. In general, these fluctuations can be recursively reduced to lower power fluctuations, namely where and , and using / , and a linear relationship between and [59,60], one may find and to be a power series in lineal energy and LET as given by Eq. (3) and (4). The constants in Eqs. (A14) and (A15) can be determined phenomenologically by fitting and to the experimental data as illustrated in the text, i.e., ∑ , , and ∑ , . As seen in these equations, the contribution from the energy loss fluctuations renormalizes the LQ biological parameters and to infinite orders in and .

Gaussian fluctuations
It is interesting to calculate , and for a widely used Gaussian distributed function as in the limit ̅ ≫ 0, the Poisson distribution can be approximated to a Gaussian (central limit theorem) with variance ̅ Similar to a general energy loss fluctuations, the contributions of the Gaussian fluctuations in energy deposition introduce correction factors to and . Because the contribution to ̅ from and beyond are identical to zero, there is no correction beyond the linear term in and hence in . This is evident from Eq. 12 as the lowest order correction to specific energy and deposition dose is quadratic. Hence the Gaussian model predicts to be only linear dependence on lineal energy and LET 1 2 . 16 Similarly for the corrections beyond fluctuations vanish exactly as the lowest order correction to the dose in Eq. 12 is cubic. Hence we find a closed form for 1 2 6 3 . 17 As seen in these equations, the contribution from the fluctuations in Gaussian model of energy deposition correct the LQ biological parameters and to scale linearly and quadratically in and . We note that the radiobiological models that start from the equations equivalent of 16 and 17 are implicitly assuming a Gaussian-type symmetry in the energy-loss processes. This includes models with linear dependence on LET in and constant (e.g., by neglecting higher order terms in and beyond quadratic terms). Such models neglect the importance of asymmetries hidden in more realistic distribution functions such as the Landau [45] and/or Vavilov [46] distribution functions, responsible for observed nonlinearities in biological responses.

Neyman's distribution of type A and DSB distribution
We devote this section to illustrate the effect of deviation from Poisson distribution on statistical fluctuations of specific energy, DSB distribution, and biological indices and . Specifically we start with construction of the Neyman's distribution function from first principles for DSB induction processes and calculate the higher order moments in specific energy needed for incorporation of the repair and misrepair mechanisms to fit globally the cell survival data.
We consider the normalized distribution function to describe stochastic energy deposition in DNA material [58,59]  where ̅ ̅ / ! exp ̅ describes the Poisson distributed events in an ensemble of single tracks. Here ̅ denotes the average number of energy deposition events and is the distribution of specific energy imparted from passage of a single track within and resulted from exactly energy deposition events. The stochastic process as such is sketched in Fig. (1A).
We denote the deposited energy resulted from exactly events in mass of DNA material corresponding to specific energy ∑ / where ∑ . The occurrence of DSBs resulted from energy deposition requires balance in energy transfer to DNA. More specifically ∑ ̅ where ̅ ∑ is the typical energy results in breaking chemical bonds for induction of a single DSB and Δ 0,1,2, … counts number of DSBs. Δ is the number of induced DSBs in an event. Hence the energy balance in an exactly events processes requires ̅ Δ/ . Further simplification can be performed by averaging over DSB population per event subsequently by averaging over events that yields Δ ̅ Δ . Here Δ is average number of DSBs per event, independent of . The double bars over denotes two independent averagings; hence the order of averaging is not an issue. Also note that and / ̅ in Eq.(A1) hence Δ / ̅ . From Eq. (A20), it is now straightforward to calculate the statistical moments of DSBs and the corrections to biological indices and by inclusion of DSB fluctuations, as discussed in the main text. Here we systematically show that the DSB partition function given by Δ ; ̅ in Eq. (A20) provides all statistical moments we need for this analysis. For clarification of notations we first check the self-consistency of equations by calculating the first moment Δ ; ̅ ̅ Δ . Here are the coefficients of expansion and 1. By further expansion of Eq. (A21) we find a power law series dependence of DSB fluctuations on Δ that we need for the expansion of Eqs. After integrating over -function we find