Simulation of depth-dose curves and water equivalent ratios of energetic proton beams in cortical bone
Ana María Zamora-Vinaroz Pablo de Vera Isabel Abril Rafael Garcia-Molina


Simulated proton depth-dose curves and WER in cortical bone]Simulation of depth-dose curves and water equivalent ratios of energetic proton beams in cortical bone




^1Departamento de Física – Centro de Investigación en Óptica y Nanofísica, Universidad de Murcia, 30100 Murcia, Spain
^2Departament de Física Aplicada, Universitat d’Alacant, 03080 Alacant, Spain
pablo.vera@um.es




  April 29, 2024


We have determined the depth-dose curve, the penetration range, and the water equivalent ratio (WER), for proton beams of clinical energies in cortical bone, by means of a detailed and accurate simulation that combines molecular dynamics and Monte Carlo techniques. The fundamental input quantities (stopping power and energy loss straggling) for the simulation were obtained from 

a reliable electronic excitation spectrum of the condensed-phase target, which takes into account the organic and mineral phases that form it. Our simulations with these inputs, that are in excellent agreement with the scarce data available for a cortical bone target, deviate from simulations performed using other stopping quantities, such as those provided in the widely used ICRU Report 49.  The results of this work emphasize the importance of an accurate determination of the stopping quantities of cortical bone in order to advance towards the millimetric precision for the proton penetration ranges and deposited dose needed in radiotherapy.





















§ INTRODUCTION
 

Nowadays, hadrontherapy  and, particularly, protontherapy, are being promoted because of their specific advantages in comparison to conventional radiotherapy. Hadrontherapy can be used to treat deep-seated tumours without damaging the surrounding healthy tissues or the surface tissues as happens with photon treatments. This is due to an increase of the radiobiological effectiveness, as well as the fact that heavy charged particles deposit their energy in the so-called Bragg peak, a narrow region at the end of their trajectories, where a pronounced dose peak is produced <cit.>.

Treatment planning in hadrontherapy includes the simulation of therapeutic ion beams impinging on human tissues with energies of the order of tens and hundreds of MeV/u and a dose assessment. The latter is performed, for simplicity, by measuring the absorbed dose and the beam range in a patient phantom made of water, which is the main constituent of the human body. Finally, the water equivalent ratio (WER), as a nice tissue equivalence estimator <cit.>, allows to compute the beam range in biological tissues <cit.>. 

As the Bragg peak is narrow, the possible uncertainties in range, which will exist in WER as well, could affect the treatment effectiveness and lead to a damage of the surrounding healthy tissues. Therefore, it is necessary to know accurately the energy deposition by swift protons in the tissues that will be sampled. This means that the electronic excitation spectrum of the target must be properly described taking into account the condensed-phase nature of the material. 

For the purpose of the present study, the biological materials are roughly classified between soft and hard (or skeletal) tissues <cit.>. 
On the one hand, the former are usually assumed to be made of  liquid water, as it has been referred to be an excellent tissue-like phantom material for determination of absorbed dose <cit.>. This material has been studied in more detail, and even though the amount of experimental data on its stopping power for protons is rather limited <cit.>,  there are several theoretical models capable of reproducing their main trends <cit.>. For other organic materials making up soft tissues, the stopping power (as well as other relevant energy loss quantities) can be calculated by a parameterisation of their excitation spectrum <cit.>.

On the other hand, there is less information available for the stopping power of heavy charged particles in hard tissues, specially in bone. As there are very few experimental measurements in hard tissues <cit.>, the main way to compute its stopping power is the additivity of atomic stopping powers using Bragg's rule <cit.>. However, 42% of bone is collagen protein and the other 58% is calcium hydroxyapatite (HAp), whose stochastic formula is Ca_10(PO_4)_6(OH)_2. The stopping power of HAp has been measured experimentally and compared with simulations based on the dielectric formalism <cit.>, finding a nice agreement between theory and experimental data. Therefore, according to Bragg's rule, the stopping power of cortical bone can be computed by the weighted sum of the stopping powers of its mineral (HAp) and organic components, taking into account its condensed-phase nature.

