Modeling of dose and sensitivity heterogeneities in radiation therapy Kristin Wiklund

Modeling of dose and sensitivity
heterogeneities in radiation therapy
Kristin Wiklund
© Kristin Wiklund, Stockholm 2012
ISBN 978-91-7447-473-2
Printed in Sweden by Universitetsservice US-AB, Stockholm 2012
Distributor: Department of Physics, Stockholm University
D edicated to Stefan and O liver,
w ho lost their battle against cancer
Sammanfattning på svenska
Den senaste tidens ökade användning av lätta joner inom strålbehandling,
företrädesvis i Japan och Europa, beror på den synnerligen precisa
dosfördelning som kan uppnås med joner tyngre än protoner och på
att jonisationerna sker tätare utmed spåret, vilket resulterar i en ökad
biologisk effekt jämfört med konventionell strålbehandling. För att kunna
utnyttja denna modalitet till fullo så måste dess inbyggda egenskaper
och dess användning i terapeutisk miljö noggrant undersökas. Ett Monte
Carlo-baserat (MC) datorprogram har därför utvecklats för simulering
av transport av primär- och sekundärelektroner i vatten. Programmet
används för olika typer av beräkningar rörande strålkarakterisering.
Vid användning av MC-teknik är beräkningarnas tillförlitlighet i stor
utsträckning beroende av tvärsnitten som är implementerade i koden.
Huruvida små förändringar i tvärsnittet påverkar elektronspårets
utbredning har därför undersökts. En semi-analytisk metod för att kunna
beräkna den radiella dosfördelningen från joner har också tagits fram
och utvärderats, som ett alternativ till MC-beräkningar. Förutom det
tidigare nämnda behovet av bättre karakterisering av strålarna så behöver
modellerna för cell- och vävnadsrespons förbättras för att bättre kunna
förutse resultatet av strålbehandlingen. Den andra delen i avhandlingen
behandlar därför just detta. De extremt heterogena och slumpmässiga
energideponeringarna från joner gör att bra modeller går att hämta från
statistiska teorier. Ett analytiskt uttryck, baserat på statistiska teorier,
för tumörkontrolls-sannolikhet för heterogena dosdeponeringar har
tagits fram och validerats med numeriska simuleringar. Dessutom kan
känslighetsvariationer i vävnader behandlas som en stokastisk effekt och
det härledda uttrycket kan användas även för detta ändamål. I praktiken
kan det förväntas bli kombinationer av dessa effekter i vävnad, så för att
undersöka detta har numeriska simuleringar gjorts.
MC-koden, KITrack, har utvecklats till att bli ett mycket användbart
verktyg för strålkarakterisering. Det framtagna analytiska uttrycket för
tumörkontroll kan hjälpa till att förbättra strålterapi. Ett nytt index, ett så
kallat anisotropi index föreslås. Det ger ett mått på frånvaron av isotropi
hos energideponeringarna och kan hjälpa till att ge en djupare förståelse
till varför en given strålmodalitet ger en viss effekt.
iv
Abstract
The increased interest in the use of light ion therapy is due to the high dose
conformity to the target and the dense energy deposition along the tracks
resulting in increased relative biological effectiveness compared to conventional radiation therapy. In spite of the good clinical experience, fundamental research on the characteristics of the ion beams is still needed in
order to be able to fully explore their use. Therefore, a Monte Carlo track
structure code, KITrack, simulating the transport of electrons in liquid water, has been developed and used for calculation of parameters of interest
for beam characterization. The influence of the choice of the cross sections
for the physical processes on the electron tracks has also been explored.
As an alternative to Monte Carlo calculations a semi-analytical approach to
calculate the radial dose distribution from ions, has been derived and validated.
In advanced radiation therapy, accurate characterization of the beams
has to be complemented by comprehensive radiobiological models, which
relate the dose deposition into the cells to the outcome of the treatment.
The second part of the study has therefore explored the influence of heterogeneity in the dose deposition into the cells as well as the heterogeneity
in the cells sensitivity to radiation on the probability of controlling the tumor. Analytical expressions for tumor control probability including heterogeneous dose depositions or variation of radiation sensitivity of cells and
tumors have been derived and validated with numerical simulations. The
more realistic case of a combination of these effects has also been explored
through numerical simulations.
The MC code KITrack has evolved into an extremely useful tool for beam
characterization. The tumor control probability, given by the analytical derived expression, can help improve radiation therapy. A novel anisotropy
index has been proposed. It is a measure of the absence of isotropy and
provides deeper understanding of the relationship between beam quality
and biological effects.
v
Contribution to papers
Paper I:
K. Wiklund derived the analytical expression in the paper. The C program
for the calculations of the radial dose distributions was developed jointly.
K. Wiklund did all the different simulations, analyzed the results and wrote
about 90% of the text in the paper.
Paper II:
K. Wiklund ran the codes to generate cross section databases. The C++
program was developed jointly. K. Wiklund analyzed the results and wrote
about 70% of the text.
Paper III
The analytical expressions were derived jointly as well as the code. K.
Wiklund ran the simulations, analyzed the results and wrote about 70% of
the text in the paper.
Paper IV
The analytical expressions were derived jointly as well as the code. K.
Wiklund ran the simulations, analyzed the results and wrote about 60% of
the text in the paper.
vi
List of Papers
This thesis is based on the following papers, which are referred to in the text
by their Roman numerals.
I
II
III
IV
Wiklund K, Olivera G H, Brahme A and Lind B K. (2008) Radial
secondary electron dose profiles and biological effects in lightion beams based on analytical and Monte Carlo Calculations
using distorted wave cross sections.
Radiation Research, 170:83-92
Wiklund K, Fernández-Varea J M and Lind B K.(2010) A Monte
Carlo program for the analysis of low-energy electron tracks in
liquid water.
Physics in Medicine and Biology, 56:1985-2003
Wiklund K, Toma-Dasu I and Lind B K. (2011) The influence
of dose heterogeneity on tumor control probability in fractionated radiation therapy.
Physics in Medicine and Biology, 56:7585-7600
Wiklund K, Toma-Dasu I and Lind B K. (2012) Impact of dose
and sensitivity heterogeneity on TCP
Manuscript
Reprints were made with permission from the publishers.
Related publications not included in this thesis
Fernandez-Varea J M, González-Muños G, Galassi M E, Wiklund K, Lind B
K., Anhesjö A and Tilly N (2011) Limitations (and merits) of PENELOPE as a
track structure code
International Journal of Radiation Biology, 88:66-70
vii
Contents
Sammanfattning på svenska . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
Contribution to papers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
List of Papers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 Development of a Monte Carlo track-structure code . . . . . . . . . . . . . . . 5
2.1 The Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 KITrack . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 MC interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
7
9
3 Cross sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1 Electronic inelastic cross sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Implemented inelastic cross section . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 The PENELOPE approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 Comparison of inelastic cross sections . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Electronic elastic cross sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Implemented elastic cross sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Ion inelastic cross sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Continuum Distorted Wave theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
16
17
18
20
21
22
25
4 Modeling of cell and tumor response to radiation . . . . . . . . . . . . . . . . . . 27
4.1 Cell survival models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Tumor response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Tumor response by simulations and analytical modeling . . . . . . . . . . . . . . . . . .
4.3.1 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Analytical modeling -approximation of an expectation value . . . . . . . . . . .
4.3.3 The influence of microscopic heterogeneity in energy deposition . . . . . . .
4.4 The Local Effect Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Three approaches to the calculation of cell survival in the framework of
LEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.2 Results of the different approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
30
33
33
38
41
44
47
49
5 Beam quality characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.1
5.2
5.3
5.4
5.5
5.6
Particle fluence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Slowing down spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Track penetration parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Specific ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Nearest neighbor of energy deposits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Radial dose profile (amorphous track modeling) . . . . . . . . . . . . . . . . . . . . . . . .
6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Summary of papers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
60
61
64
68
73
76
79
81
84
CONTENTS
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
x
1. Introduction
The main goal of radiation therapy, as of any treatment modality, is to improve the quality of life for the patient. Usually this means to eradicate the
malignant (or sometimes benign) disease, while at the same time sparing
the surrounding healthy tissue as much as possible. The end result of radiation therapy, the treatment outcome, depends on several parameters.
Some of those are the physical properties of the beam or the beam quality, dependent upon all intrinsic quantities. The beam quality determines
to a large extent the effect of the radiation in mammalian cells and tissue,
but is not alone responsible for the treatment result. In order to achieve the
best possible treatment outcome, the choice of treatment technique is essential. By accurately and carefully choosing the number of radiation fields
that should be used and the number of fractions the prescribed dose should
be delivered in, the goal of radiation therapy could be reached.
With help of new technologies such as Intensity Modulated Radiation
Therapy (IMRT) and Image Guided Radiation therapy (IGRT), a highly
conformal dose can nowadays be delivered to the targets, and changes in
patient anatomy during the course of the treatment can be monitored and
the treatment could be adapted according to them.
Radiation modality, energy, beam arrangement and fractionation are
all deterministic parameters, which are largely kept under control unlike
stochastic effects. The uncertainty associated with the delivery of the dose,
either from one fraction to another (inter-fractional) or within a fraction
(intra-fractional) is a stochastic effect that is almost always present in
radiation therapy. Those stochastic effects are hard to foresee and the
impact on the treatment outcome is not that well known. Given all the
deterministic parameters, and in rare cases even the stochastic parameters,
accurate radiobiological models of cells and tissue response to radiation
could help predict the treatment outcome.
Besides the uncertainty associated with the delivery of the dose, the
stochastic nature of all energy depositions are associated with a certain
uncertainty in the delivered dose to the target tissues. High Linear Energy
Transfer (LET) beams, such as ions with energies below 10 MeV/u, give
rise to rather heterogeneous energy deposition on microscopic scale
e.g within the cell nucleus, while conventional radiation therapy with
low-LET photons, electrons and protons gives a more homogeneous dose
1
Chapter 1. Introduction
deposition. On macroscopic scale, on the other hand, both low and high
LET beams may show rather heterogeneous dose distribution.
Radiation therapy using light ions benefits from the advantageous shape
of the ion depth dose, which presents a sharp peak, the so called Bragg Peak
(BP), and from a narrow lateral dose profile. Those favorable properties of
the ion beam makes it possible to achieve a high conformity of the dose
(and of the biological effect) to the target volume and at the same time keep
the dose low to the surrounding normal tissue. The high ionization density in and close to the BP is associated with an enhanced biological effect
(for a given absorbed dose in Gy) that often is called the Relative Biological
Effectiveness (RBE).
When a clonogenic cell is exposed to ionizing radiation it may loose
its clonogenic ability by a large number of damaging processes. Some of
the damages are lethal and others are repairable, or at least, conditionally
repairable. High LET beams cause more clustered damages than low-LET
beams. Single Strand Breaks (SSBs), where one of the DNA strands is
damaged, are efficiently repaired in cells with intact repair pathways, and
even Double Strand Breaks (DSBs), where both strands are damaged, are
correctly repaired most of the time (∼98 %). However, clustered damages
are more difficult to repair and can lead to cell death or severe damages to
the cell.
When using the Linear Quadratic model (LQ-model), cells in a tumor are
often characterized by the α/β, roughly giving the average reciprocal sublethal repair capacity of the cells in the tumor. There are other models, such
as the Repairable Conditionally Repairable (RCR) model (Lind et al., 2003),
that explicitly model the repair. Due to genetic instability there is often more
than a 30 % sensitivity variation from one cell to another within a tumor.
Furthermore, the same cell could have this high sensitivity variation from
one fraction to another due to inherent variations of the cell sensitivity in
different stages of the cell cycle, or due to the oxygen or nutrition level in the
cell.
The increased interest in the use of radiation of higher ionization densities puts even more stress on the types of models that can be used in order
to predict cell survival and or organ/tumor responses. The pioneering treatments with carbon ions in Chiba, Japan and in Darmstadt, Germany have
been made possible by the use of simplified but still useful models for the
biological effects in those beams (Scholz and Kraft, 1996; Kanai et al., 1999).
The Local Effect Model (LEM) (Scholz et al., 1997) is one of those models
and is used in Darmstadt and has been debated the recent years. Therefore,
in the summary of this thesis, a section about the LEM model has been included (section 4.4).
Since one of the most obvious characteristics of radiation with high ionization density is the extremely heterogeneous energy deposition on a microscopic scale e.g in the cell nucleus at therapeutic doses, it is very likely that
2
good model candidates can be developed from statistical theories. Variation
in the sensitivity of the targets can also be treated as stochastic effects and
be taken into account with help of statistical theories. Models estimating the
Tumor Control Probability (TCP) often do not take into account the fact that
the delivered dose is generally non uniform, either intentionally or due to
stochastics. Therefore, in this thesis an attempt have been made to include
those effects in an analytical expression based on statistics, which is validated with numerical simulations (Paper III). Sensitivity variations among
the cells in a tumor have a substantial impact on the TCP, similar to the dose
variations, and should also be included in TCP estimations. Sensitivity variations within a tumor and combinations of heterogeneous dose distributions and sensitivity variations have been investigated in Paper IV.
Due to the stochastic nature of energy deposition events in tissue, Monte
Carlo (MC) simulations are suitable for investigation of various parameters
related to RBE. Therefore, a MC track structure code has been developed in
this thesis. Paper II includes a description and benchmarking of this code.
There are mainly two different approaches developed to model the biological effects of charged particles on the basis of their specific energy deposition pattern: 1) the microdosimetric approach, which emphasizes the
stochastics of the energy distribution on the basis of single electrons; 2) the
amorphous track structure approach, which is based on the radial dose distribution as a continuous distribution, representing an average over many
tracks (Butts and Katz, 1967; Katz et al., 1971, 1994). In paper I, MC simulations have been combined with theoretical expressions to calculate the radial dose profile and verified with purely MC calculations of the radial dose.
The MC track structure calculation require accurate cross sections for all
the considered interactions in order to reach reliable results. At very low energies there is still a scarcity of reliable cross sections even though great effort has been made recently to improve existing models and to develop new
models. Paper II in this thesis concerns the impact of different cross sections
on the track structure at very low energies.
Fundamental to interpreting the results of MC calculations, is to understand some probability theory. The MC method for radiotherapy gives an
accurate and general method for dose calculations, fluence and spectrum
calculations and for calculation of dosimetric parameters such as stopping
power. It has become very efficient due to faster algorithms and faster hardware. In principle it could facilitate implementations of new treatment techniques and the technique can be used to perform ’experiments’ that can not
be done in the lab. It has to be kept in mind that for microdosimetric quantities and, to some extent, for macroscopic quantities, MC calculations only
gives an average or an expectation value of the quantity of interest, but with
an arbitrary good accuracy (within the “feasible engineering time”).
The actual radiation response of a tumor and surrounding tissues is a
rather complex function of tissue organization, cellular- and sub-cellular
3
Chapter 1. Introduction
mechanisms, time dependence, signaling pathways, bystander effects, cell
cycle sensitivity, inter cellular sensitivity together with the micro- and macroscopic ’dose’ distribution. More precise, TCP is dependent on the number of ionizations per unit volume (or per cell) or, to go one step further,
on the numbers and distribution of DSB or even more complex damages in
each cell caused by direct and indirect effects. Those measures are also dependent upon DNA density in the cells and finally, the knowledge how the
number of DSB and other complex damages correspond to different kinds of
cell death. A brief investigation of the number of Specific Ionization, SI, and
number of Specific Primary Ionization, SPI, in relation to RBE etc has been
done and can be found in section 5.4. Due to the importance of the spatial
pattern of energy deposits on the biological effect, an attempt to mathematically model the nearest neighbor distance has been made in section 5.5.
The chemical stage, i.e. indirect effects and in-vivo phenomena such as
bystander effects are disregarded. Models of the DNA and cells are not dealt
with within this thesis. Instead, focus is put on characterizing radiation
quality and find alternative measures for LET, RBE, better cell survival and
TCP models with help of approximations and simplifications.
4
2. Development of a Monte Carlo
track-structure code
2.1
The Monte Carlo method
The Monte Carlo method provides approximate solutions to problems by
performing statistical sampling experiments using random numbers, and
its origin can be traced back to G. Comte de Buffon in the 17th and 18th
century. The first to use the MC technique with the help of computers were
Neumann and Ulam in 1949. The MC method has found widespread use in
applications of imaging and high-energy physics and in medical radiation
physics, where it is mainly used for solving the Boltzmann transport equation. The stochastic nature of radiation interaction in matter makes it a good
candidate for MC technique. In this thesis, the MC method has been used
to calculate electron tracks from which a number of fundamental quantities can be deduced (Paper I, II and parts of the summary) and to calculate
Tumor Control Probability (TCP) when dose and sensitivity variations are
present (Papers III and IV and section 4.3 in the summary).
Whenever the MC method is adopted, there is always a need to verify the
results. This can be done by comparison with measured data, if available,
otherwise by comparison with the results of well-tested Monte Carlo codes
(done in paper II). The code developed in this thesis can also be tested, at
least partly, by calculation of the stopping power by a brute force method
(see section 3.1).
Simulation procedure
MC simulations are often used for complex problems and high-precision
results are associated with a considerable CPU time even on a fast computer.
Therefore, variance-reduction techniques can be employed to increase the
convergence rate. Those techniques can either be applied in the scoring step
or in the way the random walk is sampled. In many situations they are indispensable, but have to be used with care to avoid erroneous results, especially
in track-structure calculations of low-energy electrons.
All MC codes rely on good random numbers, and on accurate probabilities for the implemented physical processes. For MC track-structure simulations in the area of medical radiation physics those processes are described
by the interaction cross sections (CS) for the particles of interest (electrons,
photons and ions) in liquid water, which is being used as tissue equivalent
5
Chapter 2. Development of a Monte Carlo track-structure code
(sometimes water vapor is used instead of liquid water), or other elements
like Tungsten that can be found in the head of treatment units.
In each interaction, the particle may lose the energy, W, and/or change its
direction of movement. The angular deflection is determined by the polar
scattering angle θ, i.e. the angle between the directions of the particle before and after the interaction, and the azimuthal angle, φ. All the interaction
processes that take place when a particle passes the medium (if the target is
assumed to be free and at rest) are expressed by molecular double differential cross sections (DDCS):
∂2 σ
,
∂Ω ∂W
(2.1)
¡
¢
where ∂Ω is the element of solid angle in direction θ, φ . Due to the random
organization of the molecules in the medium, the differential cross section
(DCS) is independent of the azimuthal angle, φ, and is therefore sampled
from a uniform distribution, namely φ = 2π ξ, where ξ is a random number.
The energy-loss DCS, differential in only the energy loss, is obtained by integrating the DDCS over the solid angle and the total cross section, σ, per
water molecule is defined as the integral of the energy-loss DCS over all the
possible energy losses.
Classification of MC codes
MC particle transport codes can be classified in three groups, class I, class II
or event-by-event codes, according to the transport algorithm type (Salvat
et al., 2008). The class I and class II, condense history and mixed simulations respectively, group to different extent the effect of the various interactions into a single step that represents the change in direction and energy
loss due to the collisions. Class I simulations group together all interactions
of the primary particle, while class II only groups the soft collision and treat
the catastrophic hard interactions individually. The third approach, that has
been used almost exclusively in this thesis, is the event by event, or track
structure codes, where all the particles, primaries and their secondaries, are
followed during their passage through matter until a predetermined cutoff
energy is reached, or in absence of a cutoff, until the particle dissipates.
Many codes in the three groups described cannot follow heavier particles
like ions, but they can perform track segment simulations. Track-segment
conditions imply that, in a short segment, the energy loss of the ion is small
compared to its initial energy and can be neglected at least for some purposes. Furthermore, when performing track segment calculations, one may
neglect multiple-scattering angular deflection due to the relatively small angular deflection of the ion in that short track segment. Fragmentation of
the ion and other nuclear reactions are also disregarded by codes that don’t
transport the heavier particles. Several track-structure codes have been developed during the last two decades both for electron and ion transport. A
6
2.2. KITrack
comprehensive list of track structure codes up to the year 2006 is given in
the review article by Nikjoo et al. (2006).
2.2
KITrack
The object-oriented code developed in the present work has been written
in C++ and utilizes the GNU Scientific Library (Galassi et al., 2009) in order
to streamline it. The code follows the primary and the secondary electrons
in liquid water until their kinetic energy falls below a selected threshold energy (10 eV in the present calculations unless stated otherwise). Electrons
are then thermalized and allowed to deposit their remaining energy at a final spot at a distance from their final interaction point. This spot is chosen
to be isotropically distributed around the interaction point and the distance
is sampled from an exponential distribution. The mean value of the exponential distribution is set to the mean free path of an electron with energy of
10 eV. The energy interval presently addressed extends up to 100 keV. Consequently, the selected theoretical approaches do account for relativistic effects. No variance reduction is applied in the developed code.
Owing to the modularization of the code and the tree structure for storage
of all generations of the particles, spatial coordinates and energies together
with energy deposit (ICRU, 1998) and finally the particles generation and
descent, almost all possible quantities can be deduced or specific calculations be carried out with the code.
However, a novel way of scoring has to be introduced in order to overcome the fact that the code does not use voxel-based geometry. Applications
of this can be found in section 5.1. Any model, energy spectrum or simply
monoenergetic electrons can be adopted as a source of electrons, depending on the purpose of the MC simulations.
Basic track-segment simulations of light ions can be carried out to calculate the electron fluence spectrum from different ion species in order to
characterize beam quality in a relevant way as needed for optimized treatment planning with ions, see for example Paper I. The CDW-EIS theory, described in section 3.3.1, has been used exclusively in this thesis for all calculations regarding ions.
Random number generation
The developed track structure code has been equipped with a Linear
congruential random number generators (LCRNG), the built-in drand48
(Knuth, 1997). It has also been equipped with two additional Pseudo
Random Number Generators (PRNGs), the Mersenne twister (MT)
(Matsumoto and Nishimura, 1998) and Single Instruction Multiple Data
(SIMD)-oriented MT (Saito and Matsumoto, 2008) which are widely-used,
fast PRNGs with a long period of 219937−1(a Mersenne prime). Those
7
Chapter 2. Development of a Monte Carlo track-structure code
PRNGs were added to assure high level of randomness in simulations that
require long periods of the random numbers. The initialization seeds to the
PRNGs can be chosen to be fixed if the simulations require the same initial
condition, otherwise the seeds are set according to the value of built-in
clock in the computer and become in that case more or less random.
Sampling of interaction events
During the MC simulations, sampling has to be performed an enormous
amount of times from the Probability Density Functions (PDFs) of the interaction cross section, so it is essential that the sampling routine is fast and
reliable and PDFs of the processes invoked have to be accurate. The implemented cross section databases are recorded in histograms as PDFs by help
of the GNU scientific library (Galassi et al., 2009) and the Inverse-transform
method (Salvat et al., 2008) has been utilized as a sampling method. It is
a fast method, faster than for instance the Acceptance-Rejection method
where many more random numbers are used and thus leads to a higher time
consumption. The Inverse-transform method provides the most straightforward way of generating random samples of a random variable. The inversetransform method can only be applied to monotonous increasing functions,
i.e functions that have an inverse. A cumulative distribution always fulfill
this criteria
The inelastic cross section database was calculated with the model by
Dingfelder et al (1998; 2008) and the elastic cross section database was created by the ELSEPA code (Salvat et al., 2005). Details regarding cross sections
are presented in chapter 3.1.
While most MC codes use inelastic cross section differential in transferred
energy Eq. (2.1), an alternative and slightly different approach has been used
here. This approach may be a more pragmatic and logic approach for MC
simulations, with the cross section differential in the secondary electron energy, ε, instead of in energy transfer W :
dσ
dε
(2.2)
and the relation for the energies involved in the process:
W = E − E ′ = Ui + ε,
(2.3)
where Ui is the binding energy in the i t h molecular orbital and E ′ is the
kinetic energy of the primary electron after the collision. The angles of the
scattered electron, θprim , and recoil electron θsec , are derived from the kinematic expressions assuming that the target electron is free and at rest:
θprim =
E + 2m e c 2
E −W
E E − W + 2m e c 2
(2.4)
W E + 2m e c 2
.
E W + 2m e c 2
(2.5)
θsec =
8
2.2. KITrack
The present MC track-structure code considers eight different processes,
namely excitation, elastic scattering, dissociative attachments and ionization, where ionization is divided in five separate processes, one for each
shell of the water molecule. The total interaction cross section then becomes
σtot = σexc + σel + σdiss +
5
X
σion,i .
(2.6)
i =1
To determine the type of interaction that will take place at an interaction
point, so-called branching ratios are employed (Bielajew, 2012)
p pr ocess = σpr ocess /σt ot ,
(2.7)
by throwing a random number and follow the scheme
i f ξ < σpr oc1 /σt ot ⇒ pr oce ss 1 w i l l t ake pl ace
¢
¡
el se i f ξ < σpr oc1 + σpr oc2 /σt ot ⇒ pr oce ss 2 w i l l t ake pl ace
el se ⇒ pr oce ss 3 w i l l t ake pl ace.
(2.8)
Knowing which process that will take place, the energy deposit and/or angular deflection will be sampled from the respective cross sections. In the case
of an ionization, the energy of the secondary electron is first sampled and
the energy of the scattered electron is calculated from Eq. (2.3). The angular
deflection of the two particles is finally calculated by Eqs. (2.4)–(2.5).
Integrated and differential cross sections have been tabulated on a common logarithmic grid of ‘primary’ electron energies. Elastic differential cross
sections (DCSs) dσ/dΩ were evaluated for a grid of polar scattering angles
θ, and the various ionization and excitation DCSs dσ/dε were calculated for
a mesh of secondary electron energy ε. MC interpolation (Benedito et al.,
2001) is performed whenever the electron energy does not coincide with the
primary grid, see section 2.2.1. Notice that a single primary energy grid for
all cross sections, both total and differential, has the advantage that the index into the grid resulting from MC interpolation can be saved and reused,
thus minimizing operations and rendering the simulation faster. Tables of
differential and total cross sections were precalculated employing the theoretical formalisms and semi-empirical fits outlined in section 3.
2.2.1
MC interpolation
In MC simulations, the sampled energy of the particle seldom coincides
with an energy stored in the primary grid, which implies that an interpolation routine have to be called each time. Therefore, a fast and accurate interpolation routine has to be used to minimize the time consumption. The
most straightforward solution might be to make use of the so-called MC
interpolation, which exploits random numbers, and the fact that the MC
9
Chapter 2. Development of a Monte Carlo track-structure code
calculations are repeated until desired accuracy is achieved. Depending on
the energy dependence of the cross section, different transformations can
be used in the MC interpolation, e.g linear (lin), linear-logarithmic (lin-log)
(Benedito et al., 2001) and logarithmic-logarithmic (log-log). An investigation of the time consumption and the relative error for the three types of MC
interpolations has been made in order to conclude which one is the most
time effective in relation to the error it induces.
Interpolation error
The inelastic and elastic cross section differential in energy can be approximated (parametrized) by the function:
y = a E b,
(2.9)
where a and b are constants and E is the kinetic energy, see Figure 2.1. If a
logarithmic energy grid is used, as the present code make use of, the following relation between successive grid points holds
E i +1
= k,
Ei
(2.10)
where k is a constant. The code uses a relative fine grid with k = 101/50 (~5%
steps). The interpolation polynomial is expressed in terms of the energy and
two successive cross section values in the grid.
¡
¢
P (E ) = q (E ) y 1 − y 0 + y 0 .
(2.11)
In case of linear interpolation the expression for q(E) becomes
q (E ) =
E − E0
E 0 (k − 1)
(2.12)
and for the lin-log interpolation suggested by Benedito et al. (2001):
q (E ) =
ln (E /E 0 )
.
ln (k)
(2.13)
The resulting relative error can be expressed as
¡
P (E ) − y
y
¢
=
³
´i
E 0b h
b
1
−
q
·
k
−
1
− 1.
(E
)
Eb
(2.14)
Figure 2.2 shows the estimated relative interpolation error as a function of
the energy of the sampled energy. This error could in principle partly be
compensated for, but that would add extra time to the interpolation. The
lin-log interpolation by Benedito et al. gives a slightly better accuracy above
500 eV for excitation and ionization cross sections than the linear interpolation, below 500 eV are the opposite true. For the elastic cross section the
10
2.2. KITrack
lin-log approach gives slightly higher precision for the whole energy range
used here. The relative error above 500 eV are in the range of 0-0.0003 for
Benedito et al., 0-0.0005 for linear interpolation and below 500 eV; 0-0.001
for Benedito et al., 0-0.0005 for linear interpolation. The exact expression for
logarithmic interpolation (log-log) is given by:
i´
³
h
(ln E 1 −ln E)·ln y 0 +(lnE−ln E 0 )·lnE 1
exp
(ln E 1 −lnE 0 )
¢
¡
.
(2.15)
q (E ) =
y0 − y0
There is no need for the exact logarithmic expression for interpolation due
to the small error resulting from both lin-log and lin interpolation, but the
error depends on the density of the grid. In case of a dense grid linear interpolation can be used without loss of accuracy, while the more complex linlog or log-log interpolation routines should be used when the grid is sparse.
-18
10
Ionization
Excitation
Elastic
Cross sections / m
2
-19
10
10
10
b~-1
b~2
-20
-21
-22
10
10
-23
1
10
2
3
10
10
4
10
5
10
Energy / eV
Figure 2.1: The parametrization of the electron cross section in liquid water in different energy ranges
Time consumption
Table 2.1 shows the number of operations the different interpolation
routines uses and the resulting total time consumption for performing an
interpolation. The computing time is very high for log-log interpolation
compared to lin, lin MC, lin-log, lin-log MC interpolation. Almost no
difference is seen between lin and lin MC or lin-log and lin-log MC.
Benedito et al. claimed that a linear or linear-logarithmic interpolation on
a dense grid is only feasible for limited energy ranges. This is not true when
11
Chapter 2. Development of a Monte Carlo track-structure code
1e-03
Lin b=-1
Log b=-1
Lin b=2
Log b=2
Relative interpolation error
1e-03
8e-04
6e-04
4e-04
2e-04
0e+00
1
1.01
1.02
1.03
1.04
1.05
Fractional energy
Figure 2.2: Estimation of the interpolation error. A 101/50 grid is used. The x-axis is
the sampled energy (between to successive grid points) and y-axis the corresponding relative error for the sampled energy. b = 1 relate the part of cross section above
500 eV (se fig. 2.1) and b = 2 for below 500 eV.
indexing is used, like in the present code, where the index for the energy is
calculated and no time-consuming search function is needed.
The MC interpolation works satisfactory only in MC applications, due to
averaging over many simulations. If applied in other cases the MC interpolation will give larger deviations from the true value than, for instance, the
nearest-neighbor method. The accuracy has to be put in perspective to the
time consumption for the different interpolation routines. A dense grid and
linear interpolation are preferable as a good compromise between time consumption and accuracy.
The default MC interpolation routine in the present track code is linear
MC interpolation owing to the favorable time consumption in relation to
the relative error.
12
2.2. KITrack
Nr. of op.
Add/ Mult/
Assign
Subtr Div
ln()
PRNG if()
exp()
Time
[ns]
lin
4
2
1
0
0
0
0
8.8
lin MC
2
1
1
0
1
1
0
9.4
lin-log
4
2
1
4
0
0
0
490
lin-log MC
2
1
1
4
1
1
0
490
log-log
8
6
1
8
0
0
1
1100
Table 2.1: The number of different types of operation used by the different MC interpolation routines and the resulting total time consumption for the routines.
13
Chapter 3. Cross sections
3. Cross sections
The most fundamental quantity is the triple differential cross section, giving the probability of energy transfer, scattering angle of the projectile and
emission angle of the target. Furthermore the cross sections contain information about which orbital in the atom or molecule is involved in the interaction process. One also needs information about the rearrangement reactions of the molecular species resulting from the interaction. Although
elastic interactions do not contribute to energy dissipation of the projectile,
their proper treatment is essential because the angular deflections affect to
a large extent the spatial evolution of electron tracks. Complete and accurate data sets for inelastic and elastic scattering, either experimental or theoretical, are still missing, especially for low energy electrons. Estimations
for condensed media are usually scaled from gas-phase data. The literature list of experimental and theoretical electron cross section is extensive.
For a review of cross sections see Uehara et al. (1999) and Itikawa and Mason (2005). A very brief overview of the cross sections used in this thesis will
be presented in this section. Some discussion on the differences and similarities between the cross sections adopted in this thesis and cross sections
often used for similar calculations in the literature is also given.
3.1
Electronic inelastic cross sections
The inelastic interaction of charged particles with condensed media can
be described by means of the complex dielectric response function (DRF),
ε(W, q) = ε1 + i ε2 (as a function of energy transfer, W, and the momentum
transfer, q), with its real (ε1 ) and imaginary (ε2 ) parts related to screening
and absorption properties of matter. Within the first Born approximation,
the imaginary part of the reciprocal of the DRF in the energy-momentum
plane provides all necessary information to calculate the double differential
cross section of non-relativistic particles:
¸
·
j
ε2 (q,W )
1 X
1
1
d2 σ
=
=
Im
,
dqdW
πa 0E
ε(q,W )
πa 0 E j ε21 (q,W ) + ε22 (q,W )
(3.1)
where a 0 is the Bohr radius, E is the incident electron kinetic energy and j
corresponds to the contribution of the j t h electronic transition mode. The
imaginary part of the reciprocal of the DRF represents the so-called Energy
14
3.1. Electronic inelastic cross sections
Loss Function (ELF) which can, in the optical limit, q = 0, be expressed
e = n+iκ :
either in terms of the complex index of refraction, n
¸
·
¸
·
−1
−1
= Im
,
(3.2)
Im
ε(q ≈ 0,W )
(n − iκ)2
or in terms of the optical oscillator strength (OOS), f (q ≈ 0,W ):
¸ πE 2
·
−1
pl d f (q ≈ 0,W )
=
,
Im
ε(q ≈ 0,W )
2W
dW
(3.3)
where E pl is the plasmon energy of liquid water (E pl ≈ 20 eV). The Generalized Oscillator Strength (GOS), f (q,W ), of an atom describes completely
the effects of inelastic interaction on the projectile, within the first Born approximation.
Presently there are two experimental optical data sets for liquid water that
extend over a sufficient part of the valence excitation spectra. The first set is
from the Oak Ridge Group (Heller et al., 1974) which provides data for the
complex index of refraction and the second set of optical data is that of the
Sendai Group (Hayashi et al., 2000) who measured the GOS via inelastic Xray scattering spectroscopy.
The variation of the ELF in the (q,W ) plane, i.e. the Bethe surface, must
be known to enable calculation of the inelastic cross sections via Eq. (3.1).
Therefore various extension algorithms have been developed to model
the q-dependence. There are mainly two approaches for these extension
algorithms or dispersion relations based either on the extended Drude
function (Hamm et al., 1976; Ritchie and Howie, 1977), which is often used
by the so-called optical data models (ODM), or based on the estimation of
the GOS.
Optical data models
The deconvolution of the experimental spectrum is made by a linear combination of Drude-type functions, the so-called extended Drude function:
D i (q,W ) ≡ £
f i W γi (q)
¤2 £
¤2 ,
2
E i (q) − W 2 + W γi (q)
(3.4)
where f i , γi and E i are the oscillator strength, damping coefficient and
transition energy, respectively. The ELF is proportional to the extended
ª
©
Drude function and the model parameters f i , γi , E i at q = 0 are obtained
via a fit to the experimental optical data. In order to ensure that the optical
data model behaves reasonably well beyond the range over which the
experimental data are available, the following sum rules (Bethe sum rules)
have to be fulfilled:
ˆ ∞
π 2
(3.5)
W ε2 (q,W ) dW = E pl
2
0
ˆ ∞
¸
·
π 2
−1
dW = E pl
.
(3.6)
W Im
ε(q,W )
2
0
15
Chapter 3. Cross sections
Low energy corrections
The inadequacy of the first Born approximation for low-energy electrons
implies the need for corrections to compensate for the well-known overestimation of the cross section at low energies. As the incident electron slows
down, the strength of the interaction becomes comparable to the incident
energy, and cannot be treated as a first-order perturbation. In addition, the
exchange effects in electron-electron interaction, i.e. the indistinguishability of the scattered and ejected electron, cannot be neglected. In the opticaldata models these effects are taken into account by adding the low-energy
perturbation-corrected cross section, σp , and the exchange-corrected cross
section σexc to get the corrected cross section σ
e:
dσp dσexc
dσ
e
=
+
.
dW
dW
dW
(3.7)
Perturbation and exchange corrections are performed with slightly different
approaches in various models due to the lack of existing theory for those
corrections.
3.1.1
Implemented inelastic cross section
The present track-structure code (paper II) has been equipped with the optical data model by Dingfelder et al. (1998; 2008) as a default inelastic scattering model. In addition, the model by Dingfelder et al. with Born-Ockhur
exchange correction and the model by Dingfelder et al. with a modification
have been implemented and may be employed by the user. These two models have been calculated in this work and explained briefly below. Future
work will include implementation of the most recent ODM by Emfietzoglou
and Nikjoo (2007).
Differences in the handling of the Drude function
One of the main differences between the original model by Dingfelder et al.
(1998; 2008) and the model by Emfietzoglou and Nikjoo (2007) lies in the
way they treat the truncated Drude function. The Drude function needs to
be truncated at the threshold energy of each ionization channel. Both studies truncate the Drude function by a step function. Emfietzoglou and Nikjoo
used a smooth truncated Drude function by subtracting a simple exponential function. Dingfelder et al. employed a smeared Drude function by multiplying it with a Gaussian function. Emfietzoglou and Nikjoo redistribute
the truncated parts of the Drude function to lower level transitions and thus
obtained the full Drude function which automatically fulfills the Bethe sum
rules at all momentum transfers. Dingfelder et al. on the other hand did not
perform this redistribution and as a consequence the Bethe sum rules are
£
¤
not fulfilled. To remedy for this drawback, Im −1/ε(q,W ) has been modified by an ad hoc q-dependent renormalization factor so that the Bethe sum
16
3.1. Electronic inelastic cross sections
rules were satisfied to within 1% and the result is what is presently called the
corrected Dingfelder model.
The improved and more sophisticated dispersion relations in the model
by Emfietzoglou and Nikjoo compared to the model by Dingfelder et al.
also give rise to notable differences. In the model by Emfietzoglou and Nikjoo all three Drude terms, f i , γi , E i , are allowed to disperse, while in the
model by Dingfelder et al. only two of them, f i , E i , disperse. Moreover, the
q-dependence of E i is improved in the work by Emfietzoglou and Nikjoo by
adding a retardation effect.
The most obvious difference between the two mentioned ODMs is the
optical data sets that were used for the fitting of the Drude functions. Dingfelder et al. fit to the data of Heller et al. (1974), while Emfietzoglou and Nikjoo make use of the data set from Sendai group (Hayashi et al., 2000). However, Emfietzoglou et al. (2012) show that the results depend more sensitively upon the extension algorithm adopted than on the optical data used.
Specific low energy corrections
Dingfelder et al. used Paretzke’s low-energy correction function, Φ(E ), derived from vapor data to calculate the perturbation effect:
dσp
dW
=
X
Φ j (E )
j
dσ(E )
,
dW
(3.8)
where E is the incident energy of the electron and dσ/dW is the uncorrected differential cross section for inelastic scattering. Emfietzoglou and Nikjoo (2007) on the other hand took into account perturbation by including a
second order term:
dσp
dσ dσ(2)
=
+
.
(3.9)
dW
dW
dW
Emfietzoglou and Nikjoo applied the Mott-like modifications directly to
the cross section terms and this allows to account for exchange. Dingfelder
et al. used a slightly simpler approach where the modified Mott cross section
itself is used to model the exchange. The Born-Ochkur exchange correction
suggested by Fernández-Varea et al. (1993) has been implemented as an alternative to the exchange correction in the model by Dingfelder et al..
3.1.2
The PENELOPE approach
The PENELOPE code (Salvat et al., 2008), used in paper I for broad beam
depth-dose calculations and in paper II for depth-dose comparison with the
developed track structure code, uses an estimation of the GOS to calculate
the inelastic cross section for electron impact. The Sternheimer–Liljeqvist
model (Salvat et al., 2008) represents the generalized oscillator strength for
each shell in the atom as a δ-oscillator and is adopted in PENELOPE:
F (Wk ;Q;W ) = δ (W − Wk ) Θ (Wk −Q) + δ (W −Q) Θ (Q − Wk ) ,
17
(3.10)
Chapter 3. Cross sections
where Q is the recoil energy, Wk is the resonance energy of k t h shell, δ (x)
is the Dirac delta function and Θ (x) is the Heaviside step function. The first
term represent distant collision (low-Q interaction) and it is proportional
to the photoelectric cross section (acquired from optical data). The second
term corresponds to close collisions (large-Q interactions), in which the target electrons react as if they were free and at rest (W = Q). For more details
see Salvat et al. (2008 ).
3.1.3
Comparison of inelastic cross sections
Figure 3.1 shows cross sections calculated with the Dingfelder et al. model
(1998; 2008), the corrected Dingfelder et al. model, and the Dingfelder et
al. model with the Born–Ochkur exchange correction. The use of the Born–
Ochkur exchange correction reduces the ionization cross section at low energies compared to the exchange correction used in the original model by
Dingfelder et al. By renormalizing the ELF, a substantial decrease in the resulting cross section can be seen, most visible at intermediate energies. The
collision stopping power, i.e. the energy loss per unit length traversed by the
particle, is calculated from the inelastic cross section and is therefore a direct measure of the inelastic cross section:
S col = N
ˆ
∞
W
0
dσ
dW,
dW
(3.11)
where N is the number of molecules per unit volume. The stopping power
can also be calculated directly from the MC track structure code by allowing a particle with initial energy E to interact only once with the medium
(here liquid water) and score the two quantities, mean free path, λ, and energy transfer, W , resulting from that interaction (both elastic and inelastic)
and, therefore, this acts as a verification of the code. Comparison of stopping
power calculated by this technique with the already implemented model by
Dingfelder et al. (1998; 2008), the PENELOPE tabulated stopping power and
the stopping power calculated with the model by Emfietzoglou and Nikjoo
(2007), which is planned to be incorporated in the code, are illustrated in
Figure 3.2. It should be mentioned that below 100 eV the accuracy of the
cross sections is not well known and all the models (the models implemented in the PENELOPE code, Dingfelder et al. (1998; 2008), Emfietzoglou and
Nikjoo (2007) etc) differ largely. For more information about the influence
of slightly different cross sections on the output of MC track structure calculations, see paper II.
18
Inelastic cross section /m
2
3.1. Electronic inelastic cross sections
2×10
1×10
Corrected Dingfelder
Dingfelder
Dingfelder with Born-Ochkur
-20
-20
10
1
10
2
3
10
Energy / eV
10
4
10
5
Figure 3.1: Total inelastic cross section (e − + H2 O) implemented in the track code
as a function of energy of the incident electron calculated with the Dingfelder et al.
model, the Dingfelder et al. model with the Born–Ochkur exchange correction and
the corrected Dingfelderet al. model.
2
S / eV nm
-1
10
1
10
0
10 1
10
2
3
10
10
4
10
Energy / eV
Figure 3.2: Collision stopping power of electrons in liquid water. The solid curve is
the predictions of the default model in PENELOPE (Salvat et al., 2008), whereas the
dashed curve corresponds to the calculations by Emfietzoglou and Nikjoo (2007).
The dashed-dotted line gives the values obtained by Dingfelder et al. (1998) and the
dotted line corresponds to the corrected Dingfelder et al. model, directly calculated
with the present MC track code.
19
Chapter 3. Cross sections
3.2
Electronic elastic cross sections
The scattering of electrons (or positrons) by a central potential field V (r ) is
completely described by the direct scattering amplitude, f (θ), and the spin
flip scattering amplitude, g (θ):
dσ
= |F (θ)|2 + |G(θ)|2
dΩ
X
F (θ) = exp(iq · r i /~) f i (θ)
(3.12)
(3.13)
i
G(θ) =
X
exp(iq · r i /~)g i (θ),
(3.14)
i
where q is the momentum transfer vector and r i is the position vector of
the i t h nucleus in the molecule relative to an arbitrary origin. The direct
and spin-flip scattering amplitudes are complex functions of the polar scattering angle, θ, determined by using the large-r dependence of the Dirac
distorted plane wave, i.e. the solutions to the Dirac equation for a central
potential V (r ) which behaves asymptotically as a plane wave (alternatively
the Schrödinger equation could be used). The scattering amplitudes for the
atoms in a molecule can be summed up either coherently or incoherently.
Incoherent scattering
In case of incoherent scattering the molecule is regarded as a cluster of independent atoms and by using the additivity approximation (the Bragg sum
rule) the DCS becomes:
µ
¶
¶
µ
dσ
dσ
dσ
+2
,
(3.15)
=
dΩ
dΩ O
dΩ H
where indices O and H denote oxygen and hydrogen, respectively.
Coherent scattering
Coherent scattering takes into account the phase shift and a possible
interference of the wave from atoms in a molecule. A relatively simple
approach, which yields a good estimate of these effects is provided
by the independent-atom approximation. By the assumption of
randomly-oriented molecules the DCS is given by
dσ
dΩ
= 〈|F (θ)|2 + |G(θ)|2 〉
=
X sin(qr i j ) £
i,j
qr i j
¤
f i (θ) f j (θ) + g i (θ)g j (θ) ,
(3.16)
where r i j is the distance between i t h and j t h nucleus in the molecule.
Chemical binding and the distortion in the electron distribution for an
atom bound in a molecule has only a minor impact on the DCS and are
neglected in this approximation.
20
3.2. Electronic elastic cross sections
3.2.1
Implemented elastic cross sections
The ELSEPA package (Salvat et al., 2005) was used to create the elastic database for in the MC track structure code developed in the present work. This
package performs partial-wave analysis (PWA) of the Dirac equation for a
local central interaction potential V (r ). This is now customarily used due to
the full inclusion of relativistic effects. The molecule is regarded as consisting of spherical atoms, and each atom is regarded as a separate scattering
unit. This approach, where intramolecular multiple scattering is neglected,
is referred to as the independent atom-approximation and is valid for incident electron energies above 100 eV. The static-field approximation, where
the atom is assumed to behave like a frozen distribution of electric charge,
was adopted, and the process can be described as scattering of a projectile
by the electric field of the atom. Different corrections can be added in order
to account for deviations from that ideal situation. The resulting total scattering potential used to model the interaction between the projectile and
the target is assumed to be spherical symmetric to be manageable in the
PWA:
V (r ) = Vst (r ) + Vex (r ) + Vcp (r ) − iWabs (r ).
(3.17)
r is the distance from the center of the atomic nucleus, Vst (r ) is the static
potential, Vex (r ) is an exchange potential, Vcp (r ) is the so-called correlationpolarization potential and Wabs (r ) is the magnitude of an absorption potential. To account for exchange, i.e. the fact that the incident electron can be
interchanged with an atomic electron, the Furness and McCarthy exchange
potential was utilized (Furness and McCarthy, 1973). The electric field of the
incident electron can polarize the target atom/molecule and the polarization effect should be taken into account for energies below 10 keV. However,
polarization effects were not accounted for here, due to lack of experimental
information of such phenomena in the liquid phase. For projectiles with
kinetic energy above the first excitation threshold, there is a loss of particles
from the elastic channel to the inelastic channels. This effect has been neglected in the generation of the databases. Again, lack of experimental information is the reason behind this decision. A point nucleus was chosen
for the charge distribution of the nucleus, which results in a pure coulombic
nuclear potential. The details of the nuclear charge distribution can affect
only the DCS for projectiles with energies larger than about 50 MeV. The numerical Dirac–Fock electron densities are default in ELSEPA. They are the
most accurate expression for free atoms and have been used here for the
database generation.
In the present track code, the calculated incoherent elastic cross section
has been set as a default database and the coherent one as optional. A comparison of the two inelastic cross sections implemented in the code is seen
in Figure 3.3. The differences between incoherent and coherent inelastic
cross sections are quite substantial.
21
Chapter 3. Cross sections
The accuracy of the elastic scattering cross sections below 100 eV should
be questioned. For instance, note that by performing density scaling the
elastic scattering cross section becomes ∼2 · 10−19 m2 for an electron with
energy 10 eV and this results in a mean free path of 0.15 nm, which is less
than the average distance between the water molecules in liquid water
(∼0.3 nm)1 . This gives a significant probability for more than one elastic
scattering within the water molecule. This is contradictiously since the
definition of the cross section reads probability per molecule. This may
affect the track parameter calculations for very low energies, but most likely
not. Calculations of the depth dose and radial dose profiles and similar
macroscopic calculations will not be affected.
-18
10
Coherent sum
Incoherent sum
Elastic cross sections / m
2
-19
10
-20
10
-21
10
-22
10
1
10
2
10
3
10
Energy / eV
4
10
5
10
Figure 3.3: Electron elastic cross sections for liquid water calculated with the
ELSEPA code with incoherently (solid curve) and coherently (dashed curve)
summed scattering amplitudes.
3.3
Ion inelastic cross sections
The electron production in ion-atom collisions can be divided in four regions (see Figure 3.4) governed by different mechanisms in the interplay
between projectile, target and the outgoing electron. Figure 3.4 shows the
present CDW-EIS double differential inelastic cross section for a 1 MeV proton in water. The low-energy region is attributed to electrons produced in
1 David Liljeqvist, Personal communications
22
3.3. Ion inelastic cross sections
SC
2
-1
Double differential cross section / cm eV sr
-1
soft collisions (SC), which take place for a large impact parameter, i.e. distant (or glancing) collisions and include a small momentum transfer. It is
worth noting that about 50% of the electrons are produced in those distant collisions. The structure labeled ECC is associated with the process of
Electron Capture to the Continuum, where the outgoing electron has nearly
the same velocity as the projectile. This process is associated with the longrange coulombic character of both the target and projectile fields. The region between the ECC peak and the SC shoulder is attributed to two-center
electron emission (TCEE). The sharp peak at the high-end of the spectrum
is due to Binary-Encounter (BE) collision that involves a large momentum
transfer to the outgoing electrons.
10
10
-19
ECC
TCEE
-20
-21
10
10
BE
-22
10
-23
1
10
2
10
Electron energy / eV
10
3
Figure 3.4: Double differential inelastic ion cross section for 1 MeV protons in liquid water and electron ejection angle of zero degree, calculated with the CDW-EIS
approximation.
In the discussion of the emission of electrons in ion-atom collisions, it
is useful to categorize the various mechanisms behind the interactions in
terms of coulombic centers. The principle of the center picture, schematically depicted in Figure 3.5, is described in terms of the final-state interaction of the collision partners with the active electron. In the zero-center
case the trajectory of the electron is not significantly affected by the target
and/or projectile because of its high velocity. This helps the electron to escape the Coulombic field quickly. In this case, the perturbation of the electron by the target and projectile becomes small, as visibly via the prominent
BE peak. Soft collisions (SC) are a single-center mechanism where one of
the collision partners is neglected and the remaining form a center (target23
Chapter 3. Cross sections
or projectile-center). Theories including a two-center wave function in the
final state refer to the two-center case.
When an ion travels through a medium it can capture one or more electrons and becomes a dressed ion. The dressed ion can be transformed back
to its original state by electron loss processes. For a review on these competitive processes see Belkic (2010). In the case of dressed ions, the electron
bound to the ion can be active or passive. The electron is passive if it remains
in its initial ground state during the interaction process. The active-electron
case where the incident electron is removed results in a two-electronic process called two-center dielectronic process. The reader is referred to the
book written by Stolterfoht et al. (1997), where a detailed review on theoretical models and experimental methods for high-energy ion-atom collisions
is given. For more than one active target electron, higher-order centers need
to be utilized in the analysis of ion-atom impact (Belkic et al., 2008).
Zero-Center
Electron Emission
Single-Center
Electron Emission
q
q
q
Two-Center
Electron Emission
T
P
e-
q
q
q
Figure 3.5: Schematic sketch of the mechanisms for electron emission associated
with different trajectories. The zero center case refers to a weak two-body interaction between the projectile and the target electron. In single-center electron emission the electron is affected by the nuclear field of the target. In two-center electron
emission the electron is affected by the field of both the target (T) and projectile
(P). The picture can be more complex if the target atom has two or more active
electrons, the projectile is dressed or, the most obvious case, if the target is a molecule rather than an atom. The picture is a replica from the book of Stolterfoht et
al. (1997).
First Born and Binary Encounter Approximations
The Binary-Encounter Approximation (BEA) (Gryzinski, 1959; Bonsen
and Vriens, 1970, 1971) is a classical theory and is associated with the
projectile-center case and neglect the collisional effect of the nucleus,
which is called the free electron approximation. However, the BEA accounts
for the initial velocity of the electron bound to the target. The First Born
24
3.3. Ion inelastic cross sections
Approximation (FBA) treats the projectile interaction in the first order and
fully describes the interaction with the nucleus and is thus associated
with the target-center case. Due to the treatment of the interaction with
the projectile as a perturbation of the first order, the FBA is only valid
for weak projectile interactions. The free electron approximation is also
adopted in the simplified version of the Born theory, known as Plane-wave
Born approximation (PWBA). For projectiles with velocities larger than the
mean velocity of the bound electron, a violent binary collision gives rise
to a pronounced peak which can be observed in the DDCS. Energy and
momentum conservation laws predict that the BE peak has a maximum
at an energy given by E B E = 4T cos2 θ where T is the projectile energy
divided by its mass, expressed in units of the electron mass, and θ is the
ejection angle of the electron. This relation applies to an electron at rest
and gives a Dirac delta function, but due to the energy distribution of the
bound electrons it will appear as a peak with nonzero width. The FBA
describes rather well the BE process. The ECC peak is a two-center effect
and, therefore, not described by the FBA and the BE theory. The two-center
effect produces an enhancement in the electron emission at forward angles
and a reduction at backward angles. The FBA is applicable for the neutral
ion case because the ion is screened by its electron and the perturbation
becomes of short range, but is not sufficient for the dressed particles.
Theories like the Rutherford formula, which is the simplest approach to
ionization, and the FBA simply scale cross sections according to Z 2 where Z
is the charge of the projectile. This approach becomes invalid for slow projectiles. There are many analytical models for example the Kim model (Rudd
et al., 1992), Rudd (1988), Miller (1983) and Hansen–Kocbach–Stolterfoth
(HKS) (Hansen and Kocbach, 1989), but none of them describes the ECC
peak. The FBA describes a heavy particle as a plane wave in the entrance and
exit channels. The electron is described by the bound and continuum wave
functions of the target and/or projectile in the entrance and exit channels,
respectively. In the PWBA the continuum wave function for the electron in
the exit channels is replaced by a simple plane wave.
3.3.1
Continuum Distorted Wave theory
Salin (1969; 1972) introduced a rigorous theoretical description of the
two-center phenomena. In order to describe the ejected electron moving in
the presence of both the target and the projectile fields Salin allowed the
final wave function to be distorted by a multiplicative projectile continuum
factor. The Born approximation can make use of this distorted final wave
function and is known as the Distorted-Wave Born Approximation (DWBA).
If the distortion from the target and projectile Coulomb potentials are
included directly in both the initial and the final wave functions, one
gets the well known Continuum Distorted Wave (CDW) approximation of
25
Chapter 3. Cross sections
Belkic (1978) which has been particularly successful in reproducing the
many quantitative features of the two-center effects. For instance, the
ECC peak shows a more complicated dependence on the momentum in
the CDW approximation than predicted by the DWBA. It is worth noting,
that the classical-trajectory Monte Carlo method (Abrines and Percival,
1966a,b) is also a two-center theory that reproduces the ECC peak rather
well. The CDW theory compares well with experiments at high energies,
but overestimates the experiments at intermediate impact velocities. This
originates from the improper normalization of the initial wave functions.
In order to solve this problem, Crothers and McMann (1983) proposed to
replace the distorted full Coulomb wave function in the entrance channel
in the CDW method by its eikonal approximation (Fainstein and Rivarola,
1987). The resulting simplification of the CDW method is known as the
CDW-EIS where EIS stands for the Eikonal Initial State.
Implemented CDW-EIS model
In equilibrium, the charge-state fractions (Dalgarna and Griffing, 1955) are
used to model the contribution of the different states in the CDW-EIS model
used in this work. The possibility of having a negative charged ion becomes
increasingly important when the projectile energy decreases, but is disregarded in this CDW-EIS model. The interaction channel with the neutral ion
is treated using the FBA, while the other channels are described by means
of the CDW-EIS approximation. Depending on the charge of the projectile,
this approximation limits the validity of the CDW-EIS model used here. For
a proton this model could be used down to 0.1-0.5 MeV depending of the
target. Higher charged projectiles will increase the lower energy limit for
the approximation. The Complete Neglect of Differential Overlap (CNDO)
model is used for the molecular state description of the water molecule as
a target. Multiple ionizations of the target are not included in the CDWEIS model used here. The dielectronic process is included in the CWD-EIS
model. A more elaborate description of the implemented CDW-EIS cross
sections can be found in the publication by Olivera et al. (1995).
The ECC peak, which is not seen in FBA calculations, is clearly seen in the
CWD-EIS cross sections for 1 MeV proton in water depicted in Figure 3.4.
The SC, TCEE and BE structures are almost equally well described by these
two models.
26
4. Modeling of cell and tumor
response to radiation
Radiation damage cells either directly by ionizing the DNA molecules in the
cell nucleus or indirectly by forming free radicals that can interact with the
DNA. The DNA is the principal target for radiation damage, but damage to
other parts of the cell such as the cell membrane might also be of importance. DNA damages such as single strand breaks (SSB), where one of the
DNA strands are damaged, double strand breaks (DSB) where both DNA
strands are damaged or base damages could, if misrepaired, lead to senescence, apoptosis, mutation, chromosomal abberation, mitotic catastrophe
and/or genetic instability. Tissues with high proliferative capacity appear to
be more sensitive to ionization radiation. This is due to a manifestation of
the unrepaired or misrepaired damage when the cell divides. Tumors generally have higher proliferative capacity compared to normal tissue and often
deficient cell cycle blocks and repair systems leading to an accumulation of
damage until the cell cannot divide. This might actually be the explanation
why radiation therapy works.
To reach the doses needed to cure tumors without causing too severe
damage in healthy tissue, the dose normally has to be delivered in multiple
fractions. One reason for this is that the normal tissue is well organized and
has intact repair systems and therefore repairs most of the damage introduced in the DNA between the fractions, while tumors have deficient repair pathways and are not able to repair the damages. The fact that tumor
cells are often hypoxic and therefore quite radioresistant is another reason.
After one or more fractions their oxygen supply may improve because the
most radiosensitive, oxic, cells are killed and the remaining hypoxic cells are
turned into oxic cells. Redistribution is another factor that makes fractionated therapy so effective. Physical properties of the radiation like the Bragg
Peak (BP) for ions and different treatment techniques and delivery techniques give the opportunity to spare normal tissue by placing most of the
dose in the tumor. By considering the biological differences between normal tissue and tumors and physical properties of different radiation qualities an optimal treatment schedule in terms of number of fractions, time
between fractions and the dose per fraction can be determined. The dose
distribution in cells and normal tissue and the biological knowledge have to
be fused by mathematical models to improve the clinical outcome. Those
models continuously undergo improvements and novel models are being
27
Chapter 4. Modeling of cell and tumor response to radiation
developed at a fast rate. Two connected levels of radiation response are generally modeled, namely survival of cells and response of organs. Many organ
and tissue models incorporate an expression of cell survival, but there also
exists purely phenomenological models (Klepper, 2001), where no explicit
formula for cell survival is included.
The prediction of the clinical outcome of a radiotherapy treatment
in terms of tumor control and normal tissue complications is always
linked to a degree of uncertainty. Consequently, the expected outcome
of a treatment is given as a probability of having a certain treatment
effect. The uncertainty in the dose delivered to the patient stems from
the uncertainties in all steps of the treatment process. The fact that the
actual delivered dose is generally non uniform either intentionally or due
to stochastics is seldom taken into account in radiobiological models.
In the case of low-LET radiation such as photons, electrons and protons
the microscopic energy deposition is rather homogeneous on a scale
comparable to the cell nucleus. High-LET radiation, on the other hand,
deposit the energy in fewer and more densely ionizing tracks compared to
low-LET radiation, which results in an heterogeneous ’dose’ deposition
on subcellular level, i.e. ’dose’ variations from site to site within a cell. In a
similar manner, variations in sensitivity of cells within a tumor, intra tumor
variations, and variations from one tumor to another, intra patient/inter
tumor variations, can be defined and treated as stochastic effects.
Analytical expressions that consider dose inhomogeneities between cells
and between fractions were proposed in paper III and validated with numerical simulations. The concept of inhomogeneities in dose at different levels
was further explored in paper IV including repopulation and combined with
different kind of sensitivity variations.
Furthermore, a model used clinically for cell survival after high-LET irradiation, namely the LEM-model has been critically analyzed in this thesis
since it has been debated. This model only considers the difference in the
radial dose distribution for various beam qualities and does not account for
stochastic effects at cellular and sub-cellular level and is in that respect a
quite simple model.
4.1
Cell survival models
A crude approximation of cell survival after irradiation is given by the assumption that one ”hit” by radiation on a single sensitive target would lead
to cell death. Thus for each cell, by applying Poisson statistics
P (sur vi v al ) = P (0 hi t s) = e −
D/D
0
,
(4.1)
where D 0 is defined as the dose that gives on average one hit per
target. This gives an exponential cell survival curve (i.e a straight line in
28
4.1. Cell survival models
semi-logarithmic plot of survival against dose). For mammalian cells in
general, their response to radiation is usually described as ’shouldered’.
This means that at low doses, typical for fractionated treatment, the
response of cells is no longer a pure exponential but rather shows a
curvature. Historically, there have been two rather different explanations
of this phenomenon. One interpretation is that, owing to repair processes,
the cell survival curve will be pushed upwards to higher survival for low
and intermediate doses while at high doses the repair becomes less
efficient, which gives a shoulder. Alternatively, the curvature can arise
from interactions of initial lesions causing the survival curve to bend
towards lower survival. Many cell survival models have been proposed
from the sixties until today (Zimmer, 1961; Hug and Kellerer, 1966; Tobias,
1985; Curtis, 1986; Sontag, 1997). The most commonly used is the Linear
quadratic model (LQ-model) (Sinclair, 1966):
2
S(D) = e −αD−βD ,
(4.2)
where α is the initial slope of the curve and β is associated with the bending of the survival curve. The α/β quantifies the sensitivity of a tumor or a
organ to changes in the dose per fraction or dose rate. Prostate tumors have
a low α/β around 1.5 Gy and are thus sensitive to changes in the dose per
fraction, while early reacting normal tissues and many tumors have an α/β
of 10 Gy which makes them almost unaffected by varying dose per fraction.
The investigation of the effect of inhomogeneities in dose delivered for different fractionation schedules in paper III, was therefore performed for a
high and a low α/β respectively. At relatively high doses a purely exponential cell survival is expected, but due to the quadratic term in the LQ model
the resultant cell survival curve continues bending. A modified version of
the LQ model (Scholz and Kraft, 1994) that fulfills the exponential behavior
has been introduced:
(
2
e −αD−βD
D < Dt
S LQ (D) =
,
(4.3)
−αD t −βD 2t −(α+2βD t )(D−D t )
e
D ≥ Dt
where D t is the so-called threshold dose were the cell survival is transformed
into a purely exponential function. The repairable- conditionally repairable
(RCR) model (Lind et al., 2003) has the advantage of automatically transforming to an exponential function at sufficient high doses and it can also
model the low-dose hypersensitivity (Singh et al., 1994; Joiner et al., 2001)
that none of the other mentioned models can:
S RC R (D) = e −aD + bDe −cD .
(4.4)
where a, b and c are fitting parameters. In cases where no experimental data
are available or when the the cell survival needs to be characterized in terms
29
Chapter 4. Modeling of cell and tumor response to radiation
of α and β, the parameters can be calculated from following relations:
S RC R (D) = S LQ (D)
′
′
S RC
R (D) = S LQ (D)
′′
′′
S RC
R (D) = S LQ (D),
(4.5)
′
′′
′
′′
where S RC R , S LQ are the cell survival probability, S RC
R , S RC R and S LQ , S LQ
are its first and second derivative, with respect to the dose, for the RCR
and LQ model. An exact solution of α and β in terms of a, b and c can
be found from the first two relations. The inverse problem, to derive a, b
and c from α and β, can only be solved numerically and involve all three
relations. The increased interest in the use of radiation of higher ionization
densities has enforced development of models that can be used in order
to predict cell survival and organ or tumor responses for those radiation
qualities. One such model is the Microdosimetric-Kinetic-model (MKM) by
Hawkins (1994), another is the Local-Effect-Model (LEM) (Scholz et al.,
1997; Elsässer and Scholz, 2007), which is used in treatment planning for
carbon ion therapy at HIT (Heidelberg Ion Therapy Center). In fact, LEM is
the only high-LET model used in treatment planning today and has shown
good agreement with both in vitro and in vivo animal experiments with
carbon ions.
4.2
Tumor response
Clinical radiobiology is focused around the relationship between a given
physical absorbed dose and the resulting biological response and the factors
influencing the response. This is often expressed with dose-response plots,
either made from clinical data or predicted by models. Despite the fact that
the cell survival in tumors and normal tissues after irradiation is binomial
and obeys binomial statistics and approximately Poisson statistics, quite different approaches have to be used to calculate the probability of tumor control and normal tissue complications since the response of an organ is highly
dependent on its internal or functional structure. A tumor is assumed to always have a parallel structural organization since all clonogenic cells need
to be eradicated in order to control the tumor. Consequently, the control of
a tumor is expressed as a product of the probability of killing the individual
cells. Introducing the concept of Functional Sub Units (FSUs), introduced
by Withers (1992), and adding a factor that handles the functional structure
as well, such as amount of seriality, the tumor response models can also be
used for normal tissue response. The clonogenic cells are then replaced with
the FSUs and relations that handle the functional structure are added. The
probability of cell damage as a function of delivered dose at high doses often follows a sigmoid curve, and the dose response of a certain tumor- or
30
4.2. Tumor response
tissue type can be quantified by the normalized dose response gradient
at ´a
³
(D)
,
dose, D̃, where the dose response gradient is steepest, γ̃ ≡ D̃maxD dPdD
of the sigmoid curve and D 50 . The normalized dose response gradient is a
dimensionless number which describes how much the probability for a certain dose response endpoint changes per relative change in the delivered
dose, i.e. a change of 1% in delivered dose will results in an change of the
endpoint probability of γ̃ · 1%, at least in the steep quasilinear part of the response curve (Brahme, 1984). The γ̃ value is usually between 1 and 5. Thus,
the steepness of the slope of the sigmoid shaped dose response curve will
enhance the effect by any uncertainties in the delivered dose on the outcome of the treatment with a factor somewhere between 1 and 5. Therefore, high precision and accuracy in the dose delivery processes and dedicated quality assurance of the used radiation beams are extremely important
for providing well controlled radiotherapy. Generally, dose response curves
for a population of tumors are considerably less steep due to inter fraction
and/or intra fraction and/or intra cell variation in dose, see Figure 4.1, or
due to inter tumor/inter patient heterogeneity in sensitivity to radiation, see
Figure 4.2.
The response of the tumor can be calculated based on the number of
clonogenic cells, N0 , before the irradiation, the surviving fraction, S, of cells
after the irradiation and number of fractions, n, by either using Poisson statistics
n
P (D) = e −N0 S(d) ,
(4.6)
or binomial statistics:
¡
¢N
P (D) = 1 − S(d )n 0 .
(4.7)
In both paper III and paper IV binomial statistics are employed. D 50 , the
dose that gives 50% response probability is a descriptor for the position of
the dose response curve and the normalized dose response gradient, quantifies the steepness of it. For a Binomial distribution, by combining Eq. (4.1)
and Eq. (4.7) the expressions become
¶
µ
1
(4.8)
D 50 = −D 0 ln 1 − 1/N
2 0
µ
¶
1 N0 −1
γ̃ = ln N0 1 −
N0
(4.9)
and for the Poisson distribution:
D 50 = D 0 (ln N0 − ln ln 2)
(4.10)
ln N0
.
(4.11)
e
The parameters in the LQ model, α and β, can be calculated, when the α/β
is known, by utilizing the expression for D 50 for respective distribution and
γ̃ =
31
Chapter 4. Modeling of cell and tumor response to radiation
following relations:
α=
1
³
d
D 0 1 + α/β
β ==
³
1
D0 d + α
β
´,
(4.12)
´.
(4.13)
The parameters in the Poisson and the binomial models have a biological interpretation which made the models attractive and widely used in comparison to other phenomenological models such as Logit and Probit models.
Unfortunately, these parameter estimates will be influenced by biological
and dosimetric heterogeneity and cannot be regarded as realistic metrics
of some intrinsic biological properties of the tumor. Therefore, analytical
expressions, which preserve the biological coupling and that take into account biological or dosimetric heterogeneity have been derived in paper III.
Throughout the thesis the response of a cell is assumed to be determined
by the dose it receive and not by the dose any neighboring cells receives.
Thus any possible bystander effects are disregarded. The cells are assumed
to be uniformly distributed over the tumor volume, which is a simplification of the reality, but is a relatively good approximation when the dose and
sensitivity variations from cell to cell are modeled as stochastic effects.
32
4.3. Tumor response by simulations and analytical modeling
4.3 Tumor response by simulations and analytical
modeling
4.3.1
Numerical simulations
Two programs have been written in C++ for TCP calculations. The first one
was used for the calculations in Paper III and the second for calculations
performed in Paper IV. To account for inhomogeneity in dose or varying
sensitivity, the individual cells must be allowed to have different survival
probabilities, which means that every cell has to be followed through the
treatment and its fate has to be kept track of. Program 1 performs numerical
simulations of
Ã
!
N0
n
Y
Y
P=
1−
(4.14)
s c,i j ,
i =1
j =1
where s c,i j is the probability that cell i survives a dose d c,i j delivered in fraction j . The inner product is the product of surviving probabilities in each
fraction for cell i and the outer product calculates the probability to control
the tumor from the individual cells probabilities to die from the fractionated treatment. A different approach was used in program 2 for TCP calculations where the birth and death processes are sampled from a Bernoulli
distribution with probabilities p di v for division and 1 − s(d ) for cell death
respectively. The program returns either 0, if the number of cells left is larger than zero after the treatment, or 1, if the number of cells left is zero. This
process is repeated and an average is calculated in order to estimate the tumor control probability. In probability theory and statistics, the binomial
distribution is the discrete probability distribution of the number of successes in a sequence of n independent yes/no experiments, each of which
yields success with probability p. A single success/failure experiment is also
called a Bernoulli experiment or Bernoulli trial. Thus when n = 1, the binomial distribution is a Bernoulli distribution. In order to perform numerical
simulations of Eq. (4.14) several approximations have to be made to enable
inclusion of proliferation of cells. However, when Bernoulli trials are used,
no approximation regarding the proliferation is necessary. The dose variations are sampled from normal distributions and sensitivity variations are
sampled from normal or log-normal distributions in both programs using
the GSL library (Galassi et al., 2009). In addition to the intra tumor sensitivity, an investigation of inter patient variation in sensitivity has been made
in Paper IV by varying the average sensitivity of the cells in one tumor to
another, in other words, variation from patient to patient.
33
Chapter 4. Modeling of cell and tumor response to radiation
Dose variations
One could define an intra fractional macroscopic heterogeneity in dose as
the variation in dose from one cell to another or from a voxel containing
a number of cells to another voxel during the delivery of one single fraction due to the intrinsic stochastic nature of the dose delivery or an intentional inhomogeneous dose distribution. When analyzing high-LET beams
in respect to ’dose’ deposition it is necessary to define a microscopic intra
fractional and intra cell variation. For fractionated treatment schedules one
could also define an inter fractional heterogeneity in dose determined by intentional changes in the fractionation schedule and/or technical limitations
in delivering, with high precision, the same average dose per fraction to the
target. Deterministic heterogeneity refers to cases where the dose per fraction is known, i.e. when the dose per fraction is altered intentionally, while
stochastic heterogeneity refers to the possible inherent error in dose planning or delivery techniques etc and to the stochastic nature of energy deposition. The dose is assumed to follow a normal distribution for all the different
types of dose variations.
If both the dose per fraction and the dose to the cells is set to vary the
sampling of the dose is done as follows. Assume that a tumor is treated in n
fractions so that the tumor receives the average dose d t , j in fraction j . d t can
then be regarded as a stochastic variable belonging to a normal distribution
having a mean Ed t = D/n and variance Vd t , where D is the average total
dose given to the tumor during the therapy. d t , j is thus a sample from said
distribution. Assume furthermore that the tumor consists of N0 clonogenic
cells and that cell i receives the average dose d c,i j in fraction j . d c, j can then
also be regarded as a stochastic variable belonging to a normal distribution
with mean Ed c, j = d t , j and variance Vd c, j and d c,i j is thus a sample from
said distribution. Dose response curves where the dose to the cells vary and
where the dose per fraction vary and the combination of those are presented
in Figure 4.1. The results in Figure 4.1 are discussed further in section 4.3.3.
Sensitivity variations
The sampling of the sensitivity can be made in many different ways. Here
are some examples using the LQ model:
1. by sampling the α value, whenever a cell is born, from a distribution
while keeping the β value constant, which gives a variable cellular
α/β. The distribution is kept the same throughout the whole
treatment.
2. by resampling the α value between fractions for all cells present from
a distribution while keeping the β value constant, which gives a
cellular variable α/β.
34
4.3. Tumor response by simulations and analytical modeling
3. by resampling the α value between fractions for all cells present from
a distribution and recalcalculate β according to the sampled α value,
in order to keep a constant α/β.
4. by sampling the α value, whenever a cell is born, from a distribution
and recalcalculate β according to the sampled α value, in order to
keep a constant α/β.
5. by sampling the α value once for each tumor while keeping β value
constant, which gives a variable α/β from tumor to tumor.
6. by sampling the α value once for each tumor and recalcalculate β
according to the sampled α value, in order to keep a constant α/β
from tumor to tumor.
7. by sampling the α value once for each tumor, which in turn become
the mean α value of that specific tumor and the cells sensitivity is
sampled from a distribution with that α value as the mean value. The
β value can either be kept constant or be recalculated.
The mean α value of the statistical distributions is kept constant during the
whole treatment in all cases for a given tumor. The average α/β will in case
1) increase during the fractionated treatment due to the fact that the most
sensitive cells will be killed first. In case 2) the average α/β will be constant. Case 1) could reflect the situation of a hypoxic tumor with an initial
hypoxic fraction which increases progressively during the fractionated treatment due to the preferential killing of the normal cells, i.e. chronic hypoxia
not accompanied by reoxygenation. Case 2) could reflect the case of a tumor
with dynamic hypoxia under the assumption that the hypoxic fraction remains constant between fractions. Cases 3) and 4) could represent the case
of a tumor with heterogeneous intrinsic radiosensitivity (quantified as SF 2),
but with the same fractionation sensitivity (quantified as α/β). Cases 5) and
6) reflect inter patient variation in sensitivity. Case 1) and 2), describing intra tumor variation in sensitivity, were simulated and discussed in Paper IV
together with case 5) describing inter patient variations in sensitivity. Some
results for the cases 3) and 6) are given in Figure 4.2 and case 7) is planned
to be simulated in the future. The largest possible tumor that can be simulated is a tumor with 107 cells in absence of proliferation due to limitations
in computational power. It is very time consuming to simulate that tumor
size, though, so a smaller tumor size with 104 cells is chosen for the simulations. To maintain the intrinsic sensitivity of the cells in the smaller tumor
the D 50 is chosen to be 29 Gy, which corresponds to a D 50 of 50 Gy for the
larger tumor.
There is a lack of radiobiological insight in the variation of α and β in
relation to cell sensitivity variations. What does best characterize the cell
sensitivity variation within the cell and between tumors? Should the α/β be
35
Chapter 4. Modeling of cell and tumor response to radiation
kept constant or not? With a survival model that separates damage induction and repair, like the RCR model, the repair and damage induction can
be varied separately, which gives more control and insight in how sensitivity variations affect the outcome. However, it is not enough to know which
parameters to vary, one should also know how those parameters are distributed. Tome and Fowler (2002) showed a minimal impact on the TCP when
using a Log-Normal distributed SF 2 instead of a truncated Normal distribution for the SF 2, where SF 2 is the surviving fraction of cells after a 2 Gy fraction and can be represented as a parameter describing the clonogen sensitivity. The validity of the findings of Tome and Fowler (2002) for the current
study was performed by comparing the TCP resulting from a sampling of the
α value from a log-normal and from a Normal distribution and no significant differences were seen. The Log-Normal distribution has been assumed
to be the more realistic one and therefore used in all the calculations in the
present work.
For the intra tumor case, a constant α/β, i.e. where β is recalculated according to the sampled α value, gives a slight translation of the dose response curve towards higher doses compared to the case where α/β is allowed to vary, i.e. β is kept constant. The inter patient simulations where
the β value is recalculated, to maintain a constant α/β, gives a less steep
dose response curve, i.e. a lower γ̃, compared to the case where the α/β is allowed to vary from one tumor to another. The differences for the two ways of
sampling the sensitivity are substantial especially for the inter patient case.
All the dose response curves are calculated with 2 Gy fractions. When β is
α
where α is the sampled alpha value,
recalculated by the relation β = α/β
the TCP become independent of the α/β of the tumor because the SF 2 becomes equal for different α/β not only for a constant α value (CV = 0) but
also when the α value is sampled from a distribution.
Effects on the TCP from combinations of variation in the sensitivity, either
inter tumor or intra tumor, and variation in the dose to the cells were investigated in Paper IV. The effect on TCP from inter tumor and dose variations
add up, but the result is less than the sum of the respective effects. In the
intra tumor case the effects seems to simply add up. Sensitivity variations
will never play a positive role on TCP. This can be understood from the derived analytical expressions in Eq. (4.26) for sensitivity variation, where it is
obvious that the sensitivity among cells always play a negative role.
36
4.3. Tumor response by simulations and analytical modeling
1.0
0.9
0.8
TCP
0.7
0.6
Homogeneous dose
0.5
Heterogeneous cell dose
Heterogeneous cell and fraction dose
0.4
Heterogeneous cell and site dose
0.3
Heterogeneous fraction dose
0.2
Heterogeneous site dose
0.1
Heterogeneous site and fraction dose
Heterogeneous site, fraction and cell dose
0.0
20
25
30
35
40
45
Dose / Gy
Figure 4.1: Dose response for different cases of dose variations with a coefficient
of variation CV = 20% for a tumor with α/β = 1.5 Gy, 104 cells, D 50 = 29 Gy and
p d i v = 0. Cell dose refers to heterogeneous dose distribution between the cells in
the tumor and site dose refers to heterogeneous dose distribution within the cells.
Fraction dose refers to the case where the fractional dose is varied from one fraction
to another.
1.0
0.9
0.8
0.7
TCP
0.6
0.5
0.4
Homogeneous sensitivity
0.3
Interpatient, varying α/β
Interpatient, constant α/β
0.2
Intra tumour, varying α/β
Intra tumor, constant α/β
0.1
0.0
15
20
25
30
35
Dose / Gy
40
45
50
Figure 4.2: Dose response for different cases of sensitivity variations (CV = 20%)
for tumors with α/β = 1.5 Gy, 104 cells, D 50 = 29 Gy and p d i v = 0. Interpatient,
referes to sensitivity variation from one tumor to another and intra tumor refers to
variation in sensitivity within the tumor. Sensitivity variation means that α is varied,
whereas β could either be kept constant resulting in a varying α/β or recalculated
resulting in a constant α/β.
37
Chapter 4. Modeling of cell and tumor response to radiation
Repopulation
Repopulation is often taken into account in a rough manner by using a multiplicative factor to the initial number of clonogens, N0 , in the Poisson or the
binomial expressions for TCP according to
n
P (D) = e −R N0 S (d)
¡
¢R N0
P (D) = 1 − S(d )n
,
(4.15)
N = N0 (1 + p di v ) /T ,
(4.17)
(4.16)
where R = exp(ln2t /TD ), TD is clonogen doubling time, and t is the overall
treatment time in days. This approximation is relatively good for tissues with
low proliferate rate (high TD ) but becomes invalid for a low TD , especially
when the probability for cell survival is high. In iterative calculations of the
TCP as in the work by Deasy (1996), which was employed in Paper IV, or in
MC simulations of the TCP, made in Paper IV, the proliferation can be taken
into account in a more precise way if the probability of the cells to divide,
p di v , is known. This probability can be derived in terms of the clonogen
doubling time acquired from laboratory measurements.
Assuming that the probability, for a single cell to divide during a time interval T is p di v , the total number of cells, N , after t /T time intervals is
t
where N0 is the initial number of cells at t = 0. For a known doubling time
TD , the probability for cell division is
T/T
p di v = 2
D
−1 = e
T/T
D ln2
− 1.
(4.18)
4.3.2 Analytical modeling -approximation of an expectation
value
For simple functions of a stochastic variable it is easy to directly calculate the
expectation value and the variance, especially when the variables are independent of each other. However, for many functions (non-linear functions)
an exact determination of those quantities can be cumbersome. A common
way to approximate the mean and the variance of a function of stochastic
variables emanates from Gauss. This method is especially useful in certain
areas of science such as medical radiation physics, where the need to estimate stochastic variables, such as the dose, is large. The Gauss approximation
formulas for a function, g, of a stochastic variable, X, is of first order and
given by
£
¤
E g (X ) ≈ g [E(X )] ,
(4.19)
£
¤
£ ′
¤2
V g (X ) ≈ V(X ) g (E(X )) ,
(4.20)
£
¤
£
¤
where E g (X ) is the expectation value and V g (X ) the variance. Usually
the first order is quite sufficient, but in many cases higher order terms contain information that is of importance. Mekid and Vaja (2008) and Wang and
38
4.3. Tumor response by simulations and analytical modeling
Iyer (2005) have made attempts to include higher orders for the average and
variance but address this as propagation of errors and uncertainty in measurements, which is the conventional way to entitle the problem. Actually,
the second of the Gauss approximation formulas is also known as the law of
error propagation. In this thesis, an extension of the first of Gauss approximation formulas has been used to model the effect of uncertainty in average
dose, either intentionally introduced or the inherent uncertainty in dose at
different levels of tissue structures for various beam qualities. It could also
be used to model the sensitivity variation from cell to cell within a tumor.
Assume a vector of stochastic variables ~
X and let g ( ~
X ) be a function of
the vector of stochastic variables, whose first and second derivatives exists.
Expanding g ( ~
X ) in a Taylor series around the average of ~
X up to the second
order gives
1
X ) + (~
X −~
X )T ∇g ( ~
X ) + (~
X )T ∇2 g ( ~
X )( ~
X −~
X ).
X −~
g (d~) ≈ g ( ~
2
(4.21)
The expectation value of g ( ~
X ) then becomes
¡ ¢
1
Eg ( ~
X ) + ∇T C ~
X ).
X ∇g ( ~
X ) ≈ g (~
2
(4.22)
where C( ~
X ) is the covariance matrix.
Explicit calculations in 2D for the stochastic vector ~
X = (X 1 , X 2 )T with avX = (X 1 , X 2 ) is shown as an example:
erage ~
g (~
X) =
g (X 1 , X 2 ) + (X 1 − X 1 , X 2 − X 2 ) · ∇g ( ~
X )+
"
#
X1 − X 1
1
2 ~
(X 1 − X 1 , X 2 − X 2 )∇ g ( X )
.
2
X2 − X 2
(4.23)
Performing the matrix and vector multiplications and taking the expectation
value of the resulting expression:
£
¤
E g (~
X)
=
h
i ∂g ( ~
X)
g (X 1 , X 2 ) + E (X 1 − X 1 )(X 1 − X 1 )
+
∂X 12
h
i ∂g ( ~
X)
+
∂X 1 ∂X 2
h
i ∂g ( ~
X)
E (X 2 − X 2 )(X 1 − X 1 )
+
∂X 2 ∂X 1
h
i ∂g ( ~
X)
E (X 2 − X 2 )(X 2 − X 2 )
.
∂X 22
E (X 1 − X 1 )(X 2 − X 2 )
(4.24)
Note
that the second
term in Eq. (4.23) vanishes. Identify
i
h
E (X 1 − X 1 )(X 1 − X 1 ) as the covariance matrix element C 11 etc.
39
Chapter 4. Modeling of cell and tumor response to radiation
Assume that the covariances are constants and rewrite Eq. (4.24) in terms
of covariance:
#"
#
µ
¶"
∂
£
¤
C 11 C 12
∂
1 ∂
∂X 1
~
,
g (~
E g ( X ) = g (X 1 , X 2 ) +
X)
∂
2 ∂X 1 ∂X 2
C 21 C 22
∂X 2
=
1
X )∇g ( ~
g (~
X ).
X¯ ) + ∇T C( ~
2
(4.25)
This second order of the mean is very seldom found in the literature and certainly not in the area of medical radiation physics. In Paper III it was used
to calculate the TCP with the dose in each fraction as a stochastic variable
or the dose to each cell. In Paper IV it was used to calculate the TCP with
the cell sensitivity as a stochastic variable. Starting with Eq. (4.14) and making use of Eq. (4.22), the analytical expression for TCP with a dose variation
or sensitivity variation from cell to cell, expressed in terms of the standard
deviation of either the sensitivity or the dose, σc , becomes:
" ′
#)
(
nN0 2 s 2 s 2n−2 s ′′ s n−1
σc
.
(4.26)
+
P ≈ P (D) 1 −
2
(1 − s n )2 1 − s n
where s are the cell survival probability and s ′ and s ′′ its first and second derivative, with respect to either the dose to cell or the cell sensitivity respectively, n is the number of fractions and N0 the initial number of clonogenic
cells. In the case of the dose in each fraction as a stochastic variable the analytical expression, defined in terms of the standard deviation of the dose per
fraction, σt , becomes:
µ
µ
µ ′′
¶¶¶
n
s
s ′2
P ≈ 1 − s n 1 + σ2t
− 2
,
2
s
s
(4.27)
The LQ model was used to calculate the cell survival probability in Paper III
and IV, but any suitable model could be used.
Varying dose per fraction
In Paper III the total biological effective dose was calculated in order to obtain a tumor control probability of 50 % when the dose was homogeneously
delivered to the whole tumor using following relation:
nα/β
+
D =−
2
s
µ
nα/β
2
¶2
¡
¢
+ nD 50 2 + α/β .
(4.28)
The reason for this calculation was to ease the comparison of the effect of
inhomogeneous doses as a function of number of fractions. This approach
gives a large variation in dose per fraction as the number of fractions
changes. For tissues sensitive to changes in fractionation, i.e. late reacting
tissues characterized by low α/β values and pronounced shoulder of the
40
4.3. Tumor response by simulations and analytical modeling
cell survival curve, the dose per fraction become especially important. The
dose and cell survival relation is not linear in the shoulder region, i.e. the
increased cell kill due to doses higher than the average dose can more than
compensate for the decreased cell kill due to lower doses than average.
An alternative approach to ease the comparison could be to vary D 50 and
keep the dose per fraction and N0 constant. Figure 4.3 shows a comparison
of calculations of Eq. (4.26) using the two approaches together with a
numerical simulation of Eq. (4.14) with constant D 50 and varying dose
per fraction. If D 50 is varied while keeping the dose per fraction constant,
the TCP will increase relatively to the calculations of the TCP using the
approach in Paper III. There is one exception, at 25 fractions the dose per
fraction is 2 Gy for both cases, resulting in equal TCP.
0.6
0.5
α/β =1.5 Gy
TCP
0.4
0.3
0.2
0.1
0
0
10
20
30
Number of fractions
40
50
60
Figure 4.3: TCP as a function of number of fractions for a tumor with α/β of 1.5 Gy
and 107 initial number of cells with a CV = 20% for the dose to the cells. Solid curve
is calculated with the analytical expression given by Eq. (4.26) with a constant D 50
and varying dose per fraction, dashed curve is calculated with varying D 50 , but constant dose per fraction. The squares are numerical simulations of Eq. (4.14) with
constant D 50 and varying dose per fraction.
4.3.3 The influence of microscopic heterogeneity in energy deposition
As mentioned in the section Dose variations, a high LET beam will have a
microscopic intra fractional intra cell variation in dose. Assuming that there
are several radiation sensitive sites inside the cell, this dose variation to the
sites within a cell may affect the probability of the survival of the cell. How
the site inactivation is coupled to cell inactivation is not exactly known. If
41
Chapter 4. Modeling of cell and tumor response to radiation
for example every site is assumed to be essential, inactivation of a single site
is enough to kill the cell. An alternative proposal is that at least two sites are
needed in order to inactivate the cell (Brahme and Lind, 2010). It could be
an even more complex situation, where the location of the inactivated sites
in relation to each other is important. How to model the site response as a
function of dose to the site is still an open question. Below, one relatively
simple way to model the site response is presented, which assumes that the
cell is a serial structure, i.e. the cell dies as soon as one critical site is inactivated.
Assume that there are m quasi uniformly distributed critical sites in each
cell nucleus and that site k in cell i receives the average dose d s,i j k in fraction j. d s,i j can then also be regarded as a stochastic variable belonging to a
normal distribution with mean Ed s,i j = d c,i j and variance Vd s,i j and d s,i j k
is thus a sample from this distribution.
Provided that a cell will die as soon as at least one critical site is inactivated, the probability that cell i survives a single dose delivered in fraction j
will be given by
m
Y
s c,i j =
s s,i j k ,
(4.29)
k=1
where s s,i j k , is the probability that critical site k in cell i survives a dose
d s,i j k delivered in fraction j . The probability that cell i is still alive after n
fractions is then
n
n Y
m
Y
Y
s c,i =
s c,i j =
s s,i j k .
(4.30)
j =1
j =1 k=1
Assuming that every clonogenic cell has to be eradicated in order to control
the tumor, the probability of tumor control will be
Ã
!
N0
N0 ¡
n Y
m
Y
Y
¢ Y
1−
1 − s c,i =
P=
s s,i j k .
(4.31)
i =1
i =1
j =1 k=1
For completely homogeneous dose delivery, i.e. every site and cell receives
the dose d s,i j k in every fraction, the probability of tumor response becomes
³
´N0
n·m
P hom = 1 − s s,i
.
jk
(4.32)
Solving Eq. (4.32) for s s,i j k and substituting back into Eq. (4.31) gives
P=
N0
Y
i =1
Ã
1−
n Y
m ³
Y
j =1 k=1
1/N
0
1 − P hom
´1/n·m
!
.
(4.33)
The standard Binomial expression for tumor control as a function of cell survival when the cell receives the homogeneous dose d c,i j k in every fraction
is
´
³
n
P hom = 1 − s c,i
jk
42
N0
.
(4.34)
4.3. Tumor response by simulations and analytical modeling
Inserting Eq. (4.34) in Eq. (4.33) gives
Ã
!
N0
n Y
m
Y
Y
1/m
P=
1−
s c,i j k .
i =1
(4.35)
j =1 k=1
Eq. (4.33) has been simulated numerically and the results are presented in
Figure 4.1. It could be discussed what the targets in the cells are, i.e. what is a
site. The cell nucleus contain around 30 000 genes distributed in 46 chromosomes. Either of those could be defined as a site. There could be subregions,
other than the genes or the chromosomes, that are the critical target within
a cell. In the simulations the number of sites is chosen to be 100, but 50 and
1000 were tested without any effect on the result. It would be much more interesting if the cell inactivation could be modeled to follow after inactivation
of at least two sites or even more complex relations. A more comprehensive investigation of the interplay between the intrafractional macroscopic
heterogeneity and intra fractional microscopic heterogeneity is of great importance, where the site response is accurately modeled and the relation
between site inactivation and cell death is more realistic.
The assumption made in the derivation of the response of cells including
sites, i.e. that the response of the cell depends upon the dose deposited to
subunits of the cell and that those subunits are arranged in a serial manner,
results in the LEM expression presented by Scholz et al. 1997, as shown in
section 4.4
Figure 4.1 shows that the effect on the TCP from combining different dose
variations add up. Variation in fractional dose and in site dose acts positively on the TCP, while variation in cell dose acts negatively on the TCP for
this particular case. For very resistant cells the variation in dose can in fact
increase the TCP, see Paper IV. Variations in the fractional dose can never act
negatively on TCP and can be understood by analyzing Eq. (4.27). By studying Eq. (4.31) one can observe that variation in the dose to the sites should
influence the TCP in the same manner as a variation in the dose per fraction,
i.e. only positively.
43
Chapter 4. Modeling of cell and tumor response to radiation
4.4
The Local Effect Model
The principal assumption in the LEM is that the target, i.e. the cell nucleus,
can be divided into subvolumes, and that the local biological effect depends
on the local dose, i.e. the microscopic dose deposited in a subvolume, but
not on the radiation modality giving the dose. With this assumption, the cell
survival curve achieved from X-ray irradiation can be utilized in the calculation of biological effect of ions. The modified linear quadratic (LQ) model
is used for the description of X-ray cell survival curve. A radial dose parametrization for the ions is used as input to the LEM model, which in essence makes it an amorphous track structure model. In the original derivation of the LEM, Scholz et al. (1997) assumed Poisson distributed lethal
events. Here, an alternative derivation of the LEM expression, without the
initial assumption of Poisson statistics, will be given, similar to the derivation of Beuve (Beuve, 2009). Furthermore, an investigation of the β dependence in the model and the ability of LEM to predict the cell survival using
different approaches to introduce the shouldering seen in cell survival has
been made. Henceforth, a first attempt to implement the RCR in the framework of LEM has been made.
In the case of low LET irradiation, such as bremsstrahlung or electron
beams, the microscopic energy deposition is rather homogeneous on a scale
comparable to the cell nucleus. The coefficient of variation of the dose to
the cell nucleus (∼8 µm diameter) is generally below 1% at curative doses.
If however a small non uniformity is introduced e.g. by mixing a high and
low LET beam, the dosimetric variance is increased. In order to quantify the
radiation effect on the subcellular level the following two assumptions are
made:
1. The critical targets for cellular clonogenic capacity are quasi
uniformly distributed throughout the cell nucleus.
2. The cell lose its proliferation capacity if one or more of these critical
targets are inactivated.
If the absorbed ‘dose’ to critical target i is denoted by D i and the corresponding inactivation probability by P i (D i ), then the total probability to inactivate at least one out of n such critical targets is given by
~ = 1−
P (D)
n
Y
(1 − P i (D i )) ,
(4.36)
i =1
~ is the total microscopic dose distribution vector D
~ = (D 1 , D 2 , · · · D n )
where D
describing the dose distribution inside the cell nucleus. If the cell nucleus
instead is irradiated with a homogeneous dose D also on the microscopic
level, the total inactivation probability simply becomes:
P (D) = 1 − (1 − P i (D))n ,
44
(4.37)
4.4. The Local Effect Model
since all critical targets are assumed to have the same sensitivity. Eq. (4.37)
can be solved for P i (D)
P i (D) = 1 − (1 − P (D)) /n .
1
(4.38)
Substituting Eq. (4.38) back in Eq. (4.36) gives
~ = 1−
P (D)
n
Y
(1 − P (D i )) /n ,
1
(4.39)
i =1
where P (D i ) thus is the probability to inactivate the cell when receiving the
absorbed dose D i .
Since the probability of survival is given by S(D) = 1−P (D), the total probability of cellular survival when exposed to an arbitrary non-uniform dose
distribution is given by
)
(
n
n
Y
Y
1/n
~ = 1 − P (D)
~ = 1 − 1 − (1 − (1 − S(D i )))
=
S(D i )∆νi ,
(4.40)
S(D)
i =1
i =1
where S(D i ) is the survival probability when the cell nucleus receives the
homogeneous doses D i and ∆νi is the fractional volume of the i t h voxel. By
using the equality S = exp(ln S) Eq. (4.40) can be rewritten as a sum according to
Ã
!
n
X
~ = exp
S(D)
∆νi ln S(D i ) .
(4.41)
i =1
In the limit of infinitesimally small voxels D i becomes the quasi continuous
dose distribution D(~
r ) and Eq. (4.41) can be rewritten as an integral:
~ = exp
S(D)
µ
1
V
ˆ
3
¶
ln S x (D(~
r ))d r ,
V
(4.42)
where V is the volume of the cell nucleus and S x (D(~
r )) is the cell survival
after exposure to a dose D by X-rays.
Scholz et al. (1997) presented one simplified approach for the dose distribution in the cell nucleus where only central traversals of the ions occur, and
one rather time consuming calculation of the dose distribution in the target
by superimposing the radial dose distribution from ions with randomized
impact parameters (MC calculated ion tracks). The differences between the
two approaches are negligible (Scholz et al., 1997). For simplicity, the central traversal approach has been adopted in the present study. Since, in most
cases, the energy loss during the traversal of the cell nucleus can be neglected and by using cylindrical coordinates the 3D integral in Eq. (4.42) can
be written as
1
V
ˆ
V
3
ln S x (D(~
r ))d r
=
1
V
ˆ
V
45
ln S x (D(r, θ, z))r dzdθdr
Chapter 4. Modeling of cell and tumor response to radiation
ˆrnucˆ2πˆs̄
1
V
=
0
2πs̄
V
=
0
ln S x (D(r, θ, z))r dzdθdr
0
ˆrnuc
ln S x (D(r ))r dr
0
ˆrnuc
2π
σ
=
ln S x (D(r ))r dr,
(4.43)
0
where s̄ is the mean chord length of the ion track passing the nucleus and
σ and r nuc are the nucleus cross sectional surface and radius respectively.
If the nucleus is assumed to have a cylindrical shape, s̄ will simply be the
height of the cylinder and r nuc the radius of the cylinder.
The radial dose distribution from an ion can be approximated with the
following expression
D(r ) =

2


 λ/r mi n



λ/r
r < r mi n
2
r mi n ≤ r ≤ r max ,
0
(4.44)
r > r max
where r mi n = 10 nm and r max = γE δ , γ = 6.2 · 10−8 , is in m/MeV δ = 1.7,
where r max is in m and E in MeV/u. All the parameter values were taken
from Scholz et al. (1997).
For consistence, and by the very definition of the quantity L, the Linear
Energy Transfer, the following relation can be used to calculate the normalization constant λ
L
= 2π
ρ
ˆ∞
r D(r )dr.
(4.45)
0
The radial dose distribution in Eq. (4.44) is the average microscopic absorbed dose distribution for a single ion. If the medium is irradiated with
an average (macroscopic) dose D̄ the resultant number of particles per unit
surface area (the fluence Φ) becomes Φ = D̄ρ/L. The average number of ions
passing the cell nucleus, ν, is thus ν = σD̄ρ/L, where σ is the cross sectional
surface of the cell nucleus. Since the radial dose in Eq. (4.44) should be multiplied with ν in order to achieve the correct average dose to the medium
the L/ρ terms cancel out and the final expression will be independent of the
LET (only the particle energy matters).
For a given cell line the parameters α and β have to be known in order
to calculate S x (D(~
r )). For a given average dose to the medium and particle
energy the cell survival can be calculated by the numerical solution of Eq.
(4.42) above.
46
4.4. The Local Effect Model
4.4.1 Three approaches to the calculation of cell survival in the
framework of LEM
The cell survival values for ion beams is not simply deduced from α and β
from X-ray irradiation on the cells of interest and the radial dose distribution
of the ion, in contrast to what might generally be believed. Instead a circumstantial fitting and calculation procedure is used to get the cell survival from
the ion beam. This procedure complicates the introduction of other cell survival models such as RCR model in the framework of LEM. Due to those two
quite substantial drawbacks an attempt to derive two alternative ways to
take advantage of the X-ray data has been made. In the first approach it has
been assumed that the effect of each ion traversal on cell survival is independent from other ion traversals, i.e. non overlapping tracks, whereas in
the second approach a dependency between ion traversals leading to a synergistic effect was assumed, i.e. totally overlapping tracks.
In the two approaches it is assumed that the average number of ions, ν
passing a cell nucleus of cross sectional area σ is given by ν = σΦ, where Φ
is the number of ions per unit area, i.e. the fluence. Both approaches start
with the assumption of Poisson statistics; the probability that exactly ν ions
cross the nucleus is thus
e −σΦ (σΦ)ν
Pν =
.
(4.46)
ν!
First approach
Assume that the probability of cell survival after the passage of a single ion
through the cell nucleus is S 1 , calculated by Eq. (4.42). The probability of cell
survival with a fluence Φ of ions then becomes
∞ e −σΦ (σΦS )ν
X
1
.
(4.47)
S(Φ) =
ν!
ν=0
P
x
xn
But since ∞
n=0 /n! = e , Eq. (4.47) can be written as
S(Φ) = e −σΦ(1−S 1 ) .
(4.48)
The exponent σΦ(1 − S 1 ) in Eq. (4.48) is accordingly the mean number of
lethal events from the ion with fluence Φ. Since the absorbed dose D is given
by D = ΦL/ρ, where L is the LET and ρ the density of the medium, Eq. (4.48)
becomes
ρ
S(D) = e −σ L D(1−S 1 ) .
(4.49)
Second approach
Assume that the probability of cell survival after the passage of ν ions
through the cell nucleus is S ν (νd ), where d is the average dose in the
nucleus from a single ion passage. The probability of cell survival with a
fluence Φ of ions then becomes
∞ (σΦ)ν
X
S ν (νd ).
(4.50)
S(Φ) = e −σΦ
ν=0 ν!
47
Chapter 4. Modeling of cell and tumor response to radiation
But since the absorbed dose D is given by D = ΦL/ρ, where L is the Linear
Energy Transfer and ρ the density of the medium, the average number, ν, of
ion traversals becomes
ρ
ν = σ D.
(4.51)
L
Eq. (4.50) can thus be rewritten as
ρ
∞ (σ D)ν
X
L
ρ
S(D) = e −σ L D
ν=0
= e
ρ
−σ L D
ν!
S ν (νd )
µ
¶
ρ
1 ³ ρ ´2
S ν (0) + σ DS ν (d ) + σ D S ν (2d ) + · · · . (4.52)
L
2 L
The cell survival S ν (νd ) should be calculated from the “LEM” expression
according to
µ ˆ
¶
1
3
S ν (νd ) = exp
ln S x (νD(~
r ))d r ,
(4.53)
V V
where D(~
r ) is the lateral dose distribution from a single ion. The x-ray cell
survival S x (νd ) can be e.g. the modified LQ expression
S LQ (νd ) =
(
2
e −ανd−β(νd)
e
νd < D t
−αD t −βD 2t −(α+2βD t )(νd−D t )
νd ≥ D t
,
(4.54)
where the α and β should come from low LET data such as X-rays or alternatively the RCR model can be used
S RC R (νd ) = e −aνd + bνd e −cνd .
(4.55)
Original approach
Assume only central traversals for simplicity. In the limit of small doses one
can express the cell survival for one central traversal as
S 1 = e (−
αz σDρ/L )
,
(4.56)
L
can be deduced. Next
from which the intrinsic parameter αz = −ln(S 1 ) σρ
step is to derive a phenomenologically βz , which resembles the interaction
of some fixed number of impinging ions. The maximum slope of the cell
survival curve, s max , can be found at the dose D t for all energies:
dln Sx
= s max = αx + 2D t βx ,
dD
(4.57)
where αx and βx are from the X-ray data. This relation should hold for the
βz and αz as well:
dln S
= s max = αz + 2D t βz ,
(4.58)
dD
which gives βz = (s max − αz )/D t . The macroscopically observed quantities
α and β can be calculated from the intrinsic parameters αz and βz by the
48
4.4. The Local Effect Model
introduction of a correction for hit statistics. This transformation is necessary as the effect from exactly ν particles is not the same as the effect from a
average of ν particles. In the limit of small doses D, it is sufficient to consider
the probability of no traversal and one traversal which yields the relation:
³
ρ
ρ ´
S(D) = 1 − σ D + S 1 σ D.
L
L
(4.59)
By making use of that the derivative of the logarithm of the survival is just
the derivative of the survival in the limit of small doses, the initial slope is
calculated by:
ρ
(4.60)
α = (1 − S 1 )σ .
L
This is exactly the same as the first approach described above. Finally, the
macroscopic β is calculated by:
β=
µ
α
αz
¶2
βz .
(4.61)
The cell survival is then calculated by the LQ-model with these derived α
and β. Recall that D t is a fitting parameter and s max is taken from a calculation of the cell survival for X-rays including this fitting parameter and the
macroscopic β is derived from those parameters and scaled due to hit statistics. The strategy described will introduce an artificial shoulder and the
procedure deviates from the original concept of calculating the cell survival
from the X-ray data for the specific cell type.
4.4.2
Results of the different approaches
The motivation for the derivation of the 1st and 2nd approaches was to circumvent the illogical recalculation of α and β and restore the basic thought
behind LEM model and open it up for an implementation of other cell survival models in the LEM model. The RCR model was chosen as an alternative
to the LQ model in LEM for its capability to model hypersensitivity and for
its feature of transcending the shoulder region of the survival curve into a
linear curve at high doses, which does not involve an extra parameter such
as D t
The in vitro experiments performed on Human melanoma cell line (AA)
under aerobic conditions at TSL laboratory in Uppsala and Karolinska
University Hospital in Stockholm were chosen as a point of reference for
the studies here. The details about the experiments are given in Persson
et al (2002) and Meijer et al. (2005). Figures 4.4-4.6 show the results
from the different approaches with both LQ and RCR as cell survival
models for Boron ions with LET 40, 80 and 160 eV/nm with corresponding
energies of 36 MeV/u, 16 MeV/u and 6.4 MeV/u. The LQ parameters
αx = 0.445 Gy−1 and βx = 0.444 Gy−2 for X-ray irradiation are taken from
49
Chapter 4. Modeling of cell and tumor response to radiation
1
Cell survival fraction
0.1
0.01
Experiment
Original approach
2nd approach + LQ
2nd approach + RCR
1st approach + LQ
1st approach + RCR
0.001
0
2
4
Dose / Gy
6
8
Figure 4.4: Cell survival as a function of dose for Boron ions with LET 40 eV/nm.
maximum likelihood fit of the LQ-model to the experimental data. The
same procedure is used to find a = 1.1387 Gy−1 , b = 1.2317 Gy−1 and
c = 0.6617 Gy−1 for the RCR model. to be 3.24 µm , but a small change
in this parameter will not affect the outcome of the calculations. After a
logarithmic transformation, the total least square deviation for all three
data sets was minimized with D t as a a free parameter. The original and
1st approach got the optimum accordance with the experiments by a
D t of 35 Gy and 2nd approach by a D t of 25 Gy. The 1st approach with
LQ-model was found to have best accordance with experimental data
followed by original approach. Worst accordance gave the 1st approach
with RCR-model followed by 2nd approach with RCR-model. The advantage
of the LEM model in spite of its poor agreement with the experimental
data, at least with the present data set, is the opportunity to predict the cell
survival from ions with other LETs than those measured by interpolating
the results. The main reason for the poor agreement with experimental
data for the approaches using RCR-model results from the fact that the
extra fitting parameter D t was not used. By calculating a macroscopic α
and β, as described above, the introduction of other cell survival models
is almost impossible. It is also unnesessary as the success of the LEM
model with LQ as the cell survial model is just a coincidence and a result
50
4.4. The Local Effect Model
1
Cell survival fraction
0.1
0.01
Experiment
Original approach + LQ
2nd approach + LQ
2nd approach + RCR
1st approach + LQ
1st approach + RCR
0.001
0
2
4
Dose / Gy
6
8
Figure 4.5: Cell survival as a function of dose for Boron ions with LET 80 eV/nm.
of the fitting parameter D t that happens to lower the cell survival at doses
typical for ion impact, in such a way that it becomes in balance with the
macroscopic α and β. Figures 4.4 and 4.5 show a small shoulder for the
original approach and for the 2nd approach with LQ-model. This shoulder
is introduced artificially in some sense in the original approach by the
calculation of a macroscopic β, while in the 2nd approach its a synergistic
effect of overlapping ion tracks. Probably the LET parametrized RCR model
(Wedenberg et al., 2010) or other LET parametrized cell survival models are
a better choice than the LEM model.
51
Chapter 4. Modeling of cell and tumor response to radiation
1
Cell survival fraction
0.1
0.01
Experiment
Original approach + LQ
2nd approach + LQ
2nd approach + RCR
1st approach + LQ
1st approach + RCR
0.001
0
2
4
Dose / Gy
6
8
Figure 4.6: Cell survival as a function of dose for Boron ions with LET 160 eV/nm.
52
5. Beam quality characterization
Ions show an increased biological effectiveness compared to conventional
radiotherapy such as X-rays and electrons, but exactly how the biological effect, for various endpoints, varies with different radiation qualities has not
been uniquely defined. Radiation quality is dependent upon all intrinsic
quantities such as particle type and energy, LET, spatial pattern of energy
depositions and lateral and longitudinal distribution of dose, that affect the
biological effect.
The concept of absorbed dose in radiation therapy as defined by ICRU,
i.e. the mean energy imparted per unit mass, is a sufficient measure of the
effect of ionizing radiation as long as the ionization density is low enough.
For higher ionizing densities, this is not the case. Thus, additional parameters are needed to characterize the beam quality. The absorbed dose is
a macroscopic concept and therefore not suitable when energy depositions
on cellular and sub cellular levels are of importance. Historically, LET has
been favored as such a parameter for the characterization of beam qualities. There are however several problems with LET in this context. Perhaps
the main limitation is that different charged particles will not have the same
biological effect even though they have the same LET. A substantial amount
of the particle energy is carried away by delta rays and therefore should perhaps other parameters such as SPI, the number of primary ionization’s per
path length (in essence the inelastic cross section), (Cannell and Watt, 1985)
or SI, the number of ionizations per path length, or even more refined parameters be used. The latter was actually used before LET, but was discarded
in favor of LET. Therefore, when studying the biological outcome after irradiation with ions, not only the quantity absorbed dose but also the local
energy spectrum of secondary particles, mainly electrons, have to be considered. In this chapter an ongoing investigation of alternative measures to
LET, suitable for estimation of biological effect, will be presented with preliminary results and discussion.
Since the pioneering work by Lea (1962), which focused on the energy
depositions in specific cellular targets even before the characteristics and
the role of the DNA were known, there has been considerable interest in
relating the results of track structure calculations for different radiations and
doses to corresponding measured biological effects.
The improvements in both modeling the pattern of energy depositions
by analytical calculations, or MC track structure calculation, and the target
53
Chapter 5. Beam quality characterization
DNA have been many during the years, but still there is lot of work to be
done. Radiobiological models, like the LQ-model, use the dose as the free
parameter of radiation quality, other use both dose and LET, for instance
the MKM model (Hawkins, 1994) and the LET-parametrized RCR model
(Wedenberg et al., 2010).
The overall goal is to find a measure that uniquely defines the radiation
quality and can be beneficial to use in such radiobiological models. To oversimplify, the ultimate aim is to achieve a survival curve from the relation:
S = e −Nl ethal ,
(5.1)
where Nl et hal is the mean number of lethal events produced in the cells
by a specific beam quality. The initial lesions formed in DNA by ionizing
radiation include base damage, single strand breaks (SSBs), double strand
breaks (DSBs), DNA cross links and deletions. SSBs are believed to be repaired by the cells if the cells do not have deficient repair pathways. Even
some 99 % of the DSBs are repaired, while clustered damage within a specified distance from each other will cause severe damage that the cells are
unable to repair and might lead to cell death. It is generally believed that
large parts of the DNA are unimportant (at least for the cells proliferation capacity) and is redundant. This indicates a more complicated mapping from
ionization event to DNA damages which may cause cell death. Multiploidy
cells have a duplicate of each gene, which has given rise to a hypothesis that
not only a single critical target has to be inactivated in order to kill the cell
but two critical targets, i.e. the duplicates of the genes. It is the connection of
the spatial distribution of the energy deposition and the model of the targets
and their damage induction and repair mechanisms that allows predictions
of the radiation outcome. One of the simplest example that spatial distribution and target models are correlated is the amorphous track structure
model of Katz (Butts and Katz, 1967). Two main approaches are being used
to correlate the track structure data with DNA damage.
The first approach consists of appropriate mathematical treatment of
track structure data, such as proximity functions (Kellerer and Chmelevsky,
1975) in the theory of dual radiation action (Kellerer and Rossi, 1978) or the
promising findings made by Brenner and Ward, (1992), where correlation
between the observed yields of double strand breaks and cluster sizes was
established.
Other relevant studies, correlating energy deposition events and different
kind of DNA damage are those made by Michalik (1992), Goodhead (1994)
and Ottolenghi et al. (1999).
The other and more recent approach explicitly takes into account biological targets by incorporating models of biological matter in track structure
codes, which gives a full 3-dimensional stochastic description of the tracks
in the media and the biological structures (Nikjoo et al., 1999; Friedland
et al., 2011). The primary conclusion of all the studies is that the majority
54
5.1. Particle fluence
of DNA damage are of a simple type and the complexity of DNA damages
increases with increasing LET. The level of detail of the targets vary from homogeneous distributed DNA in a form of a cylinder to highly sophisticated
models of the DNA.
A question has to be raised: Are those detailed simulations really needed
or is it enough to determine ionization cluster size distributions in nanometric volumes of liquid water, as a substitute of sub-cellular structures or
more generally, define more appropriate physical quantities (Grosswendt,
2004)? Grosswendt has proposed that the probability of cluster-size formation and its moments, calculated by MC technique, in nanometric volumes
behaves as a function of radiation quality.
Microdosimetric single event distributions have frequently been used to
successfully characterize radiation quality, especially within the radiation
protection field. The dose average mean lineal energy, y D , which is derived
from single event distributions, is on of the most common microdosimetric quantities to characterize beam quality. The microdosimetric quantities
were developed and derived from pulse height spectrum measurements in
ion chambers and are somewhat impractical, especially nowadays with the
access to MC track structure codes. y D can be calculated from the proximity
function or by MC track structure codes. The main drawback with y D , which
is the microdosimetric analogue to dose average LET, L D , is the dependence
on the scoring volume.
Nearest neighbor analysis in three-dimensions made with MC track structure codes is a good candidate for beam characterization. Analytical modeling and analysis of nearest neighbor could be an alternative to its MC counterpart. Higher moments of nearest neighbor analysis may give further information regarding the biological outcome.
An overall conclusion may be drawn; it is more relevant to study the number (or the density) of potential lethal events along the track or in volumes of
various sizes instead of the amount of energy along a track, like LET. Different such parameters have been investigated, e.g. SI (Specific Ionization) and
SPI (Specific Primary Ionization), nearest neighbor for radiation of different
ionization densities. Many of the parameters can be effectively calculated
with the KItrack MC track structure code, developed for electron transport
in water.
5.1
Particle fluence
Fluence is a macroscopic field quantity and one of the most fundamental
radiometric quantities used in dosimetry. When the differential fluence
r ), differential in energy and angle, is known, many other useful
ΦE ,~Ω (~
fundamental quantities related to the particle field can be derived. It is a
measure of the strength of the radiation field and can be defined in two
55
Chapter 5. Beam quality characterization
or more ways. ICRU defines the fluence in terms of the number of the
particles, dN , crossing a small sampling sphere with cross sectional area
da:
À
¿
△N
dN
,
(5.2)
=
Φ=
da
△a
where the second equality is a refined and mathematically correct way of
writing due to the fact that we are dealing with particles and they can not
impinge on an infinitesimal sphere. An alternative, and more general, definition was proposed by Chilton (1978; 1979) in which the total path length, dl ,
contained within any sampling volume, dV , are used to calculate the incident fluence:
dl
.
(5.3)
Φ=
dV
The latter approach is often used in Monte Carlo simulations of individual
particle tracks, allowing the fluence to be scored in small volumes of any
shape. By multiplying the numerator and denominator in Eq. (5.2) with the
mean chord length of a the sphere s = 4dV /dA, where dA is the area of the
surface of the sphere, it is easy to see that the two definitions are identical
for the special case of a spherical volume.
Calculations of the dose deposited from the differential fluence is an intricate problem but the general (but approximate) relation is however rather
simple:
ˆ
S col (E )
dE ,
(5.4)
D = ΦE
ρ
The equation above can only be applied if δ-particle equilibrium exists otherwise more complicated approximative formulas including restricted stopping power etc have to be used. The difficulty to choose an appropriate
cutoff energy for the restricted stopping power and the integration limit of
the integral and to really know what kind of equilibrium that exists, makes
dose calculations from fluence only approximative. Low-energy electrons
contribute substantially to the fluence, due to the high probability for elastic
scattering and their resulting long path length, and can cause large errors in
the dose calculation if not taken into account properly. The chosen cutoff
energy for the MC simulations can distort the fluence spectrum and the
dose if calculated from the relation above. It is important that the stopping
power and the fluence spectrum used are calculated using the same cross
sections, otherwise an error will be introduced in the calculated dose. Calculating fluence by MC class I or II codes can lead to step-size artifacts, which
is obviously due to the simplification of the step length and angular deflection. In detailed track-structure codes that make use of voxelized geometry
for the scoring procedure there is need to do boundary checks as the particle
is followed through the media. For high resolution in the scoring of the fluence many boundary checks for those volumes have to be made and is very
56
5.1. Particle fluence
time consuming. To avoid this, a larger volume can be chosen, but the spatial resolution in the scoring of the fluence is then lost.
Fluence scoring
An alternative scoring of the fluence has been derived, which is more
natural in the present track structure code without voxelized geometry. The
alternative scoring circumvents aforementioned problems. The method is
based on the definition by Chilton, i.e. track-length per unit volume and
the MC technique, i.e. the averaging over many particles to calculate an
average value. The path lengths between successive interaction events
with coordinates (x i , y i , z i ) and (x i +1 , y i +1 , z i +1 ) are scored ‘on the fly’ in
bins, which are determined slightly different for fluence as a function
of depth and for radial fluence around a particle. In the present work
fluence in depth is represented by electrons with energies 1 and 100 keV
and radial fluence by a proton with energy 0.5 MeV. For the radial fluence
the reciprocity theorem is used along the track, which means that the
z-direction is collapsed,
path length is performed in
qand scoring of the
¢ ¢2
¡¡
2
the bin closest to r = ((x i + x i +1 ) /2) + y i + y i +1 /2 and divided by
the area of that collapsed cylindrical bin. Track segment conditions are
assumed and the secondary electron production is sampled from CDW-EIS
cross section, see section 3.3. The use of track segment condition and the
reciprocity theorem implies that the track-length per bin area has to be
divided by the mean free path of the ion. In the case of fluence as a function
of depth, the track length is scored in the bin closest to (z i + z i +1 ) /2 and
divided by the length of that bin. This means that the lateral dimensions
are collapsed, and, by using the reciprocity theorem, this gives the fluence
in a broad beam.
Simulations of the fluence as a function of depth for electrons with initial energy 1 keV and 100 keV are presented in Figure 5.1 and Figure 5.2 respectively. The scoring is performed as described above with linear scoring
bins along the initial beam axis and infinite extension in the plane resulting in fluence per unit fluence. Due to an increase in the mean scattering
angle with depth, the path length increases and the fluence, as a function
of depth, shows a build up. At low initial energies the phenomena of large
scattering angles results in fluence behind the particles point of birth, i.e.
in the negative direction of the original beam direction. The total and the
secondary fluence is largely affected by the chosen cutoff. Note that the secondary fluence is higher than the primary fluence for cutoff 10 eV but with
cutoff 50 eV the situation is the opposite. This is due to the high probability for elastic scattering for low-energy electrons along with the decreasing
probability for ionization together with the definition that the electron with
lowest energy in an ionization event becomes the secondary electron. This is
a well-recognized definition in quantum mechanics and originates from the
57
Chapter 5. Beam quality characterization
indistinguishability of the two electrons (in experiments). MC calculations
follow this convention almost exclusively, due to the cross sections available.
An alternative definition could be that the primary electron is lost if the energy transfer is above a certain energy and the two resulting electrons are
both denoted as secondaries or randomize the denotation of the electrons
as primary and secondary.
Figure 5.3 shows the fluence, differential in energy, around a 0.5 MeV proton beam. The scoring is made in logarithmic energy and radial bins due to
the rapid radial fall off for the fluence and so that the interesting low energy
part of the spectra does not get lost. A well defined Binary-Encounter peak
can be seen in the differential fluence spectrum of the 0.5 MeV proton.
Fluence per unit fluence
40
Cut off 10 eV
30
20
10
0
Cut off 50 eV
-2e-08
0e+00
2e-08
Depth / m
4e-08
6e-08
Figure 5.1: Fluence per unit primary fluence in liquid water calculated with the
present track structure code. Total (solid line), primary (dashed line) and secondary (dashed-dotted line) fluence as a function of depth for an electron with initial
energy of 1 keV with cutoff 10 eV and 50 eV.
58
5.1. Particle fluence
Fluence per unit fluence
4
3
Cut off 10 eV
2
1
Cut off 50 eV
0
0e+00
1e-04
5e-05
Depth /m
Fluence
Figure 5.2: Fluence per unit primary fluence in liquid water calculated with the
present track structure code. Total (solid line), primary (dashed line) and secondary (dashed-dotted line) fluence as a function of depth in water for an electron with
initial energy of 100 keV with cutoff 10 eV and with cutoff 50 eV.
1e+19
1e+18
1e+17
1e+16
1e+15
1e+14
1e+13
1e+12
1e+11
1e+10
1e+9
1e+8
1e+7
1
10
−10
10
−8
2
10
10
10
3
−6
10
4
10
10
5
−4
10
Radius / m
Energy / eV
Figure 5.3: Electron fluence in liquid water per ion, differential in energy for 0.5 MeV
protons as a function of the radius.
59
Chapter 5. Beam quality characterization
5.2
Slowing down spectra
A major part of energy deposition in ion beams is mediated by secondary
electrons generated along the path of the beam. The difference seen in RBE
for different radiation qualities is dependent upon those secondary electron
spectra. Therefore, the secondary electron slowing down spectra (SDS) from
ions have been calculated with the MC track-structure code assuming tracksegment conditions. Similar MC calculations have been made by Lappa et
al. (1993), Tilly et al. (2002), and González-Muñoz et al. (2008) among others.
The secondary electron production is modelled by CDW-EIS cross section and the cutoff for the transport of secondaries is 10 eV. Same track
segment approach described in the section above is adopted also for slowing down spectra computations and logarithmic bins are utilized for the
same reasons as before. To calculate the SDS by the track structure code,
the path length traveled by primary and secondary electrons in the energy
interval between E and E + dE , where dE is the energy bin size, is scored
and divided by the mean initial energy of the electrons and by dE . Electron
slowing down spectra for 1, 10 and 100 keV electrons are depicted in Figure 5.4 together with protons with energies 0.5 and 1 MeV and lithium ions
with energies 0.5 MeV/u and 10.34 MeV/u. The electron slowing down spectrum increases with decreasing electron energy because of the increase in
the number of low-energy secondary electrons and the decrease in electron
stopping power at energies below 100 eV. Protons with energy of 0.5 MeV
and 10.34 MeV/u lithium ions have the same LET (~42 eV/nm) according
to SRIM tables (Ziegler et al., 2008), but their SDS show large differences
as seen in Figure 5.4. The SDS is a relatively good estimation of the beam
quality, though not that utilizable in practice, and it shows the expected differences between ion species with same LET. Protons with energy 0.5 MeV
have higher RBE than lithium ions with energy 10.43 MeV/u, which is due
to the smaller amount of low energy electrons of the lithium ions seen in
Figure 5.4. Ions with the same velocity, 0.5 MeV/u, protons and lithium have
almost the same SDS. This is to be expected due to the very similar cross sections of electron production, classically assumed to be identical for particles
with the same velocity. The small deviations seen arise from the fact that the
CDW-EIS cross sections goes beyond the Z 2 scaling seen in first Born approximation. Interestingly, even for a relative high initial energy electron, a
large amount of all ionizations are produced by secondaries with momentary kinetic energies below 1 keV.
60
5.3. Track penetration parameters
12
10
11
10
10
-2
-1
ΦE/D / m eV Gy
-1
10
1 keV e
9
_
10
10 keV e
_
_
100 keV e
8
10
3+
10.34 MeV/u Li
7
10
0.5 MeV/u Li
3+
+
+
6
1 MeV H
0.5 MeV H
10
5
10 1
10
2
10
3
10
Energy / eV
4
10
5
10
Figure 5.4: The solid lines are calculations of slowing down spectrum of 1, 10 and
100 keV electrons and 0.5 MeV, 1 Mev protons and 10.34 MeV/u Lithium ions with
present track structure code. The dashed curve labeled 0.5 MeV/u corresponds to
Lithium ions also calculated with the present code. Dashed-dotted lines are electron slowing down spectra from Tilly et al. (2002) and dashed-dotted-dotted line
are data for 1 MeV proton from González-Muñoz et al. (2008).
5.3
Track penetration parameters
The present track structure code has been used for calculation of various
track penetration parameters, which contain information on the spatial extent of the track. Characterizing the tracks with those parameters makes
the analysis and visualization of different tracks easy, and at the same time
minimizes the loss of information that occurs when the whole track with
its many coordinates, corresponding to elastic and inelastic events, is no
longer considered. In Paper II an extensive analysis of the probability density functions of the track penetration parameters was made together with a
presentation of the mean, median, standard deviation and skewness of the
parameters.
Monoenergetic electrons start their tracks at the origin, moving initially
along the z axis. Every interaction point (xi , yi , zi ), both elastic and inelastic,
61
Chapter 5. Beam quality characterization
is considered in the calculations of the track parameters. When the primary
electron energy falls below a pre-selected cutoff kinetic energy E abs it is
stopped at that position with coordinates (xabs , yabs , zabs ) . For each electron track, the following scalar quantities were defined in Paper II, beside
the already mentioned zabs :
2
2
)1/2
+ y abs
• ρ abs =(x abs
2
2
2
)1/2
+ z abs
+ y abs
• r abs =(x abs
• z max = maxi {zi ; i = 1, ..., n}
ª
©
• ρ max =maxi (xi2 + yi2 )1/2 ; i = 1, ..., n
©
ª
• r max =maxi (xi2 + yi2 + zi2 )1/2 ; i= 1, ..., n
Subscripts ’abs’ and ’max’ refer to penetration parameters obtained from
energy deposition points belonging either to only the primary electron or
to the full track it produces, i.e. including all generations of secondary electrons. The total path of a primary, i.e. all track segments between consecutive interaction points summed, is denoted l.
To go one step further, the covariance of the parameters can be calculated, which provides a measure of the strength of the correlation between
two parameters. The covariance between the random variables X and Y is
calculated by
C ov(X , Y ) = E(X Y ) − E(X )E(Y ),
(5.5)
where the mean values are averages of a large number of simulations. A
more convenient way of interpreting the correlation between parameters is
to use the correlation coefficient, ρ X Y , which is the covariance normalized
with the standard deviations, σ X and σ X , of the two studied parameters:
ρXY =
C ov(X , Y )
σ X σY
(5.6)
This coefficient takes values between -1 and 1. The correlation coefficient is
positive if x i and y i , which are the i t h measurements of X and Y, tend to be
simultaneously greater than, or less than their respective means. The correlation coefficient is negative if x i and y i tend to lie on opposite side of their
respective means. A correlation coefficient of zero indicates no correlation
between the variables.
The correlation of the parameters
The 21 correlation coefficients resulting from the defined track penetration
parameters, labeled r abs r max etc, are displayed in Table 5.1 for energies of
100 eV, 500 eV, 1 keV, 2 keV, 10 keV and 100 keV. Interestingly, z max and ρ max
are only positively correlated for energies below 500 eV, where they has medium strong correlation of 0.45. This phenomenon can be explained by the
62
5.3. Track penetration parameters
high amount of elastic scattering below 500 eV. If the cut of is set to 50 eV instead of the default cutoff (10 eV) the aforementioned correlation decreases
from 0.45 to 0.28. The probability for the electron to have a large penetration depth is high if it, at the same time, has a large lateral penetration. It
can also be expressed as if the electron undergoes only a few interaction
events before it is absorbed. This means that those events include a large
energy transfer, it will neither penetrate deep into the medium nor experience a large lateral penetration. Again at as high energy as 100 keV a small
negative correlation can be seen between z max and ρ max . There is almost
no correlation between z abs and the total path l, at least for energies below 10 keV. At energies above 10 keV a quite weak correlation around 0.3
is seen. The high amount of elastic scattering implies the that projection
of the point of absorption is not dependent upon the electrons history, i.e.
how many events it has experienced, which kind of event and their location.
The correlation between r max and z max is quite high and increases when
the energy increases. A high energy electron will penetrate deep into the
material and its path is relative straight, which makes r max and its projection down to the z-axis, z max , to almost coincide. Strong correlation is seen
between r abs and ρ abs and between r max and ρ max for electron energies
below 1 keV. A large lateral penetration will affect the distance from origin
to absorption point (or to the maximum distant point) as long as the penetration depth is not that large, which is fulfilled for quite low energies of
the initial electron energy. A correlation coefficient close to one is found for
r abs r max ,z abs z max and ρ abs ρ max at high energies. It can be understood easily by the fact the absorption point is also the maximal penetration point for
high energies and then the other quantities will follow. With a cutoff 50 eV
also lower energies experience a high correlation for those parameters. A
quite high correlation is seen between r abs and l and r max and l for low
and high energy respectively, while this correlation disappears for intermediate energies. Electrons with high energies have relatively straight tracks
with small off scattering from incident axis, which explains the high correlation for high energies. For low energies it can be argued that if an electron
undergoes a large amount of scattering the path, l, will be long and most
certain this will be valid for the r abs also. If the electron undergoes very few
but catastrophic events, i.e. ionizations were it loses about half of its energy
(the convention is that the primary can only lose up to half of its energy to be
called a primary and r abs is a quantity associated with primaries only) and
the electron stops suddenly without any elastic scattering, the two quantities l and r abs will be closely related. In an extreme situation a very low energy
electron can undergo only one inelastic event and deposit all its energy on
the spot, which makes the two quantities exactly the same. For intermediate
energies neither of the explained cases are applicable resulting in very low
correlation between those two quantities.
63
Chapter 5. Beam quality characterization
100 eV
500 eV
1 keV
10 keV
100 keV
r abs r max
0.54
0.55
0.88
0.98
0.97
r abs z abs
0.15
0.58
0.63
0.67
0.69
r abs z max
0.37
0.48
0.64
0.70
0.69
r abs ρ abs
0.88
0.67
0.59
0.49
0.39
r abs ρ max
0.50
0.34
0.52
0.48
0.36
r abs l
0.67
0.28
0.12
0.66
0.61
r max z abs
0.71
0.34
0.58
0.66
0.69
r max z max
0.63
0.67
0.72
0.73
0.72
r max ρ abs
0.48
0.35
0.51
0.47
0.37
r max ρ max
0.93
0.77
0.61
0.49
0.36
r max l
0.47
0.19
0.14
0.73
0.69
z abs z max
0.48
0.69
0.88
0.94
0.95
z abs ρ abs
0.00
-0.02
-0.08
-0.18
-0.28
z abs ρ max
0.00
-0.01
-0.07
-0.19
-0.29
z abs l
0.03
0.01
0.03
0.28
0.25
z max ρ abs
0.22
0.03
-0.04
-0.17
-0.30
z max ρ max
0.46
0.20
-0.02
0.17
-0.31
z max l
0.31
0.11
0.09
0.44
0.40
ρ abs ρ max
0.53
0.48
0.85
0.98
0.98
ρ abs l
0.61
0.29
0.11
0.49
0.41
ρ max l
0.45
0.19
0.12
0.54
0.46
Table 5.1: Correlation coefficient of the track penetration parameters for 100 eV,
500 eV, 1 keV, 10 keV and 100 keV for a cutoff 10 eV.
5.4
Specific ionization
The RBE depends on many factors, some related to the physical
characteristics of the radiation and others on the biological system in
question. For low-LET photon and electron beams absorbed dose works
well as a descriptor of the beam quality, while other radiation modalities
such as protons and ions physical properties are not uniquely defined by
the dose imparted in the biological system. For protons and ions the RBE is
also a function of the LET. The height of the RBE maximum is continuously
decreasing with increasing atomic number of the ion while the location
of the RBE maximum is shifted to greater LET values. The high RBE of
the low energy protons becomes lower in clinical beams due to the lower
LET resulting from range and energy straggling. The heavier ions are less
affected by range and energy straggling, which mantain the high LET and
64
5.4. Specific ionization
therefore is their RBE less affected. This indicate that the amount of local
energy deposited can not be the most critical factor for RBE, but rather the
proximity of the energy deposition events. Instead, counting the number
of ionization events per path length (and not the amount deposited in
each event) could be one way to find a one to one relation with RBE. This
hypothesis was tested by calculating the Specific Primary Ionization (SPI)
per path length, i.e. the number of primary ionizations, and the Specific
Ionizations (SI), i.e. the total number of ionizations, both primary and
secondary, per path length. The aim has been to find a measure or quantity
that is not proportional to LET and can reflect the different RBE values seen
for different ions with same LET.
When an ion beam decelerates, lighter fragments are produced and in
the Bragg Peak both low LET and high LET components are present (Kempe
et al., 2007). A therapeutic beam is, due to this phenomenon, hard to analyze
in terms of SPI and SI and it is also cumbersome to simulate a therapeutic
beam. So, throughout this chapter, only monoenergetic ion beams are considered and track segment calculations are performed. The RBE peak is located around the Bragg peak and its height is dependent on the biological
system, the chosen endpoint, the irradiation set up and LET (or energy) and
particle type. The exact position of the RBE peak in relation to the BP is dependent on the charge of the ion and differs between clinical beams and
monoenergetic beams, due to clinical set ups altering the beam composition.
MC calculated SPI and SI
SPI and SI are found to vary linearly with LET for the investigated particles
(1 H+ , 4 He2+ , 6 Li3+ , 12 C6+ ) and energy ranges (0.3-30 MeV/u). The SI value
is around 3 for for all investigated particle types with an LET of ~42 eV/nm
(not shown in the Tables), but the SPI values seen in Table 5.3 decrease with
increasing charge of the ion. The reason for the choice of this particular LET
value is based on the range of the electron interaction databases implemented in the present track structure code and balanced against the validity of
CDW-EIS ionization cross sections for ions. 0.5 MeV/u is the lowest acceptable energy, at lower energies the approximations made in the CDW-EIS
theory becomes invalid. This implies that a proton beam has less distance
between primary ionizations, but the resulting secondary electrons are of
lower energy than from an ion with higher charge. The RBE maxima have
in some experiments been found to appear around 25 (or at slightly higher
LET of 30 eV/nm and corresponding energy 0.8 MeV/u ) and 105 eV/nm for
protons and 4 He2+ ions respectively (Cera et al., 1997; Belli et al., 1989) and
corresponding energies are 1 MeV/u for both protons and 4 He2+ . 6 Li3+ is
expected to have its RBE maxima around 110-130 eV/nm with corresponding energies 3-2.5 MeV/u and 12 C6+ at as high LET as 170 eV/nm and energy
10 MeV/u (Weyrather et al., 1999), but there are experiments showing a peak
65
Chapter 5. Beam quality characterization
at LET of 130 eV/nm. The RBE peak position will vary slightly with cell type
mainly due to the different cell sizes, but also difference in cell sensitivity
and DNA density and, of course, of chosen end point.
Difficulties to maintain a pure ion beam with a single LET in experiments
and the different ways of defining the average LET in those beams makes it
hard to find a corresponding LET for the ion beam simulated under track
segment conditions. The chosen RBE peak position for the different particle
species are based on experiments mostly performed on V79 cells for cell inactivation, which is a pre-stage for endpoints such as different kinds of cell
deaths and senescence.
There are no consensus in the literature where the RBE peaks are located
for different ion species so no real conclusions about coupling between the
biological effect and the SI and SPI values at those LET can be drawn. SPI
and SI values for the RBE maximum can be found in Table 5.2. The LET at
the experimental measured RBE peak for helium is surprisingly high so a
lower LET value of 80 eV/nm has been added.
The only conclusion that can be drawn from Table 5.2 is that the number
of ionizations at the RBE peak increases as the ion atomic number increases
and that the SI/SPI is around 10 for RBE peak ions. The RBE maximum does
not coincide with LET maximum for the particles mentioned, which can be
explained as a consequence of too much clustering of the energy events in a
track, and the very few tracks needed for very high-LET to achieve the same
dose as for low-LET. As stated in section 5.5, the least biological effect should
be seen for a homogeneous distribution of ionization events. The optimal
biological effect is therefore seen at a specific level of inhomogeneity in ionization events. The LET maximum ranging from 80 eV/nm for protons to
820 eV/nm for carbon ions can be found in the energy range 0.1-0.4 MeV/u.
Particles with same velocity and energy 0.5 MeV/u show SPI and SI values
that almost scale with Z 2 , which is to be expected due to the nearly Z 2 scaling of the cross section, see Table 5.3 and 5.4. Those particles also show very
similar slowing down spectra, see section 5.2.
The DNA helix is approximately 2 nm in diameter and wired around a
nucleosome with diameter 11-30 nm. To induce on average one DSB in a
DNA helix a SPI value of 0.5 per nm should be aimed at. In case of the SI
such simple relation can not be expected due to the fact that the ionizations
caused by the secondaries can take place far away from the track, and their
action has to be measured per unit volume, at least. Protons have an SPI of
0.5 at very low energy, 0.38 MeV, while carbon ions reach this at a relative low
energy of 20 MeV/u, see Table 5.3. To obtain an SI of 0.5 for ions a quite high
energy is required, which implies that a cross section database for secondary
electrons ranging from 10 eV up to at least 106 eV. The present track code
has only databases including energies from 10 eV to 105 eV. Electrons with
energy 450 eV have an SPI value of 0.5 and an SI value of 0.5 at 3.2 keV.
66
5.4. Specific ionization
SI is not sufficient to characterize the beam quality as it seems from the
investigation. The ratio of SPI/LET at the RBE peak increases with decreasing Z, and for the same LET the lighter particles have higher SPI values than
the heavier ions. The tendency of SPI is therefore correct in the sense that it
manages to explain the higher RBE for lighter particles and from that point
of view seems to be slightly better than LET. The investigation includes too
few values though and the experimental LET-energy-RBE peak relations are
too inadequate for any conclusions to be drawn and the results indicate that
higher order measures, like the spatial distribution of the ionization events,
are needed.
[MeV/u]
LET [eV/nm]
SPI [1/nm]
SI [1/nm]
0.8 1 H+
30
0.26
2.3
25
0.20
1.9
105
0.84
8.0
80
0.59
5.8
1
1 H
+
4
1 He
2+
4
1.5 He
6
2+
3+
3 Li
110
0.71
8.0
3+
2.5 Li
130
0.84
9.0
14 12 C6+
130
0.70
9.5
12 6+
170
0.95
12.6
6
10
C
Table 5.2: SPI and SI for LET values corresponding to the RBE peak for various ion
species.
e_
1
Same velocity (0.5
MeV/u)
–
Same LET ~42
eV/nm
Corresponding
energy [MeV/u]
Energy that gives
SPI ~0.5 [MeV/u]
H+
4
He2+
6
Li3+
12 6+
0.36
1.55
3.40
11.6
–
0.36
0.25
0.22
–
–
0.50
3.73
10.34
–
~0.45·10−3 ~0.38
~1.7
~5.0
~19.5
C
Table 5.3: SPI values for ions with; energy 0.5 MeV/u, same LET 42 eV/nm and the
energi of ions giving a SPI of 0.5 per nm. The SPI value for a carbon ion with LET of
42 eV/nm is not possible to calculate with the present track code.
67
Chapter 5. Beam quality characterization
1
H+
4
He2+
6
Li3+
12 6+
C
LET eV/nm
~42
~164
~311
~797
SI
~3
~13
~29
~104
Table 5.4: SI values for ions with energy 0.5 MeV/u and their corresponding LET.
5.5
Nearest neighbor of energy deposits
Provided that the spatial distribution of energy deposits (disregarding
the amount of energy per deposit) is the most important characteristic, a
concept more in line with SPI than LET should be aimed at. Such measure
could be the average (and higher moments) nearest neighbor distance in
3D space. In order to account for both inter and intra track distances the
particle fluence should also be included. Following this line of thoughts,
the least biological efficiency would be expected for completely uniform
density of energy deposits. An expression for the average nearest neighbor
distance in n-dimensions is thus of interest.
Assume that the number of interactions per unit volume follows a Poisson
distribution with probability density function
p(x; ρ) =
e −ρ ρ x
,
x!
(5.7)
where ρ is the average number of interactions per unit volume. Let ρ(r ) be
the probability density function for the nearest neighbor distance. Thus the
probability to find the nearest neighbor of an arbitrary chosen interaction
point in the spherical shell (r, r +dr ) will be ρ(r )dr . ρ(r )dr can be calculated
as the product of ρ 1 (r ), the probability of having one or more interactions in
the spherical shell and ρ 0 (r ), the probability of no interactions in the sphere
4πr 3 /3. Since the expected (average) number of interactions in the sphere
is 4πρr 3 /3, ρ 0 (r ) is easily calculated from Eq. (5.7) by putting x = 0 resulting
in ρ 0 (r ) = exp(− 34 πρr 3 ). The probability of no interactions in the spherical
shell is thus similarly given by exp(−4πρr 2 dr ) and the probability of one or
more interactions in the shell ρ 1 (r ) = 1 − exp(−4πρr 2 dr ). Taylor expansion
gives ρ 1 (r ) = 4πρr 2 dr leading to
4
3
ρ(r ) = 4πρr 2 e − 3 πρr .
(5.8)
The average nearest neighbor distance is then by definition
r¯ = 4πρ
ˆ∞
4
3
r 3 e − 3 πρr dr.
0
68
(5.9)
5.5. Nearest neighbor of energy deposits
Substituting t = 43 πρr 3 , dt = 4πρr 2 dr gives
µ
4
r¯ = πρ
3
¶− 1 ˆ∞
3
t 1/3e −t dt .
(5.10)
0
Noting that the gamma-function is defined as
Γ(z) =
ˆ∞
t z−1 e −t dt ,
(5.11)
4
0.554
Γ( ) ≈ 1/3 .
3
ρ
(5.12)
0
Eq. (5.10) can be identified as
µ
4
r¯ = πρ
3
¶− 1
3
The second moment r (2) is given by
r
(2)
= 4πρ
ˆ∞
4
3
r 4 e − 3 πρr dr,
(5.13)
0
and with the same substitution as above
r
(2)
µ
4
= πρ
3
With the solution
r
(2)
¶− 2 ˆ∞
3
t 2/3 e −t dt .
(5.14)
0
µ
4
= πρ
3
¶− 2
3
5
Γ( ).
3
The coefficient of variance, CV , is thus given by
q ¡ ¢
¡ ¢
Γ 53 − Γ2 34
CV =
≈ 0.363
¡ ¢
Γ 34
(5.15)
(5.16)
Eq. (5.12) can be generalized to n-dimensions 2 :
µ
n
r=
ρk n
¶
¶1/n µ
1
,
Γ 1+
n
(5.17)
where k n is a constant given by a closed-form expression and takes the values e.g. 2, 2π and 4π for the first 3 dimensions. In 1D the average number of
interactions per unit length, ρ, is 1/λ, which give the nearest neighbor distance of λ/2. Given the expression for the average nearest neighbor distance,
2 As suggested by Prof. Tommy Elfving
69
Chapter 5. Beam quality characterization
one can calculate which fluence, Φ, is needed to achieve uniform density of
energy deposits from charged particles with inelastic mean free path λ. Let
the fluence, Φ, be equal to the mean number of interactions, ρ, in 2D. Solve
Eq. (5.17) in 2D for ρ, i.e. the fluence, with the knowledge that the nearest
neighbor distance, r , in 1D is λ/2:
¶
µ
1
2Γ (3/2) 2
/π = 2 .
(5.18)
Φ=
λ
λ
Anisotropy index
The level of uniformity of the interactions can be estimated by the ratio
between the shorter of the average distance between interaction in 1D and
2D and the nearest neighbor distance given by Eq. (5.12). This ratio is from
here on referred to as the anisotropy index. Figures 5.5-5.7 shows the level of
uniformity, expressed as the anisotropy index, as a function of dose and energy of the primary particle for electrons, protons and carbon respectively.
In the plots showed here are the average distance between interactions in 1D
always shorter than the one in 2D. Note that by using either the known average distance in 1D or 2D only primary interactions is taken into account.
It is well known that the secondary electrons and their energy deposits are
responsible for most of the damage to the cells due to high amount of clustering of the ionization events produced. The postulate is though general
and valid, i.e. that the least biological effect is expected for totally uniform
ionization density. In order to account for all events, both primaries and secondaries, an average nearest neighbor distance of interactions in 3D including all events need to be known so the anisotropy index can be formed as
the ratio of this distance and the average distance in 3D assuming Poisson
distributed events given by 5.12. Surprisingly, an anisotropy index of nearly
one is found for electrons with energy 0.1 MeV and a typical dose per fraction of 2 Gy. Knowing that the SI/SPI is around 5 (see section 5.4 for explanation) for 0.1 MeV electrons one can expect that the longitudinal distance
between interaction decreases with less than 5. This has the effect that an
anisotropy index of one is found at slightly higher doses. The dose at which
an electron of 1 MeV have an anisotropy index of one, including secondaries, could then be estimated to be at approximately 2 Gy. Figure 5.5 gives the
impression that a high fluence can be an alternative to high LET. For high energy electrons the collision stopping power is rather constant (2 MeVcm2 /g),
so a much lower LET can not be obtained with electrons by increasing the
energy. Maybe synchrotron light could be used. The lowest level of uniformity/anisotropy index is obtained when a dose per fraction of 0.1-8 Gy is given
by ions with energies corresponding to the Bragg Peak of the ion in question. This is in line with the high biological effect seen for those ions. By
evaluating the ratio of SI and SPI for ions in the RBE peak given in Table
5.2 in chapter 5.4 one can find that this ratio ranging from 8-13. This implies an even lower anisotropy index is expected for those ions, because the
70
5.5. Nearest neighbor of energy deposits
nearest neighbor distance of all interaction events will be smaller than the
average distance in 1D for primaries, 1/λ. The biological effectiveness will
increase with decreasing uniformity and probably come to a saturation level
and even decrease when the clustering is too high. This can be coupled to
the theory of so-called overkill effect, which means that some cells receive
unnecessary large amount of energy depositions while other cells receive
no energy deposits all all. This fact and the neglecting of secondaries could
explain why the lowest level of uniformity/smallest anisotropy index does
not coincide with the RBE peak seen in experiments. Nearest neighbor distances with different criterion and higher moments of the nearest neighbor
could be calculated with the present track code under assumption of track
segments. This approach to quantify the biological efficiency in terms of anisotropy index as a function of energy (LET) and dose, could be of much help
in for example treatment planning systems and could be one way to characterize the beam quality beyond LET etc.
1000.000
Anisotrophy index
100.000
10.000
1.000
0.100
1.E+07
1.E+03
1.E+01
0.001
1.E-01
1.E+01 1.E+02 1.E+03 1.E+04
1.E+05
Dos
e
0.010
/ Gy
1.E+05
Energy / eV
Figure 5.5: Level of uniformity of the primary interactions in a electron beam as a
function of dose and energy of the particle. Data for LET-energy relation and mean
free path of the electron beam are calculated by present track code.
71
Chapter 5. Beam quality characterization
10000
1000
10
1
1.E+07
0.1
1.E+05
0.01
0.001
1.E+02
1.E+03
1.E+01
1.E+04
1.E+06
1.E+08
Dose /
Gy
Anisotropy index
100
1.E-01
Energy / eV/n
Figure 5.6: Level of uniformity of the primary interactions in a proton beam as a
function of dose and energy of the particle. Data for LET-energy relation and mean
free path of the proton beam are taken from Watt (1996).
100
Anisotropy index
10
1
0.1
0.01
1.E+07
0.001
0.0001
1.E+02
1.E+03
1.E+01
1.E+04
1.E+06
1.E+08
Dose /
G
y
1.E+05
1.E-01
Energy / eV/n
Figure 5.7: Level of uniformity of the primary interactions in a carbon ion beam as a
function of dose and energy of the particle. Data for LET-energy relation and mean
free path of the carbon ion beam are taken from Watt (1996).
72
5.6. Radial dose profile (amorphous track modeling)
5.6
Radial dose profile (amorphous track modeling)
The average or amorphous track model (Butts and Katz, 1967) uses the
response of a system to γ-rays and the radial distribution of dose about an
ion’s path to describe survival and other cellular endpoints from protons,
heavy ions, and neutron irradiation. This model neglects the stochastic
nature of individual energy depositions, but has successfully been used for
over 30 years to fit many radiobiology data sets. More recent approaches of
using the radial dose distribution from ions and the response of a system to
γ-rays is formulated in the so-called LEM model (Scholz et al., 1997), see
section 4.4. The MC method to simulate radiation transport and calculate
quantities like absorbed dose is the most accurate way, but in practical
situations, in treatment planning for instance, an advantageous solution,
with high speed and relative high accuracy can be found by combining
analytical calculations in simple geometries with information from MC
simulations in complex geometries. Radial dose distributions are calculated
on nanometric scales where the dose concept is violated. Dose is defined at
a point, but when calculating the this quantity sufficient amount of mass
has to be averaged over. The calculation of the absorbed dose can, at least
for MC calculations, be interpreted as averaging the energy deposition over
a infinite time.
In Figure 5.8, a comparison of the radial dose distribution calculated
with an analytical method and a MC method, both described in Paper
I is shown, together with calculations performed with the developed
track structure code. The CDW-EIS cross section for secondary electron
production from ions, described briefly in appendix in Paper I and in
section 3.3, was used in all three calculations and track segment condition
is assumed, i.e. no energy degradation or deflection of the ion in the
segment where the calculations are performed. The spatial distribution
of the secondary electrons was in the analytical method treated using the
generalized pencil-beam solution of the Boltzmann equation (Brahme,
1975). This solution is based on a generalization of Fermi-Eyges solution
(Brahme et al., 1981; Hollmark et al., 2004). By assuming a slow stopping
power variation, the point monodirectional pencil beam dose distribution
separates in a radial and a depth dependent term. The electron transport,
except at the lowest energies, is dominated by many small-angle scattering
events, so the angular distribution at any depth is well approximated by a
Gaussian distribution and a MC broad beam distribution, D pp (z, ∞), here
calculated by PENELOPE, was used as the depth dependence. The resulting
expression becomes:
´
³
exp r 2 /r 2 (z)
,
(5.19)
D pm (z, r ) = D pp (z, ∞)
πr 2 (z)
where r 2 (z) is the mean square radial spread or the variance.
73
Chapter 5. Beam quality characterization
10
6
Dose / Gy
E/u=0.5 MeV/u
10
10
10
4
2
0
-1
10
10
0
1
10
Radius / nm
10
2
Figure 5.8: Radial dose profile calculated with analytical method, solid black line,
MC method based on PENELOPE kernels, black stairs, and present developed track
structure code for 0.5 MeV protons with cutoff 100 eV, red stairs, and with cutoff
10 eV, blue dashed stairs.
For a more detailed description of the computational procedure see
Paper I. In the MC method used in Paper I, a scoring and sampling routine
was written in C++ and the electron transport was performed with the
PENELOPE code (Salvat et al., 2008). The parameters in the PENELOPE
code are chosen so that detailed simulations are achieved both for broad
beam depth dose and for the ’electron point sources’ centred at the ion
beam axis. The two MC methods used here perform scoring in cylindrical
bins and take advantage of the reciprocity theorem and collapse therefore
the bins in the measuring plane. To facilitate comparison with the radial
dose, calculated with the PENELOPE based method, and the analytical
method the cutoff was set to 100 eV, but a simulation with cutoff 10 eV is
added in order to point out the relative large impact the electrons below
100 eV have on the radial dose. By setting the cutoff to 100 eV the large
amount of low energy primary electrons, below 100 eV, are forced to deposit
their energy at the beam axis, i.e. zero radius. Secondary electrons with
74
5.6. Radial dose profile (amorphous track modeling)
energy below 100 eV contribute to a large extent to the dose deposited and
make the radial dose profile somewhat wider when they are allowed to be
transported. The radial dose calculated with the MC method described in
Paper I, and with the present track code show close agreement when the
cutoff is set to 100 eV, though, slightly lower dose at small radii and at large
radii is seen for the track code calculations. High energy primary electrons
and low energy secondaries contribute to the dose both at the outer part
and at the innerpart of the radial dose profile. The differences seen in the
calculated depth dose, see Paper II, between PENELOPE and the track code,
i.e. the practical ranges at high energies are shorter respectively longer
at low energies for track code calculations compared to the PENELOPE
calculated ranges, and is therefore of no help in foreseeing the alteration
in radial dose resulting from the two MC methods. For a more elaborate
discussion about cross sections, stopping power and depth dose differences
see Paper II and section 3.
75
Chapter 6. Conclusions
6. Conclusions
• A MC track structure code, written in C++, has been developed. It
has the advantage of being modular and very flexible in its scoring
routines. It can easily be upgraded with new cross sections and/or
sampling routines etc.
• The default linear MC interpolation routine has been shown to have
favorable time consumption in relation to the relative error.
• MC codes commonly make use of voxelized geometries, which can
make them quite time consuming but in the same time open up for
scoring in small compartments used in, for example, fluence scoring
or calculations of y D . The present track code does not make use of
voxelized geometries and attempts have been made to circumvent the
procedure for classical scoring. For example, a novel way of scoring
fluence, based on the definition by Chilton (1978; 1979), has been introduced in chapter 5.1. This way of scoring the fluence is promising,
but has its limitations in applicability to systems consisting of a single
material.
• Radial dose profile calculations by an analytical model presented in
Paper I agree quite well with the radial dose profiles calculated with
PENELOPE and the present track structure code.
• Specific ionization, SI, have been calculated with the MC track code
and seems not to be superior to LET in characterizing radiation quality. Specific primary ionizations, SPI, shows a tendency to be slightly
better than LET as a measure of radiation quality, but maybe not good
enough, which indicates that higher order calculations of the energy
depositions are needed.
• In depth investigations of the contribution to the fluence from low
energy electrons, which undergo mostly elastic scattering leading to
very long path lengths, have to be done in order to fully understand
the applicability of the different relations between dose, fluence and
LET.
• Calculated slowing down spectra for protons agree well with spectra
calculated by Gonzalez-Munoz et al. (2008) and the small differences
76
seen are mostly due to different inelastic ion cross section used, but
originates also from differences in the electron transport.
• Calculated slowing down spectra for electrons agree well with spectra
calculated by Tilly et al. (2002) and the small differences seen originate from differences in the electron transport like cross sections and
cutoff energies.
• Small changes in electron inelastic cross sections introduced by the
use of different exchange corrections have minor impact on the electron tracks. The use of either coherent or incoherent summed wave
functions in the calculation of elastic cross section also have impact
on the electron track. How this alteration in spatial distribution of
electron tracks will affect the damages to a biological system has not
been dealt with here.
• An expression of the nearest neighbor distance for ionizations in 3D
for primary particles has been derived. The anisotropy index, which
is defined as the ratio between the shorter of the average distance
between interaction in 1D and 2D and the nearest neighbor distance
given by Eq. (5.12), is a good estimate of the level of anisotropy of the
energy depositions and can help to estimate the biological efficiency.
Therapeutic electrons were shown to have a anisotropy index close to
one, which means low biological efficiency. Protons and carbon ions
with energies corresponding to their respective RBE peaks showed a
very small anisotropy index indicating high biological efficiency. The
level of uniformity of the interactions can be estimated by the ratio
between the shorter of the average distance between interaction in
1D and 2D and the nearest neighbor distance given by Eq. (5.12).
• Due to the calculation of macroscopic α and β values for the
LQ-model in the LEM model, other models such as the RCR model
could not be used as a substitute for the LQ-model. The success of the
LEM model should perhaps be questioned and the LET parametrized
RCR model is probably a better choice in the future.
• If the dose to a tumor with low α/β-ratio is delivered heterogeneously
at cellular level and in many fractions, the TCP can actually increase
above the TCP for a homogeneously delivered dose.
• The derived expression for TCP for heterogeneous doses delivered at
cellular level (intra fractional variation) in more than 10 fractions has
been shown to be accurate for CV s in dose up to 25%.
• In case of inter fractional variation in doses, the derived expression
showed some deviation from the numerical simulations, but the over77
Chapter 6. Conclusions
all trend is the same, i.e. in case of 30 fractions an increase in TCP is
seen when the CV for inter fractional variation of doses is increasing
• To obtain a TCP of 50 % or higher, a large heterogeneity in dose distribution has to be compensated by a rather large number of fractions
especially for the low number of clonogens.
• For tissues sensitive to changes in fractionation, i.e. late reacting tissues characterized by low α/β values and pronounced shoulder of the
cell survival curve, the increased cell kill due to doses higher than the
average dose more than compensate for the decreased cell kill due to
lower doses than average. As the number of fractions increases this
effect will be amplified.
• The effect on TCP from different dose variations add up.
• The effect on TCP from combinations of variations in the dose to the
cell and inter tumor sensitivity variations add up, but the result is less
than the sum of the respective effects.
• The effect on TCP from combinations of variations in the dose to the
cell and intra tumor sensitivity variations add up.
• Heterogeneity in the dose per fraction and in the dose to the sites
within a cell always increase the TCP, while variations of the dose to
the cells can act both positively and negatively on TCP.
• There is a large difference in the impact on the TCP between sampling
of the variation in sensitivity by keeping the α/β constant or allow it
to vary, especially for the inter patient/inter tumor case.
78
7. Outlook
• New and better cross sections, such as the improved inelastic model
by Emfietzoglou and co-workers, will be implemented in the track
code and studies similar to those in Paper II, i.e the impact on the
track parameters by different cross sections will be performed.
• It would be of great interest to improve the derived analytical expression for nearest neighbor distance so that it more accurately describes
the distance between energy depositions by including the secondary
electrons contribution.
• It is well known that the spatial pattern of the energy depositions is
the most important characteristics of a beam in terms of biological
efficiency. The present code could be used in the future to calculate
nearest neighbor distance, and higher moments of the nearest neighbor distances that could be of help in determining the biological efficiency.
• The calculations of the response of the radiosensitivity sites within the
cells is interesting in its own, but it would be of much interest to test
a hypothesis that a certain amount of critical sites has to be damaged
in order to eradicate the cell together with dose and sensitivity variations.
• In Papers III and IV the beam characteristics were kept constant for
the whole treatment. In other words, the same LET was used for all
fractions. The impact of a mixture of LETs on the TCP is of great interest, especially if the response of sites is accurately taken into account.
• It would be interesting to use the RCR model instead of the LQ-model
in the numerical simulation of tumor control probability to be able to
introduce variances in repair capacity and damage induction separately.
• In Papers III and IV the analysis of the TCP were performed on an
increasing number of fractions where the doses per fraction were allowed to vary in order to maintain a TCP of 50%. It could be beneficial
to analyze the impact of varying cellular doses, keeping the fractional
dose constant by using different approaches in the analysis. It could
79
Chapter 7. Outlook
be worth while to optimize TCP and NTCP in terms of dose per fraction.
• It could be of interest to combine an inter patient variation in sensitivity with and intertumor variation in sensitivity were the sensitivity
of a certain tumor is set as the mean sensitivity of all the cells in that
tumor.
• A better parametrization of the radial dose, based on the radial dose
calculations see chapter 5.6, may improve the LEM ability to predict
cell survival.
80
Summary of papers
Paper I
In treatment planning there is a need for fast dose calculations. MC simulations are the most accurate way to calculate dose, but are time consuming.
Therefore, in this paper an analytical pencil beam method was developed
for radial dose calculations for ions. The fact that it is possible to separate
the dose distribution of a point mono-directional, mono-energetic pencil
beam D pm (r, z) in a purely depth dependent broad beam D bb (∞, z) and a
depth dependent radial dependent term was used in this paper. The spatial
distribution of the secondary electrons was treated using the generalized
pencil beam solution of the Boltzmann equation. This solution is based on
a generalization of the Fermi-Eyges solution. The electron production from
the ion beams was sampled from CDW-EIS cross sections. Track segment
conditions were assumed in the ion beam. To obtain the radial dose distribution, integration was carried out over all possible angles and over the
energy range 100-6000 eV and over the track segment of the ion. For comparison, radial dose calculations using pencil beam kernels, computed with
PENELOPE, were also performed. Results for H + , He 2+ and Li 3+ ions with
energies 0.5 Mev/u and 1 MeV/u calculated with both methods were presented and agreed quite well. The analytical calculations overestimate the radial dose at intermediate radii, 1 nm-10 nm, compared to the MC method,
more for the ions with energies 0.5 MeV/u. As expected, ions with the same
velocity showed very similar radial dose. Small differences are seen though
and are due to the fact that CDW-EIS cross sections used do not scale with
Z 2 (only first Born approximation does). The presented results agree well at
intermediate radii, 1-10 nm, with other calculations found in the literature
whether they are analytical or MC based. At smaller radii the presented results underestimated the dose and at larger radii an overestimation was seen.
Part of the deviation at small radii might be due to the fact that only electrons with initial energy above 100 eV were simulated.
Paper II
An object-oriented track structure code written in C++ for transport of 10100 keV electrons in liquid water was presented. Modularization of the code
81
provides benefit of easy implementation of cross sections other than those
set as default and enables the use of different programming languages for
the modules. The model by Dingfelder et al. (1998) was set as default ionization and excitation processes, and Olivero et al. (1972) as optional for the
excitation process. The adopted DCSs for elastic scattering in water liquid
were generated with the ELSEPA code. As default, the incoherently summed
atomic DCSs were chosen and the coherently summed atomic DCSs were
set as optional. The code is equipped with two random number generators namely drand48 and the SIMD-oriented Mersenne Twister. For sake of
completeness, cross sections for dissociative attachment are included in the
code. Track-segment calculations for ions can be performed with the code
when double differential inelastic cross sections for ion impact are available. Several track penetration parameters describing the geometry of the
track were defined and the probability density function of those were calculated and analyzed. Detour factors, intended to quantify the shortening
of the penetration depth, were calculated with the code. The sensitivity of
some of the penetration parameters to the choice of cross sections were
explored. It was found, as expected, that at low energies the impact of the
choice of cross section on the geometry of the track was not negligible.
Paper III:
Derivations of analytical expressions of tumor control probability (TCP) for
macroscopic inter-cell dose variations and for random inter-fractional variations in average tumor dose were presented, based on binomial statistics
for the TCP and the well known LQ model for the cell survival. Numerical
calculations were performed to validate the analytical expressions. An analysis of the influence of the deterministic and stochastic heterogeneity in
dose delivery on the TCP was performed. The precision requirements in
dose delivery were discussed briefly with support of the results presented.
The main finding in this paper was that it is primarily the shape of the cell
survival curve that governs how the response is affected by macroscopic
dose variations. The analytical expressions for TCP accounting for heterogeneity in dose can quite well describe the TCP for varying dose from cell
to cell and random dose in each fraction. An increased TCP was seen when
large number of fractions are used and the variations in dose to the cells are
rather high, for tissues with low α/β.
Paper IV:
Numerical simulations of the macroscopic inter-cell dose variations, intercell sensitivity variations and combinations of those were presented. These
82
simulations are based on Bernoulli trials, which allow to include cell repopulation during the course of treatment in an exact manner, whitout performing approximations. Analytical expressions of the TCP for heterogeneous
dose distribution derived in Paper III were further explored, now in terms of
heterogeneous cell sensitivity within a tumor. Inter patient/inter tumor heterogeneity in cell sensitivity and dose variations were investigated by performing numerical simulations. The LQ-model was used to model the cell
survival both in the numerical simultions and in the analytical expressions.
The analytical expression for TCP accounting for heterogeneity in sensitivity was validated against simulations. An analysis, based on the analytical
expression, of the dependence on SF 2 and α/β on the TCP has been performed as a complement to the simulations of the dose response for particular cases of α/β and SF 2. The conclusions are that heterogeneity in sensitivity always acts negatively on TCP, while heterogeneity in dose to the cells
can act positively when the cells are resistant.
83
Acknowledgments
I would like to express my sincere gratitude to
Bengt Lind, my supervisor and co-author. Thank you for your scientific
guidance. I appreciate your prestigeless attitude, and your open mind to
all kinds of ideas, both sane and insane. We have learned a lot during the
years.
Iuliana Toma-Dasu, my co-author and co-supervisor, for your guidance
and encouragement, especially in the last minute when I thought I would
fail.
Anders Brahme, my co-author, thank you for supervising me the first years.
You are a true visionary!
Jose-Maria Fernández-Varea, my co-author, for your encouragement and
excellent advices, both scientific and non scientific. You saved me! You are
responsible for almost all my knowledge in particle physics.
Bo Nilsson and Margaretha Edgren, for giving me basic education in
radiation physics and for being such nice persons. And thank you Bo for
proofreading in the last minute.
Irena Gudowska, for sharing your expertise in the field of nuclear physics
with me.
Dzevad Belkic, for sharing your knowledge in quantum mechanics and
particle collision physics.
Hooshang Nikjoo, for giving me a broader view of research in the field of
medical radiation physics.
Lil Engström, for your kindness and for all help regarding IMPORTANT
stuff.
Marianne Eklund, for all laughs and for your friendship. DiVA is not our
friend!
Marilyn Noz for proofreading.
Karen Belkic, for being such a nice person and for your encouragement.
Pedro Andreo, for sharing your knowledge in dosimetry.
The PhD students at Department of Medical Radiation Physics: Minna
84
Wedenberg, Sara Strååt, Thiansin Liamsuwan, Reza Taleei, Laura Antonovic,
Martha Lazzeroni, Tommy Sundström, Eleftheria Alevronta, Katarzyna
Zielinska-Chomej, Chitralehka Mohanty
My former colleagues at the Department of Medical Radiation Physics:
Martha Hultqvist, my roomie, we shared joy, desperation, and Freja! It is
fast, faster than....
Susanne Larsson, for your friendship. We had a lot of fun, mostly outside
MSF, but also during work hours.
Malin Hollmark, for your encouragement.
Annelie Meijer, for being such a good friend.
Björn Andreassen, for making a sunday at MSF really fun.
Johanna Kempe, Patrick Vreede, Magdalena Adamus-Gorka, Bartosz Gorka,
Richard Holmberg, Bahram Andisheh, Malin Siddiqi, Judith Aldana, Robert
Rosengren
Ola Wimo, you say I’m solitarily responsible for your mental health...I
would like to say it is the opposite. What is the hen and what is the egg?
My friends outside work and my family, thank you for your love, patience
and understanding.
My dad, for your encouragement and for your reminders that I’m not the
type of person that gives up.
My mum, for your endless love and for keeping me aware of other things
apart from work.
Thomas, thank you for making my life better and brighter. You complete
me! Vi är som samma!
85
BIBLIOGRAPHY
Bibliography
Abrines R and Percival IC 1966a Classical theory of charge transfer and ionization of hydrogen atoms by protons Proc. Phys. Soc. 88 861–872
Abrines R and Percival IC 1966b A generalized correspondence principle and
proton-hydrogen collisions Proc. Phys. Soc. 88 873–883
Belkic D 1978 A quantum theory of ionisation in fast collisions between ions
and atomic systems J. Phys. B: Atom. Molec. Phys. 11 3529–3552
Belkic D 2010 Review of theories on double electron capture in fast ion-atom
collisions J. Math. Chem. 47 1420–1467
Belkic D, Mancev I and Hanssen J 2008 Four-body methods for high-energy
ion-atom collisions Rev. Mod. Phys. 80 250–311
Belli M, Cherubini R, Finotto S, Moschini G, Sapora O, Simone G and Tabocchini MA 1989 RBE-LET relationship for the survival of V79 cells irradiated
with low energy protons Int. J. Radiat. Biol. 55 93–104
Benedito E, Fernández-Varea JM and Salvat F 2001 Mixed simulation of the
multiple elastic scattering of electrons and positrons using partial-wave
differential cross section Nucl. Intrum. Methods B 174 91–110
Beuve M 2009 Formalization and theoretical analysis of the Local Effect
Model Radiat. Res. 172 394–402
Bielajew A 2012 Fundamentals of the Monte Carlo method for neutral and
charged particle transport
http://www-personal.umich.edu/∼bielajew/MCBook/book.pdf
Bonsen TFM and Vriens L 1970 Angular distribution of electrons ejected by
charged particles: I. Ionization of He and H2 by protons Physica 47 307–
319
Bonsen TFM and Vriens L 1971 Angular distribution of electrons ejected
by charged particles: II. the binary-encounter polarized-orbital method
Physica 54 318–324
Brahme A 1975 Simple relations for the penetration of high en energy electron beams in matter Technical report National Institute of radiation Protection, Stockholm, Sweden, Internal Report
86
BIBLIOGRAPHY
Brahme A 1984 Dosimetric precision requirements in radiation therapy Acta
Radiol. Oncol. 23 379–391
Brahme A, Lax I and Andreo P 1981 Electron beam dose planning using discrete gaussian beams. mathematical background Acta Radiol. Oncol. 20
147–158
Brahme A and Lind BK 2010 A systems biology approach to radiation therapy optimization. Radiat Environ Biophys 49 111–124
Brenner DJ and Ward JF 1992 Constraints on energy deposition and target
size of multiply damaged sites associated with DNA double-strand breaks
Int. J. Radiat. Biol. 61 737–748
Butts JJ and Katz R 1967 Theory of RBE for heavy ion bombardment of dry
enzymes and viruses Radiat. Res. 30 855–871
Cannell RJ and Watt DE 1985 Biophysical mechanisms of damage by fast
ions to mammalian cells in vitro Phys. Med. Biol. 30 255–258
Cera F, R. C, Vecchia MD, Favaretto S, Moschini G and Tiveron P 1997 Cell inactivation, mutation and DNA damage induced by light ion: Dependence
on radiation quality Microdosimetry -An interdisiplinary approach (The
Royal Society of Chemistry, Cambridge)
Chilton AB 1978 A note on the fluence concept Health Phys. 34 715–716
Chilton AB 1979 Further comments on an alternate definition of fluence
Health Phys. 36 637–638
Crothers DSF and McMann J 1983 Ionization of atoms by ion impact J. Phys
B: Atom. Molec. Phys. 16 3229–3242
Curtis SB 1986 Lethal and potentially lethal lesions induced by radiation–A
unified repair model Radiat. Res. 106 252–270
Dalgarna A and Griffing GW 1955 Energy loss pf protons passing through
hydrogen Proc. Roy. Soc. London Ser. A 232 423–434
Deasy J 1996 Poisson formulas for tumor control probability with clonogen
proliferation Radiat. Res. 145 382–384
Dingfelder M, Hantke D, Inokuti M and Paretzke H 1998 Electron inelasticscattering cross section in liquid water Radiat. Phys. Chem. 53 1–18
Dingfelder M, Ritchie RH, Turner JE, Friedland W, Paretzke HG and Hamm
RN 2008 Comparisons of calculations with PARTRAC and NOREC: transport of electrons in liquid water Radiat. Res. 169 584–594
87
BIBLIOGRAPHY
Elsässer T and Scholz M 2007 Cluster effects within the Local Effect Model
Radiat. Res. 167 319–329
Emfietzoglou D, Kyriakou I, Abril I, Garcia-Molina R and Nikjoo H 2012 Inelastic scattering of low-energy electrons in liquid water computed from
optical-data models of the Bethe surface Int. J. Radiat. Biol. 88 22–28
Emfietzoglou D and Nikjoo H 2007 Accurate electron inelastic cross sections
and stopping powers for liquid water over the 0.1-10 keV range based on
an improved dielectric description of the Bethe surface Radiat. Res. 167
110–120
Fainstein P and Rivarola R 1987 Symmetric eikonal model for ionisation in
ion-atom collisions J. Phys. B: Atom. Molec. Phys. 20 1285–1293
Fernández-Varea J, Mayol M, Liljequist D and Salvat F 1993 Inelastic scattering of electrons in solids from a generalized oscillator strength model
using optical and photoelectric data J. Phys. Condens. Matter 5 3593–3610
Friedland W, Kundrát P and Jacob P 2011 Track structure calculations on
hypothetical subcellular targets for the release of cell-killing signals in
bystander experiments with medium transfer Radiat. Prot. Dosim. 143
325–329
Furness JB and McCarthy IE 1973 Semiphenomenological optical model for
electron scattering on atoms J. Phys B: Atom. Molec. Phys. 6 2280–2291
Galassi M, Davies J, Theiler J, Gough B, Jungman G, Booth M and Rossi F
2009 GNU Scientific Library Reference Manual Reference Manual (3rd edition) ISBN 09546120-7-8
González-Muñoz G, Tilly N, Fernández-Varea JM and Ahnesjö A 2008 Monte
Carlo simulation and analysis of proton energy-deposition patterns in the
Bragg peak Phys. Med. Biol. 53 2857–2875
Goodhead DT 1994 Initial events in the cellular effects of ionizing radiations:
clustered damage in DNA Int. J. Radiat. Biol. 65 7–17
Grosswendt B 2004 Recent advances of nanodosimetry Radiat. Prot. Dosim.
110 789–799
Gryzinski M 1959 Classical theory of electronic and ionic inelastic collisions
Phys. Rev. 115 374–383
Hamm RN, Wright HA, Ritchie RH, Turner JE and Turner TP 1976 Monte
Carlo calculations of electrons through liqid water pp 1037-1053 (5th
Symp. on Microdosimetry, Verbania, Pallanza)
88
BIBLIOGRAPHY
Hansen JP and Kocbach L 1989 Ejection angles of fast-delta electrons from
K-shells ionisation induced by energetic ions J. Phys. B: Atom. Molec. Phys.
22 L71–77
Hawkins RB 1994 A statistical theory of cell killing by radiation of varying
linear energy transfer Radiat. Res. 140 366–374
Hayashi H, Watanabe N, Udagawa Y and Kao C 2000 The complete optical
spectrum of liquid water measured by inelastic x-ray scattering Proc. Natl.
Acad. Sci. U S A 97 6264–6266
Heller JM, Hamm RN, Birkhoff R and Pinter L 1974 Collective oscillation in
liquid water J. Chem. Phys. 60 3483–3486
Hollmark M, Uhrdin J, Dz. B, Gudowska I and Brahme A 2004 Influence of
multiple scattering and energy loss straggling on the absorbed dose distributions of therapeutic light ion beams: I. Analytical pencil beam model
Phys. Med. Biol. 49 3247–3265
Hug O and Kellerer AM 1966 Stochastik der Strahlenwirkung (SpringerVerlag, Berlin Heidelberg New York)
ICRU 1998 Fundamental quantities and units for ionizing radiation (Bethesda, MD)
Itikawa Y and Mason N 2005 Cross sections for electron collisions with water
molecule J. Phys. Chem. Ref. Data 34 1–20
Joiner MC, Marples B, Lambin P, Short SC and Turesson I 2001 Low-dose
hypersensitivity: current status and possible mechanisms Int. J. Radiat.
Oncol. Biol. Phys. 49 379–389
Kanai T, Endo M, Minohara S, Miyahara N, Koyama-ito H, Tomura H, Matsufuji N, Futami Y, Fukumura A, Hiraoka T, Furusawa Y, Ando K, Suzuki M,
Soga F and Kawachi K 1999 Biophysical characteristics of HIMAC clinical
irradiation system for heavy-ion radiation therapy Int. J. Radiat. Oncol.
Biol. Phys. 44 201–210
Katz R, Ackerson B, Homayoonfar M and Sharma SC 1971 Inactivation of
cells by heavy ion bombardment Radiat. Res. 43 402–425
Katz R, Zachariah R, Cucinotta FA and Zhang C 1994 Survey of cellular radiosensitivity parameters Radiat. Res. 140 356–365
Kellerer AM and Chmelevsky D 1975 Concepts of microdosimetry. III. Mean
values of the microdosimetric distributions Radiat. Environ. Biophys. 12
321–335
89
BIBLIOGRAPHY
Kellerer AM and Rossi HH 1978 A generalized formulation of dual radiation
action Radiat. Res. 75 471–488
Kempe J, Gudowska I and Brahme A 2007 Depth absorbed dose and LET
distributions of therapeutic 1 H, 4 He, 7 Li, and 1 2C beams Med. Phys. 34
183–192
Klepper LY 2001 Methods of transition from nonuniform dose distribution
in tissue to adequate dose of uniform irradiation with the same probability of postradiation complications Biomed. Eng. 35. 65–70.
Knuth DE 1997 Seminumerical Algorithms The Art of Computer Programming Vol 2 3rd ed
Lappa AV, Bigildeev EA, Burmistrov DS and Vasilyev ON 1993 "TRION" code
for radiation action calculations and its application in microdosimetry
and radiobiology Radiat Environ Biophys 32 1–19
Lea D 1962 Actions of Radiation on Living cells (Cambridge at the university
press)
Lind BK, Persson LM, Edgren MR, Hedlöf I and Brahme A 2003 Repairableconditionally repairable damage model based on dual poisson processes
Radiat. Res. 160 366–375
Matsumoto M and Nishimura T 1998 Mersenne twister: A 623dimensionally equidistributed uniform pseudo-random generator
ACM T. Model. Comput. S. 8 3–30
Meijer AE, Jernberg ARM, Heiden T, Stenerlöw B, Persson LM, Tilly N, Lind
BK and Edgren MR 2005 Dose and time dependent apoptotic response in
a human melanoma cell line exposed to accelerated boron ions at four
different LET Int. J. Radiat. Biol. 81 261–272
Mekid S and Vaja D 2008 Propagation of uncertainty: Expression of second
and third order uncertainty with third and fourth moments Measurement
41 600–609
Michalik V 1992 Model ofDNA damage induced by radiations of various
qualities Int. J. Radiat. Biol. 62 9–20
Miller J, Toburen LH and Manson S 1983 Differential cross sections for ionization of helium, neon and argon by high-velocity ions Phys. Rev. A 27
1117–1144
Nikjoo H, O’Neill P, Terrissol M and Goodhead DT 1999 Quantitative modelling of DNA damage using Monte Carlo track structure method Radiat.
Environ. Biophys. 38 31–38
90
BIBLIOGRAPHY
Nikjoo H, Uehara S, Emfietzoglou D and Cucinotta F 2006 Track-structure
codes in radiation research Radiat. Meas. 41 1052–1074
Olivera GH, Martínez AE, Rivarola RD and Fainstein PD 1995 Theoretical
calculation of electronic stopping power of water vapor by proton impact
Radiat. Res. 144 241–247
Olivero JJ, Stagat RW and Green AES 1972 Electron deposition in water vapor, with atmospheric applications J. Geophys. Res. 77 4797–811
Ottolenghi A, Ballarini F and Merzagora M 1999 Modelling radiationinduced biological lesions: from initial energy depositions to chromosome aberrations Radiat. Environ. Biophys. 38 1–13
Persson LM, Edgren MR, Stenerlöw B, Lind BK, Hedlöf I, Jernberg ARM,
Meijer AE and Brahme A 2002 Relative biological evectiveness of boron
ions on human melanoma cells Int. J. Radiat. Biol. 78 743–748
Ritchie RH and Howie A 1977 Electron excitation and the optical potential
in electron microscopy Philos. Mag. 36 463–481
Rudd M 1988 Differential cross sections for secondary electron production
by proton impact Phys. Rev. A. 38 6129–6137
Rudd M, Kim YK, Madison D and Gay T 1992 Electron production in proton
collisions with atoms and molecules: energy distributions Rev. Mod. Phys.
64 441–490
Saito M and Matsumoto M 2008 SIMD-Oriented Fast Mersenne Twister: A
128-bit Pseudorandom Number Generator, Monte Carlo and quasi-Monte
Carlo methods pp 607-622 (Springer-Verlag, Berlin)
Salin A 1969 Ionization of atomic hydrogen by proton impact J. Phys. B.:
Atom. Molec. Phys. 2 631–639
Salin A 1972 Ionization of helium by proton impact J.Phys. B.: Atom. Molec.
Phys. 5 979–986
Salvat F, Fernández-Varea JM and Sempau J 2008 PENELOPE-2008:
a code system for monte carlo simulation of electron and photon
transport (Issy-les-Moulineaux: OECD Nuclear Energy Agency)
http://www.nea.fr/lists/penelope.html
Salvat F, Jablonski A and Powell CJ 2005 ELSEPA - dirac partial-wave calculation of elastic scattering of electrons and positrons by atoms, positive
ions and molecule Comput. Phys. Commun. 165 157–190
91
BIBLIOGRAPHY
Scholz M, Kellerer AM, Kraft-Weyrather W and Kraft G 1997 Computation of
cell survival in heavy ion beams for therapy. the model and its approximation Radiat. Environ. Biophys. 36 59–66
Scholz M and Kraft G 1994 Calculation of heavy ion inactivation probabilities based on track structure, X-ray sensitivity and target size Radiat. Prot.
Dosim. 52 29–33
Scholz M and Kraft G 1996 Track structure and the calculation of biological
effects of heavy charged particles Adv. Space. Res. 18 5–14
Sinclair WK 1966 The shape of radiation survival cuves of mammalian cells
cultured in vitro. In:Biophysical aspects of radiation quality. IAEA, Vienna,
Technical reports series no. 58
Singh B, Arrand JE and Joiner MC 1994 Hypersensitive response of normal
human lung epithelial cells at low radiation doses Int. J. Radiat. Biol. 65
457–464
Sontag W 1997 A discrete cell survival model including repair after high
dose-rate of ionizing radiation Int. J. Radiat. Biol. 71 129–144
Stolterfoht N, Dubois RD and Rivarola RD 1997 Electron Emission in Heavy
Ion-Atom Collsions (Springer-Verlag, Berlin Heidelberg New York)
Tilly N, Fernández-Varea JM, Grusell E and Brahme A 2002 Comparison of
Monte Carlo calculated electron slowing-down spectra generated by 60 Co
gamma-rays, electrons, protons and light ions Phys. Med. Biol. 47 1303–
1319
Tobias CA 1985 The repair-misrepair model in radiobiology: comparison to
other models Radiat. Res. Suppl. 8 S77–S95
Tomé WA and Fowler JF 2002 On cold spots in tumor subvolumes Med. Phys.
29 1590–1598
Uehara S, Nikjoo H and Goodhead DT 1999 Comparison and assessment of
electron cross sections for Monte Carlo track structure codes Radiat. Res.
152 202–213
Wang CM and Iyer HK 2005 On higher-order corrections for propagating uncertainties Metrologia 42 406–410
Watt DE 1996 Quantities for Dosimetry of Ionizing Radiations in Liquid Water (Taylor and Francis)
Wedenberg M, Lind BK, Toma-Dasu I, Rehbinder H and Brahme A 2010
Analytical description of the LET dependence of cell survival using the
repairable-conditionally repairable damage model Radiat. Res. 174 517–
525
92
BIBLIOGRAPHY
Weyrather WK, Ritter S, Scholz M and Kraft G 1999 RBE for carbon tracksegment irradiation in cell lines of differing repair capacity Int. J. Radiat.
Biol. 75 1357–1364
Withers H 1992 Biological basis of radiation therapy, Principles and Practice
of Radiation Oncology pp 64-96 (J.B Lippincott Company, Philadelphia)
Ziegler JF, Biersack JP and D M 2008 SRIM - The Stopping and Range of Ions
in Matter.
Zimmer KG 1961 Studies on Qantitative Radiation Biology. (Oliver and Boyd,
Edingburgh and London)
93