- Email: [email protected]
A code for hadrontherapy treatment planning with the voxelscan method
Computers in Biology and Medicine 30 (2000) 311±327
www.elsevier.com/locate/compbiomed
A code for hadrontherapy treatment planning with the voxelscan method S. Berga a, F. Bourhaleb c, R. Cirio a, J. Derkaoui c, B. Gallice a, M. Hamal c, F. Marchetto a, V. Rolando b,*, S. Viscomi a a Universita di Torino e Sezione INFN, Italy Facolta di Medicina e Chirurgia, Universita del Piemonte Orientale, 28100 Novara, Italy c Faculty of Science, Mohamed I University, Oujda, Morocco
b
Received 2 August 1999; accepted 1 May 2000
Abstract A code for the implementation of treatment plannings in hadrontherapy with an active scan beam is presented. The package can determine the ¯uence and energy of the beams for several thousand voxels in a few minutes. The performances of the program have been tested with a full simulation. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Hadrons; Voxelscan method; Inverse problem; Bragg peak; Ion fragmentation
1. Introduction There are several centers worldwide where hadrons (protons or ions) are used for radiotherapy. Essentially two methods for dose delivering are used: (a) passive, where the beam is spread out by an absorber and eventually collimated to match the tumour shape; (b) active, where a narrow beam is steered to aim at only a small fraction of the tumour volume (voxel ) at a time. To fully exploit the active system, it is necessary to control precisely the beam direction in the transverse plane as well as the beam energy for the depth modulation. It is then mandatory to couple the accelerator to a precise and fast feedback system. Indeed once a
* Corresponding author. Tel.: +39-0321-660-651; fax: +39-0321-620-421. E-mail address: [email protected] (V. Rolando). 0010-4825/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 0 - 4 8 2 5 ( 0 0 ) 0 0 0 1 7 - 2
312
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
voxel has been irradiated with the correct dose, the beam has to be stopped, steered at the next voxel and the energy changed before the treatment can continue. On the other hand, any given voxel is irradiated not only when the beam is aiming at it, but also during the treatment of voxels farther along the same direction. Thus, the treatment planning for each voxel has to provide the correct number of hadrons to be delivered, their energy and direction. Summing up the contributions of each elementary fraction of the treatment, one has to obtain the required dose and shape. Once the beam data have been calculated, the treatment planning requires a further step which is the treatment veri®cation. This is done by simulating the beam throughout the target. The results of the simulation assess the correctness of the analytical calculation. This paper presents an analytical code to compute, for both protons and ions, the beam energy, direction and the hadron yield for each voxel, given the volume to be treated. Section 2 describes the code and the basic assumptions. Section 3 presents in detail the algorithm, in particular, the method used to correct the tail due to the fragmentation in the ion case. The results for two applications: (a) a water box, (b) a cranial irradiation, are discussed in Section 4. The output of the treatment planning calculation has been fully simulated with a Monte Carlo simulation. In this way the consistency of the calculation is checked in an almost independent way. Finally, Section 5 presents the conclusions. 2. Description of the code 2.1. Generalities The code, called ANCOD (ANalytical CODe), is based on two basic assumptions: . The dose delivery to the target is done in elementary steps. For each step, a volume element (voxel) of the target is irradiated. . The beam is described as a mono-energetic pencil beam aiming at the voxel center, and the energy is computed so as to have the Bragg peak right at the center of the voxel. The voxels, dimensions and positions, coincide with the X-ray Computed Tomography (CT) voxels on which the target is de®ned.1 The irradiation starts from the voxels of the transverse section farthest from the beam source. Moving from one voxel to the other in the same section requires the change of beam direction and energy. The typical angle deviation is of the order of 1 mrad, while the typical energy step is less than 1 MeV. Both variations can be achieved with limited adjustments of the accelerator magnetic ®elds, and are consequently rather fast. Once all the voxels of a given section have been treated, the beam is aimed at the next closer section by decreasing the energy. This time the required energy step has to match the voxel size and it is typically of a few MeV. The treatment ends when all the voxels within the target have been irradiated. Two righthanded reference systems have been introduced:
1
Indeed this choice is not strictly required, but we make it for the sake of simplicity.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
313
. Target Reference System (TRS) which is related the CT volume, with x-axis parallel to the sagittal, y-axis parallel to the coronal and the z-axis parallel to the axial tomography slice.2 . Beam Reference System (BRS) which is related to the beam, with the z-axis along the beam direction and the y-axis pointing upwards. The code runs on Unix platforms and is written in Fortran 77. 2.2. Input data The input data to the code are the following: . Beam parameters: ®rst, one has to specify the number of treatment ®elds, and for each ®eld, it is necessary to de®ne the nominal beam direction and the source coordinates in the BRS. Here, we de®ne the nominal beam direction as the direction of a straight line at the beam center when there is no de¯ection due to the horizontal and vertical de¯ection magnets. The source is the point (real or virtual) about which the beam pivots during the voxelscan treatment. In the calculation, we do not assume any beam cross-section. When we verify the treatment, we use a double-Gaussian shape with the same width in both directions (in the plane transverse to the beam direction). Finally, the mass and charge of the beam particles complete the input beam parameters. . CT data ®le: the pattern of energy deposition along the voxels, as they are crossed by the beam, depends on the traversed media. Thus, together with the beam parameters, one has to input the CT volume information, i.e. the dimensions, material type and electronic density for each voxel. Since there is no unique standard for the output of CT machines, it is necessary to translate the actual output to the program input format. . The required 3D dose distribution ®le: ®nally, one has to de®ne the volume to be irradiated. Since the target is fully contained in the CT volume, all the parameters of the target voxels are passed to the program with the CT data ®le. What is necessary here is to ¯ag each voxel as belonging or not belonging to the target. This is done with a 3D matrix (D(i, j, k )) ®lled up with ¯ags: +1 for the voxels belonging to the target; 0 for the voxels outside the target. The assignment has to be done by the physician on each CT slice.
2.3. Output data The output of the code is stored in two separate ®les: . The main ®le contains the treatment planning information: for each voxel of the target, the beam energy is stored, along with the number of hadrons, and the deviation of the beam with respect to the nominal beam direction. In a simpli®ed view, where at the gantry output the beam is steered along two perpendicular directions, these deviations are necessary to 2
In human anatomy, sagittal, coronal, and axial planes are de®ned as follows: sagittal plane: every plane parallel to the plane of symmetry; transverse plane (horizontal or axial): every plane orthogonal to the axis of symmetry; frontal plane (coronal): every plane orthogonal to sagittal and axial planes.
314
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
compute the dipole magnetic ®elds to aim at a given voxel. Indeed these kinds of data are those that are needed to control the treatment planning in the accelerator and the feedback system. Furthermore, the main ®le keeps a record of the relevant geometrical information, such as the relative position and orientation of the two reference systems, and the number of voxels and their dimensions along each BRS axes. . The log®le contains the main execution information such as the total dose in each slice and the total number of hadrons.
3. Calculation 3.1. Introduction The aim of the calculation is to optimize the radiation ®elds with respect to the treatment plan requirements. On one hand, one has to deliver a uniform dose to the target, as de®ned by the physician. On the other, one has to minimize the dose deposited outside the target. The parameters are: the number of irradiation ®elds and their orientation with respect to the target, the number and the energy of the hadrons for each spot. The decision on the number and the orientation of the ®elds is outside the scope of this paper and will not be discussed. Given the ®eld parameters, we approach the problem of determining the beam parameters for the ®eld under study. This is a two-step problem: ®rst, we determine the beam energy value necessary to place the Bragg peak at the center of each spot to be irradiated. This is rather straightforward, and can be done from the CT numbers. Second, one has to compute the ¯uences to be delivered at each spot. This is a complicated task which can be attacked in several ways. Let us consider a target with I voxels along the x-axis, J along y and K along z. The total number of voxels is given by N I J K: The complete treatment plan is the linear superposition of N elementary treatments, each with a beam ¯uence of wl. We indicate with d li, j, k the dose deposited in the voxel i, j, k by the elementary treatment l. Thus, the total dose deposited in the voxel i, j, k is given by X wl d li, j, k 1 D i, j, k l
This is called the inverse problem: D i, j, k DM , where DM is the prescribed dose for i, j, k within the tumour region, and D i, j, k 0 for i, j, k outside the tumour region. The matrix d li, j, k has been derived with a complete simulation of the energy deposition in a water phantom. Indeed it depends on the beam energy and on the depth of the voxel i, j, k under examination. The ¯uences wl are the unknowns. To solve the above equation, it is necessary to invert the matrix d li, j, k , which poses several major diculties due to: (a) the number of elements, which can easily be of the order of 104; (b) the fact that one can have negative values for some beam ¯uences, wl, which are obviously non-physical. This paper presents an algorithm to overcome the inverse problem and to calculate the ¯uences with an iterative process.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
315
Table 1 Range±energy relationship for protons and carbon ions obtained with GEANT C+6 ions
Protons Energy (MeV)
Range (g/cm2)
Energy (MeV)
Range (g/cm2)
10.0 15.0 20.0 25.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 125.0 150.0 175.0 200.0 250.0 300.0
0.1230 0.2539 0.4260 0.6370 0.8853 1.4890 2.2270 3.0930 4.0800 5.1840 6.3980 7.7180 11.4600 15.7700 20.6200 25.9600 37.9400 51.4500
14.1 17.4 22.0 27.0 33.0 41.0 50.0 62.0 76.0 94.0 116.0 143.0 150.0 159.0 177.0 180.0 200.0 220.0
0.075 0.109 0.166 0.236 0.337 0.510 0.712 1.077 1.539 2.260 3.286 4.750 5.167 5.715 6.871 7.083 8.471 9.974
3.2. Beam energy calculation First, we describe the algorithm that we used for the beam energy computation. As input we used an energy±range table for each type of hadron. This has been obtained using a full Monte Carlo simulation (GEANT [1]) to compute the energy deposition in a water phantom as a function of depth. Here, we de®ne as range the depth where the maximum energy deposition occurs. An example of the range table is given in Table 1 for proton and carbon ions. The same data are represented in Fig. 1 both for protons (full dots) and for carbons ions (empty dots). A given elementary treatment operation selects the irradiation ®eld, the beam source, and the voxel under treatment. Assuming a straight line for the beam trajectory, the path length, li, in a given voxel i of the target is uniquely determined. We transformed li into the water equivalent path length, l we i , by multiplying li by ri , where ri is the voxel density. This is an approximation, which holds only for densities close to that of water. P The total path length inside the target in terms of water equivalent is given by Lwe l we i , where the sum runs over all the voxels intercepted by the beam line.3 3
In the calculation, the dose is computed normalizing the energy deposition to the length of the beam path inside any voxel. For the simulation, the geometry describes precisely the target and the beam `followed' step by step along the path to compute the energy deposition.
316
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Thus, the beam energy, Eb, is given by the following linear interpolation: Eb Enÿ1
En ÿ Enÿ1 Lwe ÿ Rnÿ1 , Rn ÿ Rnÿ1
2
where Rn and Rnÿ1 are the range values from Table 1 around Lwe, and En and Enÿ1 the corresponding energies. 3.3. Fluence calculation A few assumptions are necessary to simplify the computation, and in particular, we assume a pencil-like beam. This implies that, when irradiating a given voxel, say i, the dose deposited in one of the adjacent voxels in the transverse plane, say i ÿ 1, due to the beam dimensions, is equal to the dose deposited in i when the beam is aiming at i ÿ 1: Furthermore, the transverse spot dimensions are smaller than the beam dimensions.4 The input to the ¯uence computation is the energy deposition, De Eb , z=Dz, of a pencil beam of a given energy, Eb, as a function of the depth, z. With the pencil beam assumption, one can ignore the dependence of the energy deposition on the plane transverse to the beam. We construct a table of De Eb , z=Dz for several values of Eb ranging from 100 to 400 MeV/ u for carbon ions and from 60 to 250 MeV for protons.5 For all the energies, Dz has been chosen equal to 0.1 cm. The table has been obtained with a full simulation of the energy deposition in a water phantom. We used the GEANT package which describes energy losses of the protons including the nuclear interactions. For carbon, the fragmentation processes have been taken into account using the total fragmentation cross-section quoted by Sihver [2]. For the inclusive partial cross-section for a single-channel production, we used the parametrization of Sihver et al. [3]. A section of the energy deposition table for C ions is reported in Table 2. The comparison between the energy loss as predicted by the simulation and the experimental results [4] for a carbon ion beam of 270 MeV/u is shown in Fig. 2. The ®rst step of the computation is to determine the energy deposition in each voxel crossed by the beam in a given elementary treatment. We remark that the computation of the ¯uence for the carbon ion beam is complicated by the energy tail behind the Bragg peak due to the fragmentation debris. We ®rst start by describing the algorithm ignoring this problem. The ®rst step is to determine the energy, Et, from Table 1, closest to the beam energy, Eb, for a given voxel l. In general, Et is dierent from Eb, thus the expected range Rt associated with Et diers from Lwe. We recall that Lwe is the range required to place the Bragg peak in the middle of voxel l. We then add an amount of water equivalent equal to Rt ÿ Lwe (subtract if Rt ÿ Lwe is negative) in front of the ®rst voxel along the beam path. This is equivalent to cut or to extend 4
The uniformity of the dose distribution in the transverse plane depends on the relative dimensions of the beam width and the step size. In general, as the step size decreases, the resulting convoluted beams become ¯atter. Assuming a beam with a Gaussian shape, a beam width (sigma) greater or equal to than the step size results in a uniformity better than 1%. 5 The energy interval we chose was 10 MeV/u, which corresponds to less than 1 cm in the range variation.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
317
Fig. 1. Energy±range relationship in water for protons (full dots) and for carbon ions (empty dots).
the plateau region by an amount Rt ÿ Lwe , assuming that the distribution of dE/dx is ¯at, as a function of z, for z close to zero. Indeed Rt ÿ Lwe is typically smaller than 0.1 cm and the approximation in this range is quite good. With the above correction, we make sure that the Bragg peak is centered in the voxel under treatment. The next step is straightforward: it consists of the computation of the correct energy deposition in each voxel crossed by the beam path. As input we have the depths z1, at the beam entrance, and z2, at the exit of a voxel. The depths z1 and z2 have to be mapped into the table of the energy deposition to compute the energy deposited between z1 and z2. Let us assume that the energy deposited in the most downstream voxel, say l, is El. The dose, dl, due to a single hadron is given by dl kEl =lwe ,
3
where k 1:6 10 ÿ9 Gy/MeV cm. Thus, the ¯uence wl is given by DM/dl. As one moves away from the most downstream section, the dose previously delivered has to
318
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Fig. 2. Energy loss of carbon ions of 270 MeV/u of kinetic energy as from the simulation (empty dots) and the GSI experimental results (full dots).
be accounted for. Explicitly, when treating the voxel l ÿ 1 the ¯uence wlÿ1 is given by =dlÿ1 , wlÿ1 DM ÿ wl d lÿ1 l
4
where E llÿ1 is the dose deposited in the voxel l ÿ 1 while treating voxel l. Applying the previous equation recursively, one can easily determine all the ¯uences. Explicitely, the required ¯uence for voxel i is given by ! X i wi DM ÿ wj d j =di , 5 ji1, l
where d ij is the dose deposited in the voxel i during the irradiation of the voxel j.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
319
Table 2 Linear energy deposited vs. depth for a carbon ion with 270 MeV-initial kinetic energy Depth (cm)
DE (MeV/cm)
10.6132 11.3207 11.6038 11.8868 12.1698 12.4528 12.7358 13.0189 13.3019 13.5141 13.6556 13.7972 13.9387 14.0802 14.1085 14.1368 14.1651
176.539 188.825 195.194 202.748 211.921 223.247 237.533 256.248 282.32 310.563 337.538 377.005 445.514 642.503 721.204 836.563 787.847
The proton treatment algorithm ends here, while, as stated above, the ion treatment requires a further correction due to the energy tail behind the Bragg peak. Indeed during the treatment of voxel l ÿ 1, the tail goes through voxel l, for which the ¯uence was, in fact, already
Fig. 3. Layout of the target irradiation of the water box.
320
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Fig. 4. Layout of the target irradiation with four ®elds in the XY plane.
computed. Thus, one needs to correct it by decreasing the ¯uence by the appropriate amount. The exact solution may involve an extremely lengthy calculation. Nevertheless, we notice that the tail is characterized by a height which is of the order of 5% of the plateau. One can immediately conclude that any solution at 10% level brings down the eect to less than 1%, which is acceptable. Given the spot dimensions, the beam energy, and the target size, the ratios of the ¯uences wi =wl , to obtain a ¯at dose within the target, are totally predictable and can be stored in appropriate tables. When calculating wl, the ¯uence relative to the voxel l belonging to the most downstream section, one needs to subtract the expected dose due to later treatment of all the preceeding voxels based on the expected ¯uence ratios. The same applies to all the voxels.
4. Application of the code and Monte Carlo veri®cation 4.1. Generalities As stated in Section 1, it is necessary to verify the correctness of the calculation. We have performed the treatment plan veri®cation with a full simulation. The beam parameters
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
321
Fig. 5. Energy of carbon ions vs. depth (a) and ¯uence of carbon ions vs. depth (b) in the water box phantom.
(direction and energy) and the ¯uence, as calculated by ANCOD, are the input to the simulation which uses the standard package GEANT [1].6 The beam width (rms) has been set to 0.4 cm. We have done the check of the treatment plan for two cases: 1. A water box, where the target was a parallelepiped of dimensions 4:05 4:03 4:06 cm3, with voxel dimensions of 0:15 0:13 0:14 cm3 for 24,273 voxels. The target was placed downstream of 10.46 cm of water. We used only one ®eld of irradiation, and the beam source was set at the following coordinates: x 1:95 cm, y 1:95 cm, and z ÿ200 cm [5]. Fig. 3 shows a sketch of the set up.
6
The structure of GEANT is extremely ¯exible to describe volumes of dierent characteristics and to follow the electromagnetic and hadronic interactions of the particles, parents and daughters, along the track path.
322
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Fig. 6. Monte Carlo simulation of the depth (a) and transverse (b) dose pro®le in the water box phantom.
2. A cranial CT with 98 sections, each section split into 134 voxels along the x direction and 113 voxels along the y direction.7 Each voxel was a cube of size 0.1594 cm. The target was chosen as a cube of 21 voxels along each side, placed in the middle of the cranial CT. We used four coplanar ®elds at 908 with respect to each other. The ®eld plane was perpendicular to the axial axis and the point source was placed 04 m apart from the target. In Fig. 4, we show a sketch of the setup. So far we have discussed a situation where BRS can be obtained from TRS with a set of coordinate displacements. We remark that for ®elds not parallel to the sagittal or coronal plane or, in other words, BRS rotated with respect to TRS, it is necessary to modify consequently the CT input. This is done by rotating the voxel 7
The CT numbers were already translated into water equivalent thicknesses.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
323
Fig. 7. Dose distribution in the water box phantom.
planes to be parallel or orthogonal to BRS axes and mapping the voxel densities from the old to the new reference system.8 5. Results First, we consider case (1) (water box) for a C+6 beam. Fig. 5(a) shows the beam energy as a function of the voxel position along the z coordinate. For a given direction, the beam energy required to have the Bragg peak centered in consecutive sections increases with depth. The 8 The mapping is done as follows: in TRS, a point x, y, z within the voxel i, j, k with given CT number corresponds, in BRS, to the point x 0 , y 0 , z 0 within the voxel i 0 , j 0 , k 0 ), and with the same CT number. By cycling over enough points (i.e. 10 per voxel), one can obtain a very accurate reconstruction of the CT numbers of the voxels in the rotated BRS.
324
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Fig. 8. Cranial case: dose distribution due to a single irradiation ®eld.
average energy step is 02 MeV/u for a voxel size of 0.14 cm. In the simulation, to cope with present day accelerator performances, the energy input values were rounded to 1 MeV. Thus, for the example shown in Fig. 5(a), the required number of energy steps is about 50. Fig. 5(b) shows the ¯uence per dose unit with the same horizontal scale as in Fig. 5(a). As one moves from the farthest sections to those closest to the beam entrance, the required ¯uences decrease by a factor of 050. For small ¯uences, one can observe ¯uctuations (of the order of 20.5%) which are due to the approximations used. However, to evaluate the performance of the code, one has to compare the dose deposited in every single voxel with the required dose. With the GEANT simulation, we computed the dose deposited in the voxel i, j, k), dp i, j, k: The required dose in the target, DM, is expected to be constant over the whole volume. Fig. 6(a) shows the dose pro®le in a voxel row moving along the sections. Fig. 6(b) shows the dose pro®le in a voxel row moving across the target volume in a plane perpendicular to the beam. In Fig. 6(b), we observe a depletion of the deposited dose at the borders of the target which aects mainly the border voxels. The eect is clearly due to the pencil beam approximation in ANCOD. Indeed the contributions to the dose in a voxel in
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
325
the middle of the target are due to the direct irradiation of the voxel itself and dose spill out while irradiating the neighbor voxels. For a voxel placed at the border, the latter contribution is partially missing resulting in a lack of dose. We checked the simulation results against the beam width, and found similar results for beam width smaller than 4 mm (sigma). The pencil beam approximation does not hold above this ®gure, but typical values of the beam width which is expected to be used in a voxelscan treatment is smaller than 4 mm. While examining the dose pro®le in Fig. 6(b), one can observe that the algorithm is quite satisfactory over the whole volume. One can introduce the variable as a measure of the dose uniformity s X 2 ÿ 2 , dp i, j, k ÿ DM = N DM hdFi i, j, k
Fig. 9. Cranial CT case: dose distribution due to a 4-®eld irradiation.
6
326
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
where N i j k is the total number of voxels of the target. Here hdFi, which is a normalized rms, gives directly the average deviation in terms of fraction of the required dose DM. We found hdFi 3%: Finally, Fig. 7 shows a two-dimensional view of the deposited dose for a middle section of the target. We now examine case (2) (cranial CT). Concerning the beam energy and ¯uences for each ®eld, results are very similar to the water box. Fig. 8 shows the plot of the deposited dose (in arbitrary units) in a section at the center of the target for a single ®eld irradiation. The dose in the same section is shown in Fig. 9 after the four ®elds have been applied. It is remarkable that with four ®elds one obtains a ratio between the dose on the target and the dose at the plateau close to 05. Furthermore, the target edges are rather sharp: the dose decreases from 90 to 10% of the target value in 20.4 cm, which is equivalent to two voxel widths. Also for this case, the uniformity on the target is hdFi 3%: Thus, for homogeneous and inhomogeneous media we obtained comparable results with respect to the dose uniformity. A substantial contribution to the hdFi value is due to the Monte Carlo ¯uctuation since the fraction of simulated events is limited.
6. Conclusions In this paper, we have presented a package for hadrontherapy treatment planning with the voxelscan method. The software is conceived for both protons and ions, and it is based upon two types of tables: (a) range-energy and (b) energy loss as a function of beam energy and depth. Both set of tables have been derived with a full simulation of the energy losses in a water phantom. Furthermore, the consistency of the planning has been veri®ed with a full simulation of the treatment using the GEANT package. The dose uniformity obtained over the target is within 3% for both homogeneous and inhomogeneous media.
References [1] R. Brun et al., GEANT Ð Detector Description and Simulation, Tool version 3, User's Guide, CERN DD/EE/ 84-1, 1987. [2] L. Sihver, C.H. Tsao, R. Silberberg, T. Kanai, A.F. Barghouty, Total reaction and partial cross section calculations in proton±nucleus Zt R26 and nucleus±nucleus reactions (Zp and Zt R26), Phys. Rev. C 47 (1993) 1225. [3] L. Sihver et al., A model for calculating depth-dose and ¯uence distributions, GSI Scienti®c Report 95-1, p. 221, March 1995. [4] C. Brusasco, (GSI) private communication, December 1998. [5] B. Schaner, (PSI) private communication, May 1996. Stefano Berga received his degree in Physics in 1997 from the University of Torino, Italy. His work aimed at modelling and implementing nucleus fragmentation in the energy range from 400 MeV/u down to few MeV/u. He is presently teaching Physics at the High School ``Liceo Scienti®co dell'Istituto Sociale'', Torino. Faiza Bourhaleb received her License in Theoretical Physics in 1996 from University Mohamed First of Oujda, Morocco. In accomplishing her graduate studies, she is working in particle physics on Grand Uni®ed Theories (GUT and SUSY) and their
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
327
consequences in cosmology. Since December 1997 she has been working on Hadrontherapy and cell survival models. Her research interests include Medical imaging, simulation in Biomedical Physiscs. Roberto Cirio received his Physics degree in March 1983 at the Physics Department of the Torino University. He has worked at CERN, on the e+eÿ collider LEP with the DELPHI detector. In 1992, he joined the ZEUS collaboration, at the e+ÿ/p collider HERA at DESY. Since 1998, he has been working on the CMS detector for the pp collider LHC at CERN. He has also been working on development of detectors and related items to be used in medical physics since 1992. M. Jamal E. Derkaoui received his Ph.D. degree from the Pierre et Marie Curie University (Paris 6, France) in May 1984. He has been working as a Professor of Physics at the Mohamed 1st University of Oujda (Morocco) since October 1984. His present research activity are mainly focused on hadrontherapy, particles physics and astrophysics. Barbara Gallice received her degree in Physics in 1998 from the University of Torino, Italy. Her main work, for the thesis, was devoted to develop a code for treatment planning in the hadron therapy framework. She is presently teaching Physics at the High School ``Istituto Tecnico I. Porro'' of Pinerolo, Italy. Mohamed Hamal received his M.D. (CEA) degree from Ibn Tofail University of Kenitra (Morocco) in 1994. He is presently in the last year of his Ph.D. studies at the Mohamed 1st University of Oujda (Morocco). He is studying the interaction between hadrons particles and the biological tissues and the eectiveness of hadronthearpy with respect to various paramaters as dose, energy, cell characteristics. F. Marchetto received his degree in Physics in 1973 from the University of Torino. Since then he has been working as a Research Associate at the Istituto Nazionale di Fisica Nucleare (INFN). The activities have been primarily focused on the ®eld of High Energy Physics from basic research to detector developments. He is author of numerous publications on several international journals. Valter Rolando received his degree in Physics from the University of Torino in 1994. At present he is working as System and Network Manager at the Faculty of Medicine and Surgery, University of Piemonte Orientale in Novara, Italy. He has worked in medical imaging and simulation. Silvia Viscomi received her degree in Physics in 1998 from the University of Torino, Italy. Her work was focused on treatment planning algorithms in Hadrontherapy. In 1999, she attended the School of Specialization in Medical Physics at University of Genoa. Presently, she is awailing a fellowship at the Cancer Epidemiology Unit of the University of Torino, where she is doing studies on survival and prognostic factors for childhood cancer.
www.elsevier.com/locate/compbiomed
A code for hadrontherapy treatment planning with the voxelscan method S. Berga a, F. Bourhaleb c, R. Cirio a, J. Derkaoui c, B. Gallice a, M. Hamal c, F. Marchetto a, V. Rolando b,*, S. Viscomi a a Universita di Torino e Sezione INFN, Italy Facolta di Medicina e Chirurgia, Universita del Piemonte Orientale, 28100 Novara, Italy c Faculty of Science, Mohamed I University, Oujda, Morocco
b
Received 2 August 1999; accepted 1 May 2000
Abstract A code for the implementation of treatment plannings in hadrontherapy with an active scan beam is presented. The package can determine the ¯uence and energy of the beams for several thousand voxels in a few minutes. The performances of the program have been tested with a full simulation. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Hadrons; Voxelscan method; Inverse problem; Bragg peak; Ion fragmentation
1. Introduction There are several centers worldwide where hadrons (protons or ions) are used for radiotherapy. Essentially two methods for dose delivering are used: (a) passive, where the beam is spread out by an absorber and eventually collimated to match the tumour shape; (b) active, where a narrow beam is steered to aim at only a small fraction of the tumour volume (voxel ) at a time. To fully exploit the active system, it is necessary to control precisely the beam direction in the transverse plane as well as the beam energy for the depth modulation. It is then mandatory to couple the accelerator to a precise and fast feedback system. Indeed once a
* Corresponding author. Tel.: +39-0321-660-651; fax: +39-0321-620-421. E-mail address: [email protected] (V. Rolando). 0010-4825/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 0 - 4 8 2 5 ( 0 0 ) 0 0 0 1 7 - 2
312
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
voxel has been irradiated with the correct dose, the beam has to be stopped, steered at the next voxel and the energy changed before the treatment can continue. On the other hand, any given voxel is irradiated not only when the beam is aiming at it, but also during the treatment of voxels farther along the same direction. Thus, the treatment planning for each voxel has to provide the correct number of hadrons to be delivered, their energy and direction. Summing up the contributions of each elementary fraction of the treatment, one has to obtain the required dose and shape. Once the beam data have been calculated, the treatment planning requires a further step which is the treatment veri®cation. This is done by simulating the beam throughout the target. The results of the simulation assess the correctness of the analytical calculation. This paper presents an analytical code to compute, for both protons and ions, the beam energy, direction and the hadron yield for each voxel, given the volume to be treated. Section 2 describes the code and the basic assumptions. Section 3 presents in detail the algorithm, in particular, the method used to correct the tail due to the fragmentation in the ion case. The results for two applications: (a) a water box, (b) a cranial irradiation, are discussed in Section 4. The output of the treatment planning calculation has been fully simulated with a Monte Carlo simulation. In this way the consistency of the calculation is checked in an almost independent way. Finally, Section 5 presents the conclusions. 2. Description of the code 2.1. Generalities The code, called ANCOD (ANalytical CODe), is based on two basic assumptions: . The dose delivery to the target is done in elementary steps. For each step, a volume element (voxel) of the target is irradiated. . The beam is described as a mono-energetic pencil beam aiming at the voxel center, and the energy is computed so as to have the Bragg peak right at the center of the voxel. The voxels, dimensions and positions, coincide with the X-ray Computed Tomography (CT) voxels on which the target is de®ned.1 The irradiation starts from the voxels of the transverse section farthest from the beam source. Moving from one voxel to the other in the same section requires the change of beam direction and energy. The typical angle deviation is of the order of 1 mrad, while the typical energy step is less than 1 MeV. Both variations can be achieved with limited adjustments of the accelerator magnetic ®elds, and are consequently rather fast. Once all the voxels of a given section have been treated, the beam is aimed at the next closer section by decreasing the energy. This time the required energy step has to match the voxel size and it is typically of a few MeV. The treatment ends when all the voxels within the target have been irradiated. Two righthanded reference systems have been introduced:
1
Indeed this choice is not strictly required, but we make it for the sake of simplicity.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
313
. Target Reference System (TRS) which is related the CT volume, with x-axis parallel to the sagittal, y-axis parallel to the coronal and the z-axis parallel to the axial tomography slice.2 . Beam Reference System (BRS) which is related to the beam, with the z-axis along the beam direction and the y-axis pointing upwards. The code runs on Unix platforms and is written in Fortran 77. 2.2. Input data The input data to the code are the following: . Beam parameters: ®rst, one has to specify the number of treatment ®elds, and for each ®eld, it is necessary to de®ne the nominal beam direction and the source coordinates in the BRS. Here, we de®ne the nominal beam direction as the direction of a straight line at the beam center when there is no de¯ection due to the horizontal and vertical de¯ection magnets. The source is the point (real or virtual) about which the beam pivots during the voxelscan treatment. In the calculation, we do not assume any beam cross-section. When we verify the treatment, we use a double-Gaussian shape with the same width in both directions (in the plane transverse to the beam direction). Finally, the mass and charge of the beam particles complete the input beam parameters. . CT data ®le: the pattern of energy deposition along the voxels, as they are crossed by the beam, depends on the traversed media. Thus, together with the beam parameters, one has to input the CT volume information, i.e. the dimensions, material type and electronic density for each voxel. Since there is no unique standard for the output of CT machines, it is necessary to translate the actual output to the program input format. . The required 3D dose distribution ®le: ®nally, one has to de®ne the volume to be irradiated. Since the target is fully contained in the CT volume, all the parameters of the target voxels are passed to the program with the CT data ®le. What is necessary here is to ¯ag each voxel as belonging or not belonging to the target. This is done with a 3D matrix (D(i, j, k )) ®lled up with ¯ags: +1 for the voxels belonging to the target; 0 for the voxels outside the target. The assignment has to be done by the physician on each CT slice.
2.3. Output data The output of the code is stored in two separate ®les: . The main ®le contains the treatment planning information: for each voxel of the target, the beam energy is stored, along with the number of hadrons, and the deviation of the beam with respect to the nominal beam direction. In a simpli®ed view, where at the gantry output the beam is steered along two perpendicular directions, these deviations are necessary to 2
In human anatomy, sagittal, coronal, and axial planes are de®ned as follows: sagittal plane: every plane parallel to the plane of symmetry; transverse plane (horizontal or axial): every plane orthogonal to the axis of symmetry; frontal plane (coronal): every plane orthogonal to sagittal and axial planes.
314
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
compute the dipole magnetic ®elds to aim at a given voxel. Indeed these kinds of data are those that are needed to control the treatment planning in the accelerator and the feedback system. Furthermore, the main ®le keeps a record of the relevant geometrical information, such as the relative position and orientation of the two reference systems, and the number of voxels and their dimensions along each BRS axes. . The log®le contains the main execution information such as the total dose in each slice and the total number of hadrons.
3. Calculation 3.1. Introduction The aim of the calculation is to optimize the radiation ®elds with respect to the treatment plan requirements. On one hand, one has to deliver a uniform dose to the target, as de®ned by the physician. On the other, one has to minimize the dose deposited outside the target. The parameters are: the number of irradiation ®elds and their orientation with respect to the target, the number and the energy of the hadrons for each spot. The decision on the number and the orientation of the ®elds is outside the scope of this paper and will not be discussed. Given the ®eld parameters, we approach the problem of determining the beam parameters for the ®eld under study. This is a two-step problem: ®rst, we determine the beam energy value necessary to place the Bragg peak at the center of each spot to be irradiated. This is rather straightforward, and can be done from the CT numbers. Second, one has to compute the ¯uences to be delivered at each spot. This is a complicated task which can be attacked in several ways. Let us consider a target with I voxels along the x-axis, J along y and K along z. The total number of voxels is given by N I J K: The complete treatment plan is the linear superposition of N elementary treatments, each with a beam ¯uence of wl. We indicate with d li, j, k the dose deposited in the voxel i, j, k by the elementary treatment l. Thus, the total dose deposited in the voxel i, j, k is given by X wl d li, j, k 1 D i, j, k l
This is called the inverse problem: D i, j, k DM , where DM is the prescribed dose for i, j, k within the tumour region, and D i, j, k 0 for i, j, k outside the tumour region. The matrix d li, j, k has been derived with a complete simulation of the energy deposition in a water phantom. Indeed it depends on the beam energy and on the depth of the voxel i, j, k under examination. The ¯uences wl are the unknowns. To solve the above equation, it is necessary to invert the matrix d li, j, k , which poses several major diculties due to: (a) the number of elements, which can easily be of the order of 104; (b) the fact that one can have negative values for some beam ¯uences, wl, which are obviously non-physical. This paper presents an algorithm to overcome the inverse problem and to calculate the ¯uences with an iterative process.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
315
Table 1 Range±energy relationship for protons and carbon ions obtained with GEANT C+6 ions
Protons Energy (MeV)
Range (g/cm2)
Energy (MeV)
Range (g/cm2)
10.0 15.0 20.0 25.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 125.0 150.0 175.0 200.0 250.0 300.0
0.1230 0.2539 0.4260 0.6370 0.8853 1.4890 2.2270 3.0930 4.0800 5.1840 6.3980 7.7180 11.4600 15.7700 20.6200 25.9600 37.9400 51.4500
14.1 17.4 22.0 27.0 33.0 41.0 50.0 62.0 76.0 94.0 116.0 143.0 150.0 159.0 177.0 180.0 200.0 220.0
0.075 0.109 0.166 0.236 0.337 0.510 0.712 1.077 1.539 2.260 3.286 4.750 5.167 5.715 6.871 7.083 8.471 9.974
3.2. Beam energy calculation First, we describe the algorithm that we used for the beam energy computation. As input we used an energy±range table for each type of hadron. This has been obtained using a full Monte Carlo simulation (GEANT [1]) to compute the energy deposition in a water phantom as a function of depth. Here, we de®ne as range the depth where the maximum energy deposition occurs. An example of the range table is given in Table 1 for proton and carbon ions. The same data are represented in Fig. 1 both for protons (full dots) and for carbons ions (empty dots). A given elementary treatment operation selects the irradiation ®eld, the beam source, and the voxel under treatment. Assuming a straight line for the beam trajectory, the path length, li, in a given voxel i of the target is uniquely determined. We transformed li into the water equivalent path length, l we i , by multiplying li by ri , where ri is the voxel density. This is an approximation, which holds only for densities close to that of water. P The total path length inside the target in terms of water equivalent is given by Lwe l we i , where the sum runs over all the voxels intercepted by the beam line.3 3
In the calculation, the dose is computed normalizing the energy deposition to the length of the beam path inside any voxel. For the simulation, the geometry describes precisely the target and the beam `followed' step by step along the path to compute the energy deposition.
316
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Thus, the beam energy, Eb, is given by the following linear interpolation: Eb Enÿ1
En ÿ Enÿ1 Lwe ÿ Rnÿ1 , Rn ÿ Rnÿ1
2
where Rn and Rnÿ1 are the range values from Table 1 around Lwe, and En and Enÿ1 the corresponding energies. 3.3. Fluence calculation A few assumptions are necessary to simplify the computation, and in particular, we assume a pencil-like beam. This implies that, when irradiating a given voxel, say i, the dose deposited in one of the adjacent voxels in the transverse plane, say i ÿ 1, due to the beam dimensions, is equal to the dose deposited in i when the beam is aiming at i ÿ 1: Furthermore, the transverse spot dimensions are smaller than the beam dimensions.4 The input to the ¯uence computation is the energy deposition, De Eb , z=Dz, of a pencil beam of a given energy, Eb, as a function of the depth, z. With the pencil beam assumption, one can ignore the dependence of the energy deposition on the plane transverse to the beam. We construct a table of De Eb , z=Dz for several values of Eb ranging from 100 to 400 MeV/ u for carbon ions and from 60 to 250 MeV for protons.5 For all the energies, Dz has been chosen equal to 0.1 cm. The table has been obtained with a full simulation of the energy deposition in a water phantom. We used the GEANT package which describes energy losses of the protons including the nuclear interactions. For carbon, the fragmentation processes have been taken into account using the total fragmentation cross-section quoted by Sihver [2]. For the inclusive partial cross-section for a single-channel production, we used the parametrization of Sihver et al. [3]. A section of the energy deposition table for C ions is reported in Table 2. The comparison between the energy loss as predicted by the simulation and the experimental results [4] for a carbon ion beam of 270 MeV/u is shown in Fig. 2. The ®rst step of the computation is to determine the energy deposition in each voxel crossed by the beam in a given elementary treatment. We remark that the computation of the ¯uence for the carbon ion beam is complicated by the energy tail behind the Bragg peak due to the fragmentation debris. We ®rst start by describing the algorithm ignoring this problem. The ®rst step is to determine the energy, Et, from Table 1, closest to the beam energy, Eb, for a given voxel l. In general, Et is dierent from Eb, thus the expected range Rt associated with Et diers from Lwe. We recall that Lwe is the range required to place the Bragg peak in the middle of voxel l. We then add an amount of water equivalent equal to Rt ÿ Lwe (subtract if Rt ÿ Lwe is negative) in front of the ®rst voxel along the beam path. This is equivalent to cut or to extend 4
The uniformity of the dose distribution in the transverse plane depends on the relative dimensions of the beam width and the step size. In general, as the step size decreases, the resulting convoluted beams become ¯atter. Assuming a beam with a Gaussian shape, a beam width (sigma) greater or equal to than the step size results in a uniformity better than 1%. 5 The energy interval we chose was 10 MeV/u, which corresponds to less than 1 cm in the range variation.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
317
Fig. 1. Energy±range relationship in water for protons (full dots) and for carbon ions (empty dots).
the plateau region by an amount Rt ÿ Lwe , assuming that the distribution of dE/dx is ¯at, as a function of z, for z close to zero. Indeed Rt ÿ Lwe is typically smaller than 0.1 cm and the approximation in this range is quite good. With the above correction, we make sure that the Bragg peak is centered in the voxel under treatment. The next step is straightforward: it consists of the computation of the correct energy deposition in each voxel crossed by the beam path. As input we have the depths z1, at the beam entrance, and z2, at the exit of a voxel. The depths z1 and z2 have to be mapped into the table of the energy deposition to compute the energy deposited between z1 and z2. Let us assume that the energy deposited in the most downstream voxel, say l, is El. The dose, dl, due to a single hadron is given by dl kEl =lwe ,
3
where k 1:6 10 ÿ9 Gy/MeV cm. Thus, the ¯uence wl is given by DM/dl. As one moves away from the most downstream section, the dose previously delivered has to
318
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Fig. 2. Energy loss of carbon ions of 270 MeV/u of kinetic energy as from the simulation (empty dots) and the GSI experimental results (full dots).
be accounted for. Explicitly, when treating the voxel l ÿ 1 the ¯uence wlÿ1 is given by =dlÿ1 , wlÿ1 DM ÿ wl d lÿ1 l
4
where E llÿ1 is the dose deposited in the voxel l ÿ 1 while treating voxel l. Applying the previous equation recursively, one can easily determine all the ¯uences. Explicitely, the required ¯uence for voxel i is given by ! X i wi DM ÿ wj d j =di , 5 ji1, l
where d ij is the dose deposited in the voxel i during the irradiation of the voxel j.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
319
Table 2 Linear energy deposited vs. depth for a carbon ion with 270 MeV-initial kinetic energy Depth (cm)
DE (MeV/cm)
10.6132 11.3207 11.6038 11.8868 12.1698 12.4528 12.7358 13.0189 13.3019 13.5141 13.6556 13.7972 13.9387 14.0802 14.1085 14.1368 14.1651
176.539 188.825 195.194 202.748 211.921 223.247 237.533 256.248 282.32 310.563 337.538 377.005 445.514 642.503 721.204 836.563 787.847
The proton treatment algorithm ends here, while, as stated above, the ion treatment requires a further correction due to the energy tail behind the Bragg peak. Indeed during the treatment of voxel l ÿ 1, the tail goes through voxel l, for which the ¯uence was, in fact, already
Fig. 3. Layout of the target irradiation of the water box.
320
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Fig. 4. Layout of the target irradiation with four ®elds in the XY plane.
computed. Thus, one needs to correct it by decreasing the ¯uence by the appropriate amount. The exact solution may involve an extremely lengthy calculation. Nevertheless, we notice that the tail is characterized by a height which is of the order of 5% of the plateau. One can immediately conclude that any solution at 10% level brings down the eect to less than 1%, which is acceptable. Given the spot dimensions, the beam energy, and the target size, the ratios of the ¯uences wi =wl , to obtain a ¯at dose within the target, are totally predictable and can be stored in appropriate tables. When calculating wl, the ¯uence relative to the voxel l belonging to the most downstream section, one needs to subtract the expected dose due to later treatment of all the preceeding voxels based on the expected ¯uence ratios. The same applies to all the voxels.
4. Application of the code and Monte Carlo veri®cation 4.1. Generalities As stated in Section 1, it is necessary to verify the correctness of the calculation. We have performed the treatment plan veri®cation with a full simulation. The beam parameters
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
321
Fig. 5. Energy of carbon ions vs. depth (a) and ¯uence of carbon ions vs. depth (b) in the water box phantom.
(direction and energy) and the ¯uence, as calculated by ANCOD, are the input to the simulation which uses the standard package GEANT [1].6 The beam width (rms) has been set to 0.4 cm. We have done the check of the treatment plan for two cases: 1. A water box, where the target was a parallelepiped of dimensions 4:05 4:03 4:06 cm3, with voxel dimensions of 0:15 0:13 0:14 cm3 for 24,273 voxels. The target was placed downstream of 10.46 cm of water. We used only one ®eld of irradiation, and the beam source was set at the following coordinates: x 1:95 cm, y 1:95 cm, and z ÿ200 cm [5]. Fig. 3 shows a sketch of the set up.
6
The structure of GEANT is extremely ¯exible to describe volumes of dierent characteristics and to follow the electromagnetic and hadronic interactions of the particles, parents and daughters, along the track path.
322
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Fig. 6. Monte Carlo simulation of the depth (a) and transverse (b) dose pro®le in the water box phantom.
2. A cranial CT with 98 sections, each section split into 134 voxels along the x direction and 113 voxels along the y direction.7 Each voxel was a cube of size 0.1594 cm. The target was chosen as a cube of 21 voxels along each side, placed in the middle of the cranial CT. We used four coplanar ®elds at 908 with respect to each other. The ®eld plane was perpendicular to the axial axis and the point source was placed 04 m apart from the target. In Fig. 4, we show a sketch of the setup. So far we have discussed a situation where BRS can be obtained from TRS with a set of coordinate displacements. We remark that for ®elds not parallel to the sagittal or coronal plane or, in other words, BRS rotated with respect to TRS, it is necessary to modify consequently the CT input. This is done by rotating the voxel 7
The CT numbers were already translated into water equivalent thicknesses.
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
323
Fig. 7. Dose distribution in the water box phantom.
planes to be parallel or orthogonal to BRS axes and mapping the voxel densities from the old to the new reference system.8 5. Results First, we consider case (1) (water box) for a C+6 beam. Fig. 5(a) shows the beam energy as a function of the voxel position along the z coordinate. For a given direction, the beam energy required to have the Bragg peak centered in consecutive sections increases with depth. The 8 The mapping is done as follows: in TRS, a point x, y, z within the voxel i, j, k with given CT number corresponds, in BRS, to the point x 0 , y 0 , z 0 within the voxel i 0 , j 0 , k 0 ), and with the same CT number. By cycling over enough points (i.e. 10 per voxel), one can obtain a very accurate reconstruction of the CT numbers of the voxels in the rotated BRS.
324
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
Fig. 8. Cranial case: dose distribution due to a single irradiation ®eld.
average energy step is 02 MeV/u for a voxel size of 0.14 cm. In the simulation, to cope with present day accelerator performances, the energy input values were rounded to 1 MeV. Thus, for the example shown in Fig. 5(a), the required number of energy steps is about 50. Fig. 5(b) shows the ¯uence per dose unit with the same horizontal scale as in Fig. 5(a). As one moves from the farthest sections to those closest to the beam entrance, the required ¯uences decrease by a factor of 050. For small ¯uences, one can observe ¯uctuations (of the order of 20.5%) which are due to the approximations used. However, to evaluate the performance of the code, one has to compare the dose deposited in every single voxel with the required dose. With the GEANT simulation, we computed the dose deposited in the voxel i, j, k), dp i, j, k: The required dose in the target, DM, is expected to be constant over the whole volume. Fig. 6(a) shows the dose pro®le in a voxel row moving along the sections. Fig. 6(b) shows the dose pro®le in a voxel row moving across the target volume in a plane perpendicular to the beam. In Fig. 6(b), we observe a depletion of the deposited dose at the borders of the target which aects mainly the border voxels. The eect is clearly due to the pencil beam approximation in ANCOD. Indeed the contributions to the dose in a voxel in
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
325
the middle of the target are due to the direct irradiation of the voxel itself and dose spill out while irradiating the neighbor voxels. For a voxel placed at the border, the latter contribution is partially missing resulting in a lack of dose. We checked the simulation results against the beam width, and found similar results for beam width smaller than 4 mm (sigma). The pencil beam approximation does not hold above this ®gure, but typical values of the beam width which is expected to be used in a voxelscan treatment is smaller than 4 mm. While examining the dose pro®le in Fig. 6(b), one can observe that the algorithm is quite satisfactory over the whole volume. One can introduce the variable as a measure of the dose uniformity s X 2 ÿ 2 , dp i, j, k ÿ DM = N DM hdFi i, j, k
Fig. 9. Cranial CT case: dose distribution due to a 4-®eld irradiation.
6
326
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
where N i j k is the total number of voxels of the target. Here hdFi, which is a normalized rms, gives directly the average deviation in terms of fraction of the required dose DM. We found hdFi 3%: Finally, Fig. 7 shows a two-dimensional view of the deposited dose for a middle section of the target. We now examine case (2) (cranial CT). Concerning the beam energy and ¯uences for each ®eld, results are very similar to the water box. Fig. 8 shows the plot of the deposited dose (in arbitrary units) in a section at the center of the target for a single ®eld irradiation. The dose in the same section is shown in Fig. 9 after the four ®elds have been applied. It is remarkable that with four ®elds one obtains a ratio between the dose on the target and the dose at the plateau close to 05. Furthermore, the target edges are rather sharp: the dose decreases from 90 to 10% of the target value in 20.4 cm, which is equivalent to two voxel widths. Also for this case, the uniformity on the target is hdFi 3%: Thus, for homogeneous and inhomogeneous media we obtained comparable results with respect to the dose uniformity. A substantial contribution to the hdFi value is due to the Monte Carlo ¯uctuation since the fraction of simulated events is limited.
6. Conclusions In this paper, we have presented a package for hadrontherapy treatment planning with the voxelscan method. The software is conceived for both protons and ions, and it is based upon two types of tables: (a) range-energy and (b) energy loss as a function of beam energy and depth. Both set of tables have been derived with a full simulation of the energy losses in a water phantom. Furthermore, the consistency of the planning has been veri®ed with a full simulation of the treatment using the GEANT package. The dose uniformity obtained over the target is within 3% for both homogeneous and inhomogeneous media.
References [1] R. Brun et al., GEANT Ð Detector Description and Simulation, Tool version 3, User's Guide, CERN DD/EE/ 84-1, 1987. [2] L. Sihver, C.H. Tsao, R. Silberberg, T. Kanai, A.F. Barghouty, Total reaction and partial cross section calculations in proton±nucleus Zt R26 and nucleus±nucleus reactions (Zp and Zt R26), Phys. Rev. C 47 (1993) 1225. [3] L. Sihver et al., A model for calculating depth-dose and ¯uence distributions, GSI Scienti®c Report 95-1, p. 221, March 1995. [4] C. Brusasco, (GSI) private communication, December 1998. [5] B. Schaner, (PSI) private communication, May 1996. Stefano Berga received his degree in Physics in 1997 from the University of Torino, Italy. His work aimed at modelling and implementing nucleus fragmentation in the energy range from 400 MeV/u down to few MeV/u. He is presently teaching Physics at the High School ``Liceo Scienti®co dell'Istituto Sociale'', Torino. Faiza Bourhaleb received her License in Theoretical Physics in 1996 from University Mohamed First of Oujda, Morocco. In accomplishing her graduate studies, she is working in particle physics on Grand Uni®ed Theories (GUT and SUSY) and their
S. Berga et al. / Computers in Biology and Medicine 30 (2000) 311±327
327
consequences in cosmology. Since December 1997 she has been working on Hadrontherapy and cell survival models. Her research interests include Medical imaging, simulation in Biomedical Physiscs. Roberto Cirio received his Physics degree in March 1983 at the Physics Department of the Torino University. He has worked at CERN, on the e+eÿ collider LEP with the DELPHI detector. In 1992, he joined the ZEUS collaboration, at the e+ÿ/p collider HERA at DESY. Since 1998, he has been working on the CMS detector for the pp collider LHC at CERN. He has also been working on development of detectors and related items to be used in medical physics since 1992. M. Jamal E. Derkaoui received his Ph.D. degree from the Pierre et Marie Curie University (Paris 6, France) in May 1984. He has been working as a Professor of Physics at the Mohamed 1st University of Oujda (Morocco) since October 1984. His present research activity are mainly focused on hadrontherapy, particles physics and astrophysics. Barbara Gallice received her degree in Physics in 1998 from the University of Torino, Italy. Her main work, for the thesis, was devoted to develop a code for treatment planning in the hadron therapy framework. She is presently teaching Physics at the High School ``Istituto Tecnico I. Porro'' of Pinerolo, Italy. Mohamed Hamal received his M.D. (CEA) degree from Ibn Tofail University of Kenitra (Morocco) in 1994. He is presently in the last year of his Ph.D. studies at the Mohamed 1st University of Oujda (Morocco). He is studying the interaction between hadrons particles and the biological tissues and the eectiveness of hadronthearpy with respect to various paramaters as dose, energy, cell characteristics. F. Marchetto received his degree in Physics in 1973 from the University of Torino. Since then he has been working as a Research Associate at the Istituto Nazionale di Fisica Nucleare (INFN). The activities have been primarily focused on the ®eld of High Energy Physics from basic research to detector developments. He is author of numerous publications on several international journals. Valter Rolando received his degree in Physics from the University of Torino in 1994. At present he is working as System and Network Manager at the Faculty of Medicine and Surgery, University of Piemonte Orientale in Novara, Italy. He has worked in medical imaging and simulation. Silvia Viscomi received her degree in Physics in 1998 from the University of Torino, Italy. Her work was focused on treatment planning algorithms in Hadrontherapy. In 1999, she attended the School of Specialization in Medical Physics at University of Genoa. Presently, she is awailing a fellowship at the Cancer Epidemiology Unit of the University of Torino, where she is doing studies on survival and prognostic factors for childhood cancer.