The aim of the present work is to study the main quantities related to the energy deposition by protons, having clinical energies, in both liquid water and cortical bone. From Monte Carlo simulations with the SEICS code <cit.>, we get the depth-dose curves, the penetration range and the WER in cortical bone, in order to compare them with the corresponding fot liquid water. To account for the condensed-phase nature of the studied materials, we use the MELF-GOS method (based on the dielectric formalism) to calculate the ELF of the target as well as a multi-phase approach for cortical bone <cit.>. This leads to stopping power values showing differences of the order of 10–15% between our calculations (based in the dielectric formalism and the MELF-GOS methodology) and ICRU tabulation <cit.> below 200 keV for water and 1000 keV for bone. Despite the scarcity of experimental data of depth-dose curves in bone, our depth-dose simulation is in excellent agreement with Vasiliev et al.'s experimental measurements <cit.> , even better than using ICRU data. In addition, our simulations of WER in cortical bone show small differences between using the MELF-GOS stopping power or the ICRU one. In spite of a fair agreement with recently reported WER values by Burin et al. <cit.>, the differences between both datasets are more visible.

This work is organized in the following sections. The methodological framework is  introduced in Section 2, including the calculation of the stopping quantities (Section 2.1), obtaining the energy loss function (ELF) of hydroxyapatite and cortical bone (Section 2.2), as well as a brief description of the simulation code (Section 2.3). The results found for the stopping quantities, the depth-dose curves and ranges for protons of several incident energies in both liquid water and cortical bone, as well as the WER values for bone are presented in Section 3. Finally, the main conclusions of this paper are summarized in Section 4.




    
    
 
    
    
    
    
    





§ THEORY AND SIMULATION METHODS




 §.§ Calculation of the stopping quantities


The dielectric formalism  <cit.> constitutes an appropiate framework to study the energy loss of fast protons in condensed-phase materials, as it includes the energy transfer by the projectile to the electronic excitations of the target accounting for its condensed-phase nature. 



Within this formalism, a projectile with mass M, atomic number Z_1, kinetic energy T and charge state q travelling through a medium, characterised by a dielectric function ϵ(k, ω), will be stochastically slowed down with a stopping power S_q (mean energy loss per unit path length) given by:

    S_q (T) = Me^2/πħ^2 T∫_E_-^E_+dE   E  ∫_k_-^k_+dk/k [Z_1-ρ_q(k)]^2 [ -1/ϵ(k, E)]

where E and ħ k represent the energy and momentum of the target electronic excitations, respectively. 


The stochastic nature of the energy loss processes results in a distribution of values for the projectile energy loss processes, which is characterised by the so called energy loss straggling Ω_q^2: 

    Ω_q^2 (T) = Me^2/ħ^2 π T∫_E_-^E_+d E   E^2 ∫_k_-^k_+dk/k [Z_1-ρ_q(k)]^2 [ -1/ϵ(k, E)]

In both equations eq:Sq and eq:Stg_q, e denotes the elemental charge, whereas the electronic charge distribution of the projectiles is characterized by the Fourier transform of its electronic density, ρ_q (k). 

The integration limits E_+ and E_-, as well as k_+ and k_-, are obtained by the fulfillment of the laws of conservation of energy and momentum <cit.>.

Moreover, these formulae include the so-called energy-loss function (ELF) Im[ -1/ϵ(k,E) ]. This quantity accounts for the probability of producing an excitation or ionization in the target with a certain energy and momentum, and it is calculated by using the MELF-GOS (Mermin Energy Loss Function – Generalised Oscillator Strengths) method <cit.>, as briefly described in what follows. 




The stopping power S and the energy loss straggling Ω^2 can be calculated as the weighted sum over S_q and Ω_q^2 for the different charge states q of the projectile as it travels through the target capturing and losing electrons; ϕ_q represents the probability that the projectile has a charge state q <cit.>:

    S(T) = ∑_q=0^Z_1ϕ_q (T) S_q(T)     Ω^2 (T) = ∑_q=0^Z_1ϕ_q (T) Ω_q^2 (T)


When the projectile produces an excitation in the target, the response of the target electrons will be different depending on their binding energies, so the ELF is split in two terms:

    Im[ -1/ϵ(k,E)] =  Im[ -1/ϵ(k,E)]_out + ∑_j ν_j ∑_nl Im[ -1/ϵ(k,E)]_nl^j

The first term is due to the target outer-shell electrons, while the second one comes from the inner-shell electrons. 
In the last term, ν_j specifies the stoichiometry of the jth element in a compound target, while the indexes nl refer to the (n,l) subshells of the inner electrons <cit.>.

On the one hand, inner-shell electrons preserve their atomic character because of their large binding energies. Therefore, this kind of excitations are taken into account by using hydrogenic generalised oscillator strengths (GOS) for the K-, L- and M-shells to obtain the contribution to the ELF of the inner-shell electrons <cit.>:

    Im[ -1/ϵ(k,E)]_nl^j = 2π^2 ħ^2 e^2  N/m E  df_nl^j(k, E)/ d E Θ(E - E_ th,nl^j)

