Multiscale modeling for cancer radiotherapies

Abstract

Background

Ion-beam cancer therapy, an alternative to a common radiation therapy with X-rays, has been used clinically around the world since 1990s; the number of proton therapy centers as well as facilities using heavier ions such as α-particles and carbon ions continues to grow. A number of different methods were used by various scientific communities in order to quantitatively predict therapeutic effects of application of ion beams. A Multiscale approach (MSA) reviewed in this paper is one of these methods. Its name reflects the fact that the scenario of radiation damage following the incidence of an ion beam on tissue includes large ranges of scales in time, space, and energy.

Review

This review demonstrates the motivation and scientific justification of the MSA to the physics of ion-beam therapy and its implementation to a variety of different limits and physical conditions. A number of examples of calculations at high and low values of linear energy transfer (LET), large and small ion fluences, for a single value of LET and a combination of LETs in a spread-out Bragg peak are presented. The MSA has integrated the science involved in ion-beam therapy; in the process of the development of MSA, a new physical effect of ion-induced shock waves has been predicted. Its effect on the scenario of radiation damage is discussed in detail.

Conclusions

Multiscale approach’s predictive capabilities are based on the fundamental scientific knowledge. Their strength is in relation to actual physical, chemical, and biological processes that take place following the ions incidence on tissue. This makes the approach flexible and versatile to include various conditions, such as the degree of aeration or the presence of sensitizing nanoparticles, related to particular cases. The ideas for how the MSA can contribute to an improved optimization of therapy planning are summarized in the review.

Background: multiscale scenario of radiation damage with ions

It has been 10 years since the first paper on the Multiscale approach (MSA) to ion-beam therapy (Solov’yov et al. 2009) has been published. That paper has manifested the start of the development of a phenomenon-based approach to the assessment of radiation damage with ions, fundamentally different from other methods. The first goal was to understand the scenario of radiation damage with ions in the language of physical, chemical, and biological effects, that is, to relate initial physical effects of energy loss by projectiles to the biological effects defining cell inactivation. Thus, from the very beginning, the MSA is non-dosimetric, i.e., no damage is solely defined by the locally deposited dose.

The scenario is taking place on a number of time, space, and energy scales (so its name is perfectly justified) and features physical, chemical, and biological effects. We start with its overview and then show how its understanding can be practical for calculations of a number of important quantities. The scenario starts with the propagation of ions in tissue that is in most works replaced with liquid water (since it constitutes about 75% of tissue). This propagation is dominated by ionization of molecules of the medium by incident ions and features a Bragg peak in the depth–dose curve. The location of the Bragg peak depends on the initial energy of ions. In therapeutic applications, the initial energy of ions can be manipulated so that the Bragg peak falls into the location of the tumor. The location and shape of the Bragg peak as a function of initial energy were obtained analytically (Surdutovich et al. 2009; Scifoni et al. 2010; Surdutovich and Solov’yov 2014; Solov’yov 2017) based on the singly differentiated cross sections of ionization of water molecules with ions. Even though the depth–dose curve has been obtained and adopted for treatment planning (e.g., with Monte Carlo (MC) simulations Pshenichnov et al. 2008), a successful comparison of the depth–dose curve based on the singly differentiated cross section of ionization of molecules of the medium with ions has validated our approach on the early stage.

Further analysis of singly differentiated cross sections of ionization (Scifoni et al. 2010; de Vera et al. 2013) gave us a vital understanding of the energy spectrum of secondary electrons, ejected as a result of ionization $10^{-18}{-}10^{-17}$ s after ion’s passage. It was understood that most secondary electrons ejected from molecules in the medium by collisions with ions having energies below 50 eV. More energetic δ-electrons are kinematically suppressed in the Bragg peak and remain relatively rare in the plateau region preceding the peak. At energies of about 50 eV, electrons can be treated classically (as ballistic particles) and cross sections of their interactions with molecules of the medium are nearly isotropic (Nikjoo et al. 2006). This justified the use of the random walk approximation (i.e., diffusion mechanism) to describe their transport, and this was successfully accomplished in Solov’yov et al. (2009), Surdutovich and Solov’yov (2012, 2014, 2015) and Bug et al. (2010).

There are several consequences of features of secondary electron transport that fundamentally affect the scenario of radiation damage. First, the electrons lose most of their energy within 1–1.5 nm of the ion’s path; this happens within 50 fs of the ion’s passage through the medium (Surdutovich and Solov’yov 2015). Radiation damage, such as single- and double-strand breaks (SSBs and DSBs) in nuclear DNA can result from this energy loss (inelastic collisions of secondary electrons with DNA); these lesions can also result from interactions of slower electrons via dissociative attachment. In any case, these processes happen within 3–5 nm of the ions’ path. Second, the average energy of secondary electrons only weakly depends on the energy of projectiles and in the Bragg peak is independent of the linear energy transfer (LET) of projectiles. Most of these electrons are capable of ionizing one more molecule of the medium (Surdutovich et al. 2009). Therefore, the number of secondary electrons is roughly proportional to the LET.

Third, since most of the energy lost by secondary electrons within 50 fs stay within 1–1.5 nm of the ion’s path (the so-called “hot” cylinder) and there are no means of transport of this energy (since heat conductivity and diffusion take place slowly on the ps scale), the pressure (proportional to LET) developing within the hot cylinder during the 50–1000-fs period is expected to cause a significant collective flow associated with a shock wave, provided that the LET is sufficiently large. Ion-induced shock waves predicted by the MSA have been investigated in a series of works, both analytically and computationally (Surdutovich and Solov’yov 2010, 2014; Surdutovich et al. 2013; Yakubovich et al. 2012, 2011; de Vera et al. 2016, 2017, 2018).

Fourth, multiple reactive species are formed from the molecules ionized either by primary projectiles or secondary electrons. Their effect on DNA is deemed to be more important than the direct effect of secondary electrons; therefore, the understanding of their production and transport is vital for radiation damage assessment. The reactive species are formed within 1–2 ps of the ion’s passage, and their number densities may be large, to the first approximation linear with LET. However, their rates of recombination are proportional to the square of their number densities and at large values of LET, the recombination may dominate the transport by diffusion so the number of species that diffuse out of ion tracks is suppressed. To the contrary, a strong collective flow due to an ion-induced shock wave can propagate reactive species before they could recombine, thus changing initial conditions for the chemical phase (Surdutovich and Solov’yov 2014, 2015; de Vera et al. 2018).

