The physics of proton therapy.
Journal: 2015/December - Physics in Medicine and Biology
ISSN: 1361-6560
Abstract:
The physics of proton therapy has advanced considerably since it was proposed in 1946. Today analytical equations and numerical simulation methods are available to predict and characterize many aspects of proton therapy. This article reviews the basic aspects of the physics of proton therapy, including proton interaction mechanisms, proton transport calculations, the determination of dose from therapeutic and stray radiations, and shielding design. The article discusses underlying processes as well as selected practical experimental and theoretical methods. We conclude by briefly speculating on possible future areas of research of relevance to the physics of proton therapy.
Relations:
Content
Citations
(35)
References
(158)
Chemicals
(1)
Organisms
(1)
Affiliates
(1)
Similar articles
Articles by the same authors
Discussion board
Phys Med Biol 60(8): R155-R209

The physics of proton therapy

1. Introduction

The history of proton therapy began in 1946 when Robert Wilson published a seminar paper in which he proposed to use accelerator-produced beams of protons to treat deep-seated tumors in humans (Wilson, 1946). In that paper, he explained the biophysical rationale for proton therapy as well as the key engineering techniques of beam delivery. In 1954, the first human was treated with proton beams at the Lawrence Berkeley Laboratory (Lawrence et al, 1958). In 1962, specialized radiosurgical proton treatments commenced at the Harvard Cyclotron Laboratory (Kjellberg et al, 1962a; Kjellberg et al, 1962b), followed in the mid 1970s by treatments for ocular cancers (Gragoudaset al, 1982) and larger tumors (Koehler et al, 1977). Physicists at Harvard, collaborating with clinical colleagues at Massachusetts General Hospital, the Massachusetts Eye and Ear Infirmary, and elsewhere, developed much of the physics and technology needed to treat patients with proton beams safely and effectively. Remarkably, the research and development program at Harvard continued for more than 40 years (Wilson, 2004). During the same period, physicists elsewhere were developing other key technologies, including accelerators, magnetically scanned beams, treatment planning systems, computed tomographic imaging (CT), and magnetic resonance imaging.

The widespread adoption of proton therapy has been slow in comparison to, for example, intensity-modulated photon therapy. There are several reasons for this slow adoption of proton therapy, including technical difficulty, cost, and lack of evidence of cost-competitiveness. Commercial proton delivery systems had been contemplated for decades before they finally appeared in 2001 after overcoming considerable difficulties. The cost of proton therapy equipment remains much higher than that of comparable conventional photon therapy equipment; the long-anticipated economies of scale have not, as yet, materialized. Even in times of relative prosperity, the allocation of scarce resources to proton therapy has been constrained by relatively sparse evidence of its cost-competitiveness and cost-effectiveness (Goitein and Jermann, 2003; Peeters et al, 2010; Lievens and Pijls-Johannesma, 2013).

Despite these obstacles, much progress has been made. Today there are 10 proton therapy centers in operation in the United States and 38 centers worldwide (PTCOG, 2013). The Particle Therapy Cooperative Group (PTCOG) reported that at least 93,895 patients had been treated worldwide by March, 2013 (PTCOG, 2013). The proton therapy community has stepped up efforts to conduct clinical trials that compare outcomes after proton therapy with those after other advanced technology radiation therapies (Duttenhaver et al, 1983; Shipley et al, 1995; Desjardins et al, 2003; Zietman et al, 2010).

The central rationale for proton therapy is its superior spatial dose distribution in the patient. In recent years, the advantage of protons over photons in providing a highly conformal and uniform dose to a tumor has been largely diminished by advances in photon therapies, such as intensity-modulated photon therapy and volumetric arc therapy (Weber et al, 2009). However, the relative advantage of proton therapy in sparing normal tissues has never been more apparent or important; in the United States, approximately 65% of adults and 80% of children survive 5 years after their cancer diagnosis (Valdivieso et al, 2012). About half of cancer patients receive radiotherapy as part of their treatment. Recent studies reported that the incidence of treatment-related morbidity, including second cancers, cardiovascular disease, fertility complications, and other late effects, is alarmingly high in long-term survivors of cancer (Wilson et al, 2005; Carver et al, 2007; Armstrong et al, 2009; Merchant, 2009; Sauvat et al, 2009; Newhauser and Durante, 2011; Olch, 2013). Presently, about 3% of the U.S. population are cancer survivors, corresponding to 11 million people, a figure projected to grow to 18 million by 2022 (de Moor et al, 2013) For these reasons, there is increasing interest in exploiting the tissue-sparing capabilities inherent to proton therapy to reduce the burden of treatment-related complications on patients and the healthcare system. An understanding of the physics and biology of radiogenic late effects from proton therapy has started to emerge in the literature in the last decade.

This paper reviews the basic, essential physics underlying proton therapy. The literature includes excellent reviews of various aspects of proton therapy physics, most notably the comprehensive works of Chu et al. (1993), and also those of Bonnett (Bonnett, 1993), Pedroni et al. (1995), Brahme (2004), Lomax (2009), Coutrakon (2007), and Schardt and Elsasser (2010) In addition, several reports (ICRU, 1998, 2007) and textbooks (Scharf, 1994; Breuer and Smit, 2000; DeLaney and Kooy, 2008; Linz, 2012; Paganetti, 2012; Ma and Lomax, 2013) have covered various aspects of proton therapy physics. In recent years, many textbooks dealing with general radiation oncology have included relevant chapters on proton therapy (Van Dyk, 1999; Halperin et al, 2008; Pawlicki et al, 2011). Many of the older works on proton therapy have withstood the test of time and remain excellent literature resources of continued relevance. However, in this review we mainly focus on the well-established basic physics of proton therapy and on selected advances from the last 15 years or so that are important in clinical proton therapy. In choosing advances to cover in this review, considerable selectivity was necessary because of the huge expansion of the proton therapy literature, especially in recent years. To the authors of the many studies that we were not able to mention in this review because of space limitations, we apologize and we appreciate your understanding.

2. Proton interaction mechanisms

In this section, we briefly review the predominant types of interactions of protons in matter and why they are important. Figure 1 illustrates several mechanisms by which a proton interacts with an atom or nucleus: Coulombic interactions with atomic electrons, Coulombic interactions with the atomic nucleus, nuclear reactions, and Bremsstrahlung. To a first-order approximation, protons continuously lose kinetic energy via frequent inelastic Coulombic interactions with atomic electrons. Most protons travel in a nearly straight line because the proton mass is 1832 times greater than that of an electron. In contrast, a proton passing close to the atomic nucleus experiences a repulsive elastic Coulombic interaction which, owing to the large mass of the nucleus, deflects the proton from its original straight-line trajectory. Inelastic nuclear reactions between protons and the atomic nucleus are less frequent but, in terms of the fate of an individual proton, have a much more profound effect. In a nuclear reaction, the projectile proton enters the nucleus; the nucleus may emit a proton, deuteron, triton, or heavier ion or one or more neutrons. Finally, proton Bremsstrahlung is theoretically possible, but at therapeutic proton beam energies this effect is negligible. Table 1 summarizes the proton interaction types, interaction targets, principal ejectiles, influence on the proton beam, and dosimetric manifestations. We review these interaction mechanisms, except proton Bremsstrahlung, in the following sections.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0001.jpg

Schematic illustration of proton interaction mechanisms: (a) energy loss via Coulombic interactions, (b) deflection of proton trajectory by repulsive Coulomb scattering with nucleus, (c) removal of primary proton and creation of secondary particles via non-elastic nuclear interaction. (p: proton, e: electron, n: neutron, He: Helium, γ: gamma rays)

Table 1

Summary of proton interaction types, targets, ejectiles, influence on projectile, and selected dosimetric manifestations

Interaction TypeInteraction TargetPrincipal EjectilesInfluence on ProjectileDosimetric Manifestation
Inelastic Coulomb scatteringAtomic electronsPrimary proton, ionization electronsQuasi-continuous energy lossEnergy loss determines range in patient
Elastic Coulomb scatteringAtomic nucleusPrimary proton, recoil nucleusChange in trajectoryDetermines lateral penumbral sharpness
Non-elastic nuclear reactionsAtomic nucleusSecondary protons and heavier ions, neutrons, and gamma raysRemoval of primary proton from beamPrimary fluence, generation of stray neutrons, generation of prompt gammas for in vivo interrogation
BremsstrahlungAtomic nucleusPrimary proton, Bremsstrahlung photonEnergy loss, change in trajectoryNegligible

2.1 Energy loss rate

The energy loss rate of ions, or linear stopping power, is defined as the quotient of dE and dx, where E is the energy and x is the distance. It is frequently more convenient to express the energy loss rate in a way that is independent of the mass density; the mass stopping power is defined as

Sρ=dEρdx,
(1)

where ρ is the mass density of the absorbing material. Please note that stopping power is defined for a beam, not a particle.

The energy loss rate may be described by several mathematical formulae. The simplest, yet still remarkably accurate, formula is based on the Bragg-Kleeman (BK) Rule (Bragg and Kleeman, 1905), which was originally derived for alpha particles, and is given by

Sρ=dEρdxE1pραp,
(2)

where ρ is the mass density of the material, α is a material-dependent constant, E is the initial energy of the proton beam, and the exponent p is a constant that takes into account the dependence of the proton's energy or velocity. Values of α and p may be obtained by fitting to either ranges or stopping power data from measurements or theory.

A more physically complete theory, developed by Bohr (Bohr, 1915), is based on calculation of the momentum impulse of a stationary, unbound electron and the impact parameter. A more accurate formula, attributed to Bethe and Bloch (Bethe, 1930; Bloch, 1933), takes into account quantum mechanical effects and is given by

Sρ=dEρdx=4πNAre2mec2ZAz2β2[ln2mec2γ2β2Iβ2δ2CZ],
(3)

where NA is Avogadro's number, re is the classical electron radius, me is the mass of an electron, z is the charge of the projectile, Z is the atomic number of the absorbing material, A is the atomic weight of the absorbing material, c is speed of light, β = v/c where v is the velocity of the projectile, γ = (1- β), I is the mean excitation energy of the absorbing material, δ is the density correction item arising from the shielding of remote electrons by close electrons and will result in a reduction of energy loss for higher energies, and C is the shell correction item, which is important only for low energies where the particle velocity is near the velocity of the atomic electrons. The two correction items in the Bethe-Bloch equation involve relativistic theory and quantum mechanics and need to be considered when very high or very low proton energies are used in calculations. Figure 2 plots proton stopping power as a function of proton energy in water calculated by using Equation (3).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0002.jpg

Mass stopping power (S) versus ion energy (E) for protons in liquid water. The corresponding range (R), calculated using the plotted S values and on the assumption of the continuous slowing down approximation (CSDA), is also plotted.

It is instructive to observe in Equation (3) how the projectile's characteristics govern its energy loss rate: energy loss is proportional to the inverse square of its velocity (1/v classically and 1/β relativistically) and the square of the ion charge (z = 1 for protons), and there is no dependence on projectile mass. Similarly, Equation (3) reveals that the absorber material can also strongly influence the energy loss rate. Specifically, the linear stopping power is directly proportional to the mass density. It is equivalent, but perhaps more physically meaningful, to state that the linear stopping power is proportional to the density of electrons in the absorber (NAρ Z/A), because the energy loss occurs by Coulombic interactions between the proton and atomic electrons. Z/A varies by only about 16%, from 0.5 for biologic elements such as carbon and oxygen to 0.42 for high-Z beamline components, such as lead. Hydrogen is an obvious exception to this; fortuitously, the concentration of hydrogen in the human body is low (only about 10%) and nearly uniform throughout the body. The stopping power also depends on a material's I value, and the I value depends in a monotonic way on the Z of the absorber, varying from about 19 eV for hydrogen to about 820 eV for lead. However, the stopping power goes with the logarithm of I value, so the dependence is diminished. Hence, putting these various dependencies in perspective, it is clear that the proton energy loss rate in the human body depends most strongly on the material density, which can vary by about three orders of magnitude, from air in the lung to cortical bone, and the ion velocity, which can cause the linear stopping power in water to vary by about a factor of 60 for proton energies between 1 MeV and 250 MeV.