where N is the molecular density of the target, m is the electron mass and f_nl^j (k,E)/ E is the hydrogenic GOS for the (n,l)-subshell of the jth atomic component, whose energy of ionization is E_ th,nl^j. 

Due to the weak binding energy of the outer electrons, their excitations reflect the condensed-phase nature of the target and are better described from available data at the optical limit (k=0), using a weighted-sum of Mermin-type ELFs <cit.>, as described in Ref. <cit.>: 

    Im[ -1/ϵ(k=0,E)]_ out    =    ∑_i A_i/E_i^2 Im[ -1/ϵ(E_i,γ_i;k=0,ħω)]_i
        =    ∑_i Θ(E - E_ th,i)  A_i   E  γ_i/(E^2-E_i^2)^2+(E γ_i)^2

Here, A_i, E_i and γ_i are fitting parameters related with the intensities, positions and widths of the peaks appearing in the optical ELF. The function Θ(E - E_ th,i) denotes the step function of the transferred energy ħω, whose threshold energy is E_ th,i. 



With the MELF-GOS method, the excitation spectrum of the condensed-phase target is appropriately described, automatically extending it to finite momentum transfers by means of the analytical properties of the Mermin dielectric function and the GOS, without using ad hoc hypothesis for the k-dependence of the ELF. Actually, the MELF-GOS method has reproduced nicely <cit.> the experimental ELF of liquid water for k ≥ 0 <cit.>. 



 §.§ ELF of cortical bone


In order to characterize the targets studied in this work, we apply the MELF-GOS method to liquid water and cortical bone. The former has been obtained in previous works by our group <cit.>, by a fitting of the outer ELF to experimental data <cit.> and obtaining the inner-shell contribution to the ELF from the GOS of the oxygen K-shell <cit.>. 

Cortical bone, having a more complex composition <cit.>, can be divided into a mineral contribution and an organic contribution <cit.>. Therefore, to account for the inner-shells contribution to the ELF, we have to take into account the inner-shell electrons of all elements involved in both components. The mineral part, which is called calcium hydroxyapatite (HAp), contributes with inner-shell electrons from the K- and L-shells of calcium and phosphorus, as well as from the K-shell of oxygen <cit.>. The last ones are also present in the inner-shells of the organic part of cortical bone, as well as the K- and L-shell electrons of phosphorus and sulfur, and the K-shell electrons of carbon and nitrogen. 

For the outer-shell electrons, we need to fit Im [-1/ϵ (k, E)]_out to the optical ELF of HAp or of organic bone. Unfortunately, there are no available experimental data for the former. However, ab initio calculations of the optical properties of HAp were performed by applying density functional theory (DFT) within the local-density approximation <cit.> and fitted using equation eq:ELFout, as explained in Ref. <cit.>. 



The mass composition of cortical bone is 58% of HAp and 42% of organic material <cit.>. Since the ELF of the former component is already known <cit.>, it is needed to know the ELF of the second contribution. As no experimental measurements are available of the latter constituent of cortical bone, its optical-ELF can be described by means of a single-Drude function, whose parameters depend on the mean atomic number and the density of the target <cit.>. This procedure has been applied to several organic targets <cit.>, obtaining  a nice agreement with available experimental measurements. 

Therefore, to know the optical ELF of the organic part of cortical bone the composition of this material is needed. This composition is obtained assuming that all the calcium was in the mineral part, so the HAp composition (with density ρ = 1.85 g/cm^3) was subtracted from the cortical bone composition provided in Ref. <cit.>. The density of the organic part of cortical bone was calculated taking into account the densities and mass compositions of both cortical bone and HAp <cit.>. Useful information about each condensed-phase target is shown in tab-comp and tab-prop-part. The last two columns of tab-prop-part gather the mean excitation energies for the different target materials, as obtained in the present work (I_MELF-GOS) and those reported by ICRU Report 49 (I_ICRU 49) <cit.>. It is worth to notice that our calculated value of 79.4 eV is closer to the updated ICRU value of 78 eV <cit.>.

























Composition of the targets used in this work.
   7Composition^ a (% mass)

   7

Target   1C   1H   1N   1O   1P   1Ca   1S

Liquid water   0.00   11.19   0.00   88.81   0.00   0.00   0.00

HAp Ca_10(PO_4)_6(OH)_2   0.00   0.20   0.00   41.41   18.50   39.89   0.00 

Organic part of cortical bone   30.46   9.74   8.87   48.18   1.62   0.00   1.13

