One-dimensional simulation of proton radiography to assess the feasibility to update the CT based stopping power calibration curve.

(1)

One-dimensional simulation of proton radiography to assess the feasibility to

update the CT based stopping power calibration curve.

Bachelor Thesis by Floris Drent

First examiner:

Emiel van der Graaf Second examiner:

Peter Dendooven

Rijksuniversiteit Groningen The Netherlands

6 July 2018

(2)

1 Introduction

1.1 Proton therapy

In January 2018, the medical centre of the university of Groningen(UMCG) started treating cancer patients with proton therapy. Proton therapy is using proton radiation to bombard tumor tissue instead of the well-known photon radiation. The main benefit of using proton over photon radiation is that protons deposit the bulk of their energy in one small area of their path length whereas photons deposit their energy evenly along their path length. The position where protons deposit most of their energy in one section of their path length is at the so-called Bragg peak. Therefore proton therapy causes less collateral damage to the healthy tissue surrounding the tumor. This makes proton therapy ideal for treating children and patients with tumors in vulnerable areas.

1.2 CT and proton radiography

Proton therapy can only be successfully performed if the position of the Bragg peak within the patient is known. The Bragg peak’s position can be obtained if the energy loss of the proton traversing the patient’s body is also known. The energy a proton loses per unit length through a material is called the material’s stopping power. More generally used is the term relative stopping power or in short, RSP. The RSP of a material is its stopping power relative to that of water. The various RSP-values of the various materials that make up the patient’s body are required in order to perform proton therapy. The direct method to obtain these stopping power values is proton radiography, where protons are send through a patient and the energy losses are used to create an image similar to an image produced by X-rays. However, this method requires expensive, bulky equipment that is not available at hospitals and rarely available at medical re- search facilities.

A more indirect approach to obtain the patient’ stopping power values are through information from a CT-scan on the patient. A CT-scanner virtually di- vides the patient into little blocks called voxels and by sending photons through the patient from multiple angles a computer then applies a so-called HU-value to each voxel, where the HU stands for Hounsfield unit. The Hounsfield unit of a voxel is calculated using the voxel’s linear attenuation coefficient, which is a measurement of the amount of traversing photons that are attenuated by the voxel’s tissue. One could see this as the amount of ’shadow’ the voxel produces.

Every human tissue has its own HU-value and with that information, the computer can reconstruct a 3D image of the patient’s insides. Ideally, one can take several calibration materials that correspond to human tissue, obtain its HU- and RSP-values and plot these two quantities in a calibration curve. After performing a CT-scan on the patient and the HU-values for all the patient’s voxels are obtained, this calibration curve can then be used to assign a RSP-value for each voxel.

(4)

1.3 The calibration curve

The calibration curve is thus the HU plotted against the RSP of several calibration materials. However, when this curve would be created, it can happen that it won’t follow a noticeable trend and the data points are scattered through- out the plot. Another problem that comes with using this calibration curve for obtaining the patient’s RSP is the question whether this calibration curve is suitable for every patient or should each patient has its own specific calibration curve. This depends on the error margins of the calibration curve. The question is if proton radiography is a valid method in order to acquire these error margins. This paper is about validating proton radiography as a method to acquire calibration curve error margins by simulating multiple grids filled with human tissues. Then the energy loss of a proton traversing this grid is calculated and the outcome is used to fine-tune the slopes of the calibration curve. The variation of the slopes that come out of the random grids is then plotted into several distributions. The standard deviation of each distribution is then used to determine the error margins of the calibration curve. If these error margins are around the order of 1%, then proton radiography seems like a valid method for fine-tuning the calibration curve. However, if these margins are below the order of 1%, then proton radiography cannot be validated with this method.

2 Theory

The plotting of the calibration curve starts with the calculation of HU and the RSP of 39 calibration materials. The simulation itself requires the RSP and HU of various human tissues. The contents of these calibration materials and human tissues are elaborated on in Materials and Methods. Various data is required for the calculation of the HU and RSP of a material. For example, the elemental weight fraction of each element of the material is a required for calculating both the HU and the RSP.

2.1 HU