The stopping power theory for protons and heavier ions was reviewed by Ziegler et al. (Ziegler et al, 1985; Ziegler, 1999; Ziegler et al, 2008) and in Report 49 of the International Commission on Radiation Units and Measurement (1993). Evaluated stopping power and range tables may be conveniently calculated with the SRIM code (“Stopping and Ranges of Ions in Matter,” computer program, http://www.srim.org/.).

Thus far, we have described proton energy loss in an approximate way on the assumptions that a proton loses energy along a 2-dimensional line trajectory and that energy loss is continuous. Absorption of this same energy, however, occurs in a 3-dimensional volume. Furthermore, the ionization track of a proton has an irregular 3-dimensional structure caused by random fluctuations in the location and size of primary ionization events. This is caused mainly by proton-produced recoil electrons, some of which are sufficiently energetic to create small spur tracks of ionization emanating from the main track. Because the electrons are very much lighter than the protons, each interaction can reduce the proton energy only a little. The maximum possible energy transfer in a collision of an ion of mass m with an unbound stationary electron is

Δ1max=2mec2β21β2[1+2meM11β2+(meM)2],
(4)

where me is the mass of an electron, M is the mass of the target material, c is the speed of light, and β = v/c where v is the velocity of the projectile.

Even for very energetic protons, the secondary electrons do not acquire enough energy to travel more than a few millimeters from the proton track. For example, at 200 MeV proton energy, the maximum secondary electron energy is around 500 keV, which corresponds to an electron range of approximately 2 mm in water. The probability of producing secondary electrons may be calculated with various total or differential cross-sections; these were reviewed in ICRU Report 55 (1995). Track structure models may be used to estimate the radial properties of ions (Kraft et al, 1999), although this has not found common application in clinical proton therapy. Regardless of the calculation method used, the spatial characteristics of secondary charged particles should, in principle, be taken into account near material interfaces (e.g., buildup effects in transmission beam monitoring instruments, skin, air-tumor interfaces in the lung) and in cases where the radiation quality is of interest (e.g., microdosimetric and nanodosimetric characterization of individual dose deposition events).

Finally, we note that in proton therapy water is considered an excellent tissue substitute because of its similar density, effective Z/A, and other properties. Furthermore, proton energy loss and residual range in various materials are often expressed in terms of their water-equivalent values. We discuss this in detail in section 3.1, but until then, for simplicity, we consider water as being perfectly equivalent to tissue.

2.2 Range

Range is defined as the depth at which half of protons in the medium have come to rest, as shown in the range-number curve plotted in Figure 3. There are small variations in the energy loss of individual protons (an effect called range straggling, which is discussed in the following section). Consequently, the range is inherently an average quantity, defined for a beam and not for individual particles. By convention, this means half of the protons incident on the absorber are stopped, although in some cases this is taken instead to mean that half the protons survive to near the end of range (i.e., neglecting protons removed by nuclear reactions).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0003.jpg

Relative fraction of the fluence Φ in a broad beam of protons remaining as a function of depth z in water. The gradual depletion of protons from entrance to near the end of range is caused by removal of protons from nuclear reactions. The rapid falloff in the number of protons near the end of range is caused by ions running out of energy and being absorbed by the medium. The sigmoid shape of the distal falloff is caused by range straggling or by stochastic fluctuations in the energy loss of individual protons.

The path of most protons in matter is a nearly straight line. On average, the proton's pathlength is very nearly equal to its projected pathlength and range. This simple but important fact renders many proton range calculations tractable with relatively simple numerical or analytical approaches.

Let us first consider a simple numerical calculation of proton beam range. We use proton stopping power data and perform a one-dimensional proton pathlength transport calculation on the assumptions that the ions travel only straight ahead (negligible lateral scattering) and that the protons lose energy in a continuous matter. (These assumptions are reasonable for many clinical calculations, but we examine then relax these assumptions in a later section.) In this case, the range (R) may be calculated as

R(E)=0E(dEdx)1dE0E(dEdx)1dE,
(5)

where E is the ion's initial kinetic energy. The summation denotes that the continuous transport is approximated by calculations of discrete steps. In fact, as discussed above, this equation actually gives the pathlength, which is an excellent approximation of range in most clinical situations. Figure 2 plots proton range in water calculated by using Equation (5).

Next, we calculate the proton range using an analytical approach, which may be faster and more practical than the numerical approach for many clinical calculations. We begin by noting that the interval of proton range of interest typically extends from 1 mm (about the size of a voxel in an anatomic image of patient anatomy) to about 30 cm (about the midline of a large adult male's pelvis, the deepest site in the human body). These ranges correspond to 11 MeV and 220 MeV, respectively, as seen in Figure 2. More importantly, Figure 2 reveals that the relationship between the logarithm of range and logarithm of energy is almost linear. This is fortuitous because this means that the range follows a very simple power law, as realized by Bragg and Kleeman (1905) and others early in the last century. Thus, the range of a proton may be calculated using BK Rule, or

R(E) = αE,
(6)

where, as before, α is a material-dependent constant, E is the initial energy of the proton beam, and the exponent p takes into account the dependence of the proton's energy or velocity.

The uncertainty in proton range depends on many factors. For example, the uncertainty in a range measurement depends on the precision and accuracy of the measurement apparatus and, in some cases, on the experimenter's skill. A common concern in clinical proton therapy is the uncertainty in calculated range, e.g., in calculating the settings of the treatment machine for a patient's treatment. The uncertainty in range may depend on the knowledge of proton beam energy distribution and on properties of all range absorbing materials in the beam's path. These include elemental composition, mass density, and linear stopping power. The linear stopping powers deduced from computed tomography scans have many additional sources of uncertainty, including imperfections in the calibration of CT scanners (i.e., the method used to convert image data from Hounsfield Units to linear stopping power), partial volume effects, motion artifacts, and streaks artifacts.

2.3 Energy and range straggling

In the preceding sections, we approximated the energy loss rate by assuming that the slowing of ions occurs in a smooth and continuous manner. In fact, we considered the mean energy loss rate and neglected variations in the energy loss rates of individual protons. For many clinical calculations, these assumptions are valid and lead to a reasonably good first-order approximation. However, the accumulation of many small variations in energy loss, termed energy straggling or range straggling, is one of the physical processes that strongly governs the shape of a proton Bragg curve, a subject that is covered in section 3.2. Thus, an understanding of range straggling is key to understanding the characteristics of proton dose distributions.

Figure 4 plots the relative energy loss probability density functions (PDF) for protons transmitted through water absorbers of various thickness. The curves have been normalized to enhance visual clarity. Apparently, thick absorbers result in a symmetric distribution of energy losses, whereas thin absorbers yield curves that are asymmetric with modes less than the mean and long tails of large energy losses. In principle, straggling PDFs may be calculated numerically from first principles, but usually theoretical approaches are used. Later in the section, three of the most commonly used straggling theories are described.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0004.jpg

Energy loss probability density functions (PDF) are plotted for various thicknesses of water absorbers, where the thickness is expressed in units of mean free path (mfp). For visual clarity, the energy-loss PDFs have been scaled on both the abscissa and ordinate. The single event energy loss is expressed as a fraction of the mean energy lost in the entire absorber thickness, or (Δ-Δav)/Δav. Each PDF was scaled so that the integral over all energy-loss values yields unit value. For thin absorbers (curves a-e), the PDFs are broader and asymmetric and are modeled with the Vavilov (1957) or Landau (1944) theories. For thick absorbers (curve f), the PDFs are symmetric and well-approximated with Bohr's theory (1915), i.e., a Gaussian distribution. (ICRU, 1993)

Having examined the mean energy loss rate in modest detail (section 2.1), and having conceptually introduced energy straggling, it is instructive to understand how these effects are related mathematically before delving into straggling theory. To understand these relationships, let us consider the moments of the ion energy PDF, or

Mn=0Δ1maxΔnf(Δ)dΔ,
(7)

where Δ is the energy loss of an ion in traversing an absorber, f(Δ)dΔ is the probability of energy loss occurring in the interval from Δ to Δ + dΔ, and n is the order of the moment. The zeroth moment is the total collision cross section, the first moment is the (mean) electronic stopping power, the second moment corresponds to the width (variance) of the energy straggling distribution, the third moment to its asymmetry (skewness), and the fourth to its kurtosis. The variance, sometimes denoted by σΔ, or second moment of the straggling distribution f(Δ), is

M2(Δ)=σΔ2=v0Δ1maxΔnf1(Δ)dΔ.
(8)

A more detailed discussion of the straggling moments was given by Rossi and Zaider (1996).

Next we examine theories for calculating energy straggling proposed by Bohr (1915), Landau (1944), and Vavilov (1957). These theories are valid for thick, intermediate, and thin absorbers, respectively. The criterion for selecting a valid theory for a given absorber thickness is based on a reduced energy parameter,

k=ξΔ1max,
(9)

where ξ is the approximate mean energy loss and Δ1 is the maximum energy loss possible in a single event.

According to Bohr's theory, the energy straggling distribution behaves according to a Gaussian PDF, or

f(Δ)dΔ=1σΔ2πexp(ΔΔ)22σΔ2,
(10)

where for non-relativistic ions the variance is given by

σΔ2=2πre2mec2z2NZβ2Δ1maxρx.
(11)

Bohr's theory assumes that the absorber is thick enough that there are many individual collisions (i.e., the central limit theorem holds), that the ion velocity does not decrease much in crossing the absorber, and that the absorber is made of unbound electrons. For most applications in proton therapy, Bohr's theory provides adequate accuracy. Several authors have reported convenient power law approximations to estimate sigma as a function of the proton beam range (Chu et al, 1993), or

σΔkR0m,
(12)

where R0 is the range in water in centimeters for a mono-energetic proton beam, k is a material-independent constant of proportionality, and the exponent is empirically determined. For protons in water, k = 0.012 and m =0.935 (Bortfeld, 1997).

Landau's theory relaxed the assumption that the central limit holds, i.e., there are relatively fewer individual collisions in an intermediate thickness absorber, and used an approximate expression for Δ1.

In this case, we have

f(Δ,ρx)dΔ=1ξϕL(λL),
(13)

where the parameter φL (λL) roughly corresponds to the deviation from the mean energy loss and was defined by Landau as

ϕL(λL)=1π0ey(lnyλL)sin(πy)dy.
(14)

Evaluation of the integral in Equation (14) is straightforward and the reader is referred to the literature for details (Seltzer and Berger, 1964).

Vavilov's theory is in essence a generalization of Landau's theory that utilizes the correct expression for Δ1 and is given by

f(Δ,ρx)dΔ=1ξϕV(λV,k,β2)dλV,
(15)

where

ϕV(λV,k,β2)=1πek(1+β2γ)0ekf1cos(λVy+kf2)dy,
(16)

and

λV=ΔΔΔ1maxk(1+β2γ),
(17)

where γ is Euler's constant. The evaluation of Vavilov's theory is computationally more expensive than that of Bohr's or Landau's. For additional details, the reader is referred to Vavilov's original work (Vavilov, 1957).

2.4 Multiple Coulomb scattering

As illustrated in Figure 5, a proton passing close to the nucleus will be elastically scattered or deflected by the repulsive force from the positive charge of the nucleus. While the proton loses a negligible amount of energy in this type of scattering, even a small change in its trajectory can be of paramount importance. In fact, it is necessary to take into account Coulomb scattering when designing beamlines and treatment heads (Gottschalk, 2010b) and in calculations of dose distributions in phantoms or patients, e.g., with treatment planning systems (Hong et al, 1996; Pedroni et al, 2005; Schaffner, 2008; Koch and Newhauser, 2010).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0005.jpg

Schematic diagram showing the trajectory of a proton undergoing multiple Coulomb scattering events. θ is the root mean square (rms) space scattering angle and θx is the projected rms scattering angle. (Leo, 1994)

To characterize the amount that a beam is deflected by scattering, we use the quantity of scattering power, which is defined as

T = d〈θ〉 ∕ dx,
(18)

where <θ> is the mean squared scattering angle and x is the thickness of absorber through which the proton has traveled. Similarly, the mass scattering power is simply T/ρ, where ρ is the mass density of the absorber material. Notice that the definition of scattering power utilizes the mean square of the scattering angle; scattering is symmetric about the central axis, and therefore the mean scattering angle is zero and contains no useful information. Also, note that the value of the scattering power depends on the material thickness (i.e., it depends on the material properties and dimensions of the absorber being considered).

Predictions of elastic Coulomb scattering are commonly classified according to the number of individual scattering events (Ns) that occur in a given absorber. For single scattering (Ns=1), Rutherford scattering theory applies. Plural scattering (1 < Ns < 20) is difficult to model theoretically and is not discussed further here. Multiple Coulomb scattering (MCS; Ns ≥ 20), the combined effect of all Ns scattering events, may be modeled by using a statistical approach to predict the probability for a proton to scatter by a net angle of deflection, commonly denoted by θ (Figure 5).

It is instructive to briefly examine Rutherford's theory of single scattering (Rutherford, 1911). The differential cross section dσ for scattering into the solid angle dΩ is given by

dΩdΩ=z12z22re2(mecβp)24sin4(θ2),
(19)

where z1 is the charge of the projectile, z2 is the atomic number of the absorber material, re is the classical electron radius, me is the mass of the electron, c is the speed of light, β = v/c, p is the ion momentum, and θ is the scattering angle of the proton. The angular dependence is governed by the term 1/sin(θ/2), i.e., in individual scattering events, the proton is preferentially scattered in the forward direction, at very small angles.

In clinical proton therapy, most objects of interest are thick enough to produce a great many scattering events. Thus, for clinical proton therapy, we are usually more interested in the net effect that many small-angle scattering events have on many protons, e.g., how beamline scattering in the treatment head influences the spatial distribution of dose in a patient.

Rigorous theoretical calculations of MCS are quite complex. The most complete theory was proposed by Molière (1948). Assuming that scattered particles are emitted at small deflection angles (i.e., the small-angle approximation in which sin(θ) ≈θ),

P(θ)=ηdη(2exp(η2)+F1(η)B+F2(η)B2),
(20)

where η = θ / (θ1B), and θ1=0.3965(zZpβ)ρδxA. The functions Fk(η) are defined as

Fk(η)=1k!J0(ηy)exp(y24)[y24ln(y24)]kydy,
(21)

where Jo(ηγ) is a Bessel function. Values of Fi have been tabulated in the literature by several authors. Parameter B in Equation (20) may be found by numerical methods, by solving

g(B) = ln BB + ln γ − 0.154 = 0,
(22)

where γ = 8831 q zρδx/(βAΔ) and Δ= 1.13+3.76 (Zz/137β). The following are symbols representing properties of the absorber: Z is the atomic number, A is the atomic mass, δx is the thickness, and ρ is the mass density. The proton momentum is denoted by p, β = v/c, and z = 1 is the proton charge. Even though we have not presented all of the details, clearly the evaluation of Molière's theory is indeed complex. Consequently, considerable attention has been paid in the literature to developing simpler formulae; the simplicity usually comes at the cost of reduced accuracy in modeling scattering at large angles. Gottschalk (1993) discussed several of these methods in detail. One approximation method takes the form of a Gaussian distribution,

P(θ)2θθ2exp(θζ2θ2)dθ,
(23)

where (<θ>) is the root mean square (rms) scattering angle or the width parameter of the Gaussian distribution.

Gottschalk (2010b) proposed a differential approximation to Molière's theory to predict the scattering power according to

TdM=fdM(pv,p1v1)×(ESpv)21XS,
(24)

where ES =15 MeV, p is the momentum and v is the velocity of the proton, p1 and v1 are the initial momentum and velocity, XS is scattering length of the material, and fdM is a material-independent nonlocal correction factor given by

fdM = 0.5244 + 0.1975 lg(1 − (pvp1v1)) + 0.2320 lg(pv) − 0.0098 lg(pv)lg(1 − (pvp1v1))
(25)

The factor fdM takes into account the accumulation of scattering that occurs as the proton slows from v1 to v and is material independent. The scatter length is given by

1ρXS=αNre2Z2A(2log(33219(AZ)131),
(26)

where α is the fine structure constant, N is Avogadro's number, re is electron radius, and Z, A, and ρare the atomic number, atomic weight, and density of the target material.

It is sometimes convenient to project the expected scattering angle < θ > in three-dimensional space to the corresponding value in a plane, denoted < θ x>, according to

θx = 〈θ ∕ 2,
(27)

Also, in proton therapy the lateral displacement (r) of a proton beam is typically calculated from the scattering angle. Under the Gaussian approximation, we have

P(r)dr = 6r(〈θt) exp(−3r ∕ (〈θt)dr,
(28)

where t = x/Lrad, and Lrad is the radiation length of the material, which can be calculated easily or looked up in tables. The mean square lateral displacement is given by

r〉 = 〈θt ∕ 3,
(29)

A power law approximation for rms of protons is

rrms = aR,
(30)

where a is a unitless material constant, R is the proton beam range (in the same units as rrms), and the exponent b governs the range dependence. For protons in water, a = 0.0294 and b = 0.896 (Chu et al, 1993). Koehler and Preston derived a convenient expression to calculate rrms as a function of depth in an absorber and knowledge of its maximum values, rmax, at the end of range (from an unpublished manuscript; some portions of their work were reported by Gottschalk (2010b)).

In proton therapy, MCS in the treatment head (i.e., in the scattering foils) is helpful because it allows the beam to be spread laterally to useful dimensions, e.g., to make a beam laterally large and flat so that a tumor may be completely covered with a uniform dose. Scattering foils in the nozzle are carefully designed to utilize MCS and energy loss to produce clinically useful beams (Koehler et al, 1977; Gottschalk, 2004) However, MCS in the treatment head and in the patient blurs lateral penumbral sharpness. This is manifested as penumbral growth at the edge of collimated beams and/or the growth of the lateral spot size of a scanned beam (Figure 6). Understanding and preserving penumbral sharpness is the key to realizing the full benefit of proton therapy for sparing healthy tissue.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0006.jpg

Broadening of the beam width in water due to multiple Coulomb scattering. (Pedroni et al, 2005)

Recent studies have revealed that MCS plays an important role in proton dose distribution around small implanted metal objects. Specifically, implanted fiducial markers for image-guided patient alignment have been used in proton therapy for many years (Gall et al, 1993; Welsh et al, 2004; Newhauser et al, 2007a; Newhauser et al, 2007c; Ptaszkiewicz et al, 2010; Matsuura et al, 2012). Substantial recent improvements in on-board imaging systems, patient positioning, and patient immobilization have led to increased use of radiopaque implanted fiducial markers in proton therapy to many disease sites, with the goal of improving target coverage and/or normal tissue sparing. However, recent studies revealed that some commonly used markers, even those less than 1 mm in size, can cause severe cold spots (Figure 7), compromising target coverage (Newhauser et al, 2007a; Newhauser et al, 2007c; Carnicer et al, 2013). The severity of the cold spots varies with fiducial size, material composition, and mass density. These parameters, in turn, determine the amount of MCS and energy loss in the fiducial and, hence, perturbations to the dose distribution in surrounding tissue. In essence, MCS in the fiducial causes a downstream dose shadow that may be partially or fully filled in by MCS in the surrounding tissue. While MCS is important, energy loss in the fiducial (or its water-equivalent thickness), the proximity of the fiducial to the end of the proton beam range, and of course its size and orientation with respect to the beam also should be considered. The physics of dose perturbations is explained in detail elsewhere (Newhauser et al., (2007a), and several subsequent studies (Giebeler et al, 2009; Lim et al, 2009; Cheung et al, 2010; Huang et al, 2011) have shown that it is possible to achieve good radiographic visibility using novel markers that do not significantly perturb the therapeutic dose distribution in tissue.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0007.jpg

Small implanted fiducial markers can create clinically significant dosimetric cold spots in proton therapy beams. (Upper) two-dimensional dose distribution as a function of depth in water (z) and crossfield position (x) from a Monte Carlo simulation of range-modulated proton beam incident on a water phantom containing tantalum localization markers oriented (a, b) parallel to the beam axis and (c, d) perpendicular to the beam axis. The range and modulation width are typical for uveal melanoma treatments. (Lower) simulated absorbed dose (D) as a function of depth (z) in the water phantom at various off-axis positions. The perturbed depth dose profiles are parallel to the beam axis and pass through the center of markers a–d in the plot above. For visual clarity, portions of the perturbed dose profiles upstream of the markers are not shown. An unperturbed beam is plotted with open squares. (Newhauser et al, 2007c)

Others have examined the effects of MCS in larger metal objects on clinical proton beams and characterized the suitability of approximate methods to predict MCS in practical clinical applications (Herault et al, 2005; Stankovskiy et al, 2009; Newhauser et al, 2013).

2.5 Nuclear interactions

In addition to the mechanisms already described, protons may interact with the atomic nucleus via non-elastic nuclear reactions in which the nucleus is irreversibly transformed, e.g., a reaction in which a proton is absorbed by the nucleus and a neutron is ejected, denoted by (p,n). The main effect of nuclear reactions within a therapeutic region of a proton field is a small decrease in absorbed dose due to the removal of primary protons, which is compensated to a large extent by the liberation of secondary protons and other ions. In this section, we discuss this and several other important aspects of nuclear reactions.

Before discussing reaction mechanisms, it is instructive to examine a range-number curve (Figure 3), which plots the remaining number of protons versus depth in an absorber as a beam comes to rest. The gradual depletion of protons from entrance to near the end of range is caused by removal of protons from nuclear reactions. The rapid falloff in the number of protons near the end of range is caused by ions running out of energy and being absorbed by the medium. The sigmoid shape of the distal falloff is caused by range straggling or by stochastic fluctuations in the energy loss of individual protons.

To enter the nucleus, protons need to have sufficient energy to overcome the Coulomb barrier of the nucleus, which depends on its atomic number. The total non-elastic cross-section for proton-induced nuclear reactions has a threshold, on the order of 8 MeV in the atomic nuclei of biologically relevant elements, which rises rapidly to a maximum at around 20 MeV, then asymptotically declines to about half the maximum value by about 100 MeV (Figure 8). Tabulated and graphical nuclear data may be obtained conveniently online from the Evaluated Nuclear Data File (ENDF) (http://wwwnds.iaea.org/exfor/endf.htm). ICRU Report 63 also provides extensive nuclear data for hadron therapy and radiation protection (ICRU, 2000).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0008.jpg

The total proton-induced non-elastic nuclear reaction versus proton energy, showing a threshold corresponding to the Coulomb barrier at approximately 6 MeV. (ENDF, 2013)

Several nuclear reactions are particularly important to clinical proton therapy and proton therapy research. In a proton therapy beam, proton-induced reactions can produce energetic protons, deuterons, tritons, He, He, and other ions. Secondary protons comprise as much as about 10% of the absorbed dose in a high-energy proton treatment beam; they have a small but non-negligible impact on the spatial dose distribution in a patient (Medin and Andreo, 1997; Boon, 1998; ICRU, 1998; Paganetti, 2002; Wroe et al, 2005) Deuterons and heavier ions are present in much smaller proportions; collectively they comprise about 1% or less of the therapeutic absorbed dose (ICRU, 1998). Their energy and range are very small and they deposit their kinetic energy locally, i.e., very near their point of creation.

Relatively high-current protons beams are incident on certain beam production and delivery equipment and on some patients. These proton beams produce neutrons that create significant potential safety hazards. Great care must be taken to limit exposures of personnel to neutrons (NCRP, 1971, 1990, 1993, 2005; Newhauser and Durante, 2011). Some electronic systems must also be hardened, shielded, or located so that neutron radiation does not cause soft upsets or permanent damage to semiconductor components. Attention must also be paid to neutron activation of beamline components, air, groundwater, and other materials (IAEA, 1988).

Neutrons are produced in copious quantities: they span 10 orders of magnitude in neutron energy, their energy distributions depend strongly on the proton beam energy and direction, they are extremely penetrating, and their relative biologic effectiveness is as much as about 20 times higher than that of proton radiation (ICRP, 2007). Thus they potentially increase the risk of radiogenic late effects (Hall, 2006; Brenner and Hall, 2008; Newhauser and Durante, 2011). Great care must be taken to ensure that patients and staff are adequately protected from neutron exposures. Several specific aspects of neutron exposures are considered in a later section of this paper.

Nuclear reactions inside the patient may provide a non-invasive approach to measure a variety of beam and patient properties, such as proton beam range, elemental composition of tissues, and even intra- or inter-fraction physiology. The basic approach is to detect gamma rays from proton-induced nuclear reactions, such as neutron capture reactions, denoted by (n,γ). Approaches are under development that detect photons from positron annihilation, prompt gammas, and delayed gammas. Gamma ray detection approaches have included positron emission tomography camera (Parodi et al, 2007; Moteabbed et al, 2011; Cho et al, 2013; Min et al, 2013), Compton camera (Peterson et al, 2010; Smeets et al, 2012), one-dimensional detector arrays (Min et al, 2012), and photon counting systems (Kim et al, 2012). These techniques are in various stages of research and development; none is routinely used in clinical practice. There remain many challenges to overcome, including instrument sensitivity and calibration; interpretation of measurements, including an understanding of managing measurement artifacts; and competition from alternative methods, e.g., magnetic resonance imaging techniques (Krejcarek et al, 2007; Raaymakers et al, 2008; Gensheimer et al, 2010).

2.1 Energy loss rate

The energy loss rate of ions, or linear stopping power, is defined as the quotient of dE and dx, where E is the energy and x is the distance. It is frequently more convenient to express the energy loss rate in a way that is independent of the mass density; the mass stopping power is defined as

Sρ=dEρdx,
(1)

where ρ is the mass density of the absorbing material. Please note that stopping power is defined for a beam, not a particle.

The energy loss rate may be described by several mathematical formulae. The simplest, yet still remarkably accurate, formula is based on the Bragg-Kleeman (BK) Rule (Bragg and Kleeman, 1905), which was originally derived for alpha particles, and is given by

Sρ=dEρdxE1pραp,
(2)

where ρ is the mass density of the material, α is a material-dependent constant, E is the initial energy of the proton beam, and the exponent p is a constant that takes into account the dependence of the proton's energy or velocity. Values of α and p may be obtained by fitting to either ranges or stopping power data from measurements or theory.

A more physically complete theory, developed by Bohr (Bohr, 1915), is based on calculation of the momentum impulse of a stationary, unbound electron and the impact parameter. A more accurate formula, attributed to Bethe and Bloch (Bethe, 1930; Bloch, 1933), takes into account quantum mechanical effects and is given by

Sρ=dEρdx=4πNAre2mec2ZAz2β2[ln2mec2γ2β2Iβ2δ2CZ],
(3)

where NA is Avogadro's number, re is the classical electron radius, me is the mass of an electron, z is the charge of the projectile, Z is the atomic number of the absorbing material, A is the atomic weight of the absorbing material, c is speed of light, β = v/c where v is the velocity of the projectile, γ = (1- β), I is the mean excitation energy of the absorbing material, δ is the density correction item arising from the shielding of remote electrons by close electrons and will result in a reduction of energy loss for higher energies, and C is the shell correction item, which is important only for low energies where the particle velocity is near the velocity of the atomic electrons. The two correction items in the Bethe-Bloch equation involve relativistic theory and quantum mechanics and need to be considered when very high or very low proton energies are used in calculations. Figure 2 plots proton stopping power as a function of proton energy in water calculated by using Equation (3).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0002.jpg

Mass stopping power (S) versus ion energy (E) for protons in liquid water. The corresponding range (R), calculated using the plotted S values and on the assumption of the continuous slowing down approximation (CSDA), is also plotted.

It is instructive to observe in Equation (3) how the projectile's characteristics govern its energy loss rate: energy loss is proportional to the inverse square of its velocity (1/v classically and 1/β relativistically) and the square of the ion charge (z = 1 for protons), and there is no dependence on projectile mass. Similarly, Equation (3) reveals that the absorber material can also strongly influence the energy loss rate. Specifically, the linear stopping power is directly proportional to the mass density. It is equivalent, but perhaps more physically meaningful, to state that the linear stopping power is proportional to the density of electrons in the absorber (NAρ Z/A), because the energy loss occurs by Coulombic interactions between the proton and atomic electrons. Z/A varies by only about 16%, from 0.5 for biologic elements such as carbon and oxygen to 0.42 for high-Z beamline components, such as lead. Hydrogen is an obvious exception to this; fortuitously, the concentration of hydrogen in the human body is low (only about 10%) and nearly uniform throughout the body. The stopping power also depends on a material's I value, and the I value depends in a monotonic way on the Z of the absorber, varying from about 19 eV for hydrogen to about 820 eV for lead. However, the stopping power goes with the logarithm of I value, so the dependence is diminished. Hence, putting these various dependencies in perspective, it is clear that the proton energy loss rate in the human body depends most strongly on the material density, which can vary by about three orders of magnitude, from air in the lung to cortical bone, and the ion velocity, which can cause the linear stopping power in water to vary by about a factor of 60 for proton energies between 1 MeV and 250 MeV.

The stopping power theory for protons and heavier ions was reviewed by Ziegler et al. (Ziegler et al, 1985; Ziegler, 1999; Ziegler et al, 2008) and in Report 49 of the International Commission on Radiation Units and Measurement (1993). Evaluated stopping power and range tables may be conveniently calculated with the SRIM code (“Stopping and Ranges of Ions in Matter,” computer program, http://www.srim.org/.).

Thus far, we have described proton energy loss in an approximate way on the assumptions that a proton loses energy along a 2-dimensional line trajectory and that energy loss is continuous. Absorption of this same energy, however, occurs in a 3-dimensional volume. Furthermore, the ionization track of a proton has an irregular 3-dimensional structure caused by random fluctuations in the location and size of primary ionization events. This is caused mainly by proton-produced recoil electrons, some of which are sufficiently energetic to create small spur tracks of ionization emanating from the main track. Because the electrons are very much lighter than the protons, each interaction can reduce the proton energy only a little. The maximum possible energy transfer in a collision of an ion of mass m with an unbound stationary electron is

Δ1max=2mec2β21β2[1+2meM11β2+(meM)2],
(4)

where me is the mass of an electron, M is the mass of the target material, c is the speed of light, and β = v/c where v is the velocity of the projectile.

Even for very energetic protons, the secondary electrons do not acquire enough energy to travel more than a few millimeters from the proton track. For example, at 200 MeV proton energy, the maximum secondary electron energy is around 500 keV, which corresponds to an electron range of approximately 2 mm in water. The probability of producing secondary electrons may be calculated with various total or differential cross-sections; these were reviewed in ICRU Report 55 (1995). Track structure models may be used to estimate the radial properties of ions (Kraft et al, 1999), although this has not found common application in clinical proton therapy. Regardless of the calculation method used, the spatial characteristics of secondary charged particles should, in principle, be taken into account near material interfaces (e.g., buildup effects in transmission beam monitoring instruments, skin, air-tumor interfaces in the lung) and in cases where the radiation quality is of interest (e.g., microdosimetric and nanodosimetric characterization of individual dose deposition events).

Finally, we note that in proton therapy water is considered an excellent tissue substitute because of its similar density, effective Z/A, and other properties. Furthermore, proton energy loss and residual range in various materials are often expressed in terms of their water-equivalent values. We discuss this in detail in section 3.1, but until then, for simplicity, we consider water as being perfectly equivalent to tissue.

2.2 Range

Range is defined as the depth at which half of protons in the medium have come to rest, as shown in the range-number curve plotted in Figure 3. There are small variations in the energy loss of individual protons (an effect called range straggling, which is discussed in the following section). Consequently, the range is inherently an average quantity, defined for a beam and not for individual particles. By convention, this means half of the protons incident on the absorber are stopped, although in some cases this is taken instead to mean that half the protons survive to near the end of range (i.e., neglecting protons removed by nuclear reactions).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0003.jpg

Relative fraction of the fluence Φ in a broad beam of protons remaining as a function of depth z in water. The gradual depletion of protons from entrance to near the end of range is caused by removal of protons from nuclear reactions. The rapid falloff in the number of protons near the end of range is caused by ions running out of energy and being absorbed by the medium. The sigmoid shape of the distal falloff is caused by range straggling or by stochastic fluctuations in the energy loss of individual protons.

The path of most protons in matter is a nearly straight line. On average, the proton's pathlength is very nearly equal to its projected pathlength and range. This simple but important fact renders many proton range calculations tractable with relatively simple numerical or analytical approaches.

Let us first consider a simple numerical calculation of proton beam range. We use proton stopping power data and perform a one-dimensional proton pathlength transport calculation on the assumptions that the ions travel only straight ahead (negligible lateral scattering) and that the protons lose energy in a continuous matter. (These assumptions are reasonable for many clinical calculations, but we examine then relax these assumptions in a later section.) In this case, the range (R) may be calculated as

R(E)=0E(dEdx)1dE0E(dEdx)1dE,
(5)

where E is the ion's initial kinetic energy. The summation denotes that the continuous transport is approximated by calculations of discrete steps. In fact, as discussed above, this equation actually gives the pathlength, which is an excellent approximation of range in most clinical situations. Figure 2 plots proton range in water calculated by using Equation (5).

Next, we calculate the proton range using an analytical approach, which may be faster and more practical than the numerical approach for many clinical calculations. We begin by noting that the interval of proton range of interest typically extends from 1 mm (about the size of a voxel in an anatomic image of patient anatomy) to about 30 cm (about the midline of a large adult male's pelvis, the deepest site in the human body). These ranges correspond to 11 MeV and 220 MeV, respectively, as seen in Figure 2. More importantly, Figure 2 reveals that the relationship between the logarithm of range and logarithm of energy is almost linear. This is fortuitous because this means that the range follows a very simple power law, as realized by Bragg and Kleeman (1905) and others early in the last century. Thus, the range of a proton may be calculated using BK Rule, or

R(E) = αE,
(6)

where, as before, α is a material-dependent constant, E is the initial energy of the proton beam, and the exponent p takes into account the dependence of the proton's energy or velocity.

The uncertainty in proton range depends on many factors. For example, the uncertainty in a range measurement depends on the precision and accuracy of the measurement apparatus and, in some cases, on the experimenter's skill. A common concern in clinical proton therapy is the uncertainty in calculated range, e.g., in calculating the settings of the treatment machine for a patient's treatment. The uncertainty in range may depend on the knowledge of proton beam energy distribution and on properties of all range absorbing materials in the beam's path. These include elemental composition, mass density, and linear stopping power. The linear stopping powers deduced from computed tomography scans have many additional sources of uncertainty, including imperfections in the calibration of CT scanners (i.e., the method used to convert image data from Hounsfield Units to linear stopping power), partial volume effects, motion artifacts, and streaks artifacts.

2.3 Energy and range straggling

In the preceding sections, we approximated the energy loss rate by assuming that the slowing of ions occurs in a smooth and continuous manner. In fact, we considered the mean energy loss rate and neglected variations in the energy loss rates of individual protons. For many clinical calculations, these assumptions are valid and lead to a reasonably good first-order approximation. However, the accumulation of many small variations in energy loss, termed energy straggling or range straggling, is one of the physical processes that strongly governs the shape of a proton Bragg curve, a subject that is covered in section 3.2. Thus, an understanding of range straggling is key to understanding the characteristics of proton dose distributions.

Figure 4 plots the relative energy loss probability density functions (PDF) for protons transmitted through water absorbers of various thickness. The curves have been normalized to enhance visual clarity. Apparently, thick absorbers result in a symmetric distribution of energy losses, whereas thin absorbers yield curves that are asymmetric with modes less than the mean and long tails of large energy losses. In principle, straggling PDFs may be calculated numerically from first principles, but usually theoretical approaches are used. Later in the section, three of the most commonly used straggling theories are described.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0004.jpg

Energy loss probability density functions (PDF) are plotted for various thicknesses of water absorbers, where the thickness is expressed in units of mean free path (mfp). For visual clarity, the energy-loss PDFs have been scaled on both the abscissa and ordinate. The single event energy loss is expressed as a fraction of the mean energy lost in the entire absorber thickness, or (Δ-Δav)/Δav. Each PDF was scaled so that the integral over all energy-loss values yields unit value. For thin absorbers (curves a-e), the PDFs are broader and asymmetric and are modeled with the Vavilov (1957) or Landau (1944) theories. For thick absorbers (curve f), the PDFs are symmetric and well-approximated with Bohr's theory (1915), i.e., a Gaussian distribution. (ICRU, 1993)

Having examined the mean energy loss rate in modest detail (section 2.1), and having conceptually introduced energy straggling, it is instructive to understand how these effects are related mathematically before delving into straggling theory. To understand these relationships, let us consider the moments of the ion energy PDF, or

Mn=0Δ1maxΔnf(Δ)dΔ,
(7)

where Δ is the energy loss of an ion in traversing an absorber, f(Δ)dΔ is the probability of energy loss occurring in the interval from Δ to Δ + dΔ, and n is the order of the moment. The zeroth moment is the total collision cross section, the first moment is the (mean) electronic stopping power, the second moment corresponds to the width (variance) of the energy straggling distribution, the third moment to its asymmetry (skewness), and the fourth to its kurtosis. The variance, sometimes denoted by σΔ, or second moment of the straggling distribution f(Δ), is

M2(Δ)=σΔ2=v0Δ1maxΔnf1(Δ)dΔ.
(8)

A more detailed discussion of the straggling moments was given by Rossi and Zaider (1996).

Next we examine theories for calculating energy straggling proposed by Bohr (1915), Landau (1944), and Vavilov (1957). These theories are valid for thick, intermediate, and thin absorbers, respectively. The criterion for selecting a valid theory for a given absorber thickness is based on a reduced energy parameter,

k=ξΔ1max,
(9)

where ξ is the approximate mean energy loss and Δ1 is the maximum energy loss possible in a single event.

According to Bohr's theory, the energy straggling distribution behaves according to a Gaussian PDF, or

f(Δ)dΔ=1σΔ2πexp(ΔΔ)22σΔ2,
(10)

where for non-relativistic ions the variance is given by

σΔ2=2πre2mec2z2NZβ2Δ1maxρx.
(11)

Bohr's theory assumes that the absorber is thick enough that there are many individual collisions (i.e., the central limit theorem holds), that the ion velocity does not decrease much in crossing the absorber, and that the absorber is made of unbound electrons. For most applications in proton therapy, Bohr's theory provides adequate accuracy. Several authors have reported convenient power law approximations to estimate sigma as a function of the proton beam range (Chu et al, 1993), or

σΔkR0m,
(12)

where R0 is the range in water in centimeters for a mono-energetic proton beam, k is a material-independent constant of proportionality, and the exponent is empirically determined. For protons in water, k = 0.012 and m =0.935 (Bortfeld, 1997).

Landau's theory relaxed the assumption that the central limit holds, i.e., there are relatively fewer individual collisions in an intermediate thickness absorber, and used an approximate expression for Δ1.

In this case, we have

f(Δ,ρx)dΔ=1ξϕL(λL),
(13)

where the parameter φL (λL) roughly corresponds to the deviation from the mean energy loss and was defined by Landau as

ϕL(λL)=1π0ey(lnyλL)sin(πy)dy.
(14)

Evaluation of the integral in Equation (14) is straightforward and the reader is referred to the literature for details (Seltzer and Berger, 1964).

Vavilov's theory is in essence a generalization of Landau's theory that utilizes the correct expression for Δ1 and is given by

f(Δ,ρx)dΔ=1ξϕV(λV,k,β2)dλV,
(15)

where

ϕV(λV,k,β2)=1πek(1+β2γ)0ekf1cos(λVy+kf2)dy,
(16)

and

λV=ΔΔΔ1maxk(1+β2γ),
(17)

where γ is Euler's constant. The evaluation of Vavilov's theory is computationally more expensive than that of Bohr's or Landau's. For additional details, the reader is referred to Vavilov's original work (Vavilov, 1957).

2.4 Multiple Coulomb scattering

As illustrated in Figure 5, a proton passing close to the nucleus will be elastically scattered or deflected by the repulsive force from the positive charge of the nucleus. While the proton loses a negligible amount of energy in this type of scattering, even a small change in its trajectory can be of paramount importance. In fact, it is necessary to take into account Coulomb scattering when designing beamlines and treatment heads (Gottschalk, 2010b) and in calculations of dose distributions in phantoms or patients, e.g., with treatment planning systems (Hong et al, 1996; Pedroni et al, 2005; Schaffner, 2008; Koch and Newhauser, 2010).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0005.jpg

Schematic diagram showing the trajectory of a proton undergoing multiple Coulomb scattering events. θ is the root mean square (rms) space scattering angle and θx is the projected rms scattering angle. (Leo, 1994)

To characterize the amount that a beam is deflected by scattering, we use the quantity of scattering power, which is defined as

T = d〈θ〉 ∕ dx,
(18)

where <θ> is the mean squared scattering angle and x is the thickness of absorber through which the proton has traveled. Similarly, the mass scattering power is simply T/ρ, where ρ is the mass density of the absorber material. Notice that the definition of scattering power utilizes the mean square of the scattering angle; scattering is symmetric about the central axis, and therefore the mean scattering angle is zero and contains no useful information. Also, note that the value of the scattering power depends on the material thickness (i.e., it depends on the material properties and dimensions of the absorber being considered).

Predictions of elastic Coulomb scattering are commonly classified according to the number of individual scattering events (Ns) that occur in a given absorber. For single scattering (Ns=1), Rutherford scattering theory applies. Plural scattering (1 < Ns < 20) is difficult to model theoretically and is not discussed further here. Multiple Coulomb scattering (MCS; Ns ≥ 20), the combined effect of all Ns scattering events, may be modeled by using a statistical approach to predict the probability for a proton to scatter by a net angle of deflection, commonly denoted by θ (Figure 5).

It is instructive to briefly examine Rutherford's theory of single scattering (Rutherford, 1911). The differential cross section dσ for scattering into the solid angle dΩ is given by

dΩdΩ=z12z22re2(mecβp)24sin4(θ2),
(19)

where z1 is the charge of the projectile, z2 is the atomic number of the absorber material, re is the classical electron radius, me is the mass of the electron, c is the speed of light, β = v/c, p is the ion momentum, and θ is the scattering angle of the proton. The angular dependence is governed by the term 1/sin(θ/2), i.e., in individual scattering events, the proton is preferentially scattered in the forward direction, at very small angles.

In clinical proton therapy, most objects of interest are thick enough to produce a great many scattering events. Thus, for clinical proton therapy, we are usually more interested in the net effect that many small-angle scattering events have on many protons, e.g., how beamline scattering in the treatment head influences the spatial distribution of dose in a patient.

Rigorous theoretical calculations of MCS are quite complex. The most complete theory was proposed by Molière (1948). Assuming that scattered particles are emitted at small deflection angles (i.e., the small-angle approximation in which sin(θ) ≈θ),

P(θ)=ηdη(2exp(η2)+F1(η)B+F2(η)B2),
(20)

where η = θ / (θ1B), and θ1=0.3965(zZpβ)ρδxA. The functions Fk(η) are defined as

Fk(η)=1k!J0(ηy)exp(y24)[y24ln(y24)]kydy,
(21)

where Jo(ηγ) is a Bessel function. Values of Fi have been tabulated in the literature by several authors. Parameter B in Equation (20) may be found by numerical methods, by solving

g(B) = ln BB + ln γ − 0.154 = 0,
(22)

where γ = 8831 q zρδx/(βAΔ) and Δ= 1.13+3.76 (Zz/137β). The following are symbols representing properties of the absorber: Z is the atomic number, A is the atomic mass, δx is the thickness, and ρ is the mass density. The proton momentum is denoted by p, β = v/c, and z = 1 is the proton charge. Even though we have not presented all of the details, clearly the evaluation of Molière's theory is indeed complex. Consequently, considerable attention has been paid in the literature to developing simpler formulae; the simplicity usually comes at the cost of reduced accuracy in modeling scattering at large angles. Gottschalk (1993) discussed several of these methods in detail. One approximation method takes the form of a Gaussian distribution,

P(θ)2θθ2exp(θζ2θ2)dθ,
(23)

where (<θ>) is the root mean square (rms) scattering angle or the width parameter of the Gaussian distribution.

Gottschalk (2010b) proposed a differential approximation to Molière's theory to predict the scattering power according to

TdM=fdM(pv,p1v1)×(ESpv)21XS,
(24)

where ES =15 MeV, p is the momentum and v is the velocity of the proton, p1 and v1 are the initial momentum and velocity, XS is scattering length of the material, and fdM is a material-independent nonlocal correction factor given by

fdM = 0.5244 + 0.1975 lg(1 − (pvp1v1)) + 0.2320 lg(pv) − 0.0098 lg(pv)lg(1 − (pvp1v1))
(25)

The factor fdM takes into account the accumulation of scattering that occurs as the proton slows from v1 to v and is material independent. The scatter length is given by

1ρXS=αNre2Z2A(2log(33219(AZ)131),
(26)

where α is the fine structure constant, N is Avogadro's number, re is electron radius, and Z, A, and ρare the atomic number, atomic weight, and density of the target material.

It is sometimes convenient to project the expected scattering angle < θ > in three-dimensional space to the corresponding value in a plane, denoted < θ x>, according to

θx = 〈θ ∕ 2,
(27)

Also, in proton therapy the lateral displacement (r) of a proton beam is typically calculated from the scattering angle. Under the Gaussian approximation, we have

P(r)dr = 6r(〈θt) exp(−3r ∕ (〈θt)dr,
(28)

where t = x/Lrad, and Lrad is the radiation length of the material, which can be calculated easily or looked up in tables. The mean square lateral displacement is given by

r〉 = 〈θt ∕ 3,
(29)

A power law approximation for rms of protons is

rrms = aR,
(30)

where a is a unitless material constant, R is the proton beam range (in the same units as rrms), and the exponent b governs the range dependence. For protons in water, a = 0.0294 and b = 0.896 (Chu et al, 1993). Koehler and Preston derived a convenient expression to calculate rrms as a function of depth in an absorber and knowledge of its maximum values, rmax, at the end of range (from an unpublished manuscript; some portions of their work were reported by Gottschalk (2010b)).

In proton therapy, MCS in the treatment head (i.e., in the scattering foils) is helpful because it allows the beam to be spread laterally to useful dimensions, e.g., to make a beam laterally large and flat so that a tumor may be completely covered with a uniform dose. Scattering foils in the nozzle are carefully designed to utilize MCS and energy loss to produce clinically useful beams (Koehler et al, 1977; Gottschalk, 2004) However, MCS in the treatment head and in the patient blurs lateral penumbral sharpness. This is manifested as penumbral growth at the edge of collimated beams and/or the growth of the lateral spot size of a scanned beam (Figure 6). Understanding and preserving penumbral sharpness is the key to realizing the full benefit of proton therapy for sparing healthy tissue.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0006.jpg

Broadening of the beam width in water due to multiple Coulomb scattering. (Pedroni et al, 2005)

Recent studies have revealed that MCS plays an important role in proton dose distribution around small implanted metal objects. Specifically, implanted fiducial markers for image-guided patient alignment have been used in proton therapy for many years (Gall et al, 1993; Welsh et al, 2004; Newhauser et al, 2007a; Newhauser et al, 2007c; Ptaszkiewicz et al, 2010; Matsuura et al, 2012). Substantial recent improvements in on-board imaging systems, patient positioning, and patient immobilization have led to increased use of radiopaque implanted fiducial markers in proton therapy to many disease sites, with the goal of improving target coverage and/or normal tissue sparing. However, recent studies revealed that some commonly used markers, even those less than 1 mm in size, can cause severe cold spots (Figure 7), compromising target coverage (Newhauser et al, 2007a; Newhauser et al, 2007c; Carnicer et al, 2013). The severity of the cold spots varies with fiducial size, material composition, and mass density. These parameters, in turn, determine the amount of MCS and energy loss in the fiducial and, hence, perturbations to the dose distribution in surrounding tissue. In essence, MCS in the fiducial causes a downstream dose shadow that may be partially or fully filled in by MCS in the surrounding tissue. While MCS is important, energy loss in the fiducial (or its water-equivalent thickness), the proximity of the fiducial to the end of the proton beam range, and of course its size and orientation with respect to the beam also should be considered. The physics of dose perturbations is explained in detail elsewhere (Newhauser et al., (2007a), and several subsequent studies (Giebeler et al, 2009; Lim et al, 2009; Cheung et al, 2010; Huang et al, 2011) have shown that it is possible to achieve good radiographic visibility using novel markers that do not significantly perturb the therapeutic dose distribution in tissue.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0007.jpg

Small implanted fiducial markers can create clinically significant dosimetric cold spots in proton therapy beams. (Upper) two-dimensional dose distribution as a function of depth in water (z) and crossfield position (x) from a Monte Carlo simulation of range-modulated proton beam incident on a water phantom containing tantalum localization markers oriented (a, b) parallel to the beam axis and (c, d) perpendicular to the beam axis. The range and modulation width are typical for uveal melanoma treatments. (Lower) simulated absorbed dose (D) as a function of depth (z) in the water phantom at various off-axis positions. The perturbed depth dose profiles are parallel to the beam axis and pass through the center of markers a–d in the plot above. For visual clarity, portions of the perturbed dose profiles upstream of the markers are not shown. An unperturbed beam is plotted with open squares. (Newhauser et al, 2007c)

Others have examined the effects of MCS in larger metal objects on clinical proton beams and characterized the suitability of approximate methods to predict MCS in practical clinical applications (Herault et al, 2005; Stankovskiy et al, 2009; Newhauser et al, 2013).

2.5 Nuclear interactions

In addition to the mechanisms already described, protons may interact with the atomic nucleus via non-elastic nuclear reactions in which the nucleus is irreversibly transformed, e.g., a reaction in which a proton is absorbed by the nucleus and a neutron is ejected, denoted by (p,n). The main effect of nuclear reactions within a therapeutic region of a proton field is a small decrease in absorbed dose due to the removal of primary protons, which is compensated to a large extent by the liberation of secondary protons and other ions. In this section, we discuss this and several other important aspects of nuclear reactions.

Before discussing reaction mechanisms, it is instructive to examine a range-number curve (Figure 3), which plots the remaining number of protons versus depth in an absorber as a beam comes to rest. The gradual depletion of protons from entrance to near the end of range is caused by removal of protons from nuclear reactions. The rapid falloff in the number of protons near the end of range is caused by ions running out of energy and being absorbed by the medium. The sigmoid shape of the distal falloff is caused by range straggling or by stochastic fluctuations in the energy loss of individual protons.

To enter the nucleus, protons need to have sufficient energy to overcome the Coulomb barrier of the nucleus, which depends on its atomic number. The total non-elastic cross-section for proton-induced nuclear reactions has a threshold, on the order of 8 MeV in the atomic nuclei of biologically relevant elements, which rises rapidly to a maximum at around 20 MeV, then asymptotically declines to about half the maximum value by about 100 MeV (Figure 8). Tabulated and graphical nuclear data may be obtained conveniently online from the Evaluated Nuclear Data File (ENDF) (http://wwwnds.iaea.org/exfor/endf.htm). ICRU Report 63 also provides extensive nuclear data for hadron therapy and radiation protection (ICRU, 2000).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0008.jpg

The total proton-induced non-elastic nuclear reaction versus proton energy, showing a threshold corresponding to the Coulomb barrier at approximately 6 MeV. (ENDF, 2013)

Several nuclear reactions are particularly important to clinical proton therapy and proton therapy research. In a proton therapy beam, proton-induced reactions can produce energetic protons, deuterons, tritons, He, He, and other ions. Secondary protons comprise as much as about 10% of the absorbed dose in a high-energy proton treatment beam; they have a small but non-negligible impact on the spatial dose distribution in a patient (Medin and Andreo, 1997; Boon, 1998; ICRU, 1998; Paganetti, 2002; Wroe et al, 2005) Deuterons and heavier ions are present in much smaller proportions; collectively they comprise about 1% or less of the therapeutic absorbed dose (ICRU, 1998). Their energy and range are very small and they deposit their kinetic energy locally, i.e., very near their point of creation.

Relatively high-current protons beams are incident on certain beam production and delivery equipment and on some patients. These proton beams produce neutrons that create significant potential safety hazards. Great care must be taken to limit exposures of personnel to neutrons (NCRP, 1971, 1990, 1993, 2005; Newhauser and Durante, 2011). Some electronic systems must also be hardened, shielded, or located so that neutron radiation does not cause soft upsets or permanent damage to semiconductor components. Attention must also be paid to neutron activation of beamline components, air, groundwater, and other materials (IAEA, 1988).

Neutrons are produced in copious quantities: they span 10 orders of magnitude in neutron energy, their energy distributions depend strongly on the proton beam energy and direction, they are extremely penetrating, and their relative biologic effectiveness is as much as about 20 times higher than that of proton radiation (ICRP, 2007). Thus they potentially increase the risk of radiogenic late effects (Hall, 2006; Brenner and Hall, 2008; Newhauser and Durante, 2011). Great care must be taken to ensure that patients and staff are adequately protected from neutron exposures. Several specific aspects of neutron exposures are considered in a later section of this paper.

Nuclear reactions inside the patient may provide a non-invasive approach to measure a variety of beam and patient properties, such as proton beam range, elemental composition of tissues, and even intra- or inter-fraction physiology. The basic approach is to detect gamma rays from proton-induced nuclear reactions, such as neutron capture reactions, denoted by (n,γ). Approaches are under development that detect photons from positron annihilation, prompt gammas, and delayed gammas. Gamma ray detection approaches have included positron emission tomography camera (Parodi et al, 2007; Moteabbed et al, 2011; Cho et al, 2013; Min et al, 2013), Compton camera (Peterson et al, 2010; Smeets et al, 2012), one-dimensional detector arrays (Min et al, 2012), and photon counting systems (Kim et al, 2012). These techniques are in various stages of research and development; none is routinely used in clinical practice. There remain many challenges to overcome, including instrument sensitivity and calibration; interpretation of measurements, including an understanding of managing measurement artifacts; and competition from alternative methods, e.g., magnetic resonance imaging techniques (Krejcarek et al, 2007; Raaymakers et al, 2008; Gensheimer et al, 2010).

3. Proton transport calculations

In this section, we review several aspects of proton transport physics that are encountered frequently in clinical and research situations. We describe the one-dimensional water-equivalent thickness of an arbitrary material, the shapes of a pristine Bragg curve and a spread out Bragg peak (SOBP) curve, and stray radiation exposures.

3.1 Water-equivalent thickness

As we mentioned previously, in proton therapy water closely mimics the properties of human tissues in terms of energy loss, MCS, and nuclear interactions. As such, water is a recommended phantom material for dose and range measurements and reference material for reporting corresponding calculated quantities (ICRU, 1998; IAEA, 2000). For example, it is a common and convenient clinical practice to specify the penetrating power of a proton beam by its range in water (ICRU, 1998, 2007). In this way, range losses in various beamline objects and the patient may be easily added or subtracted from one another in a physically consistent way. Viewed another way, it is also convenient to specify the range-absorbing power of various objects in the beam path, e.g., beam transmission monitors and immobilization devices, by their equivalent thickness if they were made of water.

Water-equivalent thickness (WET) is often used to characterize the beam penetration range; Figure 9 schematically illustrates the concept of WET and how it can be calculated or measured. For treatment sites with nearby critical structures, e.g., an optic nerve, the range of the planned and delivered beams must agree within a few millimeters. To accomplish this, treatment planning systems are commonly configured with the WET values of all items not included in the planning CT images, such as components in the treatment head, immobilization devices not present during the CT scan, or a treatment couch (Newhauser et al, 2007b). Similarly, to determine the measurement geometry for patient-specific clinical quality assurance measurements, the WET of measurement instruments and possibly other devices must be determined (Newhauser, 2001b, a). Thus, it is important to have methods to calculate and measure WET. In this section, we emphasize recently developed calculation methods that are convenient and suitable for clinical calculations, using the energy loss theories presented in section 2. Our discussion of WET measurement methods is very brief, mainly because they are relatively simple and obvious. In practice, however, WET measurements remain very important (Andreo, 2009; Gottschalk, 2010a; Newhauser and Zhang, 2010; Zhang et al, 2010b; Besemer et al, 2013).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0009.jpg

Schematic illustration of the concept of water equivalent thickness (WET) and how it can be calculated or measured by observing changes in the depth of a proton Bragg peak in a water tank. (Zhang and Newhauser, 2009)

The IAEA (2000) proposed that WET can be approximated by

tw=tmρmρwcm,
(31)

where the depth scaling factor, cm, can be calculated, to a good approximation, as the ratio of the continuous slowing down approximation (CSDA) range (in g cm) in water to that in the target:

cm=ρwRwρmRm.
(32)

Because the ranges in Equation (32) correspond to complete loss of ion energy, this approach is strictly valid only for stopping-length targets. An exact equation for WET that is applicable to thin targets was reported by Zhang and Newhauser (2009):

tw=tmρmSm¯ρwSw¯,
(33)

where ρw and ρm are the mass densities of water and material, respectively, and S̄m and S̄w are the mean proton mass stopping power values for the material and water, respectively, defined by

S=ESdEEdE.
(34)

For thin targets, where the proton loses a negligible fraction of its energy in the absorber material, we have

twtm(ρmρw)(SmSw)=tmαwpwαmpmEpwpm,
(35)

where the reader will recognize α and p from the discussion of the BK Rule (see section 2.1). Values of α and p for commonly encountered materials in proton therapy are provided in Table 2. Zhang and Newhauser (2009) also reported a slightly more complex analytical formula to calculate WET for targets of arbitrary thickness.

Table 2

Common materials (lung substitute plastic, high-density polyethylene (HDPE), water, polystyrene, polymethylmethacrylate (PMMA), polycarbonate resin (Lexan), bone substitute plastic, aluminum, titanium, stainless steel, lead and gold) used in heavy charged particle beams, with their mass densities ρ, values of (Z/A)eff, mole fractions and fitting parameters of α and p for these materials when applying the BK rule (The energies used in the fit were from 10 to 250 MeV).

Materialρ (g/cm)Mole fraction (%)αp
Lung substitute0.30.537H 55.577, C 32.738, N 0.927, O 7.508, Cl 0.019, Si 0.184, Mg 3.0488.994×10−31.735
HDPE0.9640.570H 66.717, C 33.2832.541×10−31.737
Water1.00.555H 66.667, O 33.3332.633×10−31.735
Polystyrene1.060.538H 49.851, C 50.1492.545×10−31.735
PMMA1.1850.539H 53.333, C 33.333, O 13.333,2.271×10−31.735
Lexan1.200.527H 42.424, C 48.485, O 9.0912.310×10−31.735
Bone substitute1.8290.516H 35.215, C 29.592, N 0.803, O 26.695, Cl 0.16, Ca 7.6791.666×10−31.730
Aluminum2.6980.482Al 1001.364×10−31.719
Titanium4.5190.459Ti 1009.430×10−41.710
Stainless steel7.850.466C 0.045, N 0.045, Si 0.450, Cr 18.150, Mn 1.250, Fe 71.460, Ni 8.550, Mo 0.0505.659×10−41.706
Lead11.3220.396Pb 1006.505×10−41.676
Gold19.3110.401Au 1003.705×10−41.677

As can be seen in the curves of the water-equivalent ratio (WER = tw/tm) plotted in Figure 10, taking into account the target thickness in calculating WER is most important for absorbers that are thick and made of high Z materials, e.g., lead scattering foils, and for protons that are of comparatively low energy when impinging on the target. For low-Z materials, such as tissue and plastic, WER depends only very weakly on the target thickness and initial proton beam energy, and the approximate (thin and stopping length) analytical methods provide sufficient accuracy for most clinical applications.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0010.jpg

Calculated water-equivalent ratio (WER) values as a function of proton beam energy. This plot illustrates the dependence of WER value on the target material, the beam energy, and the target thickness. (Zhang and Newhauser, 2009)

3.2 Features of pristine and spread-out Bragg curves

The spatial dose distribution from clinical proton therapy beams is quite similar to those from photon and electron beams. The lateral profiles are generally quite flat in the central high-dose region, then fall off rapidly in the penumbral regions, where the penumbra width increases with depth in the patient. On the other hand, the central-axis depth-dose curve from protons is somewhat similar to that from electrons, but with a sharper distal falloff. Figures 11 and and1212 compare the central-axis depth-dose curves from several radiation therapy beams, revealing the main dosimetric properties that are clinically advantageous in many cases, namely, relatively low entrance dose, large and uniform dose to cover the tumor, and rapid falloff of dose near the end of range to spare normal tissues. These properties, together with a uniform lateral dose profile and a sharp lateral penumbral width, allow proton beams to treat a wide variety to tumor sizes and locations while providing superior sparing of normal tissues in many cases.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0011.jpg

Central axis depth dose profiles from several particle beams. Note that these distributions are from solitary beams in order to clearly compare the differences in the physical properties of various radiations. The important features are that proton beams offer relatively low entrance dose and virtually no exit dose. However, many clinical treatment techniques exploit multiple field directions to enhance the uniformity of tumor coverage and to spare sensitive healthy tissues. In fact, in some cases proton treatments provide inferior skin sparing to photons and/or inferior target coverage, e.g., because of proton beams’ sensitivity to range errors. Nonetheless, beam for beam, proton beams provide excellent tissue sparing, especially beyond the end of range. (Larsson, 1993)

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0012.jpg

Comparison of depth-dose curves from proton SOBPs and from electron beams. Because the proton mass is nearly 2000 times that of an electron, proton scattering interactions (individual angular deflections and variations in collisional energy losses) are much smaller, leading to sharper lateral and distal falloff distances. (Koehler and Preston, 1972)

Having casually inspected the shape of proton depth-dose curves, we next examine their structure in greater detail, pointing out nomenclature and the physical processes that govern the shape of various regions. Figure 13 shows a pristine proton peak along with labels identifying several regions. In order of increasing depth, these are the regions of electronic buildup, protonic buildup, sub-peak, peak, and distal falloff (same as for peak). The figure also shows several characteristic depths (e.g., the depth zBP at which the Bragg peak occurs) and various characteristic lengths (e.g., the 80%-to-20% distal-falloff length l80-20 and the proximal-80%-to-distal-80% pristine-peak width).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0013.jpg

Absorbed dose D as a function of depth z in water from an unmodulated (pristine) proton Bragg peak produced by a broad proton beam with an initial energy of 154 MeV. The various regions, depths, and lengths that are labeled are defined in the text. (The electronic buildup is not visible in this plot.) This type of dose distribution is useful clinically because of the relatively low doses delivered to normal tissues in the sub-peak and distal-falloff regions relative to the target dose delivered by the peak.

The anatomic definitions of an SOBP are, in many ways, similar to those of a pristine Bragg curve, as seen in Figure 14. However, there are several unique difficulties in characterizing SOBPs because of their sometimes unusual shape. For example, SOBPs with two or more discrete pristine Bragg curves may have multiple dose maxima in the modulated-peak region (e.g., the ripple shown in Figure 15). Because of such problems we introduce a few additional quantities that are defined only for SOBPs. To a large extent, however, we have defined quantities and terminology that are common to both modulated and pristine Bragg curves.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0014.jpg

Absorbed dose D as a function of depth z in water from a spread-out proton Bragg peak (SOBP). Various locations and regions that are indicated on the plot are defined in the text. This peak was measured with a Markus-type parallel-plate ionization chamber in the Northeast Proton Therapy Center (NPTC) gantry. The measured data are plotted with open circles and the model-fit as a solid line. Note that the electronic buildup region, which spans only a few millimeters, is not visible in this plot.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0015.jpg

Absorbed dose D as a function of depth z in water from a spread-out Bragg peak (SOBP) (uppermost curve) and its constituent pristine Bragg peaks (lower curves; for clarity, all but the deepest pristine Bragg peak are only partly drawn). In many cases, the clinical target volume is larger than the width of a pristine Bragg peak. By appropriately modulating the proton range and fluence of pristine peaks, the extent of the high-dose region can be widened to cover the target volume with a uniform dose.

We have not yet mentioned how MSC affects the shape of the depth-dose curves. In fact, near the central region of a laterally “large” beam, or more correctly well inside the periphery of a large beam, there is an equilibrium in which lateral scattering away from the central axis is exactly compensated by scattering towards it. This effect is described in Figure 16, which is adapted from Koehler and Preston (Koehler and Preston, 1972). As the field size shrinks to the dimension of the rms lateral displacement due to MCS, lateral equilibrium is lost and MCS progressively depletes the proton fluence and dose along the central axis. Small proton beams have been investigated in several studies, including those by Takada (Takada, 1996), Moyers et al. (1999), Vatnitsky et al. (1999b), Bednarz et al. (2010), and Gottschalk (1999), as well as others, especially in the context of scanned beams and pencil beam dose algorithms.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0016.jpg

(Left) Proton fluence I(0; x) along the beam central axis vs. depth x in water. Curves are shown for beams with circular cross sections and radii of 1 mm to 4 mm. Some of the protons are lost because of scattering events that deflect them from the central axis. This is increasingly observed for small beams and at large depths. (Right) The corresponding central-axis absorbed-dose curves. Note how the fluence depletion reduces the absorbed dose at the peak relative to the entrance dose. (Preston and Koehler, 1998)

Here, we use a Cartesian coordinate system with the z axis parallel to and centered about the proton beam central axis. The x and y axis are mutually orthogonal and perpendicular to the z axis. The coordinate system origin is located at the front face of the absorber, e.g., the extended medium in which we consider the absorbed dose distribution.

Pristine Bragg curve

A depth-dose distribution in an absorber irradiated with a monoenergetic or nearly-monoenergetic proton beam. In other words, no device or technique has been intentionally deployed for modulating the proton fluence or spectral fluence.

Spread-out Bragg curve

A depth-dose distribution in an absorber irradiated with a beam that has been intentionally modified to increase the axial dimension of the peak region. This is accomplished by modulating the range and the fluence of the beam. Clinical systems accomplish this by combining multiple quasi-monoenergetic beams or with a continuously modulated beam.

Electronic buildup region

A small region near the surface of the absorber where the proton beam is incident. As discussed in section 2.1, high-energy proton beams liberate delta rays with sufficient kinetic energy to travel several millimeters in tissue. Under some circumstances, this region exhibits an increase of dose with increasing depth, asymptotically approaching absorbed dose in the sub-peak region within the depth corresponding to the range of the most penetrating recoil electron. In some cases, electronic buildup is not observed. There are several possible reasons for this: the presence of some material just upstream of the surface (e.g., an immobilization device or a range compensator) may provide partial or full electronic charged particle equilibrium in the absorber, it may occur in combination with protonic buildup, it may be masked by changes in the proton energy loss rate near the end of range, or the wall of a cavity dosimeter may be sufficiently thick to present electronic equilibrium to the dosimeter's sensitive volume.

Protonic buildup region

A region near the surface of the absorber where the absorbed dose increases with depth because of the buildup of secondary protons that are attributable to proton-induced non-elastic nuclear interactions (e.g., O(p,xp) reactions). As with electronic buildup, the protonic buildup may not be observed in some cases, particularly at low incident proton beam energies.

Sub-peak region

The region extending from the surface of the absorber to the depth just proximal of the peak. The physical processes involved here are, in decreasing order of importance, the stopping power's dependence on the inverse-square of the proton velocity, the removal of some protons by nuclear reactions, the liberation of secondary particles from nuclear reactions, and for very small fields, the accumulation of lateral deflections from MCS leading to lateral protonic disequilibrium and reduction of the proton fluence on the central axis. The distal extent of the sub-peak region can be calculated from zm- 2σ, where zm is the depth at the pristine Bragg peak and σ is the width of the peak. The width parameter σ can be estimated from the incident proton beam spectral fluence and the range straggling accumulated in the absorber, as discussed in section 2.3.

Pristine Bragg peak

The pristine Bragg peak is simply the maximum (or mode) dose near the end of range, and is located at zBP, which is defined next. The physical processes governing the location and/or height of the peak are mainly the proton stopping power and energy straggling, nuclear reactions to a much lesser extent and, for very small fields, MCS..

Pristine Bragg peak depth

The depth near the end of range of the primary protons at which the protons produce the maximum dose rate, denoted by zBP. Although small proton beams are not yet widely used, it is helpful to define the location of zBP in a way that is compatible with large and small beams. Figure 16 shows that the maximum dose for beams of diameter larger than 6 mm is clearly single valued and located near the end of range. For smaller beams, however, the dose at the peak near the end of range may be less than the dose in the proximal regions, creating multiple maxima to choose from. Hence, the definition of zBP restricts it to exist in the region of the R±4σ, where sigma is the distal falloff width, thereby preventing possible ambiguities, and makes zBP conceptually independent of the beam cross-sectional area.

Distal falloff region

This region extends from depths greater than that of the pristine Bragg peak depth, zBP. The width of this region is not restricted. In many practical situations, however, the distal falloff region terminates at a depth where the dose falls below a threshold value, e.g., 1% of the dose at the Bragg peak, D(zBP).

Distal-50% depth

The distalmost depth, denoted by zd50, at which the absorbed dose is equal to half of the absorbed dose at the pristine Bragg peak depth, or D(zBP)/2. For an SOBP, it is defined as the distalmost depth at which the absorbed dose is equal to half of the absorbed dose at the SOBP dose in the modulated peak region. In many cases, the dose in the modulated peak region varies with depth (perhaps by design), making selection of the absorbed dose in the modulated peak region arbitrary. Since this clearly hinders an unambiguous definition of zd50, we instead define zd50 for an SOBP as being equal to zb. If the dose in the modulated peak region varies with depth, zb may have to be determined with numerical methods. Physically, the value of zb is closely related to the RCSDA value of the mean proton energy corresponding to the most energetic peak in the SOBP. Definitions for various other distal depths are similarly defined, e.g., the distal-90% depth zd90 , zd50, and zd20,.

Proximal-50% depth

The second most distal depth, denoted by zp50, at which the absorbed dose is equal to half of the absorbed dose at the pristine Bragg peak depth, or D(zBP)/2, provided that occurs within the absorber. (For very low-energy pristine Bragg curves, the entrance dose may be greater than half the peak dose). In such cases, and for SOBPs, it is defined as the most proximal depth at which the absorbed dose is equal to half of the absorbed dose in the modulated peak region. Because of problems that are analogous to those described in the zd50 definition above, zp50 is defined as the location at which the dose is equal to half the value at depth za + Δ/2, where Δ is the width of the transition from the peak region to the sub-peak region. Δ is not a critical parameter and it may be estimated as the amount of range straggling at za or it can be determined from measured Bragg curves. Conceptually, za typically corresponds to the shallowest location that is expected to receive the maximum dose. Physically, the value of za is closely related to the RCSDA value of the least energetic pristine peak in the SOBP curve (or the only peak in the Bragg curve for a pristine peak). In many cases, particularly for SOBPs with large modulated-peak regions, the absorbed dose throughout the sub-peak region exceeds 50% of the value at za + Δ, in which case proximal 50% depth is undefined. In such cases, proximal depths at higher dose percentages, such as zp95 and zp90, may be used. If the Bragg peak occurs at zero depth, as is commonly the case for treating superficial tumors, one may simply use zp100 = 0.

80%-to-20% distal-falloff length

The distance between the distal-80% and distal-20% depths, denoted by ld80-d20. Other distal-falloff lengths are similarly defined, e.g., the 90%-to-10% distal-falloff length ld90-d10.

Proximal-80%-to-distal-80% pristine-peak width

The distance between the proximal-80% depth and the distal-80% depth, denoted by lpristine_d80-p80, where the proximal-80% and distal-80% depths are defined analogously to the distal-50% depth, zd50, already defined. Other pristine peak widths are similarly defined, e.g., the 80%-to-80% pristine peak width. In cases where the Bragg curve does not include a proximal 90 dose point, the methods described in “Proximal-50% depth” may be used.

Proximal-80%-to-distal-80% modulated peak width

The width of the modulated peak region is defined as lmod_d80_p80 = zd80 - zp80. In cases where the Bragg curve does not include a proximal 80 dose point, the methods described in “Proximal-50% depth” may be used.

Modulated-peak region

The region extending from za to zb. In general, the values of za and zb are most reliably determined using iterative numerical fitting methods. Conceptually, they are closely related to the proton ranges of the most and least penetrating pristine peaks in the SOBP.

These definitions may initially appear pedantic and overly precise for a conceptual understanding of proton Bragg curves. However, experience has shown that these definitions facilitate quantitative analysis and reporting of the characteristics of a wide variety of Bragg peaks in clinical and research settings. The definitions were developed from experience in manual and algorithmic analyses of measured clinical pristine and spread-out Bragg curves (Newhauser, 2001b, a; Newhauser et al, 2002a; Newhauser et al, 2002b) and for developing and commissioning proton dose algorithms for treatment planning purposes (Koch and Newhauser, 2005; Newhauser et al, 2007b; Koch and Newhauser, 2010). Finally, we note that in some clinical situations, the practicing medical physicist may need to define and use additional parameters, e.g., zp98 and zd98, as appropriate to a particular clinical situation or protocol. Regardless of the particular parameters chosen, it is difficult to overstate the paramount importance of using the parameters in a consistent way and being clear about their meaning when reporting results.

3.3 Model of pristine Bragg curves

In the preceding sections, we described all the major physical processes that govern the shape of proton dose distributions. Pristine Bragg curves may be calculated using a wide variety of techniques, from look-up tables of measured data to Monte Carlo simulations to analytical models. In the early days of proton therapy, indeed through the 1990s, dose algorithms in most proton treatment planning systems included very few, if any, physical models in their dose calculation algorithms. However, in the last 15 years of so, much progress has been made in modeling proton dose distributions with increasing physical completeness, realism, and accuracy. In this section, we review a method to calculate a pristine Bragg curve using a physics-based analytical model. For brevity, we present only one model, a model that describes many of the important physical processes, is computationally straightforward and fast, and is of considerable practical value in a clinical proton setting.

Bortfeld (1997) proposed an analytical equation to calculate the Bragg curve for proton energies between 10 and 200 MeV, as follows:

D(z)=Φ0eζ24σ1pΓ(1p)2πρpα1p(1+βR0)[1σD1p(ζ)+(βp+γβ+εR0)D1p1(ζ)],
(36)

where D(z) is the depth dose, z is the depth, Φ0 is the primary fluence, R0 is the range of the proton beam, σ is the standard deviation of the Gaussian distribution of the proton depth, ζ = (R0z)/σ, α and p are material-dependent constants (defined in Equation (2)), ε is the fraction of low-energy proton fluence to total proton fluence, Γ(x) is the gamma function, and Dy(x) is the parabolic cylinder function. Very good agreement is found between measured and calculated curves, and this approach has been used for the analysis and characterization of clinical proton beams (Figure 17) (Newhauser et al, 2002b) and for treatment planning dose calculations (Szymanowski and Oelfke, 2003; Soukup et al, 2005; Rethfeldt et al, 2006; Li et al, 2008).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0017.jpg

Absorbed dose D as a function of depth z in water from unmodulated (pristine) proton Bragg peaks produced by a broad proton beam with an initial energy of up to 235 MeV. With increasing depth, the accumulation of range straggling tends to broaden the peak. Beams that are more penetrating therefore have larger distal falloff distances.

3.4 Model of spread-out Bragg curves

Spread-out Bragg curves may be designed by combining multiple pristine Bragg curves. This approach was described in Wilson's seminal paper and its implementation was described by Koehler et al. (Koehler et al, 1977). Figure 15 plots several pristine Bragg curves and their resultant sum, which was similar to the SOBPs used at Harvard from the early 1970s onwards. The fluence is also modulated, which can be seen in the same figure as variations in the relative peak dose of each pristine Bragg curve. The range is typically modulated in steps that are small compared with the width of the Bragg peak, so that the SOBP contains little or no ripple. The modulation step size may be fixed or variable; generally the smallest step increments are needed for the deepest peaks, e.g., the amplitude of the ripple decreases with depth because an increasing fraction of dose is from the sub-peak region of the Bragg curves. Figure 18 plots several typical clinical SOBPs from a contemporary proton therapy system.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0018.jpg

Absorbed dose D as a function of depth z in water from various spread out Bragg peaks (SOBP) from the Northeast Proton Therapy Center (NPTC) nozzle. Open- beam central-axis depth-dose curves are plotted for each of the eight range modulator wheel tracks. The experimental values (open circles) are from charge measurements with an air-filled ionization chamber, and the model fits are shown with solid lines. Each range modulation track was designed to work over an interval of proton energies. The flat modulated-peak region was achieved by modulating the proton-beam current synchronously with the modulator wheel rotation.

SOBPs can also be modeled with analytical methods. Bortfeld and Schlegel (1996) proposed a model of the form

D(z)={34+34πln1+z^21z^+z^232πarctan2z^13:forz<za1:forzazzb0:forzzb},
(37)

where the dimensionless depth is given by

z^=zazzbza3,
(38)

and where za and zb denote the depths of the proximal and distal extents of the SOBP (see preceding section for details of their definitions) and z is the depth in water. The model was derived with several simplifying assumptions and approximations that limit its ability to predict real SOBPs. In particular, the proton energy loss rate is based on the CSDA (i.e., energy straggling is neglected) and the converted energy per unit mass approximation (where the proton energy loss and the absorbed dose in the water are assumed proportional), and non-elastic interactions of the protons with the absorber nuclei were not considered. Consequently, it does not reproduce real SOBPs well.

To attain a more realistic model, one of us (WDN) introduced terms for the finite distal-falloff distance associated with range straggling, protonic buildup effects, arbitrary slope of the modulated-peak region, and a term for the ripple in the modulated peak. In addition, the piecewise use of the function for various regions was modified to allow for a transition region between the sub-peak and modulated-peak regions, which eliminated a pronounced discontinuity there. With these additional terms and modified regions, the model becomes

D(z)=c0ΛβBUγ{c1+c24πln1+z^21z^+z^2c32πarctanc4z^1c5:forzzaΔmtz+bt:forzaΔ<z<za+Δ(mz+b)+Y:forzza+Δ},
(39)

where is defined as before, c0 is a proportionality constant, Λ takes into account the dependence on the dose due to inverse-square-law reductions in the proton fluence, βBU is the secondary-particle buildup term, c1 ~ c5 are coefficients that influence the shape of the sub-peak region, and m and b are the slope and intercept the of the modulated-peak region, respectively. The ripple in the modulated-peak region is described by the term ϒ. The distal falloff term is denoted by γ and is modeled with a cumulative normal PDF. The transition region is centered about za and is of width Δ, which is a non-critical parameter that can be estimated from the amount of range straggling expected at depth za (see previous section for detailed definition) or it may be empirically deduced from measured Bragg curves. The extended model agrees well with measured SOBPs, several of which are shown in Figure 18. This work was previously unpublished.

3.5 Analytical transport of therapeutic proton beams

Analytical models to predict therapeutic radiation are generally simple, fast, and accurate. Early treatment planning systems, many of which were still in use in the 1990s, used dose algorithms that looked up measured profiles of absorbed dose and/or used empirical formulas. Over the last 15 years, advances in our understanding of proton therapy physics, as well as huge gains in computer speed and memory capacity, have spurred the development and clinical use of increasingly realistic and complete physical models to predict absorbed dose in patients.

The basic approach underlying most analytical dose algorithms is that the proton beam's energy loss and lateral scattering may be modeled independently from one another. The physical basis for this assumption lies in the principle of conservation of momentum and energy. Lateral scattering, which is predominantly caused by MCS, occurs at small angles (see section 2.4). Hence, the energy loss associated with the vast majority of these small scattering events is negligible. Conversely, the predominant proton energy-loss mechanism, namely inelastic Coulomb interactions with atomic electrons, only negligibly deflects the protons because the proton mass is 1832 times greater than that of an electron (see section 2.1). Consequently, for an excellent approximation, the three-dimensional distribution of absorbed dose in a phantom or patient may be cast in the form of

D(x, y, z) = D(z) ⋅ OAR(x, y, z),
(40)

where x and y are lateral coordinates, z is the depth coordinate along the beam axis, D(z) is the large-field Bragg curve (i.e., with lateral protonic equilibrium, described in section 3.2), and OAR(x, y, z) is a term that takes into account the lateral properties of the dose distribution. Broad beam algorithms use Equation (40) with simple ray-casting techniques (Siddon, 1985) to determine the energy loss and range along each ray. Similarly, to determine the lateral extent of the therapeutic field, the rays are terminated in target-shaped final collimators by approximating them as black-body absorbers (Hong et al, 1996). Carefully designed broad beam algorithms are simple, fast, and accurate in media that are of relatively uniform mass density and material composition. Recently, Koch and Newhauser (2010) developed an analytical broad beam algorithm of the form

D(x,y,z)MU=i=1NDDi(zeff)MUωiSAFωiRMW(SSD+zeffSSD+z)2OARiΓi,
(41)

where N is the number of range modulation steps, ωSAF is the scatter and attenuation fluence weighting factor, ωRMW is the range modulator wheel (RMW) fluence weighting factor, OAR is off-axis ratio, and Γ is collimation modifier. This algorithm explicitly took into account edge-scattered protons and range straggling, predicted absolute dose per monitor unit (D/MU) values, and provided excellent agreement with measurements in homogeneous media, e.g., for ocular treatments. In highly heterogeneous media, however, the dose at some locations may depend strongly on the material properties of the heterogeneities in the overlying tissue, and therefore the approach of casting a single ray does not, in general, provide sufficient dosimetric accuracy.

To improve the dosimetric accuracy in heterogeneous media, the dose may be calculated at a point as the superposition of multiple pencil beams. A pencil beam is, in essence, a mathematical construct that has no exact and direct analogy in the physical world, but it is resembles a beam of infinitesimally small lateral size. By superposing enough pencil beams, the physical dose distribution can be modeled. Fortunately, as we shall see, many of the physical principles and algorithmic aspects from the broad beam approach are directly applicable to the pencil beam algorithm.

Hogstrom et al. (1981) originally proposed a pencil beam algorithm for electron therapy. Basically, an ion beam is divided into a number of finite pencil beams and the dose distribution of each pencil beam is described by the multiplication of a central-axis term and an off-axis term. The final relative absorbed dose to any point in the patient (Dp (x, y, z) ) is the summation of the dose contributions from all the pencil beams:

Dp(x,y,z)=dxdyS(x,y)PDD(zeff)(SSD+zeffSSD+z)212πσtot2(z)exp((xx)2+(yy)22σtot2(z)),
(42)

where S(x′, y′) is the relative intensity of the pencil beam at x′, y′, PDD(z) is the central-axis percentage depth dose, SSD is the source-to-surface distance, zeff is the effective depth of the point, and σtot2describes the final lateral distribution of a pencil beam by taking into account the MCS contributions from each beam-modifying device and from the patient. Many other groups have subsequently adapted and extended this concept for application in proton therapy (Petti, 1992; Hong et al, 1996; Schaffner et al, 1999; Szymanowski and Oelfke, 2002; Ciangaru et al, 2005; Schaffner, 2008; Westerly et al, 2013). Contemporary proton pencil beam algorithms provide excellent accuracy (Newhauser et al, 2007b), especially in homogeneous media, superior to that of broad beam algorithms in heterogeneous media. The improvement in accuracy comes mostly at the cost of greater computation times. Because of approximations in the pencil beam algorithm, it may not provide sufficient accuracy in extremely heterogeneous media, at material and/or density interfaces, or in other complex situations.

Currently, most proton therapy clinics utilize commercial treatment planning systems that contain pencil beam algorithms for calculation of therapeutic dose distributions. Computation speeds are generally sufficiently fast for most routine planning.

3.6 Monte Carlo transport of therapeutic proton beams

Monte Carlo transport technique has been used increasingly in radiation therapy (Seco and Verhaegen, 2013). It is more accurate than analytical models in that it takes into account the physical processes during particle transport. The most commonly used general purpose Monte Carlo codes for proton therapy research are MCNPX (Pelowitz, 2011), Geant 4 (Agostinelli et al, 2003), and FLUKA (Ferrari et al, 2005). Several groups have developed special purpose Monte Carlo codes, e.g., a treatment planning dose engine (Tourovsky et al., (2005); Yepes et al., (2009a)), typically with faster computation speed but fewer features than general purpose codes. Here we focus mainly on applications that use the MCNPX code because it is representative of most general purpose codes, we are most familiar with it of any of the codes, and constraints on journal space do not allow discussion of all codes.

The main components of a Monte Carlo simulation model are a source of protons, a treatment head, and a patient or phantom (figure 19). The proton source is typically taken as being incident on the treatment head (beamline transport is not usually simulated), where the source parameters are deduced from dose profile measurements (Newhauser et al, 2005) or taken from other beam properties, e.g., from the manufacturer (Newhauser et al, 2007b). The model of a treatment head typically includes all major mechanical components, e.g., a vacuum window, a beam profile monitor, a range modulator wheel, a second scatter, shielding plates, a range shifter assembly, backup and primary monitors, the snout and the brass collimator, and a patient-specific range compensator (Zheng et al, 2007b; Perez-Andujar et al, 2009). Models of scanned beams may model the magnetic field and transport the proton beam deflection (Peterson et al, 2009; Dowdell et al, 2012), or the proton source may be pointed, e.g., using the “thin lens” approximation to increase computation speed (Lee, 2004). The patient may be represented in voxels whose material composition and mass density are deduced from three-dimensional CT images of the patient (Taddei et al, 2009). The matrix of Hounsfield unit values in these voxels may be converted into corresponding matrices of material composition and mass density values based on a tri-linear calibration curve (Newhauser et al, 2008a). A simulation capability to accommodate time-dependent geometry was reported by Paganetti et al. (2005). Some useful applications of the Monte Carlo method are described here:

  • 1)

    Alternative or complimentary source of dosimetric data for developing, configuring, and validating analytical dose algorithms in clinical treatment planning systems. At least two studies revealed that, after a significant development effort, a Monte Carlo model of a proton therapy apparatus is sufficiently accurate and fast for prospective commissioning of an ocular treatment planning system (Koch and Newhauser, 2005; Newhauser et al, 2005) and a large-field treatment planning system (Newhauser et al, 2007b). This prospective capability reduced the time required for preclinical work by several months, resulting in substantial cost savings. Studies about scanning proton beam configuration have also been reported (Tourovsky et al, 2005; Sawakuchi et al, 2010; Dowdell et al, 2012).

  • 2)

    Alternative dose engine to calculate absorbed dose in phantoms or patient. In clinical practice, there frequently arise challenging cases where the accuracy of analytical models is unknown or suspected to be insufficient. In many such cases, measurements are not possible or practical, e.g., intracranial in vivo range verification. Multiple groups have developed Monte Carlo treatment planning engines (Tourovsky et al, 2005; Newhauser et al, 2008b; Perl et al, 2012; Mairani et al, 2013) based on various Monte Carlo codes, independently proving the technical feasibility of this approach. The physics of general purpose Monte Carlo codes for this purpose has been exhaustively validated against benchmark measurements in therapeutic proton beam data (Polf JC, 2007; Kimstrand et al, 2008; Titt et al, 2008a; Randeniya et al, 2009). Just 10 years ago, the slow execution times, memory constraints, unknown accuracy, and unknown technical feasibility of a full-blown Monte Carlo engine were serious research challenges. Today, the only major remaining obstacle to routine use of Monte Carlo engines for treatment planning is their slow execution times. For example, a proton craniospinal treatment simulation took more than 10 CPU hours (Taddei et al, 2009). Encouragingly, however, algorithmic improvements have increased speeds dramatically; Zhang et al. (2013a) found that using lattice tally in MCNPX can speed up dose simulations by one order of magnitude compared to mesh tally, especially for dose reconstructions involving large numbers of voxels. The physics of the radiation transport are identical for the lattice and mesh tallies, i.e., there was no loss of accuracy. Yepes et al. (2009b) recently reported on a fast dose calculator, a Monte Carlo track-repeating algorithm that uses a database of precomputed proton trajectories in water and applies these trajectories to heterogeneous media by scaling the proton range and scattering angles. This approach reduced the computation time for dose reconstruction in voxelized patient geometry by more than two orders of magnitude compared to MCNPX and GEANT4 (Yepes et al, 2009a; Yepes et al, 2010a). The use of low-cost alternative parallel computer architectures also appears promising, including the use of graphics processing units (Yepes et al, 2010b; Jia et al, 2012) and grid technologies (Vadapalli et al, 2011).

  • 3)

    Assessment of beam perturbation due to objects implanted in patients and development of mitigating strategies. The proportion of cancer patients with implanted devices is probably less than 10%, but it is increasing along with the aging of the population, increasing life expectancies, and increasing use of implanted devices. Examples include cardiac pacemakers, defibrillators, fiducial markers for radiotherapy localization, permanent radioactive seed implants, stents, prostheses, surgical reconstructive devices, and foreign bodies such as gunshot and shrapnel. Analytical dose algorithms in contemporary clinical treatment planning systems lack the capability to reliably assess dose perturbations from most implants. Several studies reported that implanted metallic fiducials can cause dose shadows that may compromise local tumor control for certain diseases, e.g., prostate carcinoma and uveal melanoma (Newhauser et al, 2007a; Newhauser et al, 2007c; Giebeler et al, 2009; Cheung et al, 2010; Huang et al, 2011; Matsuura et al, 2012; Zhang et al, 2012; Carnicer et al, 2013).

  • 4)

    Investigation of characteristics of therapeutic beams and beamlines. For example, edge-scattered protons can degrade the dose distribution, and the effect is difficult to model with current analytical models. Accurate predictions of D / MU require taking into account the dose from protons scattered from the edge of the patient-specific collimator, particularly for small field sizes and large depths (Titt et al, 2008b). Currently available spot scanning system only offer few options for adjusting beam spot properties like lateral and longitudinal size, and Monte Carlo simulations was used to optimize scanned beam spot (Titt et al, 2010). Because the Monte Carlo simulation technique is inherently general, it can be applied to virtually any problem involving the transport of protons or secondary radiation. Space constraints prevent us from reviewing more than a few illustrative examples. In later sections we briefly mention roles for Monte Carlo methods in research on stray radiation exposures to patients (section 5) and shielding barriers to protect staff (section 6).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0019.jpg

(Upper) Geometric model of proton therapy unit and the voxelized phantom oriented for the superior spinal proton field. The beam delivery system includes a vacuum window (A), a beam profile monitor (B), a range modulator wheel (C), a second scatterer (D), a range shifter assembly (E), backup and primary monitors (F), the snout (G), the range compensator (H), the treatment couch (I), and the patient (J). (Zhang et al, 2013a) (Lower) Simulated proton transport through the proton therapy unit and the patient.

3.1 Water-equivalent thickness

As we mentioned previously, in proton therapy water closely mimics the properties of human tissues in terms of energy loss, MCS, and nuclear interactions. As such, water is a recommended phantom material for dose and range measurements and reference material for reporting corresponding calculated quantities (ICRU, 1998; IAEA, 2000). For example, it is a common and convenient clinical practice to specify the penetrating power of a proton beam by its range in water (ICRU, 1998, 2007). In this way, range losses in various beamline objects and the patient may be easily added or subtracted from one another in a physically consistent way. Viewed another way, it is also convenient to specify the range-absorbing power of various objects in the beam path, e.g., beam transmission monitors and immobilization devices, by their equivalent thickness if they were made of water.

Water-equivalent thickness (WET) is often used to characterize the beam penetration range; Figure 9 schematically illustrates the concept of WET and how it can be calculated or measured. For treatment sites with nearby critical structures, e.g., an optic nerve, the range of the planned and delivered beams must agree within a few millimeters. To accomplish this, treatment planning systems are commonly configured with the WET values of all items not included in the planning CT images, such as components in the treatment head, immobilization devices not present during the CT scan, or a treatment couch (Newhauser et al, 2007b). Similarly, to determine the measurement geometry for patient-specific clinical quality assurance measurements, the WET of measurement instruments and possibly other devices must be determined (Newhauser, 2001b, a). Thus, it is important to have methods to calculate and measure WET. In this section, we emphasize recently developed calculation methods that are convenient and suitable for clinical calculations, using the energy loss theories presented in section 2. Our discussion of WET measurement methods is very brief, mainly because they are relatively simple and obvious. In practice, however, WET measurements remain very important (Andreo, 2009; Gottschalk, 2010a; Newhauser and Zhang, 2010; Zhang et al, 2010b; Besemer et al, 2013).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0009.jpg

Schematic illustration of the concept of water equivalent thickness (WET) and how it can be calculated or measured by observing changes in the depth of a proton Bragg peak in a water tank. (Zhang and Newhauser, 2009)

The IAEA (2000) proposed that WET can be approximated by

tw=tmρmρwcm,
(31)

where the depth scaling factor, cm, can be calculated, to a good approximation, as the ratio of the continuous slowing down approximation (CSDA) range (in g cm) in water to that in the target:

cm=ρwRwρmRm.
(32)

Because the ranges in Equation (32) correspond to complete loss of ion energy, this approach is strictly valid only for stopping-length targets. An exact equation for WET that is applicable to thin targets was reported by Zhang and Newhauser (2009):

tw=tmρmSm¯ρwSw¯,
(33)

where ρw and ρm are the mass densities of water and material, respectively, and S̄m and S̄w are the mean proton mass stopping power values for the material and water, respectively, defined by

S=ESdEEdE.
(34)

For thin targets, where the proton loses a negligible fraction of its energy in the absorber material, we have

twtm(ρmρw)(SmSw)=tmαwpwαmpmEpwpm,
(35)

where the reader will recognize α and p from the discussion of the BK Rule (see section 2.1). Values of α and p for commonly encountered materials in proton therapy are provided in Table 2. Zhang and Newhauser (2009) also reported a slightly more complex analytical formula to calculate WET for targets of arbitrary thickness.

Table 2

Common materials (lung substitute plastic, high-density polyethylene (HDPE), water, polystyrene, polymethylmethacrylate (PMMA), polycarbonate resin (Lexan), bone substitute plastic, aluminum, titanium, stainless steel, lead and gold) used in heavy charged particle beams, with their mass densities ρ, values of (Z/A)eff, mole fractions and fitting parameters of α and p for these materials when applying the BK rule (The energies used in the fit were from 10 to 250 MeV).

Materialρ (g/cm)Mole fraction (%)αp
Lung substitute0.30.537H 55.577, C 32.738, N 0.927, O 7.508, Cl 0.019, Si 0.184, Mg 3.0488.994×10−31.735
HDPE0.9640.570H 66.717, C 33.2832.541×10−31.737
Water1.00.555H 66.667, O 33.3332.633×10−31.735
Polystyrene1.060.538H 49.851, C 50.1492.545×10−31.735
PMMA1.1850.539H 53.333, C 33.333, O 13.333,2.271×10−31.735
Lexan1.200.527H 42.424, C 48.485, O 9.0912.310×10−31.735
Bone substitute1.8290.516H 35.215, C 29.592, N 0.803, O 26.695, Cl 0.16, Ca 7.6791.666×10−31.730
Aluminum2.6980.482Al 1001.364×10−31.719
Titanium4.5190.459Ti 1009.430×10−41.710
Stainless steel7.850.466C 0.045, N 0.045, Si 0.450, Cr 18.150, Mn 1.250, Fe 71.460, Ni 8.550, Mo 0.0505.659×10−41.706
Lead11.3220.396Pb 1006.505×10−41.676
Gold19.3110.401Au 1003.705×10−41.677

As can be seen in the curves of the water-equivalent ratio (WER = tw/tm) plotted in Figure 10, taking into account the target thickness in calculating WER is most important for absorbers that are thick and made of high Z materials, e.g., lead scattering foils, and for protons that are of comparatively low energy when impinging on the target. For low-Z materials, such as tissue and plastic, WER depends only very weakly on the target thickness and initial proton beam energy, and the approximate (thin and stopping length) analytical methods provide sufficient accuracy for most clinical applications.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0010.jpg

Calculated water-equivalent ratio (WER) values as a function of proton beam energy. This plot illustrates the dependence of WER value on the target material, the beam energy, and the target thickness. (Zhang and Newhauser, 2009)

3.2 Features of pristine and spread-out Bragg curves

The spatial dose distribution from clinical proton therapy beams is quite similar to those from photon and electron beams. The lateral profiles are generally quite flat in the central high-dose region, then fall off rapidly in the penumbral regions, where the penumbra width increases with depth in the patient. On the other hand, the central-axis depth-dose curve from protons is somewhat similar to that from electrons, but with a sharper distal falloff. Figures 11 and and1212 compare the central-axis depth-dose curves from several radiation therapy beams, revealing the main dosimetric properties that are clinically advantageous in many cases, namely, relatively low entrance dose, large and uniform dose to cover the tumor, and rapid falloff of dose near the end of range to spare normal tissues. These properties, together with a uniform lateral dose profile and a sharp lateral penumbral width, allow proton beams to treat a wide variety to tumor sizes and locations while providing superior sparing of normal tissues in many cases.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0011.jpg

Central axis depth dose profiles from several particle beams. Note that these distributions are from solitary beams in order to clearly compare the differences in the physical properties of various radiations. The important features are that proton beams offer relatively low entrance dose and virtually no exit dose. However, many clinical treatment techniques exploit multiple field directions to enhance the uniformity of tumor coverage and to spare sensitive healthy tissues. In fact, in some cases proton treatments provide inferior skin sparing to photons and/or inferior target coverage, e.g., because of proton beams’ sensitivity to range errors. Nonetheless, beam for beam, proton beams provide excellent tissue sparing, especially beyond the end of range. (Larsson, 1993)

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0012.jpg

Comparison of depth-dose curves from proton SOBPs and from electron beams. Because the proton mass is nearly 2000 times that of an electron, proton scattering interactions (individual angular deflections and variations in collisional energy losses) are much smaller, leading to sharper lateral and distal falloff distances. (Koehler and Preston, 1972)

Having casually inspected the shape of proton depth-dose curves, we next examine their structure in greater detail, pointing out nomenclature and the physical processes that govern the shape of various regions. Figure 13 shows a pristine proton peak along with labels identifying several regions. In order of increasing depth, these are the regions of electronic buildup, protonic buildup, sub-peak, peak, and distal falloff (same as for peak). The figure also shows several characteristic depths (e.g., the depth zBP at which the Bragg peak occurs) and various characteristic lengths (e.g., the 80%-to-20% distal-falloff length l80-20 and the proximal-80%-to-distal-80% pristine-peak width).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0013.jpg

Absorbed dose D as a function of depth z in water from an unmodulated (pristine) proton Bragg peak produced by a broad proton beam with an initial energy of 154 MeV. The various regions, depths, and lengths that are labeled are defined in the text. (The electronic buildup is not visible in this plot.) This type of dose distribution is useful clinically because of the relatively low doses delivered to normal tissues in the sub-peak and distal-falloff regions relative to the target dose delivered by the peak.

The anatomic definitions of an SOBP are, in many ways, similar to those of a pristine Bragg curve, as seen in Figure 14. However, there are several unique difficulties in characterizing SOBPs because of their sometimes unusual shape. For example, SOBPs with two or more discrete pristine Bragg curves may have multiple dose maxima in the modulated-peak region (e.g., the ripple shown in Figure 15). Because of such problems we introduce a few additional quantities that are defined only for SOBPs. To a large extent, however, we have defined quantities and terminology that are common to both modulated and pristine Bragg curves.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0014.jpg

Absorbed dose D as a function of depth z in water from a spread-out proton Bragg peak (SOBP). Various locations and regions that are indicated on the plot are defined in the text. This peak was measured with a Markus-type parallel-plate ionization chamber in the Northeast Proton Therapy Center (NPTC) gantry. The measured data are plotted with open circles and the model-fit as a solid line. Note that the electronic buildup region, which spans only a few millimeters, is not visible in this plot.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0015.jpg

Absorbed dose D as a function of depth z in water from a spread-out Bragg peak (SOBP) (uppermost curve) and its constituent pristine Bragg peaks (lower curves; for clarity, all but the deepest pristine Bragg peak are only partly drawn). In many cases, the clinical target volume is larger than the width of a pristine Bragg peak. By appropriately modulating the proton range and fluence of pristine peaks, the extent of the high-dose region can be widened to cover the target volume with a uniform dose.

We have not yet mentioned how MSC affects the shape of the depth-dose curves. In fact, near the central region of a laterally “large” beam, or more correctly well inside the periphery of a large beam, there is an equilibrium in which lateral scattering away from the central axis is exactly compensated by scattering towards it. This effect is described in Figure 16, which is adapted from Koehler and Preston (Koehler and Preston, 1972). As the field size shrinks to the dimension of the rms lateral displacement due to MCS, lateral equilibrium is lost and MCS progressively depletes the proton fluence and dose along the central axis. Small proton beams have been investigated in several studies, including those by Takada (Takada, 1996), Moyers et al. (1999), Vatnitsky et al. (1999b), Bednarz et al. (2010), and Gottschalk (1999), as well as others, especially in the context of scanned beams and pencil beam dose algorithms.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0016.jpg

(Left) Proton fluence I(0; x) along the beam central axis vs. depth x in water. Curves are shown for beams with circular cross sections and radii of 1 mm to 4 mm. Some of the protons are lost because of scattering events that deflect them from the central axis. This is increasingly observed for small beams and at large depths. (Right) The corresponding central-axis absorbed-dose curves. Note how the fluence depletion reduces the absorbed dose at the peak relative to the entrance dose. (Preston and Koehler, 1998)

Here, we use a Cartesian coordinate system with the z axis parallel to and centered about the proton beam central axis. The x and y axis are mutually orthogonal and perpendicular to the z axis. The coordinate system origin is located at the front face of the absorber, e.g., the extended medium in which we consider the absorbed dose distribution.

Pristine Bragg curve

A depth-dose distribution in an absorber irradiated with a monoenergetic or nearly-monoenergetic proton beam. In other words, no device or technique has been intentionally deployed for modulating the proton fluence or spectral fluence.

Spread-out Bragg curve

A depth-dose distribution in an absorber irradiated with a beam that has been intentionally modified to increase the axial dimension of the peak region. This is accomplished by modulating the range and the fluence of the beam. Clinical systems accomplish this by combining multiple quasi-monoenergetic beams or with a continuously modulated beam.

Electronic buildup region

A small region near the surface of the absorber where the proton beam is incident. As discussed in section 2.1, high-energy proton beams liberate delta rays with sufficient kinetic energy to travel several millimeters in tissue. Under some circumstances, this region exhibits an increase of dose with increasing depth, asymptotically approaching absorbed dose in the sub-peak region within the depth corresponding to the range of the most penetrating recoil electron. In some cases, electronic buildup is not observed. There are several possible reasons for this: the presence of some material just upstream of the surface (e.g., an immobilization device or a range compensator) may provide partial or full electronic charged particle equilibrium in the absorber, it may occur in combination with protonic buildup, it may be masked by changes in the proton energy loss rate near the end of range, or the wall of a cavity dosimeter may be sufficiently thick to present electronic equilibrium to the dosimeter's sensitive volume.

Protonic buildup region

A region near the surface of the absorber where the absorbed dose increases with depth because of the buildup of secondary protons that are attributable to proton-induced non-elastic nuclear interactions (e.g., O(p,xp) reactions). As with electronic buildup, the protonic buildup may not be observed in some cases, particularly at low incident proton beam energies.

Sub-peak region

The region extending from the surface of the absorber to the depth just proximal of the peak. The physical processes involved here are, in decreasing order of importance, the stopping power's dependence on the inverse-square of the proton velocity, the removal of some protons by nuclear reactions, the liberation of secondary particles from nuclear reactions, and for very small fields, the accumulation of lateral deflections from MCS leading to lateral protonic disequilibrium and reduction of the proton fluence on the central axis. The distal extent of the sub-peak region can be calculated from zm- 2σ, where zm is the depth at the pristine Bragg peak and σ is the width of the peak. The width parameter σ can be estimated from the incident proton beam spectral fluence and the range straggling accumulated in the absorber, as discussed in section 2.3.

Pristine Bragg peak

The pristine Bragg peak is simply the maximum (or mode) dose near the end of range, and is located at zBP, which is defined next. The physical processes governing the location and/or height of the peak are mainly the proton stopping power and energy straggling, nuclear reactions to a much lesser extent and, for very small fields, MCS..

Pristine Bragg peak depth

The depth near the end of range of the primary protons at which the protons produce the maximum dose rate, denoted by zBP. Although small proton beams are not yet widely used, it is helpful to define the location of zBP in a way that is compatible with large and small beams. Figure 16 shows that the maximum dose for beams of diameter larger than 6 mm is clearly single valued and located near the end of range. For smaller beams, however, the dose at the peak near the end of range may be less than the dose in the proximal regions, creating multiple maxima to choose from. Hence, the definition of zBP restricts it to exist in the region of the R±4σ, where sigma is the distal falloff width, thereby preventing possible ambiguities, and makes zBP conceptually independent of the beam cross-sectional area.

Distal falloff region

This region extends from depths greater than that of the pristine Bragg peak depth, zBP. The width of this region is not restricted. In many practical situations, however, the distal falloff region terminates at a depth where the dose falls below a threshold value, e.g., 1% of the dose at the Bragg peak, D(zBP).

Distal-50% depth

The distalmost depth, denoted by zd50, at which the absorbed dose is equal to half of the absorbed dose at the pristine Bragg peak depth, or D(zBP)/2. For an SOBP, it is defined as the distalmost depth at which the absorbed dose is equal to half of the absorbed dose at the SOBP dose in the modulated peak region. In many cases, the dose in the modulated peak region varies with depth (perhaps by design), making selection of the absorbed dose in the modulated peak region arbitrary. Since this clearly hinders an unambiguous definition of zd50, we instead define zd50 for an SOBP as being equal to zb. If the dose in the modulated peak region varies with depth, zb may have to be determined with numerical methods. Physically, the value of zb is closely related to the RCSDA value of the mean proton energy corresponding to the most energetic peak in the SOBP. Definitions for various other distal depths are similarly defined, e.g., the distal-90% depth zd90 , zd50, and zd20,.

Proximal-50% depth

The second most distal depth, denoted by zp50, at which the absorbed dose is equal to half of the absorbed dose at the pristine Bragg peak depth, or D(zBP)/2, provided that occurs within the absorber. (For very low-energy pristine Bragg curves, the entrance dose may be greater than half the peak dose). In such cases, and for SOBPs, it is defined as the most proximal depth at which the absorbed dose is equal to half of the absorbed dose in the modulated peak region. Because of problems that are analogous to those described in the zd50 definition above, zp50 is defined as the location at which the dose is equal to half the value at depth za + Δ/2, where Δ is the width of the transition from the peak region to the sub-peak region. Δ is not a critical parameter and it may be estimated as the amount of range straggling at za or it can be determined from measured Bragg curves. Conceptually, za typically corresponds to the shallowest location that is expected to receive the maximum dose. Physically, the value of za is closely related to the RCSDA value of the least energetic pristine peak in the SOBP curve (or the only peak in the Bragg curve for a pristine peak). In many cases, particularly for SOBPs with large modulated-peak regions, the absorbed dose throughout the sub-peak region exceeds 50% of the value at za + Δ, in which case proximal 50% depth is undefined. In such cases, proximal depths at higher dose percentages, such as zp95 and zp90, may be used. If the Bragg peak occurs at zero depth, as is commonly the case for treating superficial tumors, one may simply use zp100 = 0.

80%-to-20% distal-falloff length

The distance between the distal-80% and distal-20% depths, denoted by ld80-d20. Other distal-falloff lengths are similarly defined, e.g., the 90%-to-10% distal-falloff length ld90-d10.

Proximal-80%-to-distal-80% pristine-peak width

The distance between the proximal-80% depth and the distal-80% depth, denoted by lpristine_d80-p80, where the proximal-80% and distal-80% depths are defined analogously to the distal-50% depth, zd50, already defined. Other pristine peak widths are similarly defined, e.g., the 80%-to-80% pristine peak width. In cases where the Bragg curve does not include a proximal 90 dose point, the methods described in “Proximal-50% depth” may be used.

Proximal-80%-to-distal-80% modulated peak width

The width of the modulated peak region is defined as lmod_d80_p80 = zd80 - zp80. In cases where the Bragg curve does not include a proximal 80 dose point, the methods described in “Proximal-50% depth” may be used.

Modulated-peak region

The region extending from za to zb. In general, the values of za and zb are most reliably determined using iterative numerical fitting methods. Conceptually, they are closely related to the proton ranges of the most and least penetrating pristine peaks in the SOBP.

These definitions may initially appear pedantic and overly precise for a conceptual understanding of proton Bragg curves. However, experience has shown that these definitions facilitate quantitative analysis and reporting of the characteristics of a wide variety of Bragg peaks in clinical and research settings. The definitions were developed from experience in manual and algorithmic analyses of measured clinical pristine and spread-out Bragg curves (Newhauser, 2001b, a; Newhauser et al, 2002a; Newhauser et al, 2002b) and for developing and commissioning proton dose algorithms for treatment planning purposes (Koch and Newhauser, 2005; Newhauser et al, 2007b; Koch and Newhauser, 2010). Finally, we note that in some clinical situations, the practicing medical physicist may need to define and use additional parameters, e.g., zp98 and zd98, as appropriate to a particular clinical situation or protocol. Regardless of the particular parameters chosen, it is difficult to overstate the paramount importance of using the parameters in a consistent way and being clear about their meaning when reporting results.

Pristine Bragg curve

A depth-dose distribution in an absorber irradiated with a monoenergetic or nearly-monoenergetic proton beam. In other words, no device or technique has been intentionally deployed for modulating the proton fluence or spectral fluence.

Spread-out Bragg curve

A depth-dose distribution in an absorber irradiated with a beam that has been intentionally modified to increase the axial dimension of the peak region. This is accomplished by modulating the range and the fluence of the beam. Clinical systems accomplish this by combining multiple quasi-monoenergetic beams or with a continuously modulated beam.

Electronic buildup region

A small region near the surface of the absorber where the proton beam is incident. As discussed in section 2.1, high-energy proton beams liberate delta rays with sufficient kinetic energy to travel several millimeters in tissue. Under some circumstances, this region exhibits an increase of dose with increasing depth, asymptotically approaching absorbed dose in the sub-peak region within the depth corresponding to the range of the most penetrating recoil electron. In some cases, electronic buildup is not observed. There are several possible reasons for this: the presence of some material just upstream of the surface (e.g., an immobilization device or a range compensator) may provide partial or full electronic charged particle equilibrium in the absorber, it may occur in combination with protonic buildup, it may be masked by changes in the proton energy loss rate near the end of range, or the wall of a cavity dosimeter may be sufficiently thick to present electronic equilibrium to the dosimeter's sensitive volume.

Protonic buildup region

A region near the surface of the absorber where the absorbed dose increases with depth because of the buildup of secondary protons that are attributable to proton-induced non-elastic nuclear interactions (e.g., O(p,xp) reactions). As with electronic buildup, the protonic buildup may not be observed in some cases, particularly at low incident proton beam energies.

Sub-peak region

The region extending from the surface of the absorber to the depth just proximal of the peak. The physical processes involved here are, in decreasing order of importance, the stopping power's dependence on the inverse-square of the proton velocity, the removal of some protons by nuclear reactions, the liberation of secondary particles from nuclear reactions, and for very small fields, the accumulation of lateral deflections from MCS leading to lateral protonic disequilibrium and reduction of the proton fluence on the central axis. The distal extent of the sub-peak region can be calculated from zm- 2σ, where zm is the depth at the pristine Bragg peak and σ is the width of the peak. The width parameter σ can be estimated from the incident proton beam spectral fluence and the range straggling accumulated in the absorber, as discussed in section 2.3.

Pristine Bragg peak

The pristine Bragg peak is simply the maximum (or mode) dose near the end of range, and is located at zBP, which is defined next. The physical processes governing the location and/or height of the peak are mainly the proton stopping power and energy straggling, nuclear reactions to a much lesser extent and, for very small fields, MCS..

Pristine Bragg peak depth

The depth near the end of range of the primary protons at which the protons produce the maximum dose rate, denoted by zBP. Although small proton beams are not yet widely used, it is helpful to define the location of zBP in a way that is compatible with large and small beams. Figure 16 shows that the maximum dose for beams of diameter larger than 6 mm is clearly single valued and located near the end of range. For smaller beams, however, the dose at the peak near the end of range may be less than the dose in the proximal regions, creating multiple maxima to choose from. Hence, the definition of zBP restricts it to exist in the region of the R±4σ, where sigma is the distal falloff width, thereby preventing possible ambiguities, and makes zBP conceptually independent of the beam cross-sectional area.

Distal falloff region

This region extends from depths greater than that of the pristine Bragg peak depth, zBP. The width of this region is not restricted. In many practical situations, however, the distal falloff region terminates at a depth where the dose falls below a threshold value, e.g., 1% of the dose at the Bragg peak, D(zBP).

Distal-50% depth

The distalmost depth, denoted by zd50, at which the absorbed dose is equal to half of the absorbed dose at the pristine Bragg peak depth, or D(zBP)/2. For an SOBP, it is defined as the distalmost depth at which the absorbed dose is equal to half of the absorbed dose at the SOBP dose in the modulated peak region. In many cases, the dose in the modulated peak region varies with depth (perhaps by design), making selection of the absorbed dose in the modulated peak region arbitrary. Since this clearly hinders an unambiguous definition of zd50, we instead define zd50 for an SOBP as being equal to zb. If the dose in the modulated peak region varies with depth, zb may have to be determined with numerical methods. Physically, the value of zb is closely related to the RCSDA value of the mean proton energy corresponding to the most energetic peak in the SOBP. Definitions for various other distal depths are similarly defined, e.g., the distal-90% depth zd90 , zd50, and zd20,.

Proximal-50% depth

The second most distal depth, denoted by zp50, at which the absorbed dose is equal to half of the absorbed dose at the pristine Bragg peak depth, or D(zBP)/2, provided that occurs within the absorber. (For very low-energy pristine Bragg curves, the entrance dose may be greater than half the peak dose). In such cases, and for SOBPs, it is defined as the most proximal depth at which the absorbed dose is equal to half of the absorbed dose in the modulated peak region. Because of problems that are analogous to those described in the zd50 definition above, zp50 is defined as the location at which the dose is equal to half the value at depth za + Δ/2, where Δ is the width of the transition from the peak region to the sub-peak region. Δ is not a critical parameter and it may be estimated as the amount of range straggling at za or it can be determined from measured Bragg curves. Conceptually, za typically corresponds to the shallowest location that is expected to receive the maximum dose. Physically, the value of za is closely related to the RCSDA value of the least energetic pristine peak in the SOBP curve (or the only peak in the Bragg curve for a pristine peak). In many cases, particularly for SOBPs with large modulated-peak regions, the absorbed dose throughout the sub-peak region exceeds 50% of the value at za + Δ, in which case proximal 50% depth is undefined. In such cases, proximal depths at higher dose percentages, such as zp95 and zp90, may be used. If the Bragg peak occurs at zero depth, as is commonly the case for treating superficial tumors, one may simply use zp100 = 0.

80%-to-20% distal-falloff length

The distance between the distal-80% and distal-20% depths, denoted by ld80-d20. Other distal-falloff lengths are similarly defined, e.g., the 90%-to-10% distal-falloff length ld90-d10.

Proximal-80%-to-distal-80% pristine-peak width

The distance between the proximal-80% depth and the distal-80% depth, denoted by lpristine_d80-p80, where the proximal-80% and distal-80% depths are defined analogously to the distal-50% depth, zd50, already defined. Other pristine peak widths are similarly defined, e.g., the 80%-to-80% pristine peak width. In cases where the Bragg curve does not include a proximal 90 dose point, the methods described in “Proximal-50% depth” may be used.

Proximal-80%-to-distal-80% modulated peak width

The width of the modulated peak region is defined as lmod_d80_p80 = zd80 - zp80. In cases where the Bragg curve does not include a proximal 80 dose point, the methods described in “Proximal-50% depth” may be used.

Modulated-peak region

The region extending from za to zb. In general, the values of za and zb are most reliably determined using iterative numerical fitting methods. Conceptually, they are closely related to the proton ranges of the most and least penetrating pristine peaks in the SOBP.

These definitions may initially appear pedantic and overly precise for a conceptual understanding of proton Bragg curves. However, experience has shown that these definitions facilitate quantitative analysis and reporting of the characteristics of a wide variety of Bragg peaks in clinical and research settings. The definitions were developed from experience in manual and algorithmic analyses of measured clinical pristine and spread-out Bragg curves (Newhauser, 2001b, a; Newhauser et al, 2002a; Newhauser et al, 2002b) and for developing and commissioning proton dose algorithms for treatment planning purposes (Koch and Newhauser, 2005; Newhauser et al, 2007b; Koch and Newhauser, 2010). Finally, we note that in some clinical situations, the practicing medical physicist may need to define and use additional parameters, e.g., zp98 and zd98, as appropriate to a particular clinical situation or protocol. Regardless of the particular parameters chosen, it is difficult to overstate the paramount importance of using the parameters in a consistent way and being clear about their meaning when reporting results.

3.3 Model of pristine Bragg curves

In the preceding sections, we described all the major physical processes that govern the shape of proton dose distributions. Pristine Bragg curves may be calculated using a wide variety of techniques, from look-up tables of measured data to Monte Carlo simulations to analytical models. In the early days of proton therapy, indeed through the 1990s, dose algorithms in most proton treatment planning systems included very few, if any, physical models in their dose calculation algorithms. However, in the last 15 years of so, much progress has been made in modeling proton dose distributions with increasing physical completeness, realism, and accuracy. In this section, we review a method to calculate a pristine Bragg curve using a physics-based analytical model. For brevity, we present only one model, a model that describes many of the important physical processes, is computationally straightforward and fast, and is of considerable practical value in a clinical proton setting.

Bortfeld (1997) proposed an analytical equation to calculate the Bragg curve for proton energies between 10 and 200 MeV, as follows:

D(z)=Φ0eζ24σ1pΓ(1p)2πρpα1p(1+βR0)[1σD1p(ζ)+(βp+γβ+εR0)D1p1(ζ)],
(36)

where D(z) is the depth dose, z is the depth, Φ0 is the primary fluence, R0 is the range of the proton beam, σ is the standard deviation of the Gaussian distribution of the proton depth, ζ = (R0z)/σ, α and p are material-dependent constants (defined in Equation (2)), ε is the fraction of low-energy proton fluence to total proton fluence, Γ(x) is the gamma function, and Dy(x) is the parabolic cylinder function. Very good agreement is found between measured and calculated curves, and this approach has been used for the analysis and characterization of clinical proton beams (Figure 17) (Newhauser et al, 2002b) and for treatment planning dose calculations (Szymanowski and Oelfke, 2003; Soukup et al, 2005; Rethfeldt et al, 2006; Li et al, 2008).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0017.jpg

Absorbed dose D as a function of depth z in water from unmodulated (pristine) proton Bragg peaks produced by a broad proton beam with an initial energy of up to 235 MeV. With increasing depth, the accumulation of range straggling tends to broaden the peak. Beams that are more penetrating therefore have larger distal falloff distances.

3.4 Model of spread-out Bragg curves

Spread-out Bragg curves may be designed by combining multiple pristine Bragg curves. This approach was described in Wilson's seminal paper and its implementation was described by Koehler et al. (Koehler et al, 1977). Figure 15 plots several pristine Bragg curves and their resultant sum, which was similar to the SOBPs used at Harvard from the early 1970s onwards. The fluence is also modulated, which can be seen in the same figure as variations in the relative peak dose of each pristine Bragg curve. The range is typically modulated in steps that are small compared with the width of the Bragg peak, so that the SOBP contains little or no ripple. The modulation step size may be fixed or variable; generally the smallest step increments are needed for the deepest peaks, e.g., the amplitude of the ripple decreases with depth because an increasing fraction of dose is from the sub-peak region of the Bragg curves. Figure 18 plots several typical clinical SOBPs from a contemporary proton therapy system.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0018.jpg

Absorbed dose D as a function of depth z in water from various spread out Bragg peaks (SOBP) from the Northeast Proton Therapy Center (NPTC) nozzle. Open- beam central-axis depth-dose curves are plotted for each of the eight range modulator wheel tracks. The experimental values (open circles) are from charge measurements with an air-filled ionization chamber, and the model fits are shown with solid lines. Each range modulation track was designed to work over an interval of proton energies. The flat modulated-peak region was achieved by modulating the proton-beam current synchronously with the modulator wheel rotation.

SOBPs can also be modeled with analytical methods. Bortfeld and Schlegel (1996) proposed a model of the form

D(z)={34+34πln1+z^21z^+z^232πarctan2z^13:forz<za1:forzazzb0:forzzb},
(37)

where the dimensionless depth is given by

z^=zazzbza3,
(38)

and where za and zb denote the depths of the proximal and distal extents of the SOBP (see preceding section for details of their definitions) and z is the depth in water. The model was derived with several simplifying assumptions and approximations that limit its ability to predict real SOBPs. In particular, the proton energy loss rate is based on the CSDA (i.e., energy straggling is neglected) and the converted energy per unit mass approximation (where the proton energy loss and the absorbed dose in the water are assumed proportional), and non-elastic interactions of the protons with the absorber nuclei were not considered. Consequently, it does not reproduce real SOBPs well.

To attain a more realistic model, one of us (WDN) introduced terms for the finite distal-falloff distance associated with range straggling, protonic buildup effects, arbitrary slope of the modulated-peak region, and a term for the ripple in the modulated peak. In addition, the piecewise use of the function for various regions was modified to allow for a transition region between the sub-peak and modulated-peak regions, which eliminated a pronounced discontinuity there. With these additional terms and modified regions, the model becomes

D(z)=c0ΛβBUγ{c1+c24πln1+z^21z^+z^2c32πarctanc4z^1c5:forzzaΔmtz+bt:forzaΔ<z<za+Δ(mz+b)+Y:forzza+Δ},
(39)

where is defined as before, c0 is a proportionality constant, Λ takes into account the dependence on the dose due to inverse-square-law reductions in the proton fluence, βBU is the secondary-particle buildup term, c1 ~ c5 are coefficients that influence the shape of the sub-peak region, and m and b are the slope and intercept the of the modulated-peak region, respectively. The ripple in the modulated-peak region is described by the term ϒ. The distal falloff term is denoted by γ and is modeled with a cumulative normal PDF. The transition region is centered about za and is of width Δ, which is a non-critical parameter that can be estimated from the amount of range straggling expected at depth za (see previous section for detailed definition) or it may be empirically deduced from measured Bragg curves. The extended model agrees well with measured SOBPs, several of which are shown in Figure 18. This work was previously unpublished.

3.5 Analytical transport of therapeutic proton beams

Analytical models to predict therapeutic radiation are generally simple, fast, and accurate. Early treatment planning systems, many of which were still in use in the 1990s, used dose algorithms that looked up measured profiles of absorbed dose and/or used empirical formulas. Over the last 15 years, advances in our understanding of proton therapy physics, as well as huge gains in computer speed and memory capacity, have spurred the development and clinical use of increasingly realistic and complete physical models to predict absorbed dose in patients.

The basic approach underlying most analytical dose algorithms is that the proton beam's energy loss and lateral scattering may be modeled independently from one another. The physical basis for this assumption lies in the principle of conservation of momentum and energy. Lateral scattering, which is predominantly caused by MCS, occurs at small angles (see section 2.4). Hence, the energy loss associated with the vast majority of these small scattering events is negligible. Conversely, the predominant proton energy-loss mechanism, namely inelastic Coulomb interactions with atomic electrons, only negligibly deflects the protons because the proton mass is 1832 times greater than that of an electron (see section 2.1). Consequently, for an excellent approximation, the three-dimensional distribution of absorbed dose in a phantom or patient may be cast in the form of

D(x, y, z) = D(z) ⋅ OAR(x, y, z),
(40)

where x and y are lateral coordinates, z is the depth coordinate along the beam axis, D(z) is the large-field Bragg curve (i.e., with lateral protonic equilibrium, described in section 3.2), and OAR(x, y, z) is a term that takes into account the lateral properties of the dose distribution. Broad beam algorithms use Equation (40) with simple ray-casting techniques (Siddon, 1985) to determine the energy loss and range along each ray. Similarly, to determine the lateral extent of the therapeutic field, the rays are terminated in target-shaped final collimators by approximating them as black-body absorbers (Hong et al, 1996). Carefully designed broad beam algorithms are simple, fast, and accurate in media that are of relatively uniform mass density and material composition. Recently, Koch and Newhauser (2010) developed an analytical broad beam algorithm of the form

D(x,y,z)MU=i=1NDDi(zeff)MUωiSAFωiRMW(SSD+zeffSSD+z)2OARiΓi,
(41)

where N is the number of range modulation steps, ωSAF is the scatter and attenuation fluence weighting factor, ωRMW is the range modulator wheel (RMW) fluence weighting factor, OAR is off-axis ratio, and Γ is collimation modifier. This algorithm explicitly took into account edge-scattered protons and range straggling, predicted absolute dose per monitor unit (D/MU) values, and provided excellent agreement with measurements in homogeneous media, e.g., for ocular treatments. In highly heterogeneous media, however, the dose at some locations may depend strongly on the material properties of the heterogeneities in the overlying tissue, and therefore the approach of casting a single ray does not, in general, provide sufficient dosimetric accuracy.

To improve the dosimetric accuracy in heterogeneous media, the dose may be calculated at a point as the superposition of multiple pencil beams. A pencil beam is, in essence, a mathematical construct that has no exact and direct analogy in the physical world, but it is resembles a beam of infinitesimally small lateral size. By superposing enough pencil beams, the physical dose distribution can be modeled. Fortunately, as we shall see, many of the physical principles and algorithmic aspects from the broad beam approach are directly applicable to the pencil beam algorithm.

Hogstrom et al. (1981) originally proposed a pencil beam algorithm for electron therapy. Basically, an ion beam is divided into a number of finite pencil beams and the dose distribution of each pencil beam is described by the multiplication of a central-axis term and an off-axis term. The final relative absorbed dose to any point in the patient (Dp (x, y, z) ) is the summation of the dose contributions from all the pencil beams:

Dp(x,y,z)=dxdyS(x,y)PDD(zeff)(SSD+zeffSSD+z)212πσtot2(z)exp((xx)2+(yy)22σtot2(z)),
(42)

where S(x′, y′) is the relative intensity of the pencil beam at x′, y′, PDD(z) is the central-axis percentage depth dose, SSD is the source-to-surface distance, zeff is the effective depth of the point, and σtot2describes the final lateral distribution of a pencil beam by taking into account the MCS contributions from each beam-modifying device and from the patient. Many other groups have subsequently adapted and extended this concept for application in proton therapy (Petti, 1992; Hong et al, 1996; Schaffner et al, 1999; Szymanowski and Oelfke, 2002; Ciangaru et al, 2005; Schaffner, 2008; Westerly et al, 2013). Contemporary proton pencil beam algorithms provide excellent accuracy (Newhauser et al, 2007b), especially in homogeneous media, superior to that of broad beam algorithms in heterogeneous media. The improvement in accuracy comes mostly at the cost of greater computation times. Because of approximations in the pencil beam algorithm, it may not provide sufficient accuracy in extremely heterogeneous media, at material and/or density interfaces, or in other complex situations.

Currently, most proton therapy clinics utilize commercial treatment planning systems that contain pencil beam algorithms for calculation of therapeutic dose distributions. Computation speeds are generally sufficiently fast for most routine planning.

3.6 Monte Carlo transport of therapeutic proton beams

Monte Carlo transport technique has been used increasingly in radiation therapy (Seco and Verhaegen, 2013). It is more accurate than analytical models in that it takes into account the physical processes during particle transport. The most commonly used general purpose Monte Carlo codes for proton therapy research are MCNPX (Pelowitz, 2011), Geant 4 (Agostinelli et al, 2003), and FLUKA (Ferrari et al, 2005). Several groups have developed special purpose Monte Carlo codes, e.g., a treatment planning dose engine (Tourovsky et al., (2005); Yepes et al., (2009a)), typically with faster computation speed but fewer features than general purpose codes. Here we focus mainly on applications that use the MCNPX code because it is representative of most general purpose codes, we are most familiar with it of any of the codes, and constraints on journal space do not allow discussion of all codes.

The main components of a Monte Carlo simulation model are a source of protons, a treatment head, and a patient or phantom (figure 19). The proton source is typically taken as being incident on the treatment head (beamline transport is not usually simulated), where the source parameters are deduced from dose profile measurements (Newhauser et al, 2005) or taken from other beam properties, e.g., from the manufacturer (Newhauser et al, 2007b). The model of a treatment head typically includes all major mechanical components, e.g., a vacuum window, a beam profile monitor, a range modulator wheel, a second scatter, shielding plates, a range shifter assembly, backup and primary monitors, the snout and the brass collimator, and a patient-specific range compensator (Zheng et al, 2007b; Perez-Andujar et al, 2009). Models of scanned beams may model the magnetic field and transport the proton beam deflection (Peterson et al, 2009; Dowdell et al, 2012), or the proton source may be pointed, e.g., using the “thin lens” approximation to increase computation speed (Lee, 2004). The patient may be represented in voxels whose material composition and mass density are deduced from three-dimensional CT images of the patient (Taddei et al, 2009). The matrix of Hounsfield unit values in these voxels may be converted into corresponding matrices of material composition and mass density values based on a tri-linear calibration curve (Newhauser et al, 2008a). A simulation capability to accommodate time-dependent geometry was reported by Paganetti et al. (2005). Some useful applications of the Monte Carlo method are described here:

  • 1)

    Alternative or complimentary source of dosimetric data for developing, configuring, and validating analytical dose algorithms in clinical treatment planning systems. At least two studies revealed that, after a significant development effort, a Monte Carlo model of a proton therapy apparatus is sufficiently accurate and fast for prospective commissioning of an ocular treatment planning system (Koch and Newhauser, 2005; Newhauser et al, 2005) and a large-field treatment planning system (Newhauser et al, 2007b). This prospective capability reduced the time required for preclinical work by several months, resulting in substantial cost savings. Studies about scanning proton beam configuration have also been reported (Tourovsky et al, 2005; Sawakuchi et al, 2010; Dowdell et al, 2012).

  • 2)

    Alternative dose engine to calculate absorbed dose in phantoms or patient. In clinical practice, there frequently arise challenging cases where the accuracy of analytical models is unknown or suspected to be insufficient. In many such cases, measurements are not possible or practical, e.g., intracranial in vivo range verification. Multiple groups have developed Monte Carlo treatment planning engines (Tourovsky et al, 2005; Newhauser et al, 2008b; Perl et al, 2012; Mairani et al, 2013) based on various Monte Carlo codes, independently proving the technical feasibility of this approach. The physics of general purpose Monte Carlo codes for this purpose has been exhaustively validated against benchmark measurements in therapeutic proton beam data (Polf JC, 2007; Kimstrand et al, 2008; Titt et al, 2008a; Randeniya et al, 2009). Just 10 years ago, the slow execution times, memory constraints, unknown accuracy, and unknown technical feasibility of a full-blown Monte Carlo engine were serious research challenges. Today, the only major remaining obstacle to routine use of Monte Carlo engines for treatment planning is their slow execution times. For example, a proton craniospinal treatment simulation took more than 10 CPU hours (Taddei et al, 2009). Encouragingly, however, algorithmic improvements have increased speeds dramatically; Zhang et al. (2013a) found that using lattice tally in MCNPX can speed up dose simulations by one order of magnitude compared to mesh tally, especially for dose reconstructions involving large numbers of voxels. The physics of the radiation transport are identical for the lattice and mesh tallies, i.e., there was no loss of accuracy. Yepes et al. (2009b) recently reported on a fast dose calculator, a Monte Carlo track-repeating algorithm that uses a database of precomputed proton trajectories in water and applies these trajectories to heterogeneous media by scaling the proton range and scattering angles. This approach reduced the computation time for dose reconstruction in voxelized patient geometry by more than two orders of magnitude compared to MCNPX and GEANT4 (Yepes et al, 2009a; Yepes et al, 2010a). The use of low-cost alternative parallel computer architectures also appears promising, including the use of graphics processing units (Yepes et al, 2010b; Jia et al, 2012) and grid technologies (Vadapalli et al, 2011).

  • 3)

    Assessment of beam perturbation due to objects implanted in patients and development of mitigating strategies. The proportion of cancer patients with implanted devices is probably less than 10%, but it is increasing along with the aging of the population, increasing life expectancies, and increasing use of implanted devices. Examples include cardiac pacemakers, defibrillators, fiducial markers for radiotherapy localization, permanent radioactive seed implants, stents, prostheses, surgical reconstructive devices, and foreign bodies such as gunshot and shrapnel. Analytical dose algorithms in contemporary clinical treatment planning systems lack the capability to reliably assess dose perturbations from most implants. Several studies reported that implanted metallic fiducials can cause dose shadows that may compromise local tumor control for certain diseases, e.g., prostate carcinoma and uveal melanoma (Newhauser et al, 2007a; Newhauser et al, 2007c; Giebeler et al, 2009; Cheung et al, 2010; Huang et al, 2011; Matsuura et al, 2012; Zhang et al, 2012; Carnicer et al, 2013).

  • 4)

    Investigation of characteristics of therapeutic beams and beamlines. For example, edge-scattered protons can degrade the dose distribution, and the effect is difficult to model with current analytical models. Accurate predictions of D / MU require taking into account the dose from protons scattered from the edge of the patient-specific collimator, particularly for small field sizes and large depths (Titt et al, 2008b). Currently available spot scanning system only offer few options for adjusting beam spot properties like lateral and longitudinal size, and Monte Carlo simulations was used to optimize scanned beam spot (Titt et al, 2010). Because the Monte Carlo simulation technique is inherently general, it can be applied to virtually any problem involving the transport of protons or secondary radiation. Space constraints prevent us from reviewing more than a few illustrative examples. In later sections we briefly mention roles for Monte Carlo methods in research on stray radiation exposures to patients (section 5) and shielding barriers to protect staff (section 6).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0019.jpg

(Upper) Geometric model of proton therapy unit and the voxelized phantom oriented for the superior spinal proton field. The beam delivery system includes a vacuum window (A), a beam profile monitor (B), a range modulator wheel (C), a second scatterer (D), a range shifter assembly (E), backup and primary monitors (F), the snout (G), the range compensator (H), the treatment couch (I), and the patient (J). (Zhang et al, 2013a) (Lower) Simulated proton transport through the proton therapy unit and the patient.

4. Therapeutic absorbed dose determination

4.1 Reference dosimetry

By reference dosimetry, we mean the determination of absorbed dose in a manner that allows it to be directly related or referred to an accurate and uniform standard of absorbed dose. Clinical reference dosimetry comprises the measurement of absorbed dose in a clinic, which is related to the absorbed dose at a primary or secondary standards laboratory. This approach ensures that clinical reference dosimetry is accurate and uniform across participating institutions. Typically, clinical reference dosimetry is established by first calibrating a clinic's dosimeter at a standards laboratory, then “transferring” the calibration to the clinic's treatment beams. To minimize systematic errors introduced by this transfer process, both irradiations are made with the same dosimeter and under identical (or very similar) “reference conditions”. Therefore, reference conditions must be reproducible and clinically relevant. However, national or international calibration laboratories do not yet produce proton calibration beams of relevance to proton therapy (presumably due to prohibitively high costs, low demand, and limited resources). Because of this limitation, the proton therapy community developed alternative methods for proton reference dosimetry.

For many decades, the reference fields for calibrating proton dosimeters were characterized using a Faraday cup (Verhey et al, 1979) to measure the proton fluence in air. Today, most proton therapy institutions implement reference proton dosimetry utilizing an ionization chamber technique to measure the absorbed dose to water. With the latter technique, an ionization chamber is calibrated using reference conditions for photon therapy (i.e., Co radiation fields that are widely available at calibration laboratories) and a correction factor is applied that corrects the differences in the chamber's response to Co and proton beams. When implemented properly, the techniques agree within uncertainties (Newhauser et al, 2002a; Newhauser et al, 2002b).

Several advisory bodies have published dosimetry protocols for reference dosimetry, such as the American Association of Physicists in Medicine (AAPM, 1986), the European Clinical Heavy Particle Dosimetry Group (ECHED) (Vynckier et al, 1991, 1994), the International Commission on Radiation Units and Measurements (ICRU, 1998, 2007), and the International Atomic Energy Agency (IAEA, 2000). Vatnitsky et al. (1999a) perfomed an international proton dosimetry inter-comparison based on the ICRU Report 59 protocol (ICRU, 1998) and reported that absorbed dose to water can be delivered within 3%.

4.2 Patient field–specific dosimetry

The determination of absolute absorbed dose per monitor unit for each patient and treatment field, denoted by D/MU, rests upon the foundation of reference dosimetry. Recall that reference dosimetry is performed under simple and reproducible conditions, e.g., using a dosimeter in a water-box phantom. In contrast, dosimetry in a patient should take into account the full complexity of a patient's anatomy, e.g., irregular surface shapes, heterogeneous elemental composition and mass density, and treatment field size and shape. These additional complexities call for a conceptual framework and measurement techniques for performing routine D/MU determinations for patient treatment fields.

For several decades, D/MU values for passively scattered treatment beams were typically measured at depth in a water-box phantom the range compensator present. However, a systematic study by Fontenot et al. (2007) that measuring D/MU without the range compensator present provided more reliable results because use of the range compensator increased the severity of dose gradients near the calibration point and this resulted in larger overall uncertainties in D/MU. It is not known if this finding holds for special cases, e.g., small diameter treatment beams, large air gaps, and moving targets. In many practical situations, it is convenient or necessary to predict the absolute dose in a patient based on the corresponding dose under reference conditions (e.g., a water phantom). Because an exact theoretical approach is not possible, all approaches have utilized approximate methods, numerical calculations, measurements, or a combination of these. On the basis of previous studies, Newhauser et al. (2014a, b) reported the following method to calculate proton dose in a patient (D / MU)p:

(DMU)p=(DMU)wFCSPS=(DMU)wrefFOFFRSFSOBPFInvSqFFSFColSFCSPS,
(43)

where (D / MU)w is the D / MU in a water phantom under patient-specific field conditions; FCSPS is the compensator scatter and patient scatter factor, which takes into account differences in the scattering and attenuation of the patient and range compensator together relative to that of water and no range compensator; (D / MU)w is the D / MU in a water phantom under reference conditions and is 1 cGy/MU; FOF corrects for proton beam energy spectrum relative to the reference condition; FRS corrects for range shifter effect; FSOBP corrects for SOBP relative to the reference field; FInvSq corrects for beam divergence relative to the reference condition; FFS corrects for field size effect; and FColS corrects for differences in scatter from the reference aperture (10 cm × 10 cm) to the patient-specific aperture. They found that for prostate treatment fields, found that uncertainty in FCSPS, which is the least well understood factor, was clinically acceptable for prostate treatment (Newhauser et al, 2014b). In a companion study, the same group reported the variability in FCSPS was clinically significant for some lung treatment fields (Newhauser et al, 2014a).

For scanning proton beams, it is important to calibrate individual spot beam because patient specific collimating apertures or range compensators are not typically used in the beam line. First the number of protons per MU for each beam energy is calibrated by Faraday-cup measurements or Monte Carlo simulations. Then, either physical model or Monte Carlo simulation can be used to calculate absolute dose normalized to the number of incident protons. A simple nuclear halo model was developed at Paul Scherrer Institute (PSI) (Pedroni et al, 2005) to include secondary dose around the primary pencil beam and can predict precisely the dose directly from treatment planning without renormalization measurements.

In addition to the physical effects already discussed above, there are many additional effects that contribute to the uncertainty in D/MU that arise from imperfections in anatomic images used for treatment planning, organ motion, and in the case of scanned beams, the interplay of organ motion and the moving beam, etc. Such uncertainties may depend strongly on anatomical site, approximations utilized by the treatment planning system, beam delivery system, and experimenter's skill. Ideally, clinical treatment planning systems would provide accurate estimates of (D/MU)w and (D/MU)p with known and small uncertainties, methods for determining and reporting (D/MU)p would be standardized, and secondary methods (e.g., hand calculations similar to Eq. 43 or fast Monte Carlo simulations) would be available to independently verify (D/MU)p values determined with the primary method. In principal, it appears to be adequate knowledge of basic proton interaction physics and transport physics to reach or approach these clinical ideals. In addition, recent research studies report encouraging results for a variety of special cases, e.g., for a particular treatment technique, treatment planning system, or treatment delivery systems. As mentioned in section 3.5, Koch and Newhauser (2010) developed an analytical broad beam algorithm to predict absolute dose per monitor unit (D / MU) values with good accuracy in water for ocular proton treatments, but the ability to predict D / MU in the heterogeneous patient body still needs improvement. To date, the methods for determining (D/MU)p in the literature have not yet been validated for general application or standardized.

4.1 Reference dosimetry

By reference dosimetry, we mean the determination of absorbed dose in a manner that allows it to be directly related or referred to an accurate and uniform standard of absorbed dose. Clinical reference dosimetry comprises the measurement of absorbed dose in a clinic, which is related to the absorbed dose at a primary or secondary standards laboratory. This approach ensures that clinical reference dosimetry is accurate and uniform across participating institutions. Typically, clinical reference dosimetry is established by first calibrating a clinic's dosimeter at a standards laboratory, then “transferring” the calibration to the clinic's treatment beams. To minimize systematic errors introduced by this transfer process, both irradiations are made with the same dosimeter and under identical (or very similar) “reference conditions”. Therefore, reference conditions must be reproducible and clinically relevant. However, national or international calibration laboratories do not yet produce proton calibration beams of relevance to proton therapy (presumably due to prohibitively high costs, low demand, and limited resources). Because of this limitation, the proton therapy community developed alternative methods for proton reference dosimetry.

For many decades, the reference fields for calibrating proton dosimeters were characterized using a Faraday cup (Verhey et al, 1979) to measure the proton fluence in air. Today, most proton therapy institutions implement reference proton dosimetry utilizing an ionization chamber technique to measure the absorbed dose to water. With the latter technique, an ionization chamber is calibrated using reference conditions for photon therapy (i.e., Co radiation fields that are widely available at calibration laboratories) and a correction factor is applied that corrects the differences in the chamber's response to Co and proton beams. When implemented properly, the techniques agree within uncertainties (Newhauser et al, 2002a; Newhauser et al, 2002b).

Several advisory bodies have published dosimetry protocols for reference dosimetry, such as the American Association of Physicists in Medicine (AAPM, 1986), the European Clinical Heavy Particle Dosimetry Group (ECHED) (Vynckier et al, 1991, 1994), the International Commission on Radiation Units and Measurements (ICRU, 1998, 2007), and the International Atomic Energy Agency (IAEA, 2000). Vatnitsky et al. (1999a) perfomed an international proton dosimetry inter-comparison based on the ICRU Report 59 protocol (ICRU, 1998) and reported that absorbed dose to water can be delivered within 3%.

4.2 Patient field–specific dosimetry

The determination of absolute absorbed dose per monitor unit for each patient and treatment field, denoted by D/MU, rests upon the foundation of reference dosimetry. Recall that reference dosimetry is performed under simple and reproducible conditions, e.g., using a dosimeter in a water-box phantom. In contrast, dosimetry in a patient should take into account the full complexity of a patient's anatomy, e.g., irregular surface shapes, heterogeneous elemental composition and mass density, and treatment field size and shape. These additional complexities call for a conceptual framework and measurement techniques for performing routine D/MU determinations for patient treatment fields.

For several decades, D/MU values for passively scattered treatment beams were typically measured at depth in a water-box phantom the range compensator present. However, a systematic study by Fontenot et al. (2007) that measuring D/MU without the range compensator present provided more reliable results because use of the range compensator increased the severity of dose gradients near the calibration point and this resulted in larger overall uncertainties in D/MU. It is not known if this finding holds for special cases, e.g., small diameter treatment beams, large air gaps, and moving targets. In many practical situations, it is convenient or necessary to predict the absolute dose in a patient based on the corresponding dose under reference conditions (e.g., a water phantom). Because an exact theoretical approach is not possible, all approaches have utilized approximate methods, numerical calculations, measurements, or a combination of these. On the basis of previous studies, Newhauser et al. (2014a, b) reported the following method to calculate proton dose in a patient (D / MU)p:

(DMU)p=(DMU)wFCSPS=(DMU)wrefFOFFRSFSOBPFInvSqFFSFColSFCSPS,
(43)

where (D / MU)w is the D / MU in a water phantom under patient-specific field conditions; FCSPS is the compensator scatter and patient scatter factor, which takes into account differences in the scattering and attenuation of the patient and range compensator together relative to that of water and no range compensator; (D / MU)w is the D / MU in a water phantom under reference conditions and is 1 cGy/MU; FOF corrects for proton beam energy spectrum relative to the reference condition; FRS corrects for range shifter effect; FSOBP corrects for SOBP relative to the reference field; FInvSq corrects for beam divergence relative to the reference condition; FFS corrects for field size effect; and FColS corrects for differences in scatter from the reference aperture (10 cm × 10 cm) to the patient-specific aperture. They found that for prostate treatment fields, found that uncertainty in FCSPS, which is the least well understood factor, was clinically acceptable for prostate treatment (Newhauser et al, 2014b). In a companion study, the same group reported the variability in FCSPS was clinically significant for some lung treatment fields (Newhauser et al, 2014a).

For scanning proton beams, it is important to calibrate individual spot beam because patient specific collimating apertures or range compensators are not typically used in the beam line. First the number of protons per MU for each beam energy is calibrated by Faraday-cup measurements or Monte Carlo simulations. Then, either physical model or Monte Carlo simulation can be used to calculate absolute dose normalized to the number of incident protons. A simple nuclear halo model was developed at Paul Scherrer Institute (PSI) (Pedroni et al, 2005) to include secondary dose around the primary pencil beam and can predict precisely the dose directly from treatment planning without renormalization measurements.

In addition to the physical effects already discussed above, there are many additional effects that contribute to the uncertainty in D/MU that arise from imperfections in anatomic images used for treatment planning, organ motion, and in the case of scanned beams, the interplay of organ motion and the moving beam, etc. Such uncertainties may depend strongly on anatomical site, approximations utilized by the treatment planning system, beam delivery system, and experimenter's skill. Ideally, clinical treatment planning systems would provide accurate estimates of (D/MU)w and (D/MU)p with known and small uncertainties, methods for determining and reporting (D/MU)p would be standardized, and secondary methods (e.g., hand calculations similar to Eq. 43 or fast Monte Carlo simulations) would be available to independently verify (D/MU)p values determined with the primary method. In principal, it appears to be adequate knowledge of basic proton interaction physics and transport physics to reach or approach these clinical ideals. In addition, recent research studies report encouraging results for a variety of special cases, e.g., for a particular treatment technique, treatment planning system, or treatment delivery systems. As mentioned in section 3.5, Koch and Newhauser (2010) developed an analytical broad beam algorithm to predict absolute dose per monitor unit (D / MU) values with good accuracy in water for ocular proton treatments, but the ability to predict D / MU in the heterogeneous patient body still needs improvement. To date, the methods for determining (D/MU)p in the literature have not yet been validated for general application or standardized.

5. Stray Radiation

As with other forms of external-beam radiation beams, proton beams unavoidably produce stray radiation that impinges on the patient's entire body (Fig. 20). In the 1990s, very little was known about the exposure of patients to stray radiation produced by therapeutic proton beams. By the mid 2000s, concerns about stray radiation had become a matter of considerable speculation and controversy, particularly regarding the suitability of proton beams for treating children with cancer because of concerns about the risks from stray neutron exposures (Hall, 2006; Brenner and Hall, 2008). By the late 2000s, research on stray radiation exposures from proton and other radiation therapies had intensified and some questions have been partially or fully answered. However, many questions remain open. In this section, we review selected developments and discuss a few open questions.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0020.jpg

Schematic diagram of proton beam irradiation of the spine. There are several distinct sources of radiation exposure, including therapeutic protons (red), stray neutrons emanating from the treatment apparatus (blue), and neutrons produced by therapeutic proton radiation inside the body. A small-diameter beam of protons enters the treatment apparatus, which spreads the beam to a clinically useful size and collimates it to spare healthy tissues. The stray neutron is created by proton-induced nuclear reactions inside the treatment unit, some of which leak out and irradiate the patient. The stray radiation exposures provide no therapeutic benefit but increase the predicted risk that a patient will develop a radiogenic side effect, such as a second cancer, later in life. (Newhauser and Durante, 2011)

In a study using Monte Carlo simulations, Agosteo et al. (1998) reported that neutron exposures were the principal concern among the various types of stray radiation from proton therapy. Comprehensive measurements in clinical proton therapy beams were reported by Yan et al. (2002), including neutron spectral fluences and dose-equivalent data for large-field, radiosurgery, and ocular beamlines. That study included measurements of neutron energy spectra and neutron dose-equivalent values. Several other groups have measured neutron exposures, including Roy and Sandison (2004), Tayama et al. (2006), Mesoloras et al. (2006), Schneider et al. (2002) and Wroe et al. (2007). The early measurements were important because they suggested that the neutron exposures, while generally small in comparison to therapeutic doses, should not be neglected.

Several groups developed neutron dose reconstruction systems based on general purpose Monte Carlo codes for a wide variety of clinical proton beamlines (Siebers, 2000; Fontenot et al, 2005b; Jiang et al, 2005; Polf and Newhauser, 2005; Newhauser et al, 2007b; Zheng et al, 2007b; Moyers et al, 2008; Newhauser et al, 2008b). A wide variety of simulation approaches have been used, with varying degrees of clinical realism, assumptions, and approximation, and not surprisingly, rather disparate results. Some of the literature on measurements and simulations of neutron exposures was reviewed in NCRP Report 170 (2011). In the remainder of this section, we emphasize progress subsequent to that report.

Research methods of relevance to stray radiation estimation have progressed dramatically since Agosteo et al (1998) published their pioneering study using Monte Carlo simulation to study neutron fluence and dose in air or in a simple water box phantom. One key advance was the validation of Monte Carlo predictions against benchmark neutron measurements, for example the studies of Tayama et al. (2006), Fontenot et al. (2005a), and Polf et al. (2005). High-performance computing techniques, such as using massively parallel computing clusters, allows fast simulations of whole-body dose reconstructions involving complex treatment techniques, e.g., craniospinal irradiation (CSI) (Newhauser et al, 2009; Taddei et al, 2009; Taddei et al, 2010b; Perez-Andujar et al, 2013a; Zhang et al, 2013c), which is widely considered one of the most challenging treatments to simulate. Another key advance was increasing the level of realism in modeling anatomy. Methods have progressed in 1998 from simple water-box phantoms to stylized anthropomorphic phantoms to generic voxelized phantoms in the mid 2000s to patient-specific personalized phantoms by the late 2000s (i.e., by retrospectively using treatment planning CT scans). In fact, this year a study was completed comprising a 17-patient in silico clinical trial comparing proton CSI and photon CSI, according to the contemporary standards of care, including full Monte Carlo simulations of stray and therapeutic exposures (Zhang et al, 2014).

Another research result of considerable clinical relevance is that it is possible and indeed feasible to dramatically reduce exposures from stray neutrons that leak out of the treatment apparatus. For example, Taddei et al. (2008) found that modest modifications of the treatment unit (adding shields near the patient, using a tungsten-alloy collimator, and adding an upstream collimator near the range shifter assembly) substantially reduced the neutron dose (by 66%) for patients receiving passively scattered proton therapy for prostate cancer. Similarly, Brenner et al. (2009) investigated various precollimator/collimator combinations with different geometries and materials and concluded that an optimized design can be achieved to significantly reduce the stray neutron dose.

Despite the rapid growth of research on stray neutron radiation exposures to patients, the literature is mainly limited to a few case studies and anatomical sites. In general, the knowledge of stray radiation from proton therapy is still largely incomplete. The gaps in knowledge received considerable attention following the cautionary paper by Hall et al (2006), who suggested that the use of advanced radiotherapy modalities like proton therapy may not be justified because of incomplete knowledge of stray neutron exposures and second cancer induction. From subsequent studies, a coherent picture is gradually emerging: second cancer risk following proton therapy is lower than that after photon therapy for patients, the risk differential depends strongly on anatomic site and other host and treatment factors, the risk posed by stray radiation is small but not negligible, and the largest reductions in risk are achieved by using scattered- or scanned-proton beams instead of photon beams. Furthermore, the majority of risk of radiogenic second cancer is from in-field radiation, not out-of-field radiation. Consequently, the difference in predicted total risk (from both in-field and out-of-field radiation) from scanned versus scattered proton treatments is small and clinically negligible. At the risk of oversimplification, the theoretical risk advantage of proton therapy derives mainly from its ability to spare healthy tissues that are distal and lateral to the target volume. With judicious treatment techniques (e.g., beam orientations, location of field junctions, etc.), it is possible in many cases to utilize the rapid distal and lateral dose falloff of proton beams to reduce the dose to sensitive organs and tissues to below levels achievable with technologically comparable photon treatments.

It is worth underscoring that our knowledge and understanding of this topic in incomplete and likely to evolve substantially in the future. Current dose and risk assessments are mostly based on standard-of-care treatment plans, i.e., using clinical planning and research systems that not attempt to algorithmically minimize risk of radiation side effects. As clinical treatment planning systems are incrementally extended with these capabilities, the standards of care for proton and photon therapies will continue to evolve. This evolution will likely render our current understanding of radiation risks obsolete.

5.1 Monte Carlo simulations

During the transport of proton beams through the treatment unit to the patient, not only is therapeutic radiation deposited in the patient but also stray radiation is generated because of nuclear interactions (as explained in section 2.4). In proton therapy, secondary neutrons are dominant among this stray radiation and are a major concern (Agosteo et al, 1998; Fontenot et al, 2008). Estimation of stray neutron dose can be challenging and time consuming, and the variations in measurement or calculation results can be large (Xu et al, 2008).

The inherent accuracy of Monte Carlo simulation makes it an irreplaceable tool for stray neutron dose estimation, and indeed this technique has been used by many investigators studying stray neutrons generated during proton therapy (Agosteo et al, 1998; Jiang et al, 2005; Zheng et al, 2007a; Moyers et al, 2008; Zacharatou Jarlskog and Paganetti, 2008; Zheng et al, 2008; Newhauser et al, 2009; Taddei et al, 2009; Athar et al, 2010; Taddei et al, 2010a; Taddei et al, 2010b; Zhang et al, 2013b). Using MCNPX, both external neutrons (neutrons generated in the treatment unit) and internal neutrons (neutrons generated within the patient's body) can be accurately simulated (Taddei et al, 2009).

Polf and Newhauser (2005) found that neutron dose equivalent per therapeutic absorbed dose (H / D) at different locations around a passively scattered proton therapy unit increased with increasing range modulation and that the major source of neutrons shifted from the final collimator to the range modulation wheel. Zheng et al. (2007a; 2007b) reported that H/D generally increased with decreasing collimating aperture size, increasing proton beam energy, and increasing SOBP width, while it decreased with snout distance from the isocenter and increasing distance from the treatment unit. More than 50% of the neutron dose at all locations was from neutrons with energies higher than 10 MeV. Zheng et al. (2008) also reported that neutron spectral fluence contained two pronounced peaks, a low-energy peak around 1 MeV and a high-energy peak that ranged from around 10 MeV up to the proton energy. The mean neutron radiation weighting factors varied only slightly, from 8.8 to 10.3, with proton energy and location for a closed-aperture configuration.

5.2 Analytical model

Because of the complexity associated with Monte Carlo simulation and measurement, a simple analytical model to predict stray neutron dose in proton therapy is desirable. Zheng et al. (2007b) proposed an exponential decay model to predict H / D in air with good accuracy. This largely empirical model was extended by Zhang et al. (2010a) to predict H / D both in air and in the patient's body by taking into account off-axis effect, phantom attenuation effect, and low/high energy peaks observed in neutron spectral fluence. Perez-Andujar et al. (2013b) subsequently enhanced the model by replacing empirical components with physics-based components, obtaining

(HD)d=(HD)iso(ddiso)pC1eα1(ddiso)e(x2+y2)diso22σ12z2+C2eα2(ddiso)e(x2+y2)diso22σ22z2C3eα3(ddiso)e(x2+y2)diso22σ32z2+C4eα4(ddiso)e(x2+y2)diso22σ42z2,
(44)

where (H/D)iso is the neutron equivalent dose per therapeutic absorbed dose at isocenter. C1, C2, C3, and C4 apportion the H/D contributions from intranuclear cascade, evaporation, epithermal, and thermal neutrons, respectively. (d / diso) represents the power law that governs the neutron dose falloff as a function of distance from the effective neutron source, which is independent of the proton beam energy, d represents the distance from the neutron source in the treatment nozzle to the detecting volume and diso is the distance from the neutron source to isocenter. α1, α2, α3, and α4 denote the attenuation coefficients of intranuclear cascade, evaporation, epithermal, and thermal neutrons, respectively. diso is the distance from the phantom surface to the isocenter, and d’ is the distance from the phantom surface to the detecting volumes. The lateral distribution of neutrons is governed by σ1 for intranuclear cascade, σ2 for evaporation, σ3 for epithermal, and σ4 for thermal neutrons. z is the vertical coordinate for the neutron dose receptor and is used to scale the width parameters.

The accuracy of this model is comparable to the accuracy of typical Monte Carlo simulations or measurements of neutron dose (Figure 21). Taddei et al. (2013) modified the model reported by Zhang et al. (2010a) and used a simplified double-Gaussian model to calculation out-of-field dose in photon therapy and they also had good agreement between measurement data and model-based values.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0021.jpg

Results from Monte Carlo simulations and analytical model predictions of equivalent per therapeutic dose (H/D) for proton beams with various energies in the water phantom as a function of vertical depth (a) and lateral distance (b) (Perez-Andujar et al, 2013b).

This model and the follow-up research from our group (to extend the model for multiple energies and to implement this model in commercial treatment planning systems) push the possibility of calculating the plan-specific out-of-field dose delivered to the patient during both proton and photon therapy.

5.1 Monte Carlo simulations

During the transport of proton beams through the treatment unit to the patient, not only is therapeutic radiation deposited in the patient but also stray radiation is generated because of nuclear interactions (as explained in section 2.4). In proton therapy, secondary neutrons are dominant among this stray radiation and are a major concern (Agosteo et al, 1998; Fontenot et al, 2008). Estimation of stray neutron dose can be challenging and time consuming, and the variations in measurement or calculation results can be large (Xu et al, 2008).

The inherent accuracy of Monte Carlo simulation makes it an irreplaceable tool for stray neutron dose estimation, and indeed this technique has been used by many investigators studying stray neutrons generated during proton therapy (Agosteo et al, 1998; Jiang et al, 2005; Zheng et al, 2007a; Moyers et al, 2008; Zacharatou Jarlskog and Paganetti, 2008; Zheng et al, 2008; Newhauser et al, 2009; Taddei et al, 2009; Athar et al, 2010; Taddei et al, 2010a; Taddei et al, 2010b; Zhang et al, 2013b). Using MCNPX, both external neutrons (neutrons generated in the treatment unit) and internal neutrons (neutrons generated within the patient's body) can be accurately simulated (Taddei et al, 2009).

Polf and Newhauser (2005) found that neutron dose equivalent per therapeutic absorbed dose (H / D) at different locations around a passively scattered proton therapy unit increased with increasing range modulation and that the major source of neutrons shifted from the final collimator to the range modulation wheel. Zheng et al. (2007a; 2007b) reported that H/D generally increased with decreasing collimating aperture size, increasing proton beam energy, and increasing SOBP width, while it decreased with snout distance from the isocenter and increasing distance from the treatment unit. More than 50% of the neutron dose at all locations was from neutrons with energies higher than 10 MeV. Zheng et al. (2008) also reported that neutron spectral fluence contained two pronounced peaks, a low-energy peak around 1 MeV and a high-energy peak that ranged from around 10 MeV up to the proton energy. The mean neutron radiation weighting factors varied only slightly, from 8.8 to 10.3, with proton energy and location for a closed-aperture configuration.

5.2 Analytical model

Because of the complexity associated with Monte Carlo simulation and measurement, a simple analytical model to predict stray neutron dose in proton therapy is desirable. Zheng et al. (2007b) proposed an exponential decay model to predict H / D in air with good accuracy. This largely empirical model was extended by Zhang et al. (2010a) to predict H / D both in air and in the patient's body by taking into account off-axis effect, phantom attenuation effect, and low/high energy peaks observed in neutron spectral fluence. Perez-Andujar et al. (2013b) subsequently enhanced the model by replacing empirical components with physics-based components, obtaining

(HD)d=(HD)iso(ddiso)pC1eα1(ddiso)e(x2+y2)diso22σ12z2+C2eα2(ddiso)e(x2+y2)diso22σ22z2C3eα3(ddiso)e(x2+y2)diso22σ32z2+C4eα4(ddiso)e(x2+y2)diso22σ42z2,
(44)

where (H/D)iso is the neutron equivalent dose per therapeutic absorbed dose at isocenter. C1, C2, C3, and C4 apportion the H/D contributions from intranuclear cascade, evaporation, epithermal, and thermal neutrons, respectively. (d / diso) represents the power law that governs the neutron dose falloff as a function of distance from the effective neutron source, which is independent of the proton beam energy, d represents the distance from the neutron source in the treatment nozzle to the detecting volume and diso is the distance from the neutron source to isocenter. α1, α2, α3, and α4 denote the attenuation coefficients of intranuclear cascade, evaporation, epithermal, and thermal neutrons, respectively. diso is the distance from the phantom surface to the isocenter, and d’ is the distance from the phantom surface to the detecting volumes. The lateral distribution of neutrons is governed by σ1 for intranuclear cascade, σ2 for evaporation, σ3 for epithermal, and σ4 for thermal neutrons. z is the vertical coordinate for the neutron dose receptor and is used to scale the width parameters.

The accuracy of this model is comparable to the accuracy of typical Monte Carlo simulations or measurements of neutron dose (Figure 21). Taddei et al. (2013) modified the model reported by Zhang et al. (2010a) and used a simplified double-Gaussian model to calculation out-of-field dose in photon therapy and they also had good agreement between measurement data and model-based values.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0021.jpg

Results from Monte Carlo simulations and analytical model predictions of equivalent per therapeutic dose (H/D) for proton beams with various energies in the water phantom as a function of vertical depth (a) and lateral distance (b) (Perez-Andujar et al, 2013b).

This model and the follow-up research from our group (to extend the model for multiple energies and to implement this model in commercial treatment planning systems) push the possibility of calculating the plan-specific out-of-field dose delivered to the patient during both proton and photon therapy.

6. Shielding Design

In this section, we will briefly review the design of shielding to protect humans from stray radiation. In particular, we will focus on the bulk shielding barriers (e.g., walls, ceilings) and mazes. Bulk shielding is an extremely important aspect of facility design because proton therapy equipment is capable of producing lethal levels of stray radiation. In addition, bulk shielding intrudes on space available for equipment and personnel and it is expensive. In many ways, the design of bulk shielding is a classical engineering problem; one develops a solution that comprises an acceptable balance of the competing attributes of safety, utility, and cost.

Broadly, to design shielding one must have knowledge of the stray radiation production, transport, and attenuation in shielding barriers. In addition, the shielding design goals, i.e., the predicted exposure rate at an occupied location, depend on the fraction of time an area is occupied, its designation as a controlled or uncontrolled area, and the type of individuals present, i.e., patients, staff (radiation workers) and the general public. There are typically multiple design goals to be satisfied, e.g., one for short term exposures (averaged over one hour) and one for long term exposures (averaged over one year). In general, one does not know a priori which design goal will ultimately govern the shielding design, necessitating shielding calculations that represent multiple scenarios.

Before discussing the basic physics of shielding design, we digress briefly to point out that shielding barriers are one of several necessary components of an effective radiation protection program. Shielding barriers alone do not ensure safety. It is impractical, not to mention prohibitively expensive, to build shields that by themselves provide adequate protection for unrestricted usage of contemporary proton therapy systems. Administrative controls on beam usage are therefore necessary and, for obvious reasons, the anticipated beam usage must be taken into account in the shielding design process. Knowledge of the proton beam usage is often difficult to predict, particularly for first-of-kind treatment units where neutron production in the beam production and delivery systems are poorly known. In addition, the future outcome of a clinical trial or a change in health care policy may dramatically change proton beam usage. Leaving these uncertainties aside, one needs to know not only the beam usage (energies, charges and currents, and directions of travel) at each of the treatment locations, but also all of the corresponding proton loss rates in the accelerator, energy selection system, beam transport system, and treatment head. The loss proton losses and neutron production vary strongly on the proton beam energy and other factors (see Fig. 22), necessitating separate calculations for several proton beam energies (Newhauser et al, 2002d). Many aspects of shielding of proton therapy facilities are similar to those for photon therapy facilities, which have be reviewed in detail elsewhere (NCRP, 2008).

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0022.jpg

(Upper) The proton therapy facility treatment level: 235-MeV cyclotron (1), variable- thickness graphite energy degrader (2), momentum analysis magnets (3), slits (4), brass aperture (5), and beamstop at isocenter (6). Plan also shows the main control room (b), treatment room maze exit (m, n, o, a), and various corridors and occupied rooms on the level above (c–k). (Newhauser et al, 2002c) (Lower) Calculated neutron dose equivalent spectra at locations l–o in the gantry room and its maze (see Fig. 1) produced by a 235-MeV proton beam. The ordinate corresponds to the product of the neutron fluence-to-dose-equivalent conversion coefficient hΦ, the neutron spectral fluence ΦE, and the neutron energy En, where the product is normalized to the total neutron dose equivalent H. These spectra reveal differences in the shape due to the relative contributions from peaks due to thermal neutrons, evaporation neutrons, and cascade neutrons. The region between 10 and 10 MeV, corresponding to 1/En behavior of the spectral fluence, contributes relatively little to the total dose equivalent. (Newhauser et al, 2002c)

On the assumption that the beam usage, occupancy factors, and other aspects are known, one may calculate the neutron dose equivalent rate at a point of interest behind a slab shielding barrier (e.g., wall, ceiling, etc.) that is parallel to and at a distance a the proton beam according to

H(Ep,θ,dλ(θ))=Ho(Ep,θ)r2exp[dλ(θ)cos(α)],
(45)

where H(Ep, θ, d/λ) is the ambient dose equivalent beyond the shield, Ep is the proton beam energy, θ is the emission angle from the proton loss point and the point of interest, r is the distance between the proton loss point and the point of interest, d is the thickness of the shield, λ(θ) is the attenuation length of the shield material, Ho(Ep, θ) is the source term, α is the angle of incidence of the radiation impinging on the shield (Agosteo, 2009). This semi-empirical approach, based on the work of Moyer (Moyer, 1962), has many attractive attributes, including the inclusion of major dependencies, simplicity, and computational speed. Five decades on, in spite of its limitations, the Moyer model still is widely used, especially in cases where speed and convenience are more important than accuracy.

Personnel must be able to ingress to and egress from access shielded vaults quickly, e.g., for efficient routine clinical operation, in emergencies such as attending to a sick patient, and to maintain and repair the accelerator and beam transport system. This necessitates large openings in the bulk shielding walls. Obviously such openings result in a reduction in attenuation that must be restored. Typically the required attenuation is restored by fitting the opening with a long and curved shielding tunnel, commonly referred to as a maze or labyrinth. The maze attenuates radiation that is incident upon it in two ways: its walls attenuate the deeply penetrating radiation that is incident on the maze walls, and it attenuates the less-penetrating radiation that enters the opening of the maze by a combination of scattering and absorption processes as the radiation propagates through of the maze.

Typically the mazes of treatment rooms are not equipped with massive shielded doors because they can significantly increase access times. On the other hand, infrequently accessed vaults (e.g., for the accelerator and beam transport system) commonly utilize a maze equipped with a massive shielding door.

The dose equivalent rate outside the maze (without a door) from scattered radiation that enters the maze entrance (inside the vault) is given by

H(r1,r2,r3,,rn)=2a2Ho(a)r12i=2nfi(ri),
(46)

where r1 is the line of sight distance from the source to the point of interest in the first (closest to the source) leg of the maze maze, Ho(a) is the dose equivalent at distance a from the source, a is distance from the source to the entrance of the first leg, and fi(ri) are attenuation factors for the subsequent legs of the maze (Tesch, 1982). These are given by

f(ri)=eri0.45+0.022Ai1.3eri2.351+0.022Ai1.3
(47)

where ri is the distance from entrance of the i leg to the point of interest in that leg, and where Ai is the cross sectional area of the i leg of the maze (Agosteo, 2009). In addition, one must still (separately) consider attenuation of radiation impinging upon the outside of the maze walls (e.g., using the Moyer model) and the attenuation in maze doors, in cases where they are used. The attenuation of the maze increases with the number and length of the legs, the thickness of the maze walls, and decreases with cross sectional area. It would be difficult to overemphasize the high relative importance of the design of the mazes. Fortunately, there are available several design methods that are well understood, accurate, and that have been validated with measurements in clinical proton therapy facilities.

Materials used for the bulk shielding barriers and mazes vary somewhat, depending on the cost and availability of the shielding materials and the cost of space occupied by the shielding barriers. Typically, ordinary concrete with steel reinforcement is utilized because of its large hydrogen content and mass density (2.35 g/cm), high strength, good availability, and low cost. Contemporary proton therapy centers use concrete shielding barriers of up to several meters thick (Newhauser et al., 2002; Titt and Newhauser, 2005). The neutron attenuation length, or the thickness required to reduce the incident fluence by a factor of 1/e, is plotted in Figure 23 for ordinary concrete. At the high-energy limit, the attenuation length is approximately 50 cm of ordinary concrete (Moritz, 2001). The density and attenuating properties of concrete can be increased substantially by utilizing high-density aggregate material, but cost prohibitive in most cases. Some facilities have utilized comparatively small “local” shields, e.g., to shield beamline components that produce copious quantities of neutrons, to lessen the attenuation requirements of larger bulk shields. In some cases, higher density materials such as alloys of iron (7.9 g/cm) or tungsten (19.25 g/cm) may be advantageous for shielding against high-energy neutrons because they can be made more compact than with concrete or other bulk shielding materials. Proton beam stoppers, which prevent the primary proton beam from impinging on the bulk shielding, have been made from a variety of materials, including plastic and steel.

An external file that holds a picture, illustration, etc.
Object name is nihms-677175-f0023.jpg

The attenuation length of neutrons in concrete versus neutron energy (Moritz, 2001).

The shielding design methods most often used for proton therapy facilities are analytical methods and Monte Carlo simulations. The knowledge of the uncertainties in predictions from both of these methods is incomplete. Newhauser et al. (2002c) measured, calculated (using semi-empirical analytical models that were developed for slab and maze shielding geometrics), and simulated neutron dose rates at a 235-MeV proton therapy center, and they found that the analytical model overestimated neutron dose at most positions compared to measurements, while Monte Carlo simulations agreed more closely with the measurements. However, the Monte Carlo method is considerably more challenging because of the need to obtain convergence of solutions; generalized and automated variance reduction techniques are lacking. Titt and Newhauser (2005) evaluated the uncertainty associated with Monte Carlo shielding technique by comparing the analytical predictions with detailed Monte Carlo simulations of neutron equivalent doses at various receptor locations. They found the optimum rejection criterion to be 10% statistical uncertainty for the Monte Carlo simulations. This rejection criterion provided additional confidence because virtually all accepted simulated results had converged.

In an unpublished study, one of us (WN) estimated the cost of concrete and steel used for shielding a typical multi-room proton therapy at approximately 6M USD, with the potential to reduce shielding costs by almost one third through improved neutron shielding designs. Evidently there is considerable potential to achieve cost savings and other benefits by developing improved shielding design methods and to optimize shield designs.

Although we have not discussed small local shields of beamline components in this section, there is considerable overlap with the material presented in the previous section, for brevity we shall not consider stray radiation inside of a treatment room here.

7. Challenges and Future of Proton Therapy

The future utilization of proton therapy is difficult to predict. There are some tentative indications that it will continue to become more widely available in developed countries, perhaps owing to the clear theoretical dosimetric advantages associated with proton beam dose distributions. Perhaps the strongest argument in favor of protons is actually an indirect counter argument: there is only sparse evidence of beneficial effect from unnecessary irradiation of healthy tissues (i.e., from the exit dose of photon beams) (Terasawa et al, 2009). Other arguments in favor include an excellent record of treatment safety and efficacy. Furthermore, over the past 15 years, many arguments against proton therapy have been largely or fully disproved, e.g., it is too complex and too difficult for most clinical organizations, manufacturers cannot be counted upon to deliver systems on schedule, the range uncertainty is too large, and the risks from stray neutron doses are too great.

Nonetheless, skepticism and controversy remain regarding the likely ultimate role of proton therapy in radiation oncology. Especially in the last few years, the debate seems to focus on cost-effectiveness and cost-competitiveness. Basically, the argument goes, the cost and value of proton therapy has not been proven with evidence of improved patient outcomes, which are presumed to offset some or all of the higher costs of proton therapy systems. If the price differential between proton and photon therapies were to substantially shrink or disappear, e.g., due to economies of scale, many clinics would replace at least some photon treatment units with proton units.

For all these reasons, there are many urgent research questions of relevance to proton therapy. Some of these are enumerated here, with an emphasis on those questions that will require physics research and development to reduce cost, improve treatment quality and efficiency, and create previously new treatment capabilities of clinical importance.

  • 1)

    Can novel techniques, such as proton arc therapy (Sandison et al, 1997; Sengbusch et al, 2009; Rechner et al, 2012), be developed to improve the quality of treatment, reduce treatment time, and increase cost-competitiveness and -effectiveness?

  • 2)

    Can cost-competitiveness or treatment capability be increased significantly through incremental improvements to existing accelerator technologies, e.g., fixed-field alternating gradient synchrotron (Johnstone et al, 1999) and superconducting cyclotron accelerators (Blosser et al, 1997), or novel linear accelerators, e.g., wakefield laser accelerators (Schwoerer et al, 2006) and dielectric wall accelerators (Caporaso et al, 2008)?

  • 3)

    Can ultra-compact low-energy proton accelerators provide hand-held or robotically mounted radiation sources for treating superficial tumors or intra-operative treatment of deep-seated tumors?

  • 4)

    How can range uncertainties be quantified and minimized? In cases where uncertainties are presently too large, what roles will prospective imaging play, including proton CT (Schulte et al, 2004) and megavoltage photon CT (Langen et al, 2005; Newhauser et al, 2008a; De Marziet al, 2013)? Similarly, what role will real-time or post-treatment imaging play, including prompt gamma imaging (Peterson et al, 2010; Gueth et al, 2013), positron emission tomography (Parodi et al, 2007; Cho et al, 2013; Min et al, 2013), proton radiography (Schneider and Pedroni, 1995), and magnetic resonance imaging (Krejcarek et al, 2007)?

  • 5)

    In the future, what out-of-field dose algorithms should be developed for treatment planning systems used for research or routine clinical practice? Techniques for the calculation and measurement of therapeutic dose are reasonably well established (with a few exceptions mentioned below). In most cases, intensity-modulation techniques have enabled the clinician to provide essentially identical dose distributions to the target volumes using either proton or photon therapy. Consequently, it appears that comparative treatment planning studies of the future will mainly focus on the radiation exposure of healthy tissues and on the feasibility of delivering novel treatments. Knowledge of the physics of stray radiation exposures is incomplete, and most commercial treatment planning systems neglect or very crudely approximate the out-of-field dose.

  • 6)

    What are the optimal strategies for treating moving tumors, especially in the thorax and abdomen? For example, will scanned beam treatments with target tracking be beneficial? How should organ motion be measured during treatment? Do the benefits of incremental margin shrinkage justify increased risk of interplay dose defects? Are there cases where passively scattered beams are superior?

  • 7)

    As the use of on-board imaging is expected to increase, how will additional radiation exposures to patients be monitored and minimized? Are the radiogenic risks of future imaging schemes justified by their benefit to the patient?

  • 8)

    What is the role of future multimodality radiation therapies that include protons? The objectives of such protocols include reducing skin dose, sharpening penumbral widths and distal falloff lengths, and reducing dose delivery errors. Important sources of errors include temporal interplay of the organ motion and beam location and range uncertainties due to image artifacts caused by metal implants. Possible implementations could include a source of protons and other particles, e.g., photons, on a single rotational gantry.

  • 9)

    Is there a need for a national or international primary standards laboratory to provide reference proton beams to support a standard for dose absorbed to water?