Cortical bone^ b, c   14.43   4.72   4.20   44.61   10.50   20.99   0.32



   ^ a For simplicity, elements heavier than S with presence less than 1% have been omitted in the composition.

   ^ b The composition is the same as the one shown in the ICRU Report 49 <cit.>.

   ^ c Mass constitution of cortical bone is made of 58% HAp and 42% organic part <cit.>.



Atomic number Z, molecular mass A, density ρ and mean excitation energy I of the targets considered in this work.
      1A   Density   1I_MELF-GOS   1I_ICRU 49

Target   Z   1(g/mol)   (g/cm^3)   1(eV)   1(eV)

Liquid water   10   18.02   1.00   79.4   75.0

HAp Ca_10(PO_4)_6(OH)_2   500   1008.83   3.22   159.5-162.5   140.2

Organic part of cortical bone    1633   2832.67   1.165   74.6   -

Cortical bone   3094   5934.94   1.85   114.9   106.4










 §.§ Simulation of transport and energy deposition by a proton beam



The SEICS code <cit.> is used in this work to simulate the transport of therapeutical proton beams, with initial energies ranging from 50 MeV to 250 MeV, through liquid water and cortical bone. This simulation code is able to follow the trajectories of energetic ions through a condensed-phase target. It uses molecular dynamics methods to describe the motion of the projectile as well as Monte Carlo techniques to account for the stochastic nature of electronic energy loss, multiple elastic scattering, the electron charge-exchange between the proton and the medium, as well as the nuclear fragmentation of the projectile caused by non-elastic nuclear scattering processes <cit.>. The code uses a (relativistic) variant of the velocity Verlet's algorithm to calculate the new position r⃗ and velocity v⃗ of a projectile of mass M after a time Δ t  <cit.>:

    r⃗_i+1    =    r⃗_i + v⃗_i Δ t+F⃗_i/2M(Δ t)^2 [ 1-(v_i/c)^2 ]^3/2 ,
    v⃗_i+1    =    v⃗_i+F⃗_i+F⃗_i+1/2MΔ t [1-(v_i/c)^2 ]^3/2 ,

where c is the velocity of light, and the subscripts i and i+1 refer to successive time steps separated by a time interval Δ t.  

In equations eq:r-Verlet-rel and eq:v-Verlet-rel, the factor in brackets is an ad hoc correction which takes into account the relativistic velocity of the projectile, whereas F⃗ is the electronic stopping force <cit.>. This force, which is mainly caused by inelastic collisions with the target electrons, acts on the proton with a charge state q. Because of its stochastic nature, it is randomly sampled by a Gaussian distribution whose mean value is the stopping power S_q and whose variance is described in terms of the energy loss straggling Ω_q^2, σ = √(Ω_q^2/Δ s), where Δ s is the distance travelled by the proton in a time Δ t, i.e., Δ s = v Δ t. Therefore, the expression of the electronic stopping force is <cit.>:

    F⃗ = -[ S_q + Ω_q/√(Δ s)√(-2lnR_1)cos(2 π R_2)]  v̂ ,

where R_1 and R_2 are random numbers uniformly distributed between 0 and 1, and v̂ is the unit velocity vector. Both energy loss quantities, S_q and Ω_q^2, are calculated using equations eq:Sq and eq:Stg_q, with the ELF obtained through the MELF-GOS method.





The simulation code SEICS considers not only the elastic interactions between the projectile and the target nuclei but also the nuclear fragmentation reactions <cit.>. These are needed to properly characterize the depth-dose curves, accounting for the interactions between the primary protons and the ones of the atoms of the medium. These interactions involve the excitation of the target nucleus, its fragmentation, the emission of secondary energetic particles (such as neutrons, photons, secondary protons or heavier particles) and the relaxation of the residual nucleus. 

Nuclear fragmentation reactions are approximately included in the SEICS code by removing the primary protons from the beam according to their total nuclear fragmentation cross section, locally depositing a part of their energy <cit.>. The information related to nuclear fragmentation cross sections is obtained from ICRU Report 63 <cit.>, a validated and comprehensive cross-section compilation based on theoretical models and available experimental data. The fragmentation mean free path λ_ fr(v) for protons, which depends on their velocity v, is given by:

    λ_ fr(v) = A/N_ A ρ σ_ fr(v) ,




where A and ρ are, respectively, the mass number and density of the target, whereas N_A refers to Avogadro's number. In addition, σ_ fr(v) is the microscopic fragmentation cross section, which is calculated as the weighted sum of the fragmentation cross sections for the elements that constitute the compound target <cit.>. 

