**2. Monte Carlo methods and nuclear medicine**

#### **2.1 Monte Carlo radiation transport simulation**

2 Will-be-set-by-IN-TECH

Numerous Monte Carlo codes have been benchmarked and used for different purposes in radiation transport including nuclear medicine imaging and dosimetry. Some examples are MCPT code of Williamson (1988), EGS4 code (Luxton and Jozsef, 1999), GEANT4 code (Agostinelli, 2003) and MCNP4 code (DeMarco et al 2002b, Bohm et al 2003). The more recently developed PENELOPE (PENetration and Energy Loss Of Positrons and Electrons; photon simulation was introduced later) Monte Carlo code 2001 (Salvat et al 2001) uses cross-sections from the up-to-date EPDL97 dataset, and thus takes advantage of the latest improvements in cross-section libraries [R.D. Stewart & Strom (2001), A. Sánchez-Crespo &

Although several MC codes are available allowing to transport photons and electrons in user-defined geometries, it should be mentioned that there two codes particularly useful for

As a particular and non common approach this chapter will discuss how to implement MC

The proposed method along with the developed computation system allow to introduce as input data different types of images -both metabolic and anatomical- commonly used for nuclear medicine diagnostic, like CT, MRI, SPECT, gamma camera and PET. This capability arises from the incorporation of dedicated routines for reading and interpreting the "Digital Imaging and Communication in Medicine" (DICOM) information code, therefore significantly simplifying input data handling processes aimed to specific-patient absorbed dose distribution calculation. In this sense, suitable combination of metabolic and anatomical imaging techniques provides relevant information helpful when attempting reliable and accurate specific-patient dosimetry. When considering all together the mentioned potential advantages of the proposed method along with the implementation of dedicated voxelization techniques, it would be expected that it may constitute an important tool for daily treatment

voxel dosimetry applications, named EGS4 and MCNP [E.B. Bolch & Watson (1999)].

planning by means of MC simulations for specific-patient dosimetry assessment.

The estimation of dose from a distributed source is an important aspect of the application of labeled monoclonal antibodies to targeted radiotherapy. A major contribution to the dose arises from beta particles. The source may be described mathematically by a function giving the activity concentration at each position in the source. The dose produced at a specific site may then be calculated as the linear superposition of contributions from each volume element in the distribution treated as a point source. The dose produced by a point source of isotropic unit activity is known as dose point kernel [Prestwich WV (1985)] (DPK). Due to its significant importance and largely proved benefits and convenience for dosimetric purposes, carefully characterization of DPK has always received significant interests and efforts, therefore extended discussion will be presented about DPK calculation by means of

In summary, it would be desirable to assess a suitable and reliable method helpful for patient-specific treatment planning calculating full stochastic radiation transport by means of MC methods. Therefore, once the calculation system is already developed, it should be carried out careful and rigorous verifications of certain parameters which can point out the

In this sense, considering that the main goal of this work consists on performing dosimetry for beta-minus radiation therapy and according to the significant importance of beta DPK in nuclear medicine dosimetry, this quantity may be taken as a suitable parameter to check the feasibility of the proposed MC code. With the aim of simplifying, it is established from here

that DPK will be referred as beta-DPK unless it will be explicitly specified.

simulation in nuclear medicine by means of the PENELOPE code.

Larsson (2004)].

MC techniques.

reliability of the proposed method.

Radiation transport, including absorption and scattering processes are determined by the Boltzmann transport equation. However, it is hard to apply and to analytically solve Boltzmann equation within non-homogeneous media or regions consisting of complex boundaries. Actually, it is well known that Boltzmann equation can be analytically solved only in few cases consisting on oversimplified situations, that would strongly differ from real clinical situations. Therefore, numerical methods have been proposed for solving Boltzmann equation in complex situations performing full radiation transport by means of stochastic methods. The MC techniques belong to the most important group within this category.

Physically, the evolution of an electron-photon shower has an intrinsic random nature, so MC methods provide a convenient alternative to deal with transport problems [Salvat (2009)].

On the MC method for radiation transport, the "life" of each particle is seen as a random sequence of free paths that end with an interaction event where the particle suffers any kind of interaction. Generally, the "interaction event" refers to any change in the particle quantum state and it can arise from the change of its movement direction, the loss of energy or the production of a secondary particles. The interaction type depends on the kind of actual particle along with the medium where it is moving by random processes because of its nature. In order to point out main characteristics about how to handle with MC codes with the aim of performing nuclear medicine dosimetry, it has been selected to work with the PENELOPE main code. The fundamental motivations for using PENELOPE relays on its suitability of electron and photon transport along with the relevant characteristic of providing open source routines (written on fortran 77 language). The general considerations described below could be similarly implement on any general purposes MC main code.

#### **2.2 The PENELOPE Monte Carlo code**

The PENELOPE v.2008 is a MC algorithm and computer code for the simulation of coupled electron-photon transport. The simulation algorithm is based on scattering models that combines numerical databases with analytical cross section models for the different interaction mechanisms and it is applicable to energies (kinetic energies in the case of electrons and positrons) from a few hundred *eV* to approximated 1*GeV*. Photon transport is simulated by means of the conventional detailed method. The simulation of charged particles (electron and positron) transport is performed by means of a mixed procedure. Hard interactions, with scattering angle *θ* or energy loss *W* greater than preselected cutoff values *θ<sup>c</sup>* and *Wc* , are simulated in detail. Soft interactions correspond to scattering angle and energy loss less than the corresponded cutoffs *θc*, *Wc* and these interactions are considered by means of suitable mechanisms for simulation condensation, mainly based on multiple scattering theories, like Mollier theory. The user can select the cutoff parameters quite large looking for speeding up the calculation considerably, but it should be preliminary carefully studied the ranges for each parameter that may produce non-negligible alterations in the final scores. A characteristic feature of the here presented simulation code is that the most delicate parts of the simulation are handled internally; electrons, positrons and photons are simulated by calling the same subroutines. Thus PENELOPE makes the practical simulation of electrons and positrons as simple as that of photons (although simulating a charged particle may take a longer time).

One of the main advantages of the latest version of PENELOPE (PENELOPE v. 2008 package) regards the improvements for inner shells ionizations by electron and positron impact, which are described by using a numerical database of total cross-sections for K, L and M

 ∞ 0

*<sup>F</sup>*(*x*, *<sup>E</sup>*0) = *RCSDAS*(*E*(*x*, *RCSDA*))

many calculations of dose from various distributions of beta sources.

inhomogeneity around the emission point.

The specific absorbed fraction associated with *i*

<sup>2</sup>Φ*i*(*r*) = *Ei*

and from the Eq. (9) it follows that

then the scaled DPK for *i*

available, the CSDA range.

4*πρr*

associated with beta-decay may in general be written as:

And therefore, according to the simplified model, the scaled DPK becomes, from Eq. (7)

*E*0