In the preceding 15 years or so, much progress has been made toward a complete understanding of the physics of proton therapy. In particular, today there exist analytical models and simulation techniques to design and model most of the major aspects of clinical proton therapy systems currently in operation, several of which we have reviewed in this paper. However, even with today's state of the art knowledge of proton therapy, the answers to many questions of scientific interest and economic importance remain tantalizingly beyond the reach of current research capabilities. As in the past, clinical needs and economic forces will likely define many of the future research frontiers in proton therapy.

To fully exploit the advantages of proton beams to improve patient outcomes, it is clear that additional research is need to optimize the current generation of proton therapy systems, to make new discoveries and translate breakthroughs into novel prototype research systems, to obtain a deeper and clinically applicable understanding of the relevant aspects of radiation biology, to improve the efficiency clinical operations (e.g., industrial and process engineering), and to generate theoretical and observational evidence to assess the comparative effectiveness and cost of proton therapy relative to other comparable treatments for a wide variety of diseases, anatomical sites, and outcome endpoints. The relative priority of these goals is a matter of subjective judgment and speculation, which we shall leave to the reader.

Much of the important research will require experts and specialists from disciplines such as accelerator physics, imaging physics, various engineering specialties, oncology, biology, biostatistics and epidemiology, and computer science. Many basic and applied research studies are well suited to purely academic or clinical research environments and research teams. Other research studies will benefit from the formation of collaborative teams comprising members drawn from the academy, medicine, and industry.