After a distance s, the number of surviving particles of the beam is N(s+ ds) = N(s) e^- ds/λ_ fr, showing an exponential behaviour. In order to take this into account, a random number is sampled in each step s. If it is equal of less than the probability of fragmentation, namely, 1-N(s+ ds)/N(s) ≤ 1, then the primary proton disappears and a part of its energy is locally deposited. 

Because of their low stopping cross sections, the energy of secondary neutrals (photons and neutrons) escapes from the studied volume and this fraction of energy is disregarded. Knowing this, as well as the fact that secondary protons have long ranges compared with those of heavier particles, the treatment is normally simplified by tracking only secondary protons and depositing locally the energy of heavier secondary particles <cit.>. Nevertheless, in the SEICS code the energy of both protons and heavier particles is locally deposited for simplicity <cit.>. It will be seen later how this approximation is enough for our purposes. 















§ RESULTS AND DISCUSSION


The stopping power and energy-loss straggling of liquid water and cortical bone for protons are calculated using equations eq:Sq and eq:Stg_q, considering the equilibrium charge fractions ϕ_q of the projectile, which depend on its energy and the target. 
We have obtained the values of ϕ_q in liquid water and in cortical bone from the CasP code <cit.>, which is based on a parameterisation to experimental data.
fig:spstr(a) shows the stopping powers of liquid water and of cortical bone as a function of the incident proton energy T. Solid lines correspond to the results obtained in the current work by means of the MELF-GOS method; calculations for cortical bone are shown by a red thick solid line, while the ones for liquid water are represented by a black thin solid line. For comparison, the stopping powers provided by the ICRU Report 49 for both targets are depicted by dashed lines <cit.>. Available experimental data for ice <cit.> and liquid water <cit.> are depicted by symbols, while for cortical bone there are no measured data. 







Our results for liquid water excellently agree with the experimental data for liquid water at high energies (star symbols) <cit.>, while slight differences are found with respect to the experiments in liquid water at intermediate energies (triangle and circle symbols) <cit.>. However, it should be taken into account that the latter experiments were performed in liquid water jets, and the experimentally measured proton energy spectra had to be interpreted by means of Monte Carlo simulations in order to infer the stopping power <cit.>; such analysis was done on the basis of a scaling of the stopping power curve given by the SRIM code <cit.>. Garcia-Molina et al. <cit.> showed how the stopping power calculated from the dielectric formalism, without any ad hoc scaling, could also reproduce in detail the measured spectra in water and other cylindrical targets. Thus, the experimentally derived stopping powers of liquid water at intermediate proton energy should be taken with some care. Our calculations also closely follow the experimental stopping power for ice <cit.> down to ∼ 200 keV, underestimating it below this energy. Such differences could be attributed to the limitations of the dielectric formalism at low proton energies and/or to phase differences between ice and liquid water.

Our calculations agree well with ICRU data <cit.> for both targets at high energies. However, the latter systematically overestimate our results at energies ≲ 200 keV for liquid water and ≲ 1000 keV for cortical bone. It should be taken into account that ICRU estimates take into account the experimental information available for each target material; as for cortical bone there are no measurements, it is natural that ICRU values present larger deviations with respect to the present calculations based on the dielectric formalism.

fig:spstr(b) depicts the calculated energy-loss straggling of both targets for protons as a function of the beam energy. Again, the red thick solid line corresponds to cortical bone, whereas the black thin solid one is for liquid water. Due to the fact that there are no available data to compare our results with, the asymptotic high energy value given by the Bohr's straggling formula is included as dotted lines, as a way to contrast our calculations. This quantity has been calculated as <cit.>:

    Ω_ Bohr^2 = 4π e^4 Z_1^2 Z_2  N ,

where Z_1 is the atomic number of the projectile, whereas Z_2 and N are the atomic number and the molecular density of the target. The calculated stragglings converge reasonably well to the corresponding asymptotic Bohr's values.







It is worth to notice that the calculated energy loss quantities for cortical bone are always larger than for liquid water, which should be taken into account if detailed simulations of energetic proton beams interacting with the human body for oncological purposes (i.e., protontherapy) are needed. 