The above consequences are substantial, and they constitute the physical part of the MSA. The analytical method based on them yields an opportunity to assess chemical effects and suggests a biological model for cell inactivation. Next, the concept of lethal DNA lesion is to be defined; after that the number of such lesions per unit length of ion’s path is calculated, and the cell survival probability is obtained. The concept of lethal damage in MSA is based on two hypotheses: (i) the inactivation of cells irradiated with ions is due to nuclear DNA damage, and (ii) a DNA lesion of a certain complexity is lethal. The second hypothesis comes from a series of works (Ward 1988, 1995; Malyarchuk et al. 2008, 2009; Sage and Harrison 2011) spanning three decades. Following these hypotheses, simple DNA lesions (such as SSB or base damage), DSBs, and complex lesions consisting of several simple lesions in addition to a DSB were considered as potentially lethal. After a series of investigations, it was postulated that complex lesions consisting of a DSB and at least two more simple lesions within a length of two DNA twists are lethal, at least for a normal cell (Surdutovich and Solov’yov 2014; Verkhovtsev et al. 2016). This is the so-called criterion for lethality, which implicitly includes the probability of enzymatic repair of DNA. This criterion may be modified for different cancerous cells and some special cell lines (Verkhovtsev et al. 2016). What is even more important (distinguishing the MSA from other approaches) is that each lesion has been associated with an action of an agent such as primary particle, secondary electron, or a reactive species. An action here means a probability that a single hit will cause a lesion, not necessarily related to a particular energy deposition. This is a significant difference from nano- and microdosimetric approaches.

After the criterion of lethality is defined, the fluence of agents on a given DNA segment (located at a distance from an ion’s path) is calculated in accordance with the transport mechanism (taking into account collective flows due to ion-induced shock waves). These fluences are weighted with probabilities of chemical processes leading to lesions. After that, the yield of lethal lesions per unit of an ion’s path length is calculated using Poisson statistics (Surdutovich and Solov’yov 2014; Verkhovtsev et al. 2016). Three quantities: ion fluence, LET, and dose deposited in the cell nucleus are related. Treating two of them as independent, e.g., the LET and the dose, the average length of all tracks through the nucleus can be calculated. Then, the product of this length and the yield of lethal lesions per unit length of the ion’s path gives the yield of lethal lesions per cell. This yield depends on the dose, LET, and the oxygen concentration in the medium. Thus, the survival curves are calculated and the relative biological effectiveness (RBE) can be calculated as well. In Verkhovtsev et al. (2016), the calculated survival curves were successfully compared with those experimentally obtained for a number of cell lines.

There have been five years since the last major review of MSA was published (Surdutovich and Solov’yov 2014). The MSA has been by and large completed in a sense of its original goal. The current review shows how the above approach was applied in different conditions and demonstrates its versatility. Different effects are discussed in relation to their influence on the shape of cell survival probability curves. It is hard to proceed without showing a figure (Fig. 1) that combines the scenario of radiation damage with ions (Surdutovich and Solov’yov 2014; Solov’yov 2017); it shows several possible pathways leading from the ion that loses energy to the cell apoptosis.

Fig. 1

Scenario of biological damage with ions. Ion propagation ends with a Bragg peak, shown in the top right corner. A segment of the track at the Bragg peak is shown in more detail. Secondary electrons and radicals propagate away from the ion’s path damaging biomolecules (central circle). They transfer the energy to the medium within the hot cylinder. This results in the rapid temperature and pressure increase inside this cylinder. The shock wave (shown in the expanding cylinder) due to this pressure increase may damage biomolecules by stress (left circle), but it also effectively propagates reactive species, such as radicals and solvated electrons to larger distances (right circle). A living cell responds to all shown DNA damage by creating foci (visible in the stained cells), in which enzymes attempt to repair the induced lesions. If these efforts are unsuccessful, the cell dies; an apoptotic cell is shown in the lower right corner

In "Calculation of the average number of lethal lesions produced by ions traversing cell nuclei" and "Calculation of lesion yields and survival curves" sections, we show and discuss the major components of the MSA, the average number of lethal lesions per unit length of ion’s path through the cell nucleus and the cell survival probability, respectively. Then, in "Calculation of lesion yields and survival curves" section, we go over a number of applications of MSA to different conditions.

Calculation of the average number of lethal lesions produced by ions traversing cell nuclei

In this section, we go over the calculation of the most important quantity, the number of lethal lesions per unit length of ion’s path through the cell nucleus. This is the most physical component of the formula for calculation of cell survival probabilities, and it heavily relies on the correct understanding of the scenario of radiation damage with ions. Consequently, all future improvements in the method will be related to the modifications of this part of MSA.

This calculation is defined by the criterion of lethality that defines the number of simple lesions sufficient for the complex one to be lethal. In the previous works (Surdutovich and Solov’yov 2014; Verkhovtsev et al. 2016, 2019), the lethal lesion was defined by at least three simple lesions one of which is a DSB on the length of two DNA twists. The probability of formation of each lesion is a product of the probability of encounter of a target site with a secondary electron, reactive species, etc., and the probability of a lesion formation as a result of this collision. The cross sections of inelastic interactions different agents with DNA are being found in different experiments or calculated (Surdutovich and Solov’yov 2014; Boudaïffa et al. 2000; Huels et al. 2003; Nikjoo et al. 2002; Kumar and Sevilla 2010; Sevilla et al. 2016), and some average numbers for probabilities of DNA strand breaks and other lesions have been successfully used in Surdutovich and Solov’yov (2014) and Verkhovtsev et al. (2016, 2019). The new knowledge about these numbers can certainly be applied in the future, but these numbers may only depend on the local conditions such as oxygen density on the site. On the other side, the number of secondary particles hitting a target, or the fluence, depends on the radiation, i.e., on the LET, dose, etc., and the quantities that can be manipulated and optimized. Besides, the fluence depends on the mechanism of transport of secondary particles and this dependence is highly emphasized in the MSA.

As an ion traverses a cell nucleus, it ionizes molecules and ejected secondary electrons (first generation) start with the average energy of $\sim 45$ eV. These electrons lose most of this energy within 1–1.5 nm of the ion’s path, ionizing more molecules (including biomolecules). The second generation of slower electrons is thus formed. These secondary electrons can cause damage only within a region of a few nm. This damage can be estimated from the calculation of the average number of secondary electrons incident on uniformly distributed targets (DNA segments) in the region (Solov’yov et al. 2009; Surdutovich and Solov’yov 2014). A quantity ${{{\mathcal {N}}}}_e(r)$, the average number of simple lesions, on a target at a distance r from the ion’s path is calculated as a result. This is the secondary electrons’ contribution.