HU of the various tissues are calculated using the linear attenuation coefficient, µ, of the tissues. The coefficient determines the amount of photons emitted by the CT upon the tissue that are absorbed by the tissues. The linear attenuation coefficient is defined as .

I = I₀exp(−

Z S 0

µ(E, l)dl) (1)

where

• I is the photon intensity coming out of the tissue path.

• I0is the starting photon intensity entering the tissue path.

(5)

• µ is the linear attenuation coefficient which depends on – E, the energy of the photon entering the voxel.

– l, the position of the voxel along the tissue path.

The linear attenuation coefficient can also be directly calculated using data known about the tissue.

µ(E) = ρNa 30

X

i=1

(wi

Ai

σi(E)) (2)

where:

• ρ is the mass density in grams per cubic centimeter(g/cm³).

• Na is Avogadro’s constant of 6.022 × 10²³Mol⁻¹.

• P30

i=1 is the summation over the thirty of the most relevant elements of which the materials consist of.

• wi is the mass fraction of element i.

• Ai is the atomic mass of element i.

• σi(E) is the cross section of element i at proton energy E.

The summation over the 30 elementals is not expended beyond the 30 elements because the elements beyond atomic number 30 almost never occur in both the calibration materials and in human tissues. If they would occur, the mass fraction is neglectable.

Due to the fact that the fact that the CT works in discrete energy channels from 1 eV to 150 eV in steps of 1 eV, this formula is first split from the summation over the cross section of the various discrete energies. This gives it a clearer overview. So for every element a D-term is calculated by

Di= ρNa

wi

A_i (3)

and these D-terms are put into a sum product with the cross sections for every discrete energy channel in order to obtain the linear attenuation coefficients for each energy.

µ(Ek) =

30

X

i=1

Diσ(Ek). (4)

where k denotes the corresponding discrete energy channel. Not every energy channel of the CT is evenly used. This is due to the way the X-ray tube in the CT produces the photons. This generates a Bremsstrahlung spectrum with the same characteristic x-ray from the anode materials imposed. The spectrum per energy channel is plotted in figure 2.(The spectra were kindly provided by

(6)

Siemens Medical Solutions, Forchheim, Germany). Linear attenuation coefficients have to be intergrated over this spectrum in order to obtain the mean linear attenuation coefficients. This is done by

¯ µ =

R

kf (E)µ(E)dE R

kf (E) (5)

Where f (E) is the x-ray spectra of the CT corresponding to the voltage used.

Since the energies are discrete and the spectrum is normalised, this integration is reduced to the sum product of the spectrum and the µ over the energy channels.

¯ µ =

k

X

0

f (Ek)µ(Ek) (6)

When the mean linear attenuation coefficient for the materials are determined, the Hounsfield units of the materials can be obtained by

HU = ( µ¯

¯

µ_H₂_O − 1) × 1000 (7)

where ¯µH₂O is the mean linear attenuation coefficient of water.

2.2 RSP

The RSP of a material is the ratio between the energy loss of the proton traversing the material and the energy loss of the proton traversing water.

RSPm= (^dE_dx)m

(^dE_dx)_H₂_O (8)

so when the RSP of a material is known, the energy loss of a proton traveling though the material is

∆E = RSP × (dE

dx)H2O× ∆x (9)

The energy loss is calculated by using the Bethe formula.

−dE

dx = 4πK₀²Z²e⁴n

mec²β² [ln(2mec²β²

I(1 − β²)) − β²] (10) where

• K₀²is 3.16 × 10⁴⁵ MeV² m²C⁻⁴.

• Z is the atomic number of protons, which is one.

• e is the magnitude of the electron charge in Coulomb.

(7)

• n is the electron density in number of electrons per cubic centimeter.

• β is the ratio between the speed of the proton traversing the material and the speed of light, β = ^v_c

• I is the mean excitation energy of the material in MeV.

The mean excitation energies of the materials are calculated by

ln(I) = P30

i=1wiZ_i A_iln(Ii) P30

i=1wiZ_i A_i

(11)

where

• Zi is the atomic number of element i.

• I_i is the mean excitation energy of element i in MeV.

2.2.1 Mean excitation energy

The mean excitation energy of each element is calculated by using the following approximate empirical formulas for each element with atomic number Z.