The previous data are necessary to simulate the propagation of proton beams in liquid water and cortical bone, together with the corresponding nuclear fragmentation information for each target. fig:lfrag(a) shows the nuclear fragmentation mean free paths for protons in liquid water (dashed black line) and cortical bone (solid red line), as a function of proton energy, obtained from the data from ICRU Report 63 <cit.>, as explained in sec:SEICS. Cortical bone presents a rather shorter fragmentation mean free path, due to its larger density and the heavier elements present in its composition (see tab-comp and tab-prop-part). The fraction f_ ejectile of proton energy transferred to different ejectiles (neutrals, protons and deuterons, and heavier particles) is depicted in fig:lfrag(b), for both liquid water (thin black lines) and cortical bone (thick red lines), as a function of the primary proton energy. In cortical bone, lower fractions of neutral projectiles and secondary protons are produced with respect to water, whereas the amount of heavier ejectiles is much larger. This will imply that simulations using the SEICS code (which locally deposits the residual energy of both short-range heavier ions –correctly– and long-range secondary protons –incorrectly–) will be more accurate for cortical bone than for liquid water, which are excellent, as will be seen when comparing simulated depth-dose curves with available experimental data in fig:BraggRanges(a).




As an illustration of the influence of nuclear fragmentation reactions, fig:FrafInDose depicts the simulated depth-dose curves for 200 MeV protons in cortical bone, including (solid line) and disregarding (dashed line) the nuclear reactions. The curves have been normalized, so as the areas under the curves have a value of one. Nuclear fragmentation reactions tend to decrease the dose around the Bragg peak maximum, due to the loss of primary ions, whose slowing down is more intense around these depths. By contrast, the dose is increased in the beam entrance region, where primary protons present larger residual energies when fragmenting, as the SEICS code deposits locally this energy. 



In Ref. <cit.>, depth-dose curves for proton beams in liquid water, as simulated by the SEICS code, were compared with available experimental data <cit.> for several initial energies. Such comparison is shown in fig:BraggRanges(a), together with current results for the same beam energies in cortical bone in fig:BraggRanges(c). Experimental depth-dose curves have also been scaled to unit area. There are no experimental depth-dose curves in bone for these energies, but comparison with measured data will be reported later. For liquid water, the agreement between simulated and experimental depth-dose curves is excellent. For the largest energy shown (221.8 MeV), the simulated dose is slightly overestimated in the entrance region, due to the lower reliability of the local-energy deposition scheme used for secondary protons at these depths. Still, the agreement is very good. 

In both fig:BraggRanges(b) and (d), simulated results are shown feeding the SEICS code either with the currently calculated stopping powers and energy loss stragglings (solid lines) or with ICRU stopping powers, supplemented with our calculated stragglings (dashed lines). For liquid water, the shift of the depth-dose curves is difficult to appreciate in the figure, but the differences in the determined proton ranges (estimated as the depths were the relative dose falls to 80% after the maximum) is shown by labels, in units of mm, below each energy label. Such differences grow with the initial beam energy, going up to 2.3 mm at 221.8 MeV. Similar absolute differences in the proton ranges are found for cortical bone using the present or the ICRU slowing down data. However, it should be noted that, since the ranges in cortical bone are considerably shorter than in liquid water, the relative differences are more significant. fig:BraggRanges(b) and fig:BraggRanges(d) show, respectively, the proton ranges in liquid water and cortical bone as a function of the beam energy. Full circles represent SEICS simulations using our calculated stopping powers and energy loss stragglings, while up triangles correspond to simulations using ICRU stopping powers. The energy-range simulations are compared with results given for these targets by NIST's PSTAR database <cit.>, as well as with other simulated results for cortical bone <cit.>. While our simulations seem to agree well with PSTAR data, differences with the results of Ref. <cit.> are more visible.

Our simulated results for the proton range (in cm) are very well reproduced by the following energy–range relation:

    Range  (cm) = a + b · T  ( MeV) + c · T  ( MeV)^2  ,

where T is the initial proton beam energy (in MeV), whereas  a, b and c are fitting parameters, whose value appear in tab-fit. Red dashed lines in fig:BraggRanges(b) and fig:BraggRanges(d) denote the best fits of this equation to the ranges of protons (for energies 50–220 MeV) in liquid water and cortical bone, respectively, determined by simulations using the MELF-GOS stopping quantities. 

 Parameters used in the parabola fitting of the energy-range relationship for liquid water and cortical bone.
Material   a (MeV)   b (cm/MeV)   c (cm/MeV^2)

Liquid water   -1.23341   4.536 · 10^-2   4.54404 · 10^-4

Cortical bone   -0.85881   2.932· 10^-2   2.6587· 10^-4