Most of the reactive species (free radicals and solvated electrons, $e^-_{\text {aq}}$) are formed at locations of ionizations described above (Surdutovich and Solov’yov 2015). If the LET is relatively small, the number of reactive species is small as well, and their interaction can be neglected. Then, they very slowly (compared to secondary electrons) diffuse, reacting with DNA targets on their way. A quantity $\mathcal{N}_r(r)$, the average number of simple lesions due to reactive species on a target at a distance r from the ion’s path is calculated as a result. This is the reactive species contribution at “low-LET”.

At a high LET, the reactive species are produced in large quantities and given an opportunity they would interact much faster than they diffuse and this would lead to their recombination (Surdutovich and Solov’yov 2015). However, at high values of LET there is another mechanism for transport of radicals, i.e., the collective flow due to ion-induced shock waves. The shock waves initiated by a large pressure difference and propagating radially from each ion’s path were predicted in Surdutovich and Solov’yov (2010) and discussed in a number of works within the MSA (Surdutovich and Solov’yov 2014; Verkhovtsev et al. 2016; Surdutovich et al. 2013, 2017; Yakubovich et al. 2011; de Vera et al. 2016, 2017, 2018); the transport of radicals with a collective flow including chemical reactions was investigated by means of molecular dynamics (MD) simulations in de Vera et al. (2018). As a result, the effective ranges of the reactive species, such as hydroxyl radicals and solvated electrons, are substantially larger than those consistent with the diffusion transport mechanism. The evidence of such large ranges, inferred from the observation of the interaction of the ion tracks at large ion fluences, can be a strong argument in favor of the existence of collective flow.

Calculation of number of secondary electrons incident on a DNA target

As shown in Surdutovich and Solov’yov (2015), the number densities of the first and second generations of secondary electrons are given by,

$$\begin{aligned} n_1(t, r)& = \frac{\text {d}N_1}{\text {d}x}\frac{1}{4\pi D_1 t}\exp \left( -\frac{r^2}{4 D_1 t}-\frac{t}{\tau _1} \right), \\ n_2(t,r)& = \frac{2}{4\pi \tau _1} \frac{\text {d}N_1}{\text {d}x}\int _0^t \frac{1}{D_1 t'+D_2(t-t')}\mathrm{e}^{-\frac{r^2}{4( D_1 t'+D_2(t-t'))}-\frac{t-t'}{\tau _{2}}-\frac{t'}{\tau _1}}\text {d}t', \end{aligned}$$

(1)

where $\frac{\text {d}N_{1}}{\text {d}x}$ is the number of ionizations taking place per unit length in the longitudinal direction, x, of the ion’s trajectory, $D_1$ and $D_2$ are the diffusion coefficients, and $\tau _1$ and $\tau _2$ are the average lifetimes of electrons of the first and second generations, respectively. Since the characteristic spatial scale in the radial direction is in nanometers and in the axial direction is micrometers, $\frac{\text {d}N_{1}}{\text {d}x}$ is assumed to be constant along the length of the target.

A target is chosen to be a rectangle of area $\xi \eta$, where $\xi =6.8$ nm and $\eta =2.3$ nm are the length of two twists and the diameter of a DNA molecule, respectively. Thus, electrons or radicals hitting such a target would be hitting two rungs of a DNA molecule masked by this target. The plane of the target is chosen to be parallel to the ion’s path with dimension $\xi$ along and $\eta$ perpendicular to the path. This can be seen in Fig. 2. Then, angle $\phi =2\arctan \frac{\eta /2}{r}$ inscribes the target in a plane perpendicular to the ion’s path, where r is the distance between the target and the path.

Fig. 2

Geometry of the problem in the plane perpendicular to ion’s path. The target cylinder that encloses a DNA twist is shown as a circle. Its diameter is $\eta$. The dimension $\xi$ is perpendicular to the plane of the figure

The number of first-generation electrons hitting the described target segment of area $r \phi \xi \approx \xi \eta$ parallel to the ion’s path per unit time is

$$\begin{aligned} \Phi _{1}=- \phi \xi r D_1 \frac{\partial n_1(r,t)}{\partial r} = \frac{\phi }{2\pi }\frac{\text {d}N_{1}}{\text {d}x} \frac{r^2 \xi }{4 D_1 t^2}\exp \left( -\frac{r^2}{4 D_1 t}-\frac{t}{\tau _1}\right) . \end{aligned}$$

(2)

Its integral over time,

$$\begin{aligned} \int _0^\infty \Phi _1 \,{\rm d}t &= \frac{\phi }{2\pi }\int _0^\infty \frac{\text {d}N_{1}}{\text {d}x} \frac{r^2 \xi }{4 D_1 t^2}\exp \left( -\frac{r^2}{4 D_1 t}-\frac{t}{\tau _1}\right) \text {d}t \\ &=\frac{\phi }{2\pi }\frac{\text {d}N_{1}}{\text {d}x}\frac{r \xi }{\sqrt{D_1 \tau _1}}K_1\left( \frac{r}{\sqrt{ D_1\tau _1}}\right) ,~~ \end{aligned}$$

(3)

where $K_1$ is the Macdonald function (modified Bessel function of the second kind) (Abramowitz 1972), gives the total number of first-generation secondary electrons that hit this area. The second-generation contribution is obtained similarly:

$$\begin{aligned} \Phi _2(t,r)& = -r \phi \xi D_2\frac{\partial n_2(r,t)}{\partial r} \\&=\, {} \phi \frac{\xi r^2 D_2}{4\pi \tau _1} \frac{\text {d}N_1}{\text {d}x}\int _0^t \frac{1}{(D_1 t'+D_2(t-t'))^2} \\&\quad \times \exp \left( -\frac{r^2}{4( D_1 t'+D_2(t-t'))}-\frac{t-t'}{\tau _{2}}-\frac{t'}{\tau _1}\right) \text {d}t', \end{aligned}$$

(4)

and then,

$$\begin{aligned} \int _0^\infty \Phi _2 \,{\rm d}t &= \phi \frac{\xi r^2 D_2}{4\pi \tau _1} \frac{\text {d}N_1}{\text {d}x}\int _0^\infty \int _0^t \frac{1}{(D_1 t'+D_2(t-t'))^2} \\&\quad \times \exp \left( -\frac{r^2}{4( D_1 t'+D_2(t-t'))}-\frac{t-t'}{\tau _{2}}-\frac{t'}{\tau _1}\right) \text {d}t'\text {d}t \end{aligned}$$

(5)

gives the number of second-generation secondary electrons that hit the same area. The average number of simple lesions due to a single ion, ${{{\mathcal {N}}}}_e(r)$, can now be obtained as the sum,