I ≈ 19.0eV f or Z = 1 (12)

I ≈ 11.2 + 11.7Z eV f or 2 ≤ Z ≤ 13 (13) I ≈ 52.8 + 8.71Z eV f or Z > 13 (14) 2.2.2 Relativistic speed

The value for β depends on the starting speed, and therefore the kinetic energy, of the proton emitted at the patient.

E_kin= m_pc²

p1 − β² − m_pc² (15)

where mpc² is the rest mass of the proton in MeV. This formula rewritten for β² is

β²= 1 − ( m_pc² mpc²+ Ekin

)² (16)

(8)

3 Materials and Methods

3.1 The calibration curve

The creation of the first calibration curve is done using 39 calibration materials.

[1] These calibration materials consist out of both liquids and solids. Examples are plastics, metals, solutions and hydrocarbons. The HU and RSP of these materials are calculated. The HU and RSP of the calibration and human tissue samples are calculated in OpenLibre. The HU of every calibration materials is calculated by first determining the D-term for each of the 30 elements by using formula (3). After this, the µ of each material for each element and energy channel of the cross section is calculated using formula (4). All these µ’s are then integrated over the spectrum in order to obtain the ¯µ. The ¯µ are filled into formula (7) to calculate the HU of the materials. The HU and RSP of the calibration materials are plotted against each other. This creates a continuous curve when ignoring all materials that do not contain hydrogen atoms. This curve is categorized into three domains: the first domain from -800 HU to -40 HU, the second domain from -40 HU to 200 HU and the third domain from 200 HU to 2000 HU. These domains are choses such that each domain corresponds with the HU and RSP of lung, blood and bone respectively. For each domain, the slope, and therefore the ration between RSP and HU, is determined by a local linear fit within the respective domain.

Table 1: Some examples of the 39 calibration materials along with their respective weight fraction per element.

(9)

Figure 1: The normalised spectrum of the photon intensity per energy channel of the CT. The linear attenuation coefficient, µ, is integrated with this spectrum in order to obtain the mean linear attenuation coefficient, ¯µ.

Figure 2: The calibration curve where the HU is plotted against the RSP for each of the 39 calibration materials. The data points in this plot are shown in table 2.

(10)

Table 2: All calibration materials with their respective HU and RSP which are used to plot the calibration curve.

(11)

Figure 3: The first three-piecewise calibration curve with the domains [-800 HU;-40 HU], [-40 HU; 200 HU] and [200 HU; 2000HU]

(12)

3.2 The proton radiography simulation

The proton radiography simulation is performed using Excel. From a data set of 81 human tissues, such as bone, organ and tissue, the HU and the RSP is calculated.[2][3][4] A grid of 50 by 50 voxels is filled with randomly generated integer numbers between 1 and 81. Each number represents one of the human tissues in the list. The grid of random numbers is applied on the list of 81 names, calculated RSP and calibrated RSP so three more grids are generated.

The calculated RSP are the values for stopping power of the human tissues which are calculated by using formulas (8) and (10). The calibrated RSP are the values for stopping power which are obtained by first calculating the HU of the materials and then obtaining the RSP that corresponds with the HU in the calibration curve. These three grids are filled with the corresponding names and the two different types of RSP, calculated RSP and calibration RSP. These two grids represent the tissue of the patient. Fifty rows with a length of fifty blocks of various, randomised human tissues are simulated. The grid can be considered to be 50 separate rows stacked on each other since there will be no interaction between them. The width of these blocks are one millimeter in the simulation. The following procedure is created for the calculated RSP and then copied using the calibrated RSP instead.

1. A grid is created where the first column has the value of 200 MeV. This represents the kinetic energy of the protons that enters their respective row.

2. When a proton traverses through a voxel filled with material, it loses an amount of energy and thus also speed. Therefore the β of the proton will change. For every time a proton travels thought a voxel in the simulation, the new β has to be recalculated in order to obtain the new dE/dx of water. For the first voxel in the row, the β is calculated by filling in 200 MeV in formula (16).

3. Using the newly calculated β in the voxel, the energy loss of water per unit length is calculated for the first voxel. This is done by filling in the new β into the Bethe-formula.

