

ORIGINAL ARTICLE 

Year : 2017  Volume
: 13
 Issue : 1  Page : 6979 

Tabulated squareshaped source model for linear accelerator electron beam simulation
Navid Khaledi^{1}, Mahmood Reza Aghamiri^{2}, Hossein Aslian^{3}, Ahmad Ameri^{1}
^{1} Department of Clinical Oncology, Imam Hossein Hospital, Shahid Beheshti University of Medical Sciences, Tehran, Iran ^{2} Department of Radiation Medicine Engineering, Shahid Beheshti University, Tehran, Iran ^{3} Department of Physics, University of Trieste, Trieste, Italy
Date of Web Publication  16May2017 
Correspondence Address: Ahmad Ameri Department of Radiation Oncology, Imam Hossein Hospital, Shahid Beheshti University of Medical Sciences, Tehran Iran
Source of Support: None, Conflict of Interest: None  Check 
DOI: 10.4103/09731482.206235
Context: Using this source model, the Monte Carlo (MC) computation becomes much faster for electron beams. Aims: The aim of this study was to present a source model that makes linear accelerator (LINAC) electron beam geometry simulation less complex. Settings and Design: In this study, a tabulated squareshaped source with transversal and axial distribution biasing and semiGaussian spectrum was investigated. Subjects and Methods: A low energy photon spectrum was added to the semiGaussian beam to correct the bremsstrahlung Xray contamination. After running the MC code multiple times and optimizing all spectrums for four electron energies in three different medical LINACs (Elekta, Siemens, and Varian), the characteristics of a beam passing through a 10 cm × 10 cm applicator were obtained. The percentage depth dose and dose profiles at two different depths were measured and simulated. Results: The maximum difference between simulated and measured percentage of depth doses and dose profiles was 1.8% and 4%, respectively. The low energy electron and photon spectrum and the Gaussian spectrum peak energy and associated full width at half of maximum and transversal distribution weightings were obtained for each electron beam. The proposed method yielded a maximum computation time 702 times faster than a complete head simulation. Conclusions: Our study demonstrates that there was an excellent agreement between the results of our proposed model and measured data; furthermore, an optimum calculation speed was achieved because there was no need to define geometry and materials in the LINAC head. Keywords: Electron beam, linear accelerator, Monte Carlo NParticle, source model
How to cite this article: Khaledi N, Aghamiri MR, Aslian H, Ameri A. Tabulated squareshaped source model for linear accelerator electron beam simulation. J Can Res Ther 2017;13:6979 
How to cite this URL: Khaledi N, Aghamiri MR, Aslian H, Ameri A. Tabulated squareshaped source model for linear accelerator electron beam simulation. J Can Res Ther [serial online] 2017 [cited 2021 Sep 20];13:6979. Available from: https://www.cancerjournal.net/text.asp?2017/13/1/69/206235 
> Introduction   
Monte Carlo (MC) calculation methods are widely used in radiation therapy.^{[1],[2],[3],[4],[5],[6],[7],[8]} New treatment planning systems (TPSs) for external radiation therapy are utilizing MC algorithms for more accurate dose calculation and distribution by transporting all (or most) the interactions between particles and matter.^{[9]} This pointbypoint transport is a timeconsuming process. The less material we define, the less particle interaction we have, and therefore, calculation time could be dramatically decreased.
Due to the fact that there are several materials in a linear accelerator (LINAC) head, namely, primary and secondary scattering foils, jaws, and applicator, several interactions usually occur, and as a result, the typical running time of the MC calculation is noticeably long. Accordingly, employing the final exited beam from LINAC specifications can decrease the computation time.^{[4]}
Many studies on beam modeling of photon and electrons produced by LINACs have been performed, and this subject is still being investigated. Jiang et al. worked on modeling of the LINAC head (including the jaws and primary collimator) for electron beams, in which upper head simulation data were separately stored and used.^{[2],[6],[7],[10],[11],[12],[13],[14],[15],[16],[17],[18],[19]} Then, the lower head for the electron beam transport was simulated. Their model was based on electron and photon point sources and considered scattering and photon contamination in the upper and lower head. Other models in previous methods have also been purposed as an alternative method to reduce calculation time and stored data. For example, the virtual source model does this by parameterization of the beam path and neglecting the full simulation of the LINAC head and applicator.^{[12]} In this method, the source is often divided into two sources with different axial distances and propagation angles to reproduce the scattering of electrons in the applicator body and cutout edges.
Meanwhile, several studies on simulation of the whole LINAC head and applicator have been performed. In some studies, Monte Carlo NParticle (MCNP) or EGSnrc MC codes for electron beams have been used for modeling of all the materials in the LINAC head components to obtain their specific aims.^{[5],[20]}
Ali et al. modeled the Elekta Precise LINAC electron source by assuming two point sources. They validated this model by measuring electron beams ranging from 4 to 15 MeV.^{[4]} Chang et al. investigated the effect of mean energy and the radial intensity on percentage depth dose (PDD) and the dose profile of a photon beam using an EGSnrcbased DOSXYZnrc code by generating the phase space.^{[21]}
In our study, we simulated the semiGaussian spectrum of the final beam that passed through the applicator. It was composed of a lowenergy electron beam, Xray contamination, and a Gaussian spectrum. The applied phase space was positioned at 10 cm from the phantom surface with 1 cm axial and different lateral weightings that were a mix of photon and electron beams. It was also benchmarked for a set of cutouts.
> Subjects and Methods   
Measurements
To ensure the accuracy of the MC results, the quality and distribution of the beam were benchmarked using PDD and profile measurement in a water phantom. The measurements were performed using calibrated dosimetry equipment including planeparallel and thimble chambers in a 50 cm × 50 cm × 50 cm MP3 (PTW Freiburg) water phantom. The LINAC machines were Elekta Precise, Siemens Primus, and Varian Clinac 2100CD. The standard applicator field size was set to 10 cm × 10 cm at 100 cm sourcetosurface distance for the Varian and Siemens machines and 95 cm for the Elekta. The PTW parallelplane model was the TM34045 with a sensitive volume of 0.02 cm ^{3} for obtaining the PDD and profile curves according to the IAEA TRS398 dosimetry protocol.^{[22]} A thimble chamber model TM31010 with a sensitive volume of 0.125 cm ^{3} was used as the reference chamber, which was placed outside of the water phantom at the corner of the radiation field.
The electron nominal energies studied for Elekta were 6, 10, 12, and 14 MeV. These numbers were 6, 9, 12, and 16 MeV for the Varian and 5, 10, 12, and 14 MeV for the Siemens machines. Beam profile dosimetry was performed at two different depths for each energy in the crossplane direction.
Monte Carlo simulation
Source and geometry specification
MCNP4C Monte Carlo code was used for simulation. For electron transport, the standard library EL03 was utilized.^{[23]} Considering bremsstrahlung photon contamination generated in the LINAC head components such as the scattering foils, jaws, and applicator, the MCPLIB04 photon crosssection library was used. The electron (CUT: E) and photon (CUT: P) energy cutoffs were 0.521 and 0.02 MeV, respectively.
Due to multiple scattering of the electrons, the electrons that have exited from the applicator have multiple direction angles, and that is why defining a point source is irrational. Hence, we considered a tabulated squareshaped source in MCNP code with transversal (X and Y) and axial (Z) weight samplings [Figure 1]. The source had a directional (DIR) distribution to consider the different direction of the electrons, in which the solid angle was 6° for Siemens and Varian and 7° for Elekta.^{[24]} The place (thickness) of the tabulated source was 9–10 cm (1 cm thickness) from the simulated phantom surface.  Figure 1: Distribution of the tabulated squareshaped source used in the Monte Carlo simulation (a) lateral view, (b) top view
Click here to view 
The source particle (PAR) is composed of three discrete beam types: (a) lowenergy electrons contamination spectrum, (b) Gaussian electrons spectrum, and (c) photon contamination spectrum. The probability (weight) of the photon spectrum was 50% for all electron nominal energies, and the remaining 50% was related to the lowenergy and Gaussian spectrums. The lowenergy spectrum forms the lower energy tail of the semiGaussian spectrum.^{[5],[24],[25]}
A 50 cm × 50 cm × 50 cm cube of H_{2}O was designed as a water phantom in the simulation. The central axis (CAX) was divided into 0.2 cm × 0.2 cm × 0.2 cm cubes (cells) to obtain the PDD. The depth of maximum dose (d_{max}) and a deeper depth for each electron beam energy were divided into 0.2 cm × 0.2 cm × 0.2 cm cubes to obtain the profile curves. Dose calculations were performed for each cube using tally type 8 (*F8) and pulse height distribution in a cell.^{[26],[27]} The simulated PDD and profile curves points were obtained by dividing the simulation (tally) results by the mass of a tally cell. The maximum number of histories was 10^{9}.
Percentage depth dose and profile simulation
Electron beam quality is a function of R_{50} (depth in water, in which the dose falls off to 50% of maximum dose) and the practical range, R_{p}, so that:^{[22],[28],[29],[30]}
Where E_{p,0} and E_{0(mean)} are the most probable energy at a phantom surface and the mean energy at a phantom surface, respectively, in MeV. R_{p} and R_{50} are in centimeters.
In the first step, the obtained E_{p,0} from measurements was considered as the peak energy of the Gaussian spectrum in the beam simulation. Consequently, by changing the peak energy of the Gaussian spectrum and its full width at half of maximum (FWHM), the desired beam quality was obtained. A lowenergy electron spectrum was added to the main Gaussian spectrum to constitute a semiGaussian spectrum. In addition, a lowenergy photon spectrum separately was added to take into account the Xray background. This procedure was done by changing the peak energy, FWHM, and low energy electron and photons to obtain the best matching between the simulation and measurement.
To increase the buildup dose of the PDD curve, the weight of the low energy electron spectrum was increased. The superficial parts of the PDD buildup region were affected by lesser energies of the low energy electron spectrum tail, and deeper parts of the PDD buildup were controlled by changing the higher energies of the tail. Furthermore, the peak energy of the Gaussian electron spectrum was increased by increasing the E_{p0}. The FWHM variation of the electron spectrum had a direct relation with the E_{0(mean)} and R_{p}. The weight of the low energy photon spectrum was related to the level of background Xrays. According to these headlines, the simulated PDDs were tuned and matched with the measured PDD's R_{p}, E_{p0}, and E_{0(mean)}. Benchmarking of simulation and measurement of profile flatness and penumbra were performed by changing the transversal and axial weight biasing of the tabulated squareshaped source as well as the distribution angle. For this purpose, a maximum of eight strength points (SPs) were considered in the + x and + y directions. For example, the used SPs for the 6 MeV electron beam of Elekta machine in the MC simulation are shown in [Figure 2]. More details for all other electron beam nominal energies are presented in [Table 1]. The method was the same for all radiation fields, and only the position and strength of the probability points were different. As shown in [Figure 2], the strength probability on the CAX was less than the adjacent point. This arises from the fact that any point on the curve line affects nearby points in the phantom due to the angle of incident particles and some scattering.  Figure 2: The strength probability curve used for a 6 MeV beam of the Elekta linear accelerator
Click here to view 
The field size at each depth was defined as the distance between two 50% dose points of the profile.^{[29]} Finally, the model was verified for 3 standard cutouts with widths of 3, 5, and 6 cm for Elekta, Siemens, and Varian LINACs.
The computational time required for running each code on a 3.0 GHz Core i5 Laptop, with 4 GB of RAM, ranged from ~10 to ~15 min for lowest to highest electron energy, respectively, by setting 3 × 10^{6} histories to reach an acceptable relative error.
Electron spectrum parameterization
To achieve the best result for the beam, the lower energy tail of the Gaussian electron beam equation was obtained by dynamic fitting of the electron energy probability on a curve with a specific type using SigmaPlot version 11 (Systat Software Inc., San Jose, California, USA).
The normalized Gaussian peak probability equation is presented as Eq. (3). According to this equation, the Gaussian spectrum was defined by placing the FWHM (a) and the peak energy (b).
In this equation E, b, and a are in MeV.
> Results   
The source distributions in the transversal and axial directions affect the flatness, field size, and penumbra of the profile. Whereas, changing the position and probability of the transversal SPs can alter the field size and penumbra of the beam profile, and owing to this fact, the flat plateau region width will also change. In addition, to decrease the volume of the simulated penumbra, similar to a profile dose falloff, the probability of SPs should be rapidly decreased. The penumbra and field size can be changed by changing the axial biasing weight and position so that there is a rough decrease of the axial sampling of the source, which increases the penumbra and decreases the field size smoothly. In this case, the penumbra region was changed, and the relative dose for the tail of the profile remained almost unchanged. However, we tuned it by changing the position and probability of the axial weighting probability and transversal SPs. The axial weighting probability was chosen to be from 0.1 to 10. Furthermore, by increasing the solid propagation angle, both field size and penumbra were increased.
To consider the plateau region of the beam profile, the SPs had a large amount [points 1–4 in [Table 1]. In comparison, the 6^{th} to 8^{th} SPs had small values, which were proportional to the penumbrae region of the dose profile. As shown in [Table 1], there was no sharp falloff for SPs of the lowenergy beam. The penumbra values for Elekta were also higher than the two other machines. Therefore, SP probabilities were decreased more smoothly for the Elekta (point 3, 4 and 5) in comparison to the Siemens and Varian machines. In addition, a severe decrease in SP probabilities was observed for Varian's 9, 12, and 16 MeV electrons (5^{th} and 6^{th} points), and the distance between these points was 0.5 mm. For Elekta, the SPs in the edge of the profile were increased with the energy increscent.
Due to the larger field size for the Elekta at the phantom surface in comparison with the two other LINACs, the distribution angle of the electrons was higher for this machine (above 6.2°).
As shown in [Figure 3],[Figure 4],[Figure 5] the relative dose differences between MC and measurement were under 4% for the profiles. However, the spatial difference was smaller than this value. As illustrated in [Table 2], the maximum difference between the MC and measurement field size was <1%, except for the second depth profile of the 9 and 16 MeV electron beam produced by the Varian machine, which had a difference of about 1.9%.  Figure 3: Monte Carlo and measurement profile results at two different depths for Elekta machine and their pointbypoint relative dose differences
Click here to view 
 Figure 4: Monte Carlo and measurement profile results at two different depths for Siemens machine and their pointbypoint relative dose differences
Click here to view 
 Figure 5: Monte Carlo and measurement profile results at two different depths for Varian machine and their pointbypoint relative dose differences
Click here to view 
According to [Table 2], the difference between MC and measurement penumbra value was smaller than 0.1 cm. The maximum difference of flatness at the first depth between MC and measurement was 2.2% and for the second depth was 3.4% for all machines.
The differences between computed and measured profiles are most palpable in the penumbra region, which could be due to the steepness of this region of the profile, such that a small change in the spatial position of the profile points leads to a large change in the dose level of the penumbra region.
MCsimulated PDDs were achieved by changing the previously mentioned parameters, and the results are shown in [Table 3]. The weight of the Gaussian electron beam for low energy electron beams was slightly higher than for high energy beams. Except for 12 MeV beams, the FWHM was directly increased by a nominal energy increscent. By fitting the lowenergy electron spectrum points in the MC simulation on the curve, their associated equation was obtained and is presented in [Table 3]. [Figure 6] shows the mixed lowenergy electrons spectrum to the Gaussian electron peak. The obtained semiGaussian spectrum and the photon spectrum are shown in [Figure 7]. [Figure 7] shows that due to the increase in the total Xray background for higher electron energies, the normalized probability of photon contamination was increased in the MC simulation.  Table 3: Source energy spectrum specification used in Monte Carlo simulation
Click here to view 
 Figure 6: The lower energy electron tail for main Gaussian spectrum curves. The dots indicate the points used in the Monte Carlo NParticle code and the solid lines are the results of the curve fitting
Click here to view 
 Figure 7: The normalized photon and electron spectrum used in the Monte Carlo simulation
Click here to view 
[Table 3] shows that the larger FWHMs belong to the Elekta machine energies. Due to the lower FWHM value of Siemens electron beams in comparison with two other machines, it had high steepness in the descending part of the PDD curve.
[Table 4] represents the MC and measurement results for the important parameters of PDD curves. As shown in [Table 4], the MC values are in good agreement with the measured PDDs.  Table 4: Monte Carlo and measurement results for percentage depth dose curve
Click here to view 
The 6 MeV electron beam of Elekta, 12 MeV electron beam of Siemens, and 16 MeV electron beam of Varian LINACs showed the best results by the present method, with the difference between simulated and measured PDDs under 2.1%. The cutouts width for Elekta, Siemens, and Varian was 3, 5, and 6 cm, respectively. Simulated and measured PDDs and beam profiles are illustrated in [Figure 8] and [Figure 9].  Figure 8: Beam profiles for cutouts with a width of 3, 5, and 6 cm for Elekta, Siemens, and Varian, respectively, at two different depths
Click here to view 
 Figure 9: Percentage depth dose curves for cutouts with a width of 3, 5, and 6 cm for Elekta, Siemens, and Varian, respectively
Click here to view 
> Discussion   
Compared to a full simulation of the materials in the LINAC head, an upper and lower head of the Siemens Primus LINAC was simulated in MCNP code.^{[5]} The calculation time and the fluence on the phantom surface were compared with the tabulated squareshaped source. The calculation time for the tabulated squareshaped and 700,000 histories was 180 s, and the full LINAC head simulation took 21,000 s. The fluence of a 12 MeV electron beam on the phantom surface for a full LINAC head simulation and in the same history was 0.162 per source history. For this presented source model, it was 0.997 per source particle. Therefore, the speed of the presented source model is about 117 times faster than the full simulation, and the particle number efficiency was 6 times better than the full simulation of the LINAC head. The reason why the fluence in the two simulation methods is different is due to the fact that a majority of the number of particles are absorbed by the head materials before they reach the phantom surface. It was deduced that to achieve the same fluence, the full simulation method was 702 times slower than the presented model.
> Conclusion   
In the present study, we investigated a new source model to simplify and boost the MC LINAC electron beam simulation that can be used in radiotherapy research and MCbased treatment planning algorithms. This model is hundreds of times faster than a LINAC whole head simulation. Moreover, by employing the present model, there is no need to define LINAC head components geometry. Therefore, the timeconsuming process of geometry definition can be skipped. This method can be employed by researchers in their investigations of electron beam simulation, as well as for TPS algorithm developers. Using a powerful PC with parallel computation ability, reduce the computation time and make the presented model faster than the current results.
Our work was performed with a 10 cm × 10 cm applicator size, considered as a reference applicator, but using other applicator sizes to reach more comprehensive results can be addressed in future studies.
Acknowledgment
We would like to thank Dr. Kourosh Sheibani, for helping us in language editing of the text.
Financial support and sponsorship
This study was supported by the PhD thesis funds.
Conflicts of interest
There are no conflicts of interest.
> References   
1.  Henzen D. Monte Carlo based treatment planning for modulated electron radiotherapy using a photon multileaf collimator. Diss., Eidgenössische Technische Hochschule Zurich; Nr. 21855; 2014. 
2.  SheikhBagheri D, Rogers DW. Monte Carlo calculation of nine megavoltage photon beam spectra using the BEAM code. Med Phys 2002;29:391402. 
3.  Turian JV, Smith BD, Bernard DA, Griem KL, Chu JC. Monte Carlo calculations of output factors for clinically shaped electron fields. J Appl Clin Med Phys 2004;5:4263. 
4.  Ali OA, Willemse CA, Shaw W, O'Reilly FH, du Plessis FC. Monte carlo electron source model validation for an Elekta Precise linac. Med Phys 2011;38:236673. 
5.  Khaledi N, Arbabi A, Sardari D, Mahdavi SR, Aslian H, Dabaghi M, et al. Monte Carlo investigation of the effect of small cutouts on beam profile parameters of 12 and 14 MeV electron beams. Radiat Meas 2013;5152:4854. 
6.  Khaledi N, Arbabi A, Sardari D, Mohammadi M, Ameri A. Simultaneous production of mixed electron – Photon beam in a medical LINAC: A feasibility study. Phys Med 2015;31:3917. 
7.  Jiang SB, Kapur A, Ma CM. Electron beam modeling and commissioning for Monte Carlo treatment planning. Med Phys 2000;27:18091. 
8.  Sudahar H, Gopalakrishna Kurup PG, Murali V, Velmurugan J. Influence of smoothing algorithms in Monte Carlo dose calculations of cyberknife treatment plans: A lung phantom study. J Cancer Res Ther 2012;8:36772. 
9.  Reynaert N, van der Marck S, Schaart D, van der Zee W, van VlietVroegindeweij C, Tomsej M, et al. Monte Carlo treatment planning for photon and electron beams. Radiat Phys Chem 2007;76:64386. 
10.  Lin SY, Chu TC, Lin JP. Monte Carlo simulation of a clinical linear accelerator. Appl Radiat Isot 2001;55:75965. 
11.  Tzedakis A, Damilakis JE, Mazonakis M, Stratakis J, Varveris H, Gourtsoyiannis N. Influence of initial electron beam parameters on Monte Carlo calculated absorbed dose distributions for radiotherapy photon beams. Med Phys 2004;31:90713. 
12.  Fix MK, Keall PJ, Dawson K, Siebers JV. Monte Carlo source model for photon beam radiotherapy: Photon source characteristics. Med Phys 2004;31:310621. 
13.  SheikhBagheri D, Rogers DW. Sensitivity of megavoltage photon beam Monte Carlo simulations to electron beam and other parameters. Med Phys 2002;29:37990. 
14.  Edimo P, Clermont C, Kwato MG, Vynckier S. Evaluation of a commercial VMC++ Monte Carlo based treatment planning system for electron beams using EGSnrc/BEAMnrc simulations and measurements. Phys Med 2009;25:11121. 
15.  Habib B, Poumarede B, Tola F, Barthe J. Evaluation of PENFAST – A fast Monte Carlo code for dose calculations in photon and electron radiotherapy treatment planning. Phys Med 2010;26:1725. 
16.  GlideHurst C, Bellon M, Foster R, Altunbas C, Speiser M, Altman M, et al. Commissioning of the Varian TrueBeam linear accelerator: A multiinstitutional study. Med Phys 2013;40:031719. 
17.  Dovbnya A, Nikiforov V, Uvarov V. Modeling and optimization of electron linac exit systems for nuclear technologies. Nucl Instrum Methods Phys Res A 2006;558:199201. 
18.  Makarashvili V, Chemerisov S, Micklich B. Simulations of a LINACbased photoneutron source. Nucl Instrum Methods Phys Res A 2012;696:13640. 
19.  Zheng L, Du Y, Huang W, Tang C. Simulation of dark current and dark currentinduced background photons in the Thomson scattering Xray source. Nucl Instrum Methods Phys Res A 2015;800:127. 
20.  Shimozato T, Okudaira K, Fuse H, Tabushi K. Monte Carlo simulation and measurement of radiation leakage from applicators used in external electron radiotherapy. Phys Med 2013;29:38896. 
21.  Chang K, Wang Z, Shiau A. Determining optimization of the initial parameters in Monte Carlo simulation for linear accelerator radiotherapy. Radiat Phys Chem 2014;54:1615. 
22.  IAEA. Absorbed Dose Determination in External Beam Radiotherapy: An International Code of Practice for Dosimetry Basedon Standards of Absorbed Dose to Water (TRS398). Vienna: IAEA; 2001. 
23.  Briemeister J. MCNP: A General Monte Carlo NParticle Transport Code (Version4C). Los Alamos National Laboratory, LA13709M; 2000. 
24.  Iaccarino G, Strigari L, D'Andrea M, Bellesi L, Felici G, Ciccotelli A, et al. Monte Carlo simulation of electron beams generated by a 12 MeV dedicated mobile IORT accelerator. Phys Med Biol 2011;56:457996. 
25.  Ding G, Rogers D. Energy spectra, angular spread, and dose distributions of electron beams from various accelerators used in radiotherapy. National Research Council of Canada Report. NRC Institute for National Measurement Standards; Report no. PIRS0439 1995. 
26.  Koivunoro H, Siiskonen T, Kotiluoto P, Auterinen I, Hippelainen E, Savolainen S. Accuracy of the electron transport in mcnp5 and its suitability for ionization chamber response simulations: A comparison with the egsnrc and penelope codes. Med Phys 2012;39:133544. 
27.  Goorley T, Olsher D. Using MCNP5 for medical physics applications. American Nuclear Society Topical MeetingMonte Carlo 2005; 2005. p. 1721. 
28.  Faddegon BA, Sawkey D, O'Shea T, McEwen M, Ross C. Treatment head disassembly to improve the accuracy of large electron field simulation. Med Phys 2009;36:457791. 
29.  Khan FM. The Physics of Radiation Therapy. 5 ^{th} ed. Baltimore, MD: Williams and Wilkins; 2009. 
30.  Khaledy N, Arbabi A, Sardari D. The effects of cutouts on output, mean energy and percentage depth dose of 12 and 14 MeV electrons. J Med Phys 2011;36:2139. [ PUBMED] [Full text] 
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8], [Figure 9]
[Table 1], [Table 2], [Table 3], [Table 4]