$$\begin{aligned} {{{\mathcal {N}}}}_e(r)={{{\mathcal {N}}}}_1(r)+{{{\mathcal {N}}}}_2(r)=\Gamma _e\int _0^\infty \Phi _1 \,{\rm d}t+\Gamma _e\int _0^\infty \Phi _2 \text {d}t, \end{aligned}$$

(6)

where ${{{\mathcal {N}}}}_1(r)$ and ${{{\mathcal {N}}}}_2(r)$ are the average numbers of simple lesions produces by secondary electrons of the first and second generations, respectively, and $\Gamma _e$ is the probability for an electron to induce a simple lesion on a hit. The dependencies of ${{{\mathcal {N}}}}_1(r)$ and ${{{\mathcal {N}}}}_2(r)$ are shown in Fig. 3.

Fig. 3

Average numbers of simple lesions due to a single carbon ion with a Bragg peak energy propagating through a uniform chromatin as functions of radial distance from the ion’s path. The lesions are produced by secondary electrons of the first (solid line) and second (dashed line) generations, ${{{\mathcal {N}}}}_1(r)$ and ${{{\mathcal {N}}}}_2(r)$. These dependencies are calculated using the corresponding number of hits, Eqs. (3) and (5), multiplied by the probability of the production of a simple lesion per hit, $\Gamma _e=0.03$ (used in Surdutovich and Solov’yov 2014). A straight (dotted) line is the values for reactive species, ${{{\mathcal {N}}}}_r(r)$, calculated using Eq. (13) with numbers from Verkhovtsev et al. (2016)

Equation (6) gives the average number of simple DNA lesions due to secondary electrons of the first and second generations as a function of the distance of the target DNA segment from the ion’s path. The next step is adding to this the contribution of reactive species, which is a product of the average number of hits on the chosen area by reactive species $N_r$ (this value is similar to $\int _0^\infty \Phi _1 \text {d}t$ and $\int _0^\infty \Phi _2 \text {d}t$) and the probability of lesion production per hit, $\Gamma _r$. The number $N_r$ depends on the value of the LET, since at small values of LET the transport of radicals is defined by diffusion and at high values, the collective flow is expected to dominate this process.

Calculation of the reactive species contribution for small values of LET

The number of produced reactive species, such as free radicals and solvated electrons, depends on the LET. If the LET is not very high, it is expected that the number of reactive species is proportional to the secondary electron production, $\text {d}N_1/\text {d}x$, and, therefore, increases nearly linearly with the value of LET (Surdutovich and Solov’yov 2014). At sufficiently high values of LET, extra production of radicals is possible due to water radiolysis at locations adjacent to the ion’s path. This effect has not yet been quantified and will be accounted for in future works along with the definition of the domain of the LET, where this effect becomes significant. In this work, a linear dependence between the number of reactive species and LET is assumed and the difference between high and low values of LET is defined only by the mechanism of transport of the reactive species; at low LET, this transport is defined by diffusion. Moreover, this means that chemical reactions such as $2\mathrm{OH}\rightarrow \mathrm{H}_2\mathrm{O}_2$ and $e^-_{\text {aq}}+\mathrm{OH}\rightarrow \mathrm{OH}^-$ are rare and their frequency can be neglected compared to the diffusion term in the diffusion equation (Surdutovich and Solov’yov 2015). Thus, the transport of reactive species in the low LET case can be calculated, by solving a diffusion equation,

$$\begin{aligned} \frac{\partial n_r}{\partial t}=D_r \nabla ^2 n_r, \end{aligned}$$

(7)

where $n_r$ is the number density and $D_r$ is the diffusion coefficient for reactive species.

The initial conditions for this equation can be taken from Surdutovich and Solov’yov (2015),

$$\begin{aligned} \frac{\partial n_{r}(r,t)}{\partial t}=\frac{\text {d}N_1}{\text {d} x}\delta ^{(2)}(r)\delta (t)+\frac{n_{1}({r}, t)}{\tau _{1}}+\frac{n_{2}({r}, t)}{\tau _{2}}, \end{aligned}$$

(8)

where the first term describes the species formed at sites of original ionizations by the projectile, while other two terms are due to inelastic processes involving secondary electrons of the first and secondary generations, respectively. Ionizations and excitations that lead to the production of reactive species, $n_{r}(r,t)$, through the mechanism of Eq. (8) take place by about 50 fs (Surdutovich and Solov’yov 2015). By that time, the forming reactive species are localized within 3 nm of the ion’s path. These are the initial conditions for the following propagation of reactive species by the diffusion and/or collective flow that happen on much larger scales, up to 100 ps in time and 50 nm in distance. Therefore, in this paper, a simplified initial condition is used,

$$\begin{aligned} \frac{\partial n_{r}(r,t)}{\partial t}=K\frac{\text {d}N_1}{\text {d} x}\delta ^{(2)}(r)\delta (t), \end{aligned}$$

(9)

where K is the number of reactive species produced due to each secondary electron of the first generation ejected by an ion. The value of $K \approx 6$ can be evaluated as follows. The primary ionization produces $\hbox {H}_{2}\hbox {O}^{+}$, which is likely to produce a hydroxyl radical (von Sonntag 1987). The same thing happens when the secondary electron of the first generation ionizes a water molecule (and thus becomes an electron of the second generation). Then, two electrons of the second generation (the ionizing and ejected) can produce about four reactive species, two as a result of further energy loss in inelastic processes and two more if they become solvated electrons. A more accurate number for K can be obtained if the probabilities of the above processes are combined following a comprehensive radiochemical analysis.

The solution to Eq. (7) with the initial condition (9) is given by,

$$\begin{aligned} n_r(r,t) = K \frac{\text {d} N_1}{\text {d} x}\frac{1}{4\pi D_r t}\exp {\left( -\frac{r^2}{4D_r t}\right) }~. \end{aligned}$$

(10)

The next step is to find the number of reactive species, $\Phi _r$, incident on the target at a distance r from the ion’s path per unit time. We proceed similarly to Eqs. (2) and (3).

$$\begin{aligned} \Phi _{r}=- \phi \xi r D_r \frac{\partial n_r(r,t)}{\partial r} = \frac{\phi }{2\pi }K \frac{\text {d}N_{1}}{\text {d}x} \frac{r^2 \xi }{4 D_r t^2}\exp \left( -\frac{r^2}{4 D_r t}\right) , \end{aligned}$$

(11)

and its integral over time is simply,

$$\begin{aligned} \int _0^\infty \Phi _r \text {d}t=\frac{\phi \xi }{2\pi }K \frac{\text {d}N_{1}}{\text {d}x}=K \frac{\text {d}N_{1}}{\text {d}x}\frac{\xi }{\pi }\arctan \frac{\eta /2}{r}~. \end{aligned}$$