4. The calculated energy loss for water is multiplied with the RSP of the corresponding voxel and also multiplied with a factor of 10⁻³ due to the voxel width being one millimeter. The outcome of this is the energy the proton loses when it travels though the voxel.

5. The energy for the next voxel in the row is calculated by taking the energy of the previous voxel and subtracting the energy loss which was calculated using the β corresponding to the that previous voxel. For the first step, the energy of the previous voxel is set at 200 MeV. After this, this process is repeated until all fifty voxels in the row are filled.

(13)

Multiple additional grids are generated such that each voxel is filled with a quantity which is required for the entire process. These additional grids are filled with the following quantities, each grid having its own quantity:

1. The energy of the proton entering the voxel.

2. The β of the proton entering the voxel.

3. The energy loss of the proton for water for that voxel.

4. The energy loss of the proton in the material for that voxel.

This process is done twice. Once using the calculated RSP and once using the calibrated RSP. The calibrated RSP are determined by calculating the HU of the materials and using the initial calibration curve to obtain the RSP. The final outcome of this simulation are the kinetic energies of the protons after they travel through a row of filled with fifty voxels of random human tissues. This creates two columns filled with the final energies of the protons. One for the calculated RSP and one for the calibrated RSP.

For every row r and for both types of RSP, the difference between the starting energy of 200 MeV and the final energy is noted, ∆E_r^calcand ∆E_r^cali. Calc and cali stand for the types of RSP, respectively calculated and calibrated. This is the total energy loss per row, per type of RSP. The difference in total energy loss between the types of RSP is taken. This difference is then squared and summed over separately for each row. The difference between the total energy losses of both types of RSP is the smallest if this value is at its minimum. Therefore this value is called the to-be-minimised (TBM).

T BM =

50

X

r=1

(∆E_r^calc− ∆E^cali_r )² (17)

The total energy loss using the calibrated RSP is subtracted from the total energy loss using the calculated RSP. This is done for every row and the average is taken. This produces the average difference in energy loss(ADEL). This average can be positive or negative where the TBM is always positive. If the grid with random numbers is generated numerous times, in the end results, the amount of positive ADEL and negative ADEL should be the same.

ADEL = 1 50

50

X

r=1

(∆E^calc_r − ∆E_r^cali) (18)

(14)

Figure 4: Schematic drawing of the process of filling up one row. The starting kinetic energy of the proton entering the first voxel is set at 200 MeV. At row one, the starting energy is filled into the first voxel. In row two, the β is calculated from the starting energy. This β is used to calculate the energy loss per unit length in water in the third row. This (^dE_dx)wateris multiplied with the length of the voxel and the RSP corresponding to the voxel in order to obtain the energy loss. This energy loss is filled into row four and subtracted from the starting energy in order to determine to kinetic energy of the proton entering the next voxel. This process thus repeats itself until all voxels are filled.

(15)

Figure 5: Schematic drawing of the simulation where the block in the middle represents a grid of 50 by 50 voxels. The energies are the initial and final kinetic energies of the protons entering and leaving the grid from left to right.

The solver of Excel is applied upon the TBM. The parameters that are changed are the three slopes of the three domains in the calibration curve, a₁, a₂, a₃. The solver changes these slopes such that the TBM becomes the smallest value it can find. The slopes that produces the smallest TBM, and the ADEL that comes with it, are noted. The random number grid is generated again with a new setup. The process of finding the parameters that cause the smallest TBM and noting them along side the corresponding ADEL is repeated fifty times, each with a newly generated random number grid.

This produces a distribution of values for ideal slopes for each of the three domains in the calibration curve. These values are put into a histogram in order to create a distribution curve. Also the number of positive and negative ADEL are counted and compared with each other. The histograms are used to obtain the distribution curves with the desired parameters µ and σ. The formula for a distribution curve is

f_a,µ,σ(x) = a exp(−(x − µ)²

2σ² ) (19)

A table is created for each slope where the middle point of the distribution brackets is set as x and the corresponding frequency is set as y(x). The sum of the differences between y and f_a,µ,σ(x) squared for each point x is minimised in

(16)

Figure 6: The process of obtaining the parameters and the ADEL. This process is repeated multiple times in order to obtain a distribution of slopes a₁, a₂ and a₃.