The experimental information both on the stopping power and on the depth-dose curves for proton beams in cortical bone is very scarce, despite the importance of this target in radiotherapy treatments. The only measured depth-dose curve we are aware of is that of 200 MeV protons interacting with a cortical bone phantom material reported by Vasiliev et al. <cit.>. In the experimental setup, the 200 MeV beam, having an energy width of ± 0.3% , was first passed through a liquid water energy degrader 100 mm in depth. The resulting, further energy widened, beam then bombarded the cortical bone phantom. This setup is sketched in the inset of Fig. <ref>.

The SEICS code allows to introduce an energy spectrum to the incident beam, as well as to measure the energy spectrum of the protons leaving the target slab. The spectra of the protons entering and leaving the liquid water target of 100 mm width are depicted in the inset of Fig. <ref>. The energy loss straggling of liquid water makes the beam to increase its energy width from the initial ± 0.3% to ± 0.93% when exiting the liquid water slab. This latter value has been then used to simulate the proton propagation in the cortical bone target. The red solid line in Fig. <ref> corresponds to the simulation using the stopping powers and energy loss stragglings obtained in this work (both for liquid water and cortical bone), while the black dashed line uses the ICRU stopping powers, complemented with the current energy loss straggling values. As can be clearly seen, MELF-GOS slowing down quantities yield a simulated depth-dose curve which almost perfectly agrees with the experimentally measured data from Ref. <cit.>, except for a slight overestimation of the dose around the Bragg peak maximum (again, results have been normalized, so as the areas under the depth-dose curves have a value of one). Such a good result (particularly, the excellent reproduction of the depth for maximum dose) has been obtained without the need of slight adjustments in the nominal initial beam energy, which sometimes are needed due to the initial energy uncertainties <cit.>. In contrast, the simulation using the ICRU stopping power (for both target materials) results in a Bragg curve maximum slightly shifted to shallower depths, together with a somewhat more intense dose around the maximum, as easily observed in the figure. The range variation associated with the use of one set of stopping powers or the other is of 1.3 mm. This comparison supports the accuracy of the slowing down quantities calculated by means of the MELF-GOS method, as well as the good performance of the SEICS simulation code.



Simulated Bragg curves for proton beams in cortical bone can be used to obtain the energy-range relation (Fig. <ref>(d)) and, from it, to calculate the water equivalent ratio (WER) for this material, an important quantity in protontherapy treatment planning. The WER is obtained as the ratio of the cortical bone range and the liquid water range for a given proton beam energy. WER values were recently reported for cortical and compact bone by <cit.> on the basis of Monte Carlo simulations by means of the MCNP6.1 and TOPAS codes. These data are shown as lines in Fig. <ref>, together with our calculations based on the dielectric formalism (ICRU) calculations and the SEICS code, depicted by full (open) symbols. In general terms, Burin et al.'s and our results are compatible, being both close to 1.70 in the energy range from 50 to 220 MeV. The results based on SEICS, using either the MELF-GOS or the ICRU slowing down quantities, present a similar behaviour with the beam energy, slowly decreasing from ∼ 1.70 at 50 MeV to ∼ 1.65 from around 120 MeV, with small variations among them. However, the evolution of the WER curve obtained by Burin et al. <cit.> is different, presenting a slight monotonic increase in the studied energy range, being always very close to 1.70. Further simulations, using other Monte Carlo codes and different stopping power input data, would be needed to more accurately determine the WER of protons in cortical bone, particularly at the energies > 100 MeV commonly employed in protontherapy.




§ SUMMARY AND CONCLUSIONS




Cortical bone is a very relevant material, almost unavoidably encountered by energetic proton beams during cancer treatment by means of protontherapy. However, despite the importance of  the presence of this material in treatment plans, which will introduce uncertainties in the penetration range of proton beams in the human body, there is still a current lack of accurate information on its energy loss quantities, such as the stopping power and the energy loss straggling. In fact, the assessment of the water equivalent ratio (WER) of cortical bone has only very recently been undertaken by means of general purpose Monte Carlo codes such as TOPAS or MCNP <cit.>, though using simple approximations to the stopping power, such as the Bragg's rule or ICRU estimates <cit.>. The former introduces uncertainties related to the chemical bonding and condensed-phase state of the target, while the accuracy of the latter relies on the availability of measured data, which is lacking in this case. 

The dielectric formalism represents a convenient theoretical method to calculate accurate proton stopping powers and energy loss stragglings in condensed-phase materials.