where the track average stopping power (*< S >*= *E*0/*RCSDA*)has been introduced. In this sense, it could be possible to interpret the scaled DPK as the ratio of the stopping power at a particular point along the electron track to the average stopping power over the entire track. It is also customary to use radial distributions of dose around isotropic point sources of electrons or beta-emitters in an infinite water medium, so called dose DPKs, as the basis of

Dosimetry for Beta-Emitter Radionuclides by Means of Monte Carlo Simulations 269

In fact, DPKs have proved to be adequate and remarkable useful for several dosimetric purposes in nuclear medicine. Actually, DPKs are commonly used for daily applications and calculations when some approximations may be acceptable, like non high tissue

Contrary to monoenergetic sources, for the case of beta-emitting radionuclides it will be necessary to calculate the DPK as resulting from a spectrum of electrons. The spectrum

> *N* ∑ *i*=1

which represents a decomposition into *N* groups, each of which has a branching probability

 *Ei* 0

*n*(*E*) =

of *β<sup>i</sup>* and end-point energy *Ei*. Also associated with each group is an average energy

*E*0 *< Ei >*

*E*0 *RCSDA < Ei >*

*Fi*(*x*) = 4*πρr*

A scaled DPK for a beta group may also be introduced. To this end a scaled distance is defined as the fraction of the CSDA range for an electron with energy *Ei* . Designating the latter *ri* ,

where *x* = *r*/*ri* . Scaling the distance in this manner differs from the suggested by Berger, who uses the 90% point. It has the advantage, however, of using a quantity that is more readily

As indicated before, calculation of the DPK requires the evaluation of the CSDA range. In order to accomplish this task approximations have been employed. As first approach, it may

*< Ei >*=

given that the spectral distribution *ni*(*E*) is normalized to unit area.

<sup>Φ</sup>*i*(*r*) = *Ei*

0

0

*th* beta group becomes

*F*(*x*, *E*0)*dx* = 1 (8)

*βini*(*E*) (10)

*Eni*(*E*)*dE* (11)

*ni*(*E*0)Φ(*r*, *E*0)*dE*<sup>0</sup> (12)

, *E*<sup>0</sup> 

<sup>2</sup>*ri*Φ*i*(*r*) (14)

*dE*<sup>0</sup> (13)

*th* beta-group is given by

 *r RCSDA*

*ni*(*E*0)*F*

*<sup>&</sup>lt; <sup>S</sup> <sup>&</sup>gt;* (9)

<sup>=</sup> *<sup>S</sup>*(*E*(*x*, *RCSDA*))

shells derived from the theory described by Bote and Salvat (2008). It was also included photon polarization effects in Rayleigh and Compton scattering, as well as refinements in the modeling of Rayleigh photon scattering and of inelastic collisions of electrons and positrons. These improvement would provide significantly better description of charged particle transport and energy deposition.

In this study, PENELOPE v. 2008 MC code has been considered for use in nuclear medicine dosimetry assessment as well as treatment planning approach when possible. During the last years, PENELOPE has been widely employed for different studies by physicists on the radiation transport field [J. Asenjo & Sánchez-Reyes (2002), S.J. Ye & Naqvi (2004), J. Sempau & Salvat (2004), Bielajew & Salvat (2001), J. Sempau & Fernández-Varea (2001)], therefore supporting the selected choice.

#### **2.3 Dose Point Kernel**

The Dose Point Kernel can be defined considering an extremely simplified model of a monoenergetic isotropic point source emitting electrons which move outwards within an homogeneous medium and slow down continuously according to a stopping power function *S*(*E*) [W.V. Prestwich & Kwok (1989)]. If the source energy is *E*0, therefore the distance *r* from the origin, at which the source is located, is related with the remaining energy *E*(*r*) by:

$$r = \int\_{E(r)}^{E\_0} \frac{dE}{S(E)}\tag{1}$$

The CSDA Range (*RCSDA* also called *r*<sup>0</sup> by other authors) satisfies

$$R\_{\mathbb{C}SDA} = \int\_0^{E\_0} \frac{dE}{\mathbb{S}(E)}\tag{2}$$

Then, the absorbed dose per source transformation in a spherical shell is:

$$D(r) = \frac{\mathcal{S}(E(r))}{4\pi\rho r^2} \tag{3}$$

The specific absorbed fraction is given by

$$\Phi(r, E\_0) = \frac{S(E(r))}{4\pi\rho r^2 E\_0} \tag{4}$$

It can be shown that Φ(*r*, *E*0) satisfies the constraint

$$\int\_0^\infty 4\pi \rho r^2 \Phi(r, E\_0) dr = 1\tag{5}$$

It is convenient to introduce a dimensionless quantity to represent distance as the fraction of the CSDA range, asigning:

$$\mathbf{x} = \frac{r}{R\_{\rm CSDA}}\tag{6}$$

The scaled electron DPK is defined through the relation

$$F(\mathbf{x}, E\_0) = 4\pi\rho r^2 R\_{\rm CSDA} \Phi(r, E\_0) \tag{7}$$

It will be useful to introduce a new quantity satisfying

4 Will-be-set-by-IN-TECH

shells derived from the theory described by Bote and Salvat (2008). It was also included photon polarization effects in Rayleigh and Compton scattering, as well as refinements in the modeling of Rayleigh photon scattering and of inelastic collisions of electrons and positrons. These improvement would provide significantly better description of charged

In this study, PENELOPE v. 2008 MC code has been considered for use in nuclear medicine dosimetry assessment as well as treatment planning approach when possible. During the last years, PENELOPE has been widely employed for different studies by physicists on the radiation transport field [J. Asenjo & Sánchez-Reyes (2002), S.J. Ye & Naqvi (2004), J. Sempau & Salvat (2004), Bielajew & Salvat (2001), J. Sempau & Fernández-Varea (2001)], therefore

The Dose Point Kernel can be defined considering an extremely simplified model of a monoenergetic isotropic point source emitting electrons which move outwards within an homogeneous medium and slow down continuously according to a stopping power function *S*(*E*) [W.V. Prestwich & Kwok (1989)]. If the source energy is *E*0, therefore the distance *r* from the origin, at which the source is located, is related with the remaining energy *E*(*r*) by:

*dE*

*dE*

 *E*<sup>0</sup> 0

*<sup>D</sup>*(*r*) = *<sup>S</sup>*(*E*(*r*))

<sup>Φ</sup>(*r*, *<sup>E</sup>*0) = *<sup>S</sup>*(*E*(*r*))

It is convenient to introduce a dimensionless quantity to represent distance as the fraction of

*<sup>x</sup>* <sup>=</sup> *<sup>r</sup> RCSDA*

4*πρr*2*E*<sup>0</sup>

*<sup>S</sup>*(*E*) (1)

*<sup>S</sup>*(*E*) (2)

<sup>4</sup>*πρr*<sup>2</sup> (3)

<sup>2</sup>Φ(*r*, *E*0)*dr* = 1 (5)

<sup>2</sup>*RCSDA*Φ(*r*, *E*0) (7)

(4)

(6)

*r* = *E*<sup>0</sup> *E*(*r*)

*RCSDA* =

The CSDA Range (*RCSDA* also called *r*<sup>0</sup> by other authors) satisfies

Then, the absorbed dose per source transformation in a spherical shell is:

 ∞ 0

4*πρr*

*F*(*x*, *E*0) = 4*πρr*

particle transport and energy deposition.

The specific absorbed fraction is given by

the CSDA range, asigning:

It can be shown that Φ(*r*, *E*0) satisfies the constraint

The scaled electron DPK is defined through the relation

It will be useful to introduce a new quantity satisfying

supporting the selected choice.

**2.3 Dose Point Kernel**

$$\int\_0^\infty F(\mathbf{x}, E\_0) d\mathbf{x} = 1 \tag{8}$$

And therefore, according to the simplified model, the scaled DPK becomes, from Eq. (7)

$$F(\mathbf{x}, E\_0) = \frac{R\_{\rm CSDA} S(E(\mathbf{x}, R\_{\rm CSDA}))}{E\_0} = \frac{S(E(\mathbf{x}, R\_{\rm CSDA}))}{~~} \tag{9}~~$$

where the track average stopping power (*< S >*= *E*0/*RCSDA*)has been introduced. In this sense, it could be possible to interpret the scaled DPK as the ratio of the stopping power at a particular point along the electron track to the average stopping power over the entire track.

It is also customary to use radial distributions of dose around isotropic point sources of electrons or beta-emitters in an infinite water medium, so called dose DPKs, as the basis of many calculations of dose from various distributions of beta sources.

In fact, DPKs have proved to be adequate and remarkable useful for several dosimetric purposes in nuclear medicine. Actually, DPKs are commonly used for daily applications and calculations when some approximations may be acceptable, like non high tissue inhomogeneity around the emission point.

Contrary to monoenergetic sources, for the case of beta-emitting radionuclides it will be necessary to calculate the DPK as resulting from a spectrum of electrons. The spectrum associated with beta-decay may in general be written as:

$$m(E) = \sum\_{i=1}^{N} \beta\_i n\_i(E) \tag{10}$$

which represents a decomposition into *N* groups, each of which has a branching probability of *β<sup>i</sup>* and end-point energy *Ei*. Also associated with each group is an average energy

$$ = \int\_0^{E\_l} E n\_i(E) dE \tag{11}$$

given that the spectral distribution *ni*(*E*) is normalized to unit area. The specific absorbed fraction associated with *i th* beta-group is given by

$$\Phi\_i(r) = \int\_0^{E\_l} \frac{E\_0}{} n\_i(E\_0)\Phi(r, E\_0)dE\_0\tag{12}$$

and from the Eq. (9) it follows that

$$4\pi\rho r^2\Phi\_i(r) = \int\_0^{E\_i} \frac{E\_0}{R\_{\rm CSDA} < E\_i >} n\_i(E\_0) F\left(\frac{r}{R\_{\rm CSDA}}, E\_0\right) dE\_0 \tag{13}$$

A scaled DPK for a beta group may also be introduced. To this end a scaled distance is defined as the fraction of the CSDA range for an electron with energy *Ei* . Designating the latter *ri* , then the scaled DPK for *i th* beta group becomes

$$F\_i(\mathbf{x}) = 4\pi\rho r^2 r\_i \Phi\_i(r) \tag{14}$$

where *x* = *r*/*ri* . Scaling the distance in this manner differs from the suggested by Berger, who uses the 90% point. It has the advantage, however, of using a quantity that is more readily available, the CSDA range.

As indicated before, calculation of the DPK requires the evaluation of the CSDA range. In order to accomplish this task approximations have been employed. As first approach, it may

(a) (b)

Dosimetry for Beta-Emitter Radionuclides by Means of Monte Carlo Simulations 271

The accuracy of Monte Carlo simulations strongly depends on the accuracy in the probability density functions, and thereby the material cross-section libraries used for radiation transport

The PENELOPE parameter *Eabs* determines -for each simulated material and transported particle type- the energy threshold for particle absorption, defining the energy value below which transport is no longer simulated and residual energy is assumed to be deposited locally. The lowest possible *Eabs* value is 100*eV* according to manual recommendations based on the

With the aim of performing radiation transport for dosimetric purposes, the first step may be the design and development of a code devoted to the DPK calculations, which may basically consist on a main program able to allow the user to set up the energy or the spectrum to be radiated by the point source centered in a geometry defined by concentric spheres defining shells of thickness *RCSDA*/40 (the corresponding *RCSDA* for each material is calculated for the actual or effective energy for monoenertic sources or radionuclides, respectively). Due to phisical properties of charged particles, it may be enough to tally DPK up to 1.5*RCSDA*. Finally, the material for which the DPK needs to be investigated along with the corresponding

The geometry file can be obtained from a dedicated user-designed program devoted to set up the *RCSDA* obtained from MC or reliable databases. Once the user provides *RCSDA* the program will automatically generate the corresponding file according to internal PENELOPE

Preliminary simulations are performed to achieve *Eabs* parameter values in order to optimize computation process assessing suitable balances between statistical errors and the CPU computation times. Therefore, after exhaustive examinations, absorption cutoff values of 100*eV* (minimum recommendable by PENELOPE developers) have been considered for emitting sources of energies lower than 100*keV* (*i.e.* 5*keV*, 10*keV*, 20*keV* and 50*keV*) whereas cutoff values of 1*keV* have been set for emitting sources with energies equal or higher than

As example from the preliminary simulations devoted to achieving suitable absorption energy cutoff valued, it can be mentioned that in the case of the 100*keV* monoenergetic source, it was found a maximum difference in the DPK of 0.16% when changing cutoff value from 1*keV* to the minimum possible of 100*eV*, which may be taken as negligible. However, the situation is significantly different when dealing with low energy emitting sources. DPK differences up to 5% may be found when changing absorption cutoff values from 1*keV* to 100*eV* in the case of

Fig. 1. 90Y and 177Lu spectra compared with its *EEff* calculated.

simulation parameters may be adjusted by the user in the input file.

100*keV* (*i.e.* 100*keV*, 200*keV*, 500*keV*, 1*MeV*, 2*MeV* and 5*MeV*).

lower energy (50*eV*) in PENELOPE v. 2008 database.

requirements for simulation geometries.

5*keV* monoenergetic source.

calculations.

be useful to look for databases containing relevant and required information, like *Stopping Power and Range Tables for Electrons* which may be available from the National Institute of Standards and Technology1. A good recommemdation is that database values should be compared with the corresponding data used by the MC code (it can be assessed using the Table subroutine in the case of the PENELOPE v. 2008 code). In he case of PENELOPE v. 2008, it can be seen that NIST tabulated data and PENELOPE database values present negligible differences. Therefore, CSDA ranges required for shell thickness calculations can be derived from internal MC database or NIST tables without distinction.

On the other hand, in order to establish the corresponding dose delivering shell thicknesses in the case of isotope emission spectra, it may be considered the approximation of treating an isotope spectrum as monoenergetic source with an effective energy (*EEff* ) associated with the average (properly weighted according to the actual spectrum) energy of the whole emission spectrum. Due to the fact that shell thicknesses depend upon source energy, once material is given, it becomes clear the motivations for introducing and determining *EEff* in order to assess some effective or equivalent energy value with the aim of calculating the corresponding shell thickness for the DPK calculations.

As a consequence of this assumption, the corresponding shell thickness for the k-isotope (*RCSDA*) can be calculated from the relation:

$$R\_{\rm CSDA}^k = R\_{\rm CSDA} \left( E\_{Eff}^k \right) \tag{15}$$

where *k* =90Y, 131I, 89Sr, 153Sm, 177Lu, ... and *EEff* is computed according to:

$$E\_{Eff}^k = \frac{\sum\_{i=1}^N p\_i E\_i}{\sum\_{i=1}^N E\_i} \tag{16}$$

Where *N* is the number of channels of the *k*-isotope energy spectrum and *pi* and *Ei* are its weight and energy respectively. It can be seen, as reported in Table 1, two examples of the calculated *E<sup>k</sup> Eff* and *<sup>R</sup><sup>k</sup> CSDA* along with the corresponding graphic sketch (Figure 1) outlining how *EEff* may be estimated. Radionuclide emission spectra were extracted from the NIST database tables2.


Table 1. *RCSDA* used and *E<sup>k</sup> Eff* calculated for different isotopes.

#### **2.4 Dedicated subroutines based on Monte Carlo techniques**

The Monte Carlo method has been successfully applied to radiation transport problems in clinical dosimetry. The method is particularly useful for complex problems that can not be solved analytically and/or deterministically or when measurements are not practical.

<sup>1</sup> http://www.nist.gov

<sup>2</sup> Idem 1.

6 Will-be-set-by-IN-TECH

be useful to look for databases containing relevant and required information, like *Stopping Power and Range Tables for Electrons* which may be available from the National Institute of Standards and Technology1. A good recommemdation is that database values should be compared with the corresponding data used by the MC code (it can be assessed using the Table subroutine in the case of the PENELOPE v. 2008 code). In he case of PENELOPE v. 2008, it can be seen that NIST tabulated data and PENELOPE database values present negligible differences. Therefore, CSDA ranges required for shell thickness calculations can be derived

On the other hand, in order to establish the corresponding dose delivering shell thicknesses in the case of isotope emission spectra, it may be considered the approximation of treating an isotope spectrum as monoenergetic source with an effective energy (*EEff* ) associated with the average (properly weighted according to the actual spectrum) energy of the whole emission spectrum. Due to the fact that shell thicknesses depend upon source energy, once material is given, it becomes clear the motivations for introducing and determining *EEff* in order to assess some effective or equivalent energy value with the aim of calculating the corresponding

As a consequence of this assumption, the corresponding shell thickness for the k-isotope

 *Ek Eff* 

*<sup>i</sup>*=<sup>1</sup> *piEi* ∑*<sup>N</sup> <sup>i</sup>*=<sup>1</sup> *Ei*

*CSDA* along with the corresponding graphic sketch (Figure 1) outlining

*CSDA* used [*cm*]

(15)

(16)

*CSDA* = *RCSDA*

*Eff* <sup>=</sup> <sup>∑</sup>*<sup>N</sup>*

Where *N* is the number of channels of the *k*-isotope energy spectrum and *pi* and *Ei* are its weight and energy respectively. It can be seen, as reported in Table 1, two examples of the

how *EEff* may be estimated. Radionuclide emission spectra were extracted from the NIST

*Eff* calculated [*keV*] *<sup>R</sup><sup>k</sup>*

90Y 936.9 0.4029 177Lu 138.9 0.0246 89Sr 583.6 0.2182 153Sm 231.3 0.0565 131I 188.2 0.0407

*Eff* calculated for different isotopes.

The Monte Carlo method has been successfully applied to radiation transport problems in clinical dosimetry. The method is particularly useful for complex problems that can not be

solved analytically and/or deterministically or when measurements are not practical.

from internal MC database or NIST tables without distinction.

*Rk*

where *k* =90Y, 131I, 89Sr, 153Sm, 177Lu, ... and *EEff* is computed according to:

*Ek*

shell thickness for the DPK calculations.

*Eff* and *<sup>R</sup><sup>k</sup>*

Table 1. *RCSDA* used and *E<sup>k</sup>*

<sup>1</sup> http://www.nist.gov

<sup>2</sup> Idem 1.

Isotope *E<sup>k</sup>*

**2.4 Dedicated subroutines based on Monte Carlo techniques**

calculated *E<sup>k</sup>*

database tables2.

(*RCSDA*) can be calculated from the relation:

Fig. 1. 90Y and 177Lu spectra compared with its *EEff* calculated.

The accuracy of Monte Carlo simulations strongly depends on the accuracy in the probability density functions, and thereby the material cross-section libraries used for radiation transport calculations.

The PENELOPE parameter *Eabs* determines -for each simulated material and transported particle type- the energy threshold for particle absorption, defining the energy value below which transport is no longer simulated and residual energy is assumed to be deposited locally. The lowest possible *Eabs* value is 100*eV* according to manual recommendations based on the lower energy (50*eV*) in PENELOPE v. 2008 database.

With the aim of performing radiation transport for dosimetric purposes, the first step may be the design and development of a code devoted to the DPK calculations, which may basically consist on a main program able to allow the user to set up the energy or the spectrum to be radiated by the point source centered in a geometry defined by concentric spheres defining shells of thickness *RCSDA*/40 (the corresponding *RCSDA* for each material is calculated for the actual or effective energy for monoenertic sources or radionuclides, respectively). Due to phisical properties of charged particles, it may be enough to tally DPK up to 1.5*RCSDA*. Finally, the material for which the DPK needs to be investigated along with the corresponding simulation parameters may be adjusted by the user in the input file.

The geometry file can be obtained from a dedicated user-designed program devoted to set up the *RCSDA* obtained from MC or reliable databases. Once the user provides *RCSDA* the program will automatically generate the corresponding file according to internal PENELOPE requirements for simulation geometries.

Preliminary simulations are performed to achieve *Eabs* parameter values in order to optimize computation process assessing suitable balances between statistical errors and the CPU computation times. Therefore, after exhaustive examinations, absorption cutoff values of 100*eV* (minimum recommendable by PENELOPE developers) have been considered for emitting sources of energies lower than 100*keV* (*i.e.* 5*keV*, 10*keV*, 20*keV* and 50*keV*) whereas cutoff values of 1*keV* have been set for emitting sources with energies equal or higher than 100*keV* (*i.e.* 100*keV*, 200*keV*, 500*keV*, 1*MeV*, 2*MeV* and 5*MeV*).

As example from the preliminary simulations devoted to achieving suitable absorption energy cutoff valued, it can be mentioned that in the case of the 100*keV* monoenergetic source, it was found a maximum difference in the DPK of 0.16% when changing cutoff value from 1*keV* to the minimum possible of 100*eV*, which may be taken as negligible. However, the situation is significantly different when dealing with low energy emitting sources. DPK differences up to 5% may be found when changing absorption cutoff values from 1*keV* to 100*eV* in the case of 5*keV* monoenergetic source.

Fig. 4. Convergence analysis for the 25*th* (a), 35*th* (b), 45*th* (c) and 55*th* (d) shell. Error bars

The distinction between the scattering and the primary contributions by the developed code is assessed by "marking" or "painting" each particle on the main code distinguishing energy deposition arising from primary or scattered particles. It should be emphasized that this subtlety is able because of the great -actually absolute- aperture of PENELOPE as an open source code. Technically, this aim can be accomplished by adding a special new parameter in the particle quantum state along with the inclusion of suitable flags during particle tracking. With this tool it becomes possible to obtain the dose delivered by each particle while

Dosimetry for Beta-Emitter Radionuclides by Means of Monte Carlo Simulations 273

Therefore, dose contributions calculation is on-run performed and separated files are

The first step consists on developing a code capable of reading and interpreting patient-specific medical images in order to use them as input data for further computation

According to typical medical imaging techniques, the most commonly found image storage formats consist on space voxelization. In this sense, it would be desiderable to have at disposal some mechanism devoted to interpret and handle voxelized information in order to perform patient-specific dosimetry calculations. If the MC package does not provide a voxelization toolkit, it would be useful and/or necessary that the user may implement some dedicated voxelization technique in order to allow the use of patient-specific data from DICOM images within the simulation. The role of the voxelization technique regards the interpretation by the simulation code of the patient data in order to provide the corresponding inputs for simulation geometry and mass distribution. Some MC main codes provide voxelization capability with user distribution, but unfortunately this is not the case of PENELOPE v.2008. Therefore, it was necessary to develop a dedicated program for patient data upload. This aim was accomplished by considering patient geometry as a tensor whose elements correspond to a patient image voxel and therefore its material and mass density may be inferred from CT

generated for total, primary and scattering dose delivering tallied at each shell.

calibration, which provides the Hounsfield (CT) index vs. electronic density.

correspond to 3*σ*, *i.e.* 99.9% of confidence.

discriminating each kind of contribution.

with the simulation code.

**2.4.2 Monte Carlo direct simulation on patient images**

2.4.2.1 Geometry voxelization in Monte Carlo methods

As reported in Figure 2, it can be shown that DPK curves do not exhibit strong differences when calculated for 100 or 1000eV considering emission energy of 100*keV*.

Fig. 2. Total sDPK curves for different cutoffs (100*eV*, 1*keV* and 5*keV*)

#### **2.4.1 Monte Carlo DPK calculation**

It is considered a punctual source (spectral or monoenergetic) centered in a water sphere divided in many shells. Shell thicknesses are defined according to *RCSDA*/40 and dose delivered within each shell is tallied using the developed MC code system, based on physics from PENELOPE v. 2008.

In addition, dedicated subroutines have been introduced in the MC main code with the aim of computing separately the different contributions to the total absorbed dose. Suitable distinction may be made between primary and scattering components, considering scattering contribution as arising from all energy depositions belonging to non primary particles.

Also, it is necessary to establish an optimal number of showers to be simulated in view of optimizing the compromise between short CPU calculation time and accurate enough results. In this sense, it may be useful to carry out preliminary calculations with different number of showers to assess a suitable choice, according to Figures 3 and 4, for example. In base of these figures it may be suggest that 107 primary showers would be an acceptable selection. Actually, it should be usefull to add that according to this choice the CPU time is not longer than 8 hours for a typical microprocessor on a standard PC.

Fig. 3. Convergence analysis of the 15*th* shell for different numbers of initial showers tallying deposited energy along with uncertainties calculated as three standard deviations.

8 Will-be-set-by-IN-TECH

As reported in Figure 2, it can be shown that DPK curves do not exhibit strong differences

It is considered a punctual source (spectral or monoenergetic) centered in a water sphere divided in many shells. Shell thicknesses are defined according to *RCSDA*/40 and dose delivered within each shell is tallied using the developed MC code system, based on physics

In addition, dedicated subroutines have been introduced in the MC main code with the aim of computing separately the different contributions to the total absorbed dose. Suitable distinction may be made between primary and scattering components, considering scattering contribution as arising from all energy depositions belonging to non primary particles. Also, it is necessary to establish an optimal number of showers to be simulated in view of optimizing the compromise between short CPU calculation time and accurate enough results. In this sense, it may be useful to carry out preliminary calculations with different number of showers to assess a suitable choice, according to Figures 3 and 4, for example. In base of these figures it may be suggest that 107 primary showers would be an acceptable selection. Actually, it should be usefull to add that according to this choice the CPU time is not longer

Fig. 3. Convergence analysis of the 15*th* shell for different numbers of initial showers tallying

deposited energy along with uncertainties calculated as three standard deviations.

when calculated for 100 or 1000eV considering emission energy of 100*keV*.

Fig. 2. Total sDPK curves for different cutoffs (100*eV*, 1*keV* and 5*keV*)

than 8 hours for a typical microprocessor on a standard PC.

**2.4.1 Monte Carlo DPK calculation**

from PENELOPE v. 2008.

Fig. 4. Convergence analysis for the 25*th* (a), 35*th* (b), 45*th* (c) and 55*th* (d) shell. Error bars correspond to 3*σ*, *i.e.* 99.9% of confidence.

The distinction between the scattering and the primary contributions by the developed code is assessed by "marking" or "painting" each particle on the main code distinguishing energy deposition arising from primary or scattered particles. It should be emphasized that this subtlety is able because of the great -actually absolute- aperture of PENELOPE as an open source code. Technically, this aim can be accomplished by adding a special new parameter in the particle quantum state along with the inclusion of suitable flags during particle tracking. With this tool it becomes possible to obtain the dose delivered by each particle while discriminating each kind of contribution.

Therefore, dose contributions calculation is on-run performed and separated files are generated for total, primary and scattering dose delivering tallied at each shell.

#### **2.4.2 Monte Carlo direct simulation on patient images**

The first step consists on developing a code capable of reading and interpreting patient-specific medical images in order to use them as input data for further computation with the simulation code.

#### 2.4.2.1 Geometry voxelization in Monte Carlo methods

According to typical medical imaging techniques, the most commonly found image storage formats consist on space voxelization. In this sense, it would be desiderable to have at disposal some mechanism devoted to interpret and handle voxelized information in order to perform patient-specific dosimetry calculations. If the MC package does not provide a voxelization toolkit, it would be useful and/or necessary that the user may implement some dedicated voxelization technique in order to allow the use of patient-specific data from DICOM images within the simulation. The role of the voxelization technique regards the interpretation by the simulation code of the patient data in order to provide the corresponding inputs for simulation geometry and mass distribution. Some MC main codes provide voxelization capability with user distribution, but unfortunately this is not the case of PENELOPE v.2008. Therefore, it was necessary to develop a dedicated program for patient data upload. This aim was accomplished by considering patient geometry as a tensor whose elements correspond to a patient image voxel and therefore its material and mass density may be inferred from CT calibration, which provides the Hounsfield (CT) index vs. electronic density.

(a) (b)

Dosimetry for Beta-Emitter Radionuclides by Means of Monte Carlo Simulations 275

Fig. 5. CT (a), and SPECT (b) images showing both a transversal cut of the geometry and the

Over the two past decades, the linear-quadratic (LQ) model has emerged as a convenient tool to quantify biological effects for radiotherapy. Therefore, just for didactics purposes ilustrating how to proceed for implementing radiobiological calculations, the LQ model will be taken as example. From a general point of view and without interest on going deeper, it should be mentioned that the LQ model provides a simple method to calculate the cell survival fraction (*SF*), and on its consequences it is possible to see a reliable measure of the

Assuming the Poisson statistics to survival cells and considering the control probability of the tumor as the probability of having them death, it can be assessed the LQ model, assuming a dose *D* homogeneously distributed over a homogeneous tumor, establishes that the *SF* is

where *α* and *β* are empiric biological parameters associated to the probability of breaking DNA helix. The *α* parameter represents the probability to induce a double strand break in DNA helix with a single energy deposition, whereas the *β* parameter represents the probability to have two single strand breaks with two separate energy depositions, one on the first DNA helix, the other on the second one, which are close enough -both spatially and temporally- to combine together producing the same effect of a double strand break, *i.e.* a

> *SF* <sup>=</sup> *Ns N*<sup>0</sup>

where *N*<sup>0</sup> refers to the initial quantity of tumoral cells and *Ns* to the survival tumoral cells.

<sup>−</sup>*SF* <sup>=</sup> exp

Whereas NTCP formulation is considered by the interpretation of empirical clinical data, since the organ behavior cannot be easily described by equations, and also the sparing of a normal tissue is not simply linked to the number of surviving cells, but also to the maintenance of its

The Lyman-Kutcher-Burman (LKB) [V.A. Semenenko (2008), Z. Xu & Jiang (2006)] model

−*e*

<sup>−</sup>*αD*−*βD*<sup>2</sup>

<sup>−</sup>*αD*−*βD*<sup>2</sup>

(17)

(18)

(19)

*SF* = *e*

activity information.

given by the relationship

lethal (non reparable) damage.

Then, the TCP is given by the equation

Also, eq. (17) considers *SF* and it is defined by the equation,

*TCP* = *e*

functionality, which depends on the organ architecture.

proposes that the NTCP may be assessed by

TCP.

Once patient geometry is already voxelized, 3D energy and dose deposition are tallied at each voxel.

2.4.2.2 Input data files from DICOM images

DICOM is a commonly used imaging file format in almost all medical imaging fields, particularly in nuclear medicine. For the purposes of the present work it was developed a specific code for reading these image files and further convertion to suitable ascii files able to be introduced as input for the developed radiation transport code.

One comfortable way for reading and interpreting DICOM images is available in the MATLAB® environment. The implemented method consists on exploiting the useful MATLAB package called "dicomread" that allows the user to read DICOM files as 3D tensors. The obtained tensor is used to generate a suitable array which elements refer to the information at each voxel in the DICOM file. This array is further uploaded as input ascii data providing the required patient-specific radiation transport properties for MC simulations.

2.4.2.3 Metabolic & anatomical patient-specific input data

As described above, analogous treatments apply for uploading DICOM functional or anatomical images. Patient anatomy is used to define simulation geometry and physical properties may be inferred by means of suitable calibrations or implementing some "*ad hoc*" criterion consulting expert physicians. On the other hand, functional images, like SPECT or PET are used to provide activity spatial distribution and knowing the corresponding infused radionuclide spectrum, it may be possible to use this information as simulation radiation source. Combining both of these imaging techniques as simulation inputs, it would be possible to assess reliable patient-specific dosimetry.

The DICOM functional images provide the activity information which is interpreted by the code to establish the source distribution and its activity everywhere.

2.4.2.4 Specific-patient dose distribution calculation

Once input data are already uploaded by means of the developed subroutines aimed to converting DICOM images into simulation geometry and radiation source, it becomes possible to perform the corresponding MC simulation to the actual patient-specific situation. As reported in Figure 5, the DICOM images are successfully interpreted to establish both a *N*-material (four different materials in the reported example, namely: compact bone, soft tissue, air and lung have been selected as a simplified first approximation) patient-specific geometry. In addition, activity intensity and spatial distribution are also satisfactory interpreted and combined with the corresponding anatomical data show a suitable definition of the patient body as a virtual phantom for MC calculations.

### **2.5 Radiobiological models**

Once 3D dose distribution is assessed, radiobiological considerations can be derived by applying the adequate models with this aim.

Most calculations of the biological effect of radiation on tumors assume that the clonogenic cell density is uniform even if account is taken of non-uniform dose distribution. Anyway, after implementing some methodology for to identifying tumor/normal tissues on patient images it will be necessary to introduce a model capable of assesing *Tumor Control Probability* (TCP) on the tumor [Dale (1988),D.J. Brenner & Sachs (1995),M. Guerrero (2004)] and *Normal Tissue Complications Probability* (NTCP) over the normal tissue and organs closed to the treatment region [V.A. Semenenko (2008),L.A. Dawson & Haken (2002)].

10 Will-be-set-by-IN-TECH

Once patient geometry is already voxelized, 3D energy and dose deposition are tallied at each

DICOM is a commonly used imaging file format in almost all medical imaging fields, particularly in nuclear medicine. For the purposes of the present work it was developed a specific code for reading these image files and further convertion to suitable ascii files able to

One comfortable way for reading and interpreting DICOM images is available in the MATLAB® environment. The implemented method consists on exploiting the useful MATLAB package called "dicomread" that allows the user to read DICOM files as 3D tensors. The obtained tensor is used to generate a suitable array which elements refer to the information at each voxel in the DICOM file. This array is further uploaded as input ascii data providing the required patient-specific radiation transport properties for MC simulations.

As described above, analogous treatments apply for uploading DICOM functional or anatomical images. Patient anatomy is used to define simulation geometry and physical properties may be inferred by means of suitable calibrations or implementing some "*ad hoc*" criterion consulting expert physicians. On the other hand, functional images, like SPECT or PET are used to provide activity spatial distribution and knowing the corresponding infused radionuclide spectrum, it may be possible to use this information as simulation radiation source. Combining both of these imaging techniques as simulation inputs, it would be

The DICOM functional images provide the activity information which is interpreted by the

Once input data are already uploaded by means of the developed subroutines aimed to converting DICOM images into simulation geometry and radiation source, it becomes possible to perform the corresponding MC simulation to the actual patient-specific situation. As reported in Figure 5, the DICOM images are successfully interpreted to establish both a *N*-material (four different materials in the reported example, namely: compact bone, soft tissue, air and lung have been selected as a simplified first approximation) patient-specific geometry. In addition, activity intensity and spatial distribution are also satisfactory interpreted and combined with the corresponding anatomical data show a suitable definition

Once 3D dose distribution is assessed, radiobiological considerations can be derived by

Most calculations of the biological effect of radiation on tumors assume that the clonogenic cell density is uniform even if account is taken of non-uniform dose distribution. Anyway, after implementing some methodology for to identifying tumor/normal tissues on patient images it will be necessary to introduce a model capable of assesing *Tumor Control Probability* (TCP) on the tumor [Dale (1988),D.J. Brenner & Sachs (1995),M. Guerrero (2004)] and *Normal Tissue Complications Probability* (NTCP) over the normal tissue and organs closed to the treatment

voxel.

2.4.2.2 Input data files from DICOM images

be introduced as input for the developed radiation transport code.

2.4.2.3 Metabolic & anatomical patient-specific input data

possible to assess reliable patient-specific dosimetry.

2.4.2.4 Specific-patient dose distribution calculation

code to establish the source distribution and its activity everywhere.

of the patient body as a virtual phantom for MC calculations.

region [V.A. Semenenko (2008),L.A. Dawson & Haken (2002)].

**2.5 Radiobiological models**

applying the adequate models with this aim.

Fig. 5. CT (a), and SPECT (b) images showing both a transversal cut of the geometry and the activity information.

Over the two past decades, the linear-quadratic (LQ) model has emerged as a convenient tool to quantify biological effects for radiotherapy. Therefore, just for didactics purposes ilustrating how to proceed for implementing radiobiological calculations, the LQ model will be taken as example. From a general point of view and without interest on going deeper, it should be mentioned that the LQ model provides a simple method to calculate the cell survival fraction (*SF*), and on its consequences it is possible to see a reliable measure of the TCP.

Assuming the Poisson statistics to survival cells and considering the control probability of the tumor as the probability of having them death, it can be assessed the LQ model, assuming a dose *D* homogeneously distributed over a homogeneous tumor, establishes that the *SF* is given by the relationship

$$SF = e^{-aD - \beta D^2} \tag{17}$$

where *α* and *β* are empiric biological parameters associated to the probability of breaking DNA helix. The *α* parameter represents the probability to induce a double strand break in DNA helix with a single energy deposition, whereas the *β* parameter represents the probability to have two single strand breaks with two separate energy depositions, one on the first DNA helix, the other on the second one, which are close enough -both spatially and temporally- to combine together producing the same effect of a double strand break, *i.e.* a lethal (non reparable) damage.

Also, eq. (17) considers *SF* and it is defined by the equation,

$$SF = \frac{N\_{\rm s}}{N\_0} \tag{18}$$

where *N*<sup>0</sup> refers to the initial quantity of tumoral cells and *Ns* to the survival tumoral cells. Then, the TCP is given by the equation

$$T\mathbb{C}P = e^{-SF} = \exp\left(-e^{-aD-\beta D^2}\right) \tag{19}$$

Whereas NTCP formulation is considered by the interpretation of empirical clinical data, since the organ behavior cannot be easily described by equations, and also the sparing of a normal tissue is not simply linked to the number of surviving cells, but also to the maintenance of its functionality, which depends on the organ architecture.

The Lyman-Kutcher-Burman (LKB) [V.A. Semenenko (2008), Z. Xu & Jiang (2006)] model proposes that the NTCP may be assessed by

Fig. 6. sDPK calculated for all the monoenergetic beams simulated.

Fig. 7. All monoenergetic beams simulated on log-log scale normalized to maximum,

Fig. 8. sDPK for 100*keV* showing total, primary and scattering contributions.

reported in Figures 11 and 12.

expressed in length units [cm] highlighting differences for a wide range of emission energies.

Dosimetry for Beta-Emitter Radionuclides by Means of Monte Carlo Simulations 277

commonly used radionuclides for targeted radiotherapy treatments have been studied with the developed calculation code. The obtained results for 90Y, 131I, 89Sr, 153Sm and 177Lu are

$$\text{NTCP} \equiv c(u) = \frac{1}{\sqrt{2\pi}} \int\_{\infty}^{u} e^{-t^2/2} dt \tag{20}$$

where

$$
\mu = \frac{D - TD\_{50}}{m \cdot TD\_{50}} \tag{21}
$$

when conditions are assumed as in the LQ model for the TCP calculations, *m* is a dimensionless correction factor and *TD*<sup>50</sup> is the dose, over all the tumor, for that NTCP is equal to 50%. And it is possible to characterize *c*(*u*) as

$$\begin{aligned} \mathcal{L}(u) &= \frac{1}{2} \exp\left(\kappa u - \frac{\kappa^2 u^2}{2}\right) \qquad \text{for} \quad u < 0\\ \mathcal{L}(u) &= 1 - \frac{1}{2} \exp\left(\kappa u - \frac{\kappa^2 u^2}{2}\right) \quad \text{for} \quad u \ge 0 \end{aligned} \tag{22}$$

where *κ* is a dimensionless parameter approximated by 0.8154.

#### **3. Application of patient-specific dosimetry system**

According to dosimetry calculation purposes, DPKs are all what is needed to calculate dose distributions in homogeneous media starting from activity distributions by simply means of analytical convolution. However, it can be shown that in cases of inhomogeneous media, these analytical methods may not provide accurate enough results for dose distribution calculations, as expected.

Therefore, this section will present some results about the calculation of DPK both for monoenergetic and spectral *β*− sources within homogeneous as well as inhomogeneous media with the aim of pointing out the corresponding differences.

In addition, results for preliminary dose maps along with TCP and NTCP calculations for patient-specific images will be reported, the last ones will be considered as examples of the wide range of possibilities in radiobiology offered by the versatility of the implemented calculation system based on 3D dose distribution.

#### **3.1 Monoenergetic beta-emitting sources scaled DPK**

Beta-minus monoenergetic sources were simulated and their sDPK were calculated for initial primary particles energies of 5*keV*, 10*keV*, 20*keV*, 50*keV*, 100*keV*, 200*keV*, 500*keV*, 1*MeV*, 2*MeV* and 5*MeV*. The obtained results are shown in Figures 6 and 7.

As mentioned, the implemented calculation code incorporates dedicated routines devoted to quantify the different dosimetric contributions. Figures 8 and 9 show some obtained results for the total sDPK along with the corresponding separation between the primary and scattering contributions for different monoenergetic sources.

In addition, it is interesting to report the relative contribution from scattering component to the total sDPK at different distances, as shown in Figure 10.

#### **3.2 Spectral beta-emitting sources Dose Point Kernels**

With the aim of preliminary benchmarking the capability of the developed calculation code in the framework of clinical implementation of beta-minus radiolabeled therapy, sDPK at different energies were calculated for homogeneous media. Therefore, some of the most 12 Will-be-set-by-IN-TECH

*<sup>u</sup>* <sup>=</sup> *<sup>D</sup>* <sup>−</sup> *TD*<sup>50</sup> *m* · *TD*<sup>50</sup>

when conditions are assumed as in the LQ model for the TCP calculations, *m* is a dimensionless correction factor and *TD*<sup>50</sup> is the dose, over all the tumor, for that NTCP is

*<sup>κ</sup><sup>u</sup>* <sup>−</sup> *<sup>κ</sup>*2*u*<sup>2</sup> 2

According to dosimetry calculation purposes, DPKs are all what is needed to calculate dose distributions in homogeneous media starting from activity distributions by simply means of analytical convolution. However, it can be shown that in cases of inhomogeneous media, these analytical methods may not provide accurate enough results for dose distribution

Therefore, this section will present some results about the calculation of DPK both for monoenergetic and spectral *β*− sources within homogeneous as well as inhomogeneous

In addition, results for preliminary dose maps along with TCP and NTCP calculations for patient-specific images will be reported, the last ones will be considered as examples of the wide range of possibilities in radiobiology offered by the versatility of the implemented

Beta-minus monoenergetic sources were simulated and their sDPK were calculated for initial primary particles energies of 5*keV*, 10*keV*, 20*keV*, 50*keV*, 100*keV*, 200*keV*, 500*keV*, 1*MeV*,

As mentioned, the implemented calculation code incorporates dedicated routines devoted to quantify the different dosimetric contributions. Figures 8 and 9 show some obtained results for the total sDPK along with the corresponding separation between the primary and

In addition, it is interesting to report the relative contribution from scattering component to

With the aim of preliminary benchmarking the capability of the developed calculation code in the framework of clinical implementation of beta-minus radiolabeled therapy, sDPK at different energies were calculated for homogeneous media. Therefore, some of the most

*f or u <* 0

*f or u* ≥ 0

*<sup>κ</sup><sup>u</sup>* <sup>−</sup> *<sup>κ</sup>*2*u*<sup>2</sup> 2

<sup>√</sup>2*<sup>π</sup>*

 *u* ∞ *e* −*t*

2/2*dt* (20)

(21)

(22)

*NTCP* <sup>≡</sup> *<sup>c</sup>*(*u*) = <sup>1</sup>

equal to 50%. And it is possible to characterize *c*(*u*) as

*<sup>c</sup>*(*u*) = <sup>1</sup>

*<sup>c</sup>*(*u*) = <sup>1</sup> <sup>−</sup> <sup>1</sup>

**3. Application of patient-specific dosimetry system**

calculations, as expected.

where *κ* is a dimensionless parameter approximated by 0.8154.

media with the aim of pointing out the corresponding differences.

2*MeV* and 5*MeV*. The obtained results are shown in Figures 6 and 7.

scattering contributions for different monoenergetic sources.

the total sDPK at different distances, as shown in Figure 10.

**3.2 Spectral beta-emitting sources Dose Point Kernels**

calculation system based on 3D dose distribution.

**3.1 Monoenergetic beta-emitting sources scaled DPK**

<sup>2</sup> exp

<sup>2</sup> exp

where

Fig. 6. sDPK calculated for all the monoenergetic beams simulated.

Fig. 7. All monoenergetic beams simulated on log-log scale normalized to maximum, expressed in length units [cm] highlighting differences for a wide range of emission energies.

Fig. 8. sDPK for 100*keV* showing total, primary and scattering contributions.

commonly used radionuclides for targeted radiotherapy treatments have been studied with the developed calculation code. The obtained results for 90Y, 131I, 89Sr, 153Sm and 177Lu are reported in Figures 11 and 12.

Fig. 12. Total radionuclide DPK normalized to global maximum pointing out the differences

Dosimetry for Beta-Emitter Radionuclides by Means of Monte Carlo Simulations 279

in penetration capabilities.

Fig. 13. Different contributions to the 90Y sDPK.

Fig. 14. Different contributions to the 177Lu sDPK.

**3.3 Scaled Dose Point Kernels for non water-equivalent media**

As mentioned above, it is commonly found in clinical applications the requirement of considering different materiales for patient dosimetry. In this sense, it may be useful to

Fig. 9. sDPK for 1*MeV* showing total, primary and scattering contributions.

Fig. 10. sDPK along with relative scattering contribution for 500*keV* monoenergetic source.

Fig. 11. Total sDPK for different *β*− radionuclides.

In the case of radionuclides, the different contributions to the total delivered dose have been assessed by means of an analogue procedure to that of the monoenergetic sources. Figures 13 and 14 report the obtained results for the contributions to total sDPK allowing to appreciate qualitative different behaviours according to each radionuclide.

14 Will-be-set-by-IN-TECH

Fig. 9. sDPK for 1*MeV* showing total, primary and scattering contributions.

Fig. 10. sDPK along with relative scattering contribution for 500*keV* monoenergetic source.

In the case of radionuclides, the different contributions to the total delivered dose have been

Figures 13 and 14 report the obtained results for the contributions to total sDPK allowing to

assessed by means of an analogue procedure to that of the monoenergetic sources.

appreciate qualitative different behaviours according to each radionuclide.

Fig. 11. Total sDPK for different *β*− radionuclides.

Fig. 12. Total radionuclide DPK normalized to global maximum pointing out the differences in penetration capabilities.

Fig. 13. Different contributions to the 90Y sDPK.

Fig. 14. Different contributions to the 177Lu sDPK.

#### **3.3 Scaled Dose Point Kernels for non water-equivalent media**

As mentioned above, it is commonly found in clinical applications the requirement of considering different materiales for patient dosimetry. In this sense, it may be useful to

reference medium and bone or lung inhomogeneities have been introduced considering wide

Dosimetry for Beta-Emitter Radionuclides by Means of Monte Carlo Simulations 281

Figures 16 and 17 report the obtained results for the above described cases. It can be appreciated that, at least *a priori*, it seems to be an anomalous behaviour at the interfaces between different materials. It has been obtained that this "peculiar" situation is always found for all the considered cases independently of inhomogeneity type, thickness or relative

ranges of inhomogeneity thickness and relative position within the irradiated sphere.

Fig. 16. sDPK corresponding to different relative positions of a 6-shell ( <sup>6</sup>

Fig. 17. sDPK corresponding to different bone inhomogeneity thicknesses.

radionuclides, as reported in Figure 18 for 131I.

**3.5 Patient-specific dose distribution along with TCP and NTCP calculations**

Once both the anatomical and metabolic images were already uploaded to the developed calculation code it became possible to perform the corresponding patient-specific dosimetry establishing the desidered radionuclide and corresponding activity. As example of the proposed method in practical cases, dual SPECT-CT images have been used to carry out the simulation. Dose distribution calculations have been performed for different beta-minus

Although relative dose distribution have been presented, it should be emphasized that once the source activity is established, it becomes possible to determine the corresponding number

<sup>40</sup> · *RCSDA*) thick

position.

bone inhomogeneity.

investigate the effect on sDPK due to the consideration of irradiated media constituted of materials typically used in clinics other than water or normal soft tissue.

In order to show the capability of the developed calculation code to perform radiation transport in different materials, two typical clinical tissues have been selected for this investigation. Bone (defined according to International Commission on Radiation Units and Measurements, ICRU) and lung (International Commission on Radiological Protection, ICRP) tissues have been considered because they represent somehow the extreme cases in typical clinical situations. Figures 14 and 15 report he obtained sDPK calculated over bone and lung spheres for 100*keV* and 500*keV* sources. As a consequence of the sDPK suitable dimensionless length normalization, it has been obtained that the sDPK qualitatively preserve the same shape disregarding of the specific considered medium.

Fig. 15. Contributions to the sDPK calculated for 100*keV* (a) and 500*keV* (b) centered at a bone and lung spheres, respectively.

#### **3.4 Scaled Dose Point Kernels on inhomogeneous media**

In order to highlight the effects of inserting an inhomogeneity in the calculation of the sDPK different situations have been suitably designed and investigated. Water has been selected as 16 Will-be-set-by-IN-TECH

investigate the effect on sDPK due to the consideration of irradiated media constituted of

In order to show the capability of the developed calculation code to perform radiation transport in different materials, two typical clinical tissues have been selected for this investigation. Bone (defined according to International Commission on Radiation Units and Measurements, ICRU) and lung (International Commission on Radiological Protection, ICRP) tissues have been considered because they represent somehow the extreme cases in typical clinical situations. Figures 14 and 15 report he obtained sDPK calculated over bone and lung spheres for 100*keV* and 500*keV* sources. As a consequence of the sDPK suitable dimensionless length normalization, it has been obtained that the sDPK qualitatively preserve the same

(a)

(b)

In order to highlight the effects of inserting an inhomogeneity in the calculation of the sDPK different situations have been suitably designed and investigated. Water has been selected as

Fig. 15. Contributions to the sDPK calculated for 100*keV* (a) and 500*keV* (b) centered at a

bone and lung spheres, respectively.

**3.4 Scaled Dose Point Kernels on inhomogeneous media**

materials typically used in clinics other than water or normal soft tissue.

shape disregarding of the specific considered medium.

reference medium and bone or lung inhomogeneities have been introduced considering wide ranges of inhomogeneity thickness and relative position within the irradiated sphere.

Figures 16 and 17 report the obtained results for the above described cases. It can be appreciated that, at least *a priori*, it seems to be an anomalous behaviour at the interfaces between different materials. It has been obtained that this "peculiar" situation is always found for all the considered cases independently of inhomogeneity type, thickness or relative position.

Fig. 16. sDPK corresponding to different relative positions of a 6-shell ( <sup>6</sup> <sup>40</sup> · *RCSDA*) thick bone inhomogeneity.

Fig. 17. sDPK corresponding to different bone inhomogeneity thicknesses.

#### **3.5 Patient-specific dose distribution along with TCP and NTCP calculations**

Once both the anatomical and metabolic images were already uploaded to the developed calculation code it became possible to perform the corresponding patient-specific dosimetry establishing the desidered radionuclide and corresponding activity. As example of the proposed method in practical cases, dual SPECT-CT images have been used to carry out the simulation. Dose distribution calculations have been performed for different beta-minus radionuclides, as reported in Figure 18 for 131I.

Although relative dose distribution have been presented, it should be emphasized that once the source activity is established, it becomes possible to determine the corresponding number

been extracted from literature [G. Luxton & King (2008)], namely *TD*<sup>50</sup> = 24.5 and *m* = 0.18.

Dosimetry for Beta-Emitter Radionuclides by Means of Monte Carlo Simulations 283

Fig. 20. NTCP for the slice of interest obtained for 1*mCi* activity delivered by 131I.

The adequacy of the developed radiation transport simulation code based on PENELOPE v.2008 for nuclear medicine and focused on beta-emitters dosimetry has been preliminary investigated checking the consistency of the results expected when considering electron

In this sense, dose point kernels have been calculated for several monoenergetic and spectral sources including the most commonly used radionuclides for typical nuclear medicine

Contrary to other works, it was considered the *RCSDA* range for both monoenergetic and spectrum sources, instead of the *X*<sup>90</sup> parameter used commonly on spectrum sources. Therefore it could be not straightforward to compare the obtained results for radionuclides with those reported by other authors. However, it does not constitute a significant drawback because comparisons can be performed for monoenergetic emission sources, which may help

It is interesting to comment about the obtained results for sDPK when inserting inhomogeneities. It has been found that non negligible abrupt changes appear in sDPK curves poducing remarkable discontinuity on the sDPK derivate. This behaviour may make difficult or actually avoid adequate analytical calculation of sDPK, as could be the case of realistic clinical situations. In this sense, it constitutes one of the main advantages of developing and

Technical features of simulation process have been studied in detail including convergence of representative observable mean values and standard deviations. In view of the specific requirements for most of the calculation performed in this work, it has been found that 107 primary showers seems to be a convenient and suitable compromise between calculation time

Regarding the investigations about the separation and characerization of the different contributions to the total absorbed dose, it should be mentioned that the main contribution comes from primary radiation, as expected. Although clearly non negligible, scattering contributions are always significantly lower than the primary component to absorbed dose. In addition, it has been shown how to obtain the different dose contributions at different

implementing MC methods devoted to reliable radiation transport computation.

The obtained NCTP distribution is presented in Figure 20.

**4. General overview and final remarks**

therapies.

and accurate enough results.

transport in the energy range from 5*keV* to 5*MeV*.

for assessing the reliability and accuracy of the proposed method.

Fig. 18. Relative dose distribution on virtual patient irradiated by 131I.

of primary showers and therefore obtaining straightforwardly the corresponding absolute dose distribution. This information will be necessary for the calculation of radiobiological quantities.

The 3D absolute dose distribution constitutes the basic required information for implementing suitable algorithms devoted to estimate the TCP and NTCP. As mentioned previously, the linear quadratic model has been implemented for TCP estimation. As example of how to proceed, it has been arbitrary taken a total infused activity of 1mCi and typical values for the LQ radiobiological parameters have been used according to the literature, namely *α* = 0.12*Gy*−<sup>1</sup> and *β* = 0.0137*Gy*−<sup>2</sup> given a relationship of *α*/*β* = 8.7591. Figure 19 reports the obtained TCP distribution for the slice of interest.

Fig. 19. TCP for the slice of interest obtained for 1*mCi* activity delivered by 131I.

However, it is necessary to remark that in real clinical situations TCP and NTCP estimation will require the implementation of some mechanism dedicated to the identification of normal/tumor tissue regions in order to proceed with the assigment of the corresponding radiobiological parameters. Actually, TCP evaluations have to be applied to tumor region, whereas NTCP evaluations apply to normal organs, each having different radiobiological parameter values. Once the previous step is accomplished the calculation method is performed following the procedure described above.

Similarly, it can be proposed some model for calculating the NTCP. With the aim of highlighting the required procedure, typical radiobiological values for the LKB model have been extracted from literature [G. Luxton & King (2008)], namely *TD*<sup>50</sup> = 24.5 and *m* = 0.18. The obtained NCTP distribution is presented in Figure 20.

Fig. 20. NTCP for the slice of interest obtained for 1*mCi* activity delivered by 131I.