Acknowledgements

This work was supported in part by the National Cancer Institute (award 1R01 CA131463-01A1), Department of Defense (award W81XWH-10-1-0170), and the Bella Bowman Foundation.

Medical Physics Program, Department of Physics and Astronomy, Louisiana State University, 202 Nicholson Hall, Baton Rouge, LA, 70803, USA
Mary Bird Perkins Cancer Center, 4950 Essen Lane, Baton Rouge, LA, 70809, USA
Corresponding author: Department of Physics and Astronomy, Louisiana State University, Baton Rouge, LA, USA Tel.: 225-578-2762; Fax: 225-215-1376 ude.usl@resuahwen (W. Newhauser).

Abstract

The physics of proton therapy has advanced considerably since it was proposed in 1946. Today analytical equations and numerical simulation methods are available to predict and characterize many aspects of proton therapy. This article reviews the basic aspects of the physics of proton therapy, including proton interaction mechanisms, proton transport calculations, the determination of dose from therapeutic and stray radiations, and shielding design. The article discusses underlying processes as well as selected practical experimental and theoretical methods. We conclude by briefly speculating on possible future areas of research of relevance to the physics of proton therapy.

Keywords: proton therapy, interaction mechanism, transport, Monte Carlo, dosimetry, shielding
Abstract
Collaboration tool especially designed for Life Science professionals.Drag-and-drop any entity to your messages.