(12)

Calculation of the reactive species contribution for large values of LET

If the reactive species are formed in large quantities as a result of a high-LET-ion’s traverse, the collective flow due to the shock wave is the main instrument for the transport of these species away from the ion’s path. Interestingly, ranges of propagation of radicals used to be in the realm of chemistry (von Sonntag 1987; LaVerne 1989; Alpen 1998). However, in the case of high LET, this issue is addressed by physicists; the MD simulation (with a use of MBN Explorer package (Solov’yov et al. 2012; Sushko et al. 2016a) showed that the range depends on the value of LET (de Vera et al. 2018), but a more extensive investigation is needed to obtain a more detailed dependence.

In Verkhovtsev et al. (2016), a simple model was used to describe this transport. The value of the average number of lesions at a distance r from the ion’s path, ${{{\mathcal {N}}}}_r=\Gamma _r N_r$, was considered to be a constant within a certain LET-dependent range R, i.e.,

$$\begin{aligned} {{{\mathcal {N}}}}_r(r)={{{\mathcal {N}}}}_{r}\Theta (R-r), \end{aligned}$$

(13)

where $\Theta$ is the Heaviside step function. The value $\mathcal{N}_r$ also depends on the degree of oxygenation of the medium, since the concentration of oxygen dissolved in the medium affects the number of formed radicals as well as the effectiveness of lesion repair. In principle, more information about ${{{\mathcal {N}}}}_r$ is needed. For example, at high LET, more reactive species are expected to be produced through the radiolysis of water in the cores of the ion tracks at times $\ge ~50$ fs after the energy transfer from secondary electrons to the medium has taken place. This process can now be studied by MD simulations using the MBN Explorer package (Solov’yov et al. 2012; Sushko et al. 2016b), which is capable of resolving the corresponding temporal and spatial scales.

The comprehensive picture of transport of reactive species includes diffusion (dominant at low values of LET), collective flow (dominant at high values of LET), and chemical reactions. With this understanding, as LET increasing Eq. (12) should gradually transform into Eq. (13). In addition to these equations, the effective range of reactive species is limited by the criterion of lethality that requires a minimal fluence at each site. More discussion on this topic can be found in Verkhovtsev et al. (2019).

Calculation of lesion yields and survival curves

Within the MSA, the probability of lesions is calculated using Poisson statistics and the next step is the calculation of the average number of simple lesions, ${{{\mathcal {N}}}}$.

$$\begin{aligned} {{{\mathcal {N}}}}={{{\mathcal {N}}}}_e(r)+{{{\mathcal {N}}}}_r(r)~. \end{aligned}$$

(14)

Based on this, the probability of lethal damage according to the criterion of lethality is (Surdutovich and Solov’yov 2014; Verkhovtsev et al. 2016),

$$\begin{aligned} P_l(r)=\lambda \sum _{\nu =3}^\infty \frac{{{{\mathcal {N}}}}^\nu }{\nu !}\exp {\left[ -{{{\mathcal {N}}}}\right] }, \end{aligned}$$

(15)

where $\lambda =0.15$. This criterion states that three DNA lesions, one of which is a double-strand break, have to occur within two DNA twists. The probability given by Eq. (15) is then integrated over space ($2\pi r\text {d}r$) giving the number of lethal lesions per unit segment of the ion’s path, $\text {d}N_{l}/\text {d}x$,

$$\begin{aligned} \frac{\text {d} N_{l}}{\text {d} x}=2\pi n_s \int _0^\infty P_l(r)r \text {d}r, \end{aligned}$$

(16)

where $n_s$ is the target density calculated as in Verkhovtsev et al. (2016).

At this point, the results of the previous sections can be combined in the expression for the yield of lethal lesions. Such an expression was obtained in Surdutovich and Solov’yov (2014) and Verkhovtsev et al. (2016) for the case of non-interfering ion paths as,

$$\begin{aligned} Y_{l} = \frac{{\text {d}}N_{l}}{{\text {d}}x} \, {\bar{z}} \, N_{\text {ion}}(d), \end{aligned}$$

(17)

where $N_{\text {ion}}$ is the number of ions that traverse a target and ${\bar{z}}$ is the average length of the trajectory of the ion’s traverse. This yield is a product of the yield per unit length of the ion’s path and the average length within a target passed by all ions (${{{\bar{z}}}}N_{\text {ion}}$).

Equation 17 gives the number of lethal damage sites per cell nucleus, and therefore, according to the Poisson statistics, the probability of cell deactivation is,

$$\begin{aligned} \Pi _{d} = 1-\exp {(-Y_l)}, \end{aligned}$$

(18)

i.e., unity less the probability of zero lethal lesions. Then, the probability of cell’s survival is $\Pi _{\text {surv}}=1-\Pi _{d} = \exp {(-Y_l)}$, which is usually written as

$$\begin{aligned} -\ln {\Pi _{\text {surv}}}=Y_{l}, \end{aligned}$$

(19)

i.e., the natural logarithm of the inverse cell survival probability is equal to the yield of lethal lesions in the nuclear DNA. This expression has been used since Surdutovich and Solov’yov (2014). The yield given by Eq. (17) was used in a number of applications (Verkhovtsev et al. 2016, 2019; Surdutovich and Solov’yov 2017, 2018). It can be rewritten in several ways,

$$\begin{aligned} Y_{l} = \frac{{\text {d}}N_{l}}{{\text {d}}x} \, {\bar{z}} \, N_{\text { ion}}(d)=\frac{\pi }{16}N_g\frac{\sigma (S_e)}{S_e}d=\frac{\pi }{16}N_g\sigma (S_e)F_{\text {ion}}, \end{aligned}$$

(20)

where $F_{\text {ion}}$ is the ion fluence. Now, we want to dwell on the universality and versatility of this expression. Its first representation, $\frac{{\text {d}}N_{l}}{{\text {d}}x} \, {\bar{z}} \, N_{\text {ion}}(d)$, indicates that the yield is just a product of two quantities, the number of lethal lesions per unit length of ion’s trajectory and the total length of ions path through the cell nucleus, which can be broken into average length of a traverse by the number of ions passing through the nucleus. This number depends on the dose. However, the dose in the case of ions is not an independent parameter, it is regulated by the ions fluence:

$$\begin{aligned} d=\frac{S_e {\bar{z}} \, N_{\text {ion}}}{\rho V}=\frac{S_e \, F_\text{ion}}{\rho }, \end{aligned}$$

(21)

where $\rho$ is the mass density of the nucleus. Before we analyze a number of effects and limits, we want to acknowledge a successful comparison of calculated survival curves at a range of LET values for a number of different cell lines shown in Fig. 4 (Verkhovtsev et al. 2016).

Fig. 4

Survival curves for different human cell lines: adenocarcinomic A549 cells (a), normal fibroblasts AG1522 (b), cervical cancer HeLa cells (c), normal skin fibroblasts NB1RGB (d), glioblastoma A172 cell line (e), and endothelial EAhy926 cells (f). The calculated survival probabilities are shown with lines and experimental data from Wéra et al. (2011, 2013) (A549), Raju (1991), Autsavapromporn (2011), Hamada (2006) (AG1522), Zhao (2013), Usami (2016) (HeLa), Tsuruoka (2005), Suzuki (2000) (NB1RGB), Suzuki (2000), Tsuboi (1998) (A172), and Riquier (2013) (EAhy926) are shown with symbols

First, it is interesting to analyze the limits of $N_{\text {ion}}$; the minimum (nonzero) value for it is one. Then, both the dose and yield are defined by LET, at that the former is linear with it and the latter is linear if the LET is small, but may be quadratic if the LET is larger. This enhancement is expected as a result of the transport of reactive species to larger distances by ion-induced shock waves. If the LET is too large, the lethal damage may happen already on a fraction of ${\bar{z}}$. This means that the “rest” of the dose is wasted, the relative biological effectiveness is reduced, and the so-called overkill effect is observed. On the other side, when $N_{\text {ion}}$ is very large, ion tracks are likely to overlap. This corresponds to the case of large ion fluences, which was discussed in Surdutovich and Solov’yov (2018). This limit may be important in the case of applications of laser-driven proton beams.

Second, Eq. (21) is only valid when the LET is the same for all ions; when it is not, e.g., in the case of a spread-out Bragg peak, then

$$\begin{aligned} d=\sum _j\frac{S_{ej} {\bar{z}} \, F_j}{\rho }, \end{aligned}$$

(22)

where a subscript j indicates a corresponding component of the ion beam. This dependence was exploited in Surdutovich and Solov’yov (2017), and it will be discussed below because the spread-our Bragg peak (SOBP) is used clinically and in many experiments as well.

Third, more intriguing effects are seen in the second representation, $\frac{\pi }{16}N_g\frac{\sigma (S_e)}{S_e}d$. As was mentioned, at relatively small values of LET, the $\frac{{\text {d}}N_{l}}{{\text {d}}x}$ is linear with LET, i.e., $\sigma (S_e)=\xi _1 S_e$, where $\xi _1$ is a coefficient. Then, the yield is linear with the dose. However, when LET increases, the nonlinearity of yield of dose dependence comes from an expected quadratic dependence (Verkhovtsev et al. 2019) $\sigma (S_e)=\xi _2 S_e^2$ due to the ion-induced shock wave effect. Fourth, $N_g$ in this representation is the number of base pairs in the whole cell nucleus, which gets in this formula from the expectation that the cell is in the interphase and chromatin is uniformly distributed over the nucleus. In particular, this means that the yield for all human cells would be the same. As this may be true for healthy cells of normal tissue, this may not be true for cancerous cells. More research is needed to clarify this point.

Fifth, the oxygen concentration dependence is “hidden” in the value of $\frac{{\text {d}}N_{l}}{{\text {d}}x}$. It affects the reactive species effect through the value of ${{{\mathcal {N}}}}_{r}$ that enters Eq. (14). The map of oxygen concentration automatically produces the map of the oxygen enhancement ratio (OER), which is the ratio of doses required to achieve the same biological effect with a given oxygen concentration to that with the maximum oxygen concentration. The comparison of OER calculated using the MSA to that measured experimentally is shown in Fig. 5 (Verkhovtsev et al. 2016). The map of OER is deemed to be an important component of therapy optimization.

Fig. 5

Oxygen enhancement ratio at the 10% survival level for V79 and CHO cells irradiated with carbon ions. Symbols denote the experimental data taken from Tinganelli (2015), Furusawa (2000), Hirayama et al. (2005) and Chapman et al. (1977)

Sixth, if the LET is fixed, Eq. (20) suggests that the yield and therefore the logarithm of survival probability are linear with dose, thus making survival curves in their traditional coordinates straight lines. A comparison of a number of survival curves at a range of LET values shown in Fig. 4 supports this observation; however, there are experiments that the so-called shouldered survival curves are observed. At this point, it is worthwhile to remind a reader that a vast research of X-rays survival curves (Alpen 1998) that the straight survival curves indicate a single-hit scenario of radiation damage. This means that a single hit of a target (in our case with an ion) leads to cell inactivation with a given probability. This probability includes the probability of DNA damage repair. In the framework of molecular theories developed from 1950s to 1990s (Alpen 1998), including the microdosimetric kinetic model (MKM) (Hawkins 1996, 2009), the shouldered survival curves are the result of either nonlinear damage or repair. It is interesting to place the MSA on this map.

The criterion of lethality and Eq. (20) produce linear survival curves for cells irradiated with ions. This model includes the probability of enzymatic repair, embedded into the criterion. The criterion itself can be different for different cell lines, but it will lead to straight lines nevertheless. The “shoulderness through damage” translates into MSA language as tracks overlap. In this case, the $\frac{{\text {d}}N_{l}}{{\text {d}}x}$ depends on fluence and therefore or dose and Eq. (20) becomes nonlinear with dose and predicts a shouldered survival curve (Surdutovich and Solov’yov 2018). However, this happens at very large values of fluence and dose, far larger than those used clinically. Therefore, it is more likely that a shouldered curve in ion therapy may be due to repair process. The solution to this problem was suggested in Verkhovtsev et al. (2016) and it is as follows.

This solution does not change the expression for the yield given by Eq. (20), except for a constant coefficient. What changes is the logarithm of survival probability (19); instead of being linear with the yield, it becomes a quadratic function,

$$\begin{aligned} -\log \Pi _{\text {surv}} = Y_{l} - (\chi _0 - \chi _1 Y_l) Y_{l}= (1 - \chi _0) Y_{l} + \chi _1 Y_l^2, \end{aligned}$$

(23)

where $\chi _0$ and $\chi _1$ are positive constants. The first representation can be phenomenologically interpreted in such a way that the cell lines for which the survival curves are shouldered are more resistive than those for which the survival is linear, at small values of yield the r.h.s. is linear with respect to $Y_l$ with a coefficient $1-\chi _0<1$; however, as the yield increases the resistivity decreases linearly and when $(\chi _0 - \chi _1 Y_l)$ turns to zero, the survival becomes “normal”. This is formalized as,

$$\begin{aligned} - \ln {\Pi _{\text {surv}}}=\, & {} (1 - \chi ) Y_l =\, Y_l - \Theta (\chi _0 - \chi _1 Y_l) \, (\chi _0 - \chi _1 Y_l) \, Y_l, \\ \chi=\, & {} \left( \chi _0 - \chi _1 \, Y_l \right) \, \Theta (\chi _0 - \chi _1 \, Y_l)~. \end{aligned}$$

(24)

The coefficient $\chi$ gradually approaches zero with increasing number of lesions until it becomes equal to zero at a critical value, ${\tilde{Y}}_l = \chi _0/\chi _1$, which depends, in particular, on dose and LET. Above this critical value, Eq. (19) remains valid. Thus, the critical yield ${\tilde{Y}}_l$ is the transition point in the survival curve from the linear-quadratic to the linear regime. The examples of application of this model are shown in Fig. 6 (Verkhovtsev et al. 2016).

Fig. 6

Survival curves for a repair-efficient CHO cell line. The calculated survival probabilities are shown with lines and experimental data from Weyrather et al. (1999) and Usami (2008) are shown by symbols. The survival curves are calculated using Eq. (24) with $\chi _0 = 0.35$ and $\chi _1 = 0.04$

For $Y_l < \chi _0/\chi _1$, the survival probability given by Eq. (24) can be rewritten as,

$$\begin{aligned} - \ln {\Pi _{\text {surv}}} = (1 - \chi _0) \, \frac{\pi }{16} \sigma \, N_{g} \, \frac{d}{S_e} + \chi _1 \left( \frac{\pi }{16} \sigma \, N_{g} \right) ^2 \frac{d^2}{S_e^2}~. \end{aligned}$$

(25)

At this point, the famous empirical parameters $\alpha$ and $\beta$ of the linear-quadratic model (Alpen 1998) given by

$$\begin{aligned} - \ln {\Pi _{\text {surv}}} = \alpha d + \beta d^2, \end{aligned}$$

(26)

can be introduced. Equation (25) provides the molecular-level expressions for these parameters at doses $d \le \frac{16}{\pi } \frac{S_e}{\sigma N_g} \frac{\chi _0}{\chi _1}$:

$$\begin{aligned} \alpha = (1 - \chi _0) \, \frac{\pi }{16} \sigma \, N_{g} \, \frac{1}{S_e} \ , \qquad \qquad \beta = \chi _1 \, \left( \frac{\pi }{16} \sigma \, N_{g} \right) ^2 \frac{1}{S_e^2} \ . \end{aligned}$$

(27)

At $Y_l > \chi _0/\chi _1$, i.e., for $d > \frac{16}{\pi } \frac{S_e}{\sigma N_g} \frac{\chi _0}{\chi _1}$, survival curves are linear, and the parameter $\alpha$ is given by

$$\begin{aligned} \alpha = \frac{\pi }{16} \frac{\sigma \, N_{g}}{S_e} \ . \end{aligned}$$

(28)

Thus, the MSA methodology has been discussed. The main result is given by Eq. (20), which gives the expression for the yield of lethal lesions. This expression is obtained as a result of analysis of physical, chemical, and biological effects on the corresponding scales. Each of its components can be further refined, but its scientific clarity is sound. For instance, in recent years, the product of LET and dose, i.e., $S_e d,$ is used for proton therapy optimization (Underwood and Paganetti 2016). In the Bragg peak region, $\sigma (S_e)=\xi _2 S_e^2$ and this optimization parameter is a consequence of Eq. (20). While we are leaving the outlook of what has to be done along the MSA in the future to the Conclusion section, we get to some applications of MSA promised above.

Application of MSA at different limits of LET

Survival curves along a spread-out Bragg peak

The goal of Surdutovich and Solov’yov (2017) was to suggest an algorithm for choosing the energy distribution of ion fluence at the entrance in order to achieve the uniform cell survival distribution throughout the SOBP. In the beginning, it was shown that the uniform dose distribution leads to an increase in cell inactivation along the SOBP toward a sharp maximum at its distal end. In this review, we will just show the algorithm in order to achieve the uniform cell survival at a constant oxygen concentration along the SOBP.

Let the maximum initial energy at the entrance be $E_0$ and let it change by step $\Delta E$ to construct the SOBP; the depth of each pristine Bragg peak can be denoted by $x_j$, where $j=0,1,2,\ldots, J$. According to Eqs. (20) and (22), at a given depth x, the yield is

$$\begin{aligned} Y_{l}=\frac{\pi }{16}N_g\sum _j \sigma (S_j(x))F_j=Y_0, \end{aligned}$$

(29)

where $Y_0$ is the target yield throughout the SOBP. The goal is to obtain the distribution of $F_j$. Clearly,

$$\begin{aligned} F_0=Y_0 \frac{16}{\pi N_g \sigma (S_0(x_0))}, \end{aligned}$$

(30)

the fluence at maximum energy corresponds to the desired yield at the distal end of the Bragg peak. Then,

$$\begin{aligned} \frac{\pi }{16}N_g (\sigma (S_1(x_1))F_1+\sigma (S_0(x_1))F_0)=\frac{\pi }{16}N_g \sigma (S_0(x_0))F_0, \end{aligned}$$

(31)

which gives

$$\begin{aligned} F_1=\frac{\sigma (S_0(x_0))-\sigma (S_0(x_1))}{\sigma (S_1(x_1))}F_0, \end{aligned}$$

(32)

on the next step we find $F_2$ from

$$\begin{aligned} \frac{\pi }{16}N_g (\sigma (S_2(x_2))F_2+\sigma (S_1(x_2))F_1+\sigma (S_0(x_2))F_0)=\frac{\pi }{16}N_g \sigma (S_0(x_0))F_0, \end{aligned}$$

(33)

and so on. If the oxygen concentration depends on x, this affects all $S_j(x)$ and can be easily included in the algorithm. Figure 7 (Surdutovich and Solov’yov 2017) shows the application of the algorithm for a proton SOBP example.

Fig. 7

Solid line shows the profile of dependence of yield of lethal lesions in cells along the SOBP as a function of distance on the distal end of the SOBP. The dashed line shows the profile of the depth–dose curve that produced the above result

The overkill effect at large LET

In this section, we want to briefly discuss the limit of large values of LET, so large that $N_{\text {ion}}$ is close to one. In this limit, it is important that even though $N_{\text {ion}}$ in Eqs. (20) and (21) is an average number of ions traversing the cell nucleus, in reality the number of ions is integer. Therefore, $N_{\text {ion}}$ can be redefined as the minimum number of ions required to cause the damage reflected by the survival fraction of $\Pi _0$ and the corresponding yield $Y_0$. Then, (since $F_{\text {ion}}=N_{\text {ion}}/A_n$, where $A_n$ is a cross-sectional area of a cell nucleus) Eq. (20) can be solved for $N_{ion}$ as,

$$\begin{aligned} N_{\text {ion}} =\left[ \frac{16 Y_0 A_n}{\pi N_g \sigma (S_e)} \right] + 1, \end{aligned}$$

(34)

where square brackets denote the integer part of their content. The relative biological effectiveness (RBE) is given by the ratio of dose delivered by photons, $d_\gamma$ to that delivered by ions in order to achieve the same survival fraction or yield. Then, in virtue of Eq. (21),

$$\begin{aligned} {\text{RBE}}\,=\,\frac{d_\gamma }{d}=\frac{d_\gamma \rho V}{S_e \,{\bar{z}} \, N_{\text {ion}}} =\frac{d_\gamma \rho V}{S_e \,{\bar{z}} \, \left( \left[ \frac{16 Y_0 A_n}{\pi N_g \sigma (S_e)} \right] + 1\right) }~. \end{aligned}$$

(35)

This equation explains the overkill effect. When LET is small, the integer part in the numerator is large compared to unity. In this limit, RBE is given by

$$\text{RBE}\,=\,\frac{\pi N_g d_\gamma \rho }{16\, Y_0}\,\frac{\sigma (S_e)}{S_e}.$$

(36)

Since $\sigma (S_e)\propto S_e$ in this limit, RBE is independent of LET. Then, with increasing LET, $\sigma (S_e)\propto S_e^2$ and RBE becomes linear with LET until $\frac{16\, Y_0 \,A_n}{\pi N_g \sigma (S_e)}$ becomes close to unity. This is the limit of large LET, in which RBE becomes inversely proportional to LET,

$$\text{RBE}\,=\,\frac{d_\gamma \, \rho \, V}{S_e \,{\bar{z}}}.$$

(37)

This dependence is discussed in more detail in Verkhovtsev et al. (2019), and the dependence of RBE corresponding to Eq. 35 is shown in Fig. 8 (Verkhovtsev et al. 2019). A piecewise dependence at increasing values of LET corresponding to small values of $N_{\text{ion}}$ deserves a comment. Nothing is wrong with such a dependence mathematically; physically, the uncertainty in LET leads to a continuous curve traced in figure.

Fig. 8

RBE at 10% cell survival for human normal tissue cells irradiated with carbon ions. The results are obtained using Eq. 35. In the high-LET region, the RBE becomes inversely proportional to LET, and the absolute values of RBE depend on the number of ions that traverse the cell nucleus. The values of $N_{\text {ion}}$ corresponding to different segments of the calculated curve are indicated. The dashed line is a guide to the eye connecting median points of the hyperbolas. Symbols depict experimental data from Suzuki et al. (1996), Suzuki (2000), Tsuruoka (2005) and Belli (2008)

Conclusions and outlook

We reviewed the major methodological concepts of the Multiscale approach to the physics of ion-beam therapy and demonstrated that the whole approach converges to a single formula that calculates the yield of lethal lesions in a cell irradiated with ions. This yield, equal to the logarithm of the inverse probability of survival of the cell, depends on the depth, the composition of tissue in front of the cell, oxygen concentration, and the type of the cell. It was demonstrated that the MSA allows one to calculate the probability of cell survival in a variety of conditions, such as high and low values of LET, large and small values of fluence, and aerobic and hypoxic environment. MSA generically predicts linear survival curves, but can explain shouldered curves in special cases. Thus, it is a truly universal and robust method of assessment of radiation damage with ions. Besides its effectiveness, the method answers many questions about the nature of effects that are taking place on a plethora of scales in time, space, and energy.

This review was not intended to compare the MSA with other approaches leading to calculations of survival curves, such as microdosimetric kinetic model (MKM) (Hawkins 1996, 2009) and following modified MKM (Kase et al. 2006), local effect model (LEM) (Scholz and Kraft 1996; Friedrich et al. 2012; Elsaesser and Scholz 2007), and track structure simulations (Friedland et al. 2017; Stewart 2015; McNamara et al. 2017). Such comparisons are desirable for many reasons, but it will require efforts on different sides. All of these approaches are based on dosimetry (nanodosimetry or microdosimetry), i.e., one way or another assuming that the dose per se does the damage. Other parameters and assumptions are present as well, depending on the approach. MSA is the only phenomenon-based approach, i.e., the radiation damage is deemed to be a consequence of series of effects. By design, MSA has to answer why certain effects (e.g., the decrease in cell survival probability at the distal end of the SOBP with a uniform physical dose) take place. Other methods may “include” effects (like the above mentioned) in updated versions and claim that their approach can be used for therapy optimization. However, the optimization of therapy planning deserves a solid theoretical base rather than a solution that somehow works and hopefully treats patients well. Our claim is that the MSA is uniquely designed in response to this quest; it has outstanding predictive qualities, and its reliance on the fundamental science makes it exceptionally valuable for the optimization of treatment planning as was demonstrated in a number of examples in this review. In general, if different methods containing different physics manage to predict comparable cell survival curves, it would be at least interesting to know why. The MSA was designed as an inclusive scientific approach, and so far it lives to the expectations. Its additional strength is in its capability of adjustment to changing external conditions, e.g., the presence of sensitizing nanoparticles (Haume et al. 2018) (of given composition, size, and density). In such cases, additional effects are just included in the scenario. The ion-induced shock wave phenomenon may change the initial conditions for the chemical phase of radiation damage. This prediction could be compared with the track structure simulations if the shock waves were included effectively in their scenario (e.g., by increasing diffusion coefficients for reactive species depending on their positions in the track for some time on ps scale).

Summarizing the future directions for the MSA, much more research, should be done in order to improve and deepen the understanding of the scenario of radiation damage with ions. First, the discovery of ion-induced shock waves, predicted by the authors and already included in the scenario of radiation damage, would be the most significant step toward the recognition of the MSA. Second, a more elaborated scenario of transport of reactive species including the collective flow due to the shock waves as a function of LET and its comparison with MC simulations will also be an important development. Third, a comprehensive study of survival curves for a large variety of cell lines and conditions is definitely desired. Fourth, experiments with high fluences and disabled DNA repair function could explore the effects of tracks overlap, measure the effective radii of ion tracks, which can help better understanding of transport of reactive species. Fifth, a better understanding of SOBP features will certainly improve the optimization of therapy planning and bring it to a more scientific level. Sixth, the sensitizing effect of nanoparticles should be further explored in contact with experimentalists. Finally, the MSA should be applied on the next, larger, scale to optimize the achievement of tumor control as a function of relevant external and internal conditions.