However, its application requires the knowledge of the target ELF, which is difficult to determine for cortical bone, given its multi-phase nature. In Ref. <cit.>, the results of DFT calculations of the ELF of calcium hydroxyapatite (the mineral phase of cortical bone) <cit.> were combined with an estimate of the ELF of the organic part of cortical bone, obtained from a parametric model to predict the ELF of organic materials <cit.>. Experimental measurements of the stopping power of proton and helium ions in hydroxyapatite performed by <cit.> demonstrated the good accuracy of the dielectric formalism, together with the MELF-GOS method to describe the target ELF. An appropriate combination of the stopping powers for the mineral and organic parts of cortical bone also predicted its stopping power in good agreement with the scarce available information for this target <cit.>.

Still, a full validation of the calculated slowing down quantities requires the reproduction of other experiments, of relevance for protontherapy, by means of radiation transport simulations. We sought to simulate proton depth-dose curves in cortical bone by means of the SEICS code, implementing our calculated slowing down quantities. However, to the best of our knowledge, there is only one experimental measurement <cit.>, whose data are, nonetheless, very useful: a 200 MeV proton beam (with a well characterized energy dispersion) first traversed a liquid water energy degrader, to then bombard cortical bone. Thus this experiment tests the stopping quantities of both liquid water and cortical bone.

As shown in the present work, our calculated stopping power of cortical bone presents significant deviations from the ICRU49 estimate <cit.> for energies below 1000 keV, while these differences only appear below 200 keV for liquid water (a better characterized target, for which there exist experimental stopping power data, even though scarce).

Such deviations are significant enough to produce visible differences in the simulated depth-dose curves in cortical bone: while our calculated stopping quantities lead to a simulated Bragg curve in very good agreement with the measurements by Vasiliev et al. <cit.>, the use of ICRU49's stopping power yields a depth-dose curve with the peak position slightly off, and with a somewhat worse general agreement with experiment. Thus, these results demonstrate the high accuracy of our calculated slowing down quantities, and their visible impact in clinically relevant magnitudes for both important constituents of human body, namely, liquid water and cortical bone.

Provided the excellent agreement of our calculated energy loss quantities, we used them in our SEICS simulation code to study the proton range in liquid water and cortical bone as a function of the beam initial energy. The ratio of the ranges in cortical bone and liquid water can then be used to obtain the water equivalent ratio (WER) for cortical bone as a function of the beam energy. The differences in the simulated penetration range of protons in these materials, when using our current stopping powers or those given by ICRU <cit.> progressively grow with the beam energy. Range differences are around 0.7–0.8 mm for 121.2 MeV protons, while they increase to 2.0–2.3 mm for 221.8 MeV. In relative terms, such differences are more significant for cortical bone than for liquid water, as absolute ranges are rather shorter in the former. The resulting ranges obtained by using the currently calculated stopping powers and energy loss stragglings are well fitted by parabolic functions in the energy range studied. 

The differences in penetration range reflect in differences in the WER values for cortical bone. Simulated WER with SEICS, using either our calculated stopping quantities or those given by ICRU49, present slight deviations, even though their evolution with beam energy is similar. They decrease from ∼ 1.70 at 50 MeV (a value similar to that predicted by Burin et al. <cit.>) to ∼ 1.65 at energies around 120–160 MeV, with a very slight increase afterwards. In contrast, the WER derived in Ref. <cit.>, even though close, presents a slightly different evolution with energy, with a very small but constant increase, being always around 1.70.

As current treatment plannings typically require millimetric precision, and given the fact that cortical bone is difficult to avoid in many clinical situations, the observed differences in the proton ranges and WER values may contribute to undesired uncertainties during protontherapy treatments. It is clear then the importance to know the slowing down quantities (especially the stopping power, but also the energy loss straggling) of protons in cortical bone, at least as accurately as they are known for liquid water. In this work we have presented theoretical calculations and simulations, based on the dielectric formalism, which, by comparison with the experiments performed with calcium hydroxyapatite <cit.> and the measured depth-dose curve for cortical bone phantom <cit.>, seem to offer a very good accuracy. Experimental determinations of the stopping power of protons in cortical bone would be very useful to further check the presented results, as well as to better understand the differences in WER as determined in the present work and those recently obtained in Ref. <cit.>.




This work is part of the R&D project no. PID2021-122866NB-I00 funded by the Spanish Ministerio de Ciencia e Innovación (MCIN/AEI/10.13039/501100011033/) and by the European Regional Development Fund (“ERDF A way to make Europe”), as well as of the R&D project no. 22081/PI/22 funded by the Autonomous Community of the Region of Murcia through the call “Projects for the development of scientific and technical research by competitive groups”, included in the Regional Program for the Promotion of Scientific and Technical Research (Action Plan 2022) of the Fundación Séneca – Agencia de Ciencia y Tecnología de la Región de Murcia.




§ REFERENCES





iopart-num 

