order the find the most optimal parameters a, µ, σ.

T BMa_i =

6

X

x=1

(fa,µ,σ(x)i− y(x)i)² (20)

Here T BM_a_i is the to-be-minimised of the slopes a₁, a₂ and a₃ in order to find the parameters for f_a,µ,σ(x).

(17)

4 Results

The process of minimising and noting down the slopes a1, a2, a3 is repeated 50 times. The resulting distributions are shown in tables 3, 4 and 5. Every time a value for a slope is within one of the brackets, the slope frequency of that bracket goes up by one. These three tables produce histograms and corresponding fitted normal distributions as plotted in figures 7, 8 and 9. The parameters of the normal distributions, the amplitude a, the expectation value µ and the standard deviation σ, for each of the three plots are shown in table 6. The tables for approaching the histogram in order to create the normal distribution curves are in appendix A.

Table 3: Distribution table for slope a1

Table 4: Distribution table for slope a2

(18)

Table 5: Distribution table for slope a₃

Table 6: Table of distribution parameters a, µ, σ for each slope distributions of a1, a2, a3

(19)

Figure 7: Histogram distribution and normal distribution of slope a1 plotted over each other.

(20)

(21)

(22)

5 Discussion

5.1 Relative deviation

The standard deviations of each slope are plugged back into the calibration curve to see the difference the deviation causes. This is done by determining the RSP_cali^±σ where the slope is ai± σ. This RSP_cali^±σ is then compared to the RSPcaliwhere the slope without ±σ is used.

(1 − RSP_cali^±σ RSP_cali⁽⁰⁾

) × 100% = (1 −(a_i± σ)HU + b_i aiHU + bi

×)100% (21)

The outcome of this equation is the relative RSP difference or more generally known as the calibration curve error margin. The HU that are plugged into this equation are the end points of the domains with which the first calibration curve is split into three linear pieces. The i denotes three domains. The relative RSP difference the standard deviation creates in the calibration curve is noted in table 7. There are several things to note for the figures in table 7.

• The difference between relative RSP differences using +σ and the relative RSP differences using −σ is in the order of 10⁻¹⁴ and therefore there is no significant difference between applying either +σ or −σ.

• The discontinuity due to the fact that only the slopes are changed when fitting the calibration curve and not the intersects of the linear pieces is shown at the HU -40 and 200. There are significant differences in relative RSP differences at these points. The discontinuity causes causes gaps between the three linear pieces of the three piece-wise calibration curve.

• The relative RSP difference at the HU of -800 is much larger than at the HU of 1000.

(23)

The last remark raises suspicion that the relative RSP difference is larger for negative HU than for positive HU. Therefore the relative RSP differences for HU between -1000 and 1000 in steps of 100 HU are shown in table 8. It is shown in this table that the relative RSP difference for a HU of -1000 is three times as large as that of a HU of 1000. This table is put into a histogram in figure 10. The relative RSP correspondence is 100% minus the relative RSP difference. Ignoring the dissent of the bars at the HU -100, 100 and 200, the histogram shows a much sharper decline in relative RSP correspondence in the negative HU domain compared to that in the positive HU domain. This is due to the fact that the expectation value for the slope that corresponds with the negative HU, µ1is about twice as large as that of the expectation value for the slope that corresponds with the positive HU, µ3. This, while their respective standard deviations are about the same. This shows that the error margin for calibration is larger for lung tissue, which is in the HU domain of between -1000 and -400, compared to soft tissue, which is in the middle HU domain and bone tissue, which is in the higher HU domain. This is in agreement with previous measured error margins. (G Poludniowski et al.)

Finally, the most important outcome to take note of, is that for almost each HU, the relative RSP difference is in the order of a 1-5%. This is quite significant when considering that above 1% differences could already imply that using the same calibration curve for several patients could cause the Bragg peak to be off target.[5]

5.2 Error Discussion

Although the relative RSP differences are in the order of 1%, the method of acquiring the distribution of slopes has several flaws compared to a realistic proton radiography.

• This simulation of proton radiography is severely simplistic. It is one- dimensional, done in excel with only applying the Bethe-formula. This is of course not a good representation of reality.

• Additionally, the simulated patients are randomised grids human tissue data. Only one set of human tissue data is used, so the accuracy of these grids also depend of the accuracy of this data set. These grids are not a good representation of patients partly due to the fact that if enough randomised grids are used, the total energy loss would approach an average due to every row set-up that is possible with the 81 human tissues would be used multiple times.

(24)

• The intersect of the three pieces of the three piece-wise linear fit of the calibration curve are set and do not change during the simulation or fitting.

Therefore the calibration curve can be seen as three lose lines rotating at a fixed axis when being fitted. This causes discontinuity at the end points of the three domains. This can be seen in the relative RSP differences in table 7.

Despite these errors in the simulation and acquisition of relative RSP differences, one could argue that if these decently-sized relative differences of in the order of 1% already appear in a very crude, simplistic model, it could also quite likely appear in a more realistic model.

5.3 ADEL

For every iteration of obtaining the slopes for the distribution curves, the average difference in energy loss was also noted. All the ADEL are summed over all the iterations. This gives a general ADEL of 3.26. Out of the course of 50 iterations, 6 more positive ADEL than negative ADEL are acquired. What this tells is that, on average, the total energy loss using calculated RSP is slightly larger than the total energy loss when using calibrated RSP. This shows another error in the simulation.

Table 7: The relative RSP differences when using either a_i+ σ_i or a_i− σi. The first two rows are for a₁, the middle two rows are for a₂ and the last two rows are for a₃.

(25)

Table 8: The array of RSP^±σ for HU between -1000 and 1000. The two columns on the right show the relative RSP difference for both +σ and −σ.

(26)

Figure 10: Histogram distribution showing how much the RSP^±σ corresponds with the RSP without the ±sigma. Despite the discontinuities at the HU of -100, 100 and 200, this histogram shows that the calibration error margins are bigger for human tissues such as lung than for soft tissue.

(27)

6 Conclusion

The relative RSP differences, or calibration curve error margins, are found to be around 1-5% for most human tissues . This is in correspondence with previously found values by others. Therefore it seems like proton radiography is a valid method in order to acquire margin errors in the calibration curve. However due to the simplistic nature of the simulation, a significantly more realistic approach is recommended in order to find out if the calibration curve error margins can be obtained by performing proton radiography and finally if the calibration curves have to be patient-specific or not.

(28)

7 Bibliography

7.1 References

References

[1] van Abbema J. Accurate relative stopping power prediction from dual energy CT for proton therapy Phys. Med. Biol. 2015; 60:3825-3846, 2017.

[2] D.R. White, H. Q. Woodard, and S.M. Hammond. Average soft-tissue and bone models for use in radiation dosimetry. Br. J. Radiol. 60(717), 907-913, 1987.

[3] D.R. White, E.M. Widdowson, H.Q. Woodard and J.W.T. Dickerson. The composition of body tissues (ii) fetus to young adult. Br. J. Radiol. 64(758), 149-159, 1991.

[4] H.Q. Woodard and D.R. White. The composition of body tissues Br. J. Ra- diol. 45(708), 1209-1218, 1986.

[5] Poludniowski G, Allinson NM, Evans PM. Proton radiography and tomogra- phy with application to proton therapy. Br. J. Radiol. 2015; 88:20150134.

7.2 Additional literature

• P.J. Doolan et al Patient-specific stopping power calibration for proton therapy planning based on single-detector proton radiography Phys. Med.

Biol.60 1901, 2015

• James E. Turner Atoms, Radiation, and Radiation Protection

• Wilfried Schneider et al Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions Phys. Med. Biol.45 459, 2000

(29)

Appendices

Appendix A Distribution-Approachment Tables

Table 9: The distribution-approach table for determining the distribution parameters a, µ, σ of slope a1.The T BM in the bottom right is the summation of the right column. The distribution parameters of f (x) are changed such that the TBM is the smallest.

(30)

Table 10: The distribution-approach table for determining the distribution parameters a, µ, σ of slope a2.

Table 11: The distribution-approach table for determining the distribution parameters a, µ, σ of slope a3.

One-dimensional simulation of proton radiography to assess the feasibility to update the CT based stopping power calibration curve.