Introduction
The Monte-Carlo method is a statistical random sampling technique for solving complex multi-dimensional integral equations that are difficult to solve analytically. The general method was introduced in 1949 by Metropolis and Ulam. Reference Metropolis and Ulam1 X-ray interactions at energies of interest in radiotherapy and medical imaging typically result in the production of charged particles such as electrons and positrons. The Boltzmann transport equation can be used to describe the motion of the coupled X-ray and charged particle transport through a defined geometry. Reference Duderstadt and Martin2,Reference do Amaral Rodriguez and Vilhena3 An exact deterministic solution to this complex multi-dimensional equation is difficult to obtain, providing motivation for the use of the more efficient statistical Monte-Carlo method to solve the problem. One of the early demonstrations of the value of the Monte-Carlo method for modelling X-ray interactions was published in 1957 by Bruce and Johns. Reference Bruce and Johns4 They showed how the Monte-Carlo method could be used to calculate the spectra of scattered X-rays in water for a radiotherapy X-ray beam. A solution to the problem of using Monte-Carlo methods for simulating the transport of electrons was proposed in 1963 by Berger, Reference Berger5 introducing the method known as the condensed history technique. Modelling the transport of ionising radiation using Monte-Carlo lends itself to implementation on digital computers, and Berger notes the implementation of his algorithm on an IBM 704 computer in the FORTRAN computing language. Indeed, it is the increased availability and rapid increase in computing power since the early 1990’s that has been a driver for the significant increase in the publication and citation of research into the use of Monte-Carlo for radiotherapy applications. Figure 1 shows the results of a Web of Science search (using the terms ‘Monte-Carlo’ AND ‘Radiotherapy’ OR ‘Radiation Therapy’) that reveals a total of 8419 publications which together have accumulated over 148,000 citations.
The main components of a Monte-Carlo ionising radiation transport simulation are the geometry, the physics models and the cross-section data representing the probabilities of the different physics occurring as a function of energy and material. Developing computer code that accurately encapsulates these three components, particularly for complex geometries is a serious undertaking. However, implementing a basic Monte-Carlo algorithm to model X-ray or gamma ray transport in simple geometries should be within the capabilities of most final-year physics, maths or engineering undergraduate students. The interested reader can find examples of simple implementations of the Monte-Carlo method in two interesting educational papers Reference Arqueros and Montesinos6,Reference Sharifzadeh, Afarideh, Khala and Gholipour7 . Fortunately, for those wishing to accurately model the detailed physics of the coupled X-ray and charged particle interactions in the complex geometries found in radiotherapy and medical imaging, a number of freely available codes are available that require relatively little programming experience. Additionally, Monte-Carlo algorithms are also now a common feature in most of the commercial radiotherapy treatment planning systems. Reference Chetty, Curran and Cygler8 This has made Monte-Carlo simulation accessible to those interested in the application of the technique to their own user-specific problem without having to go through the time-consuming process of developing code from scratch.
A constant caveat to the use of Monte-Carlo techniques for simulating radiotherapy treatments is the trade-off between statistical accuracy and the calculation time, known as the efficiency of the simulation. In the context of Monte-Carlo simulations, statistical accuracy is achieved by using a large number of random samples or particle ‘histories’ e.g. the number of electrons or X-rays, for radiation transport problems. As the number of histories is increased, the calculated quantities e.g. absorbed dose or fluence converge to those expected from real-world measurement. As well as the number of histories, other important factors that influence the statistical accuracy include the quality of the random number generator and the complexity of the physics and geometry models used in the simulation. Improvements in efficiency have been made through the use of variance reduction techniques and parallel processing using multiple Central Processing Unit (CPUs), and more recently, implementation on graphical processing units (GPUs). Reference Verhaegen and Seco9–Reference Jia, Gu, Graves, Folkerts and Jiang12 Variance reduction techniques are methods used to decrease the statistical variability (the variance) of calculated quantities such as absorbed dose or fluence without increasing the computational time, thus improving the efficiency of the use of computational resources.
In this first of two overview papers, we will begin by introducing the reader to the different Monte-Carlo codes that have been most widely used for modelling the different aspects of radiotherapy treatments. A companion paper will provide an overview of the main areas of application in radiotherapy, including modelling the production of beams of ionising radiation for radiotherapy and medical imaging, treatment verification, patient dosimetry and radiobiology.
Monte-Carlo codes
There are a number of Monte-Carlo codes that can be used to model the radiation transport problem. Short introductions to the codes that are freely available, still supported and most widely used in radiotherapy will now follow. It is acknowledged that there are other codes, such as VMC/VMC++ Reference Jia, Gu, Graves, Folkerts and Jiang12–Reference Sempau, Wilderman and Bielajew17 and DPM/gDPM, Reference Jia, Gu, Graves, Folkerts and Jiang12,Reference Sempau, Wilderman and Bielajew17 that have been developed and described in the literature over the years but are no longer freely available or supported as stand-alone codes. In some cases, they have been integrated into other more user-friendly or commercial treatment planning systems. The FLUKA code should also be acknowledged, as it is used extensively at CERN for modelling high energy particle physics and has also been used for simulating charged particle and heavy ion radiotherapy. Reference Böhlen, Cerutti and Chin18
EGSnrc/BEAMnrc/DOSXYZnrc
EGSnrc is one of the most widely used Monte-Carlo codes for radiotherapy and medical imaging applications. It is based on the Electron Gamma Shower (EGS) code developed at the Stanford Linear Accelerator (SLAC), Reference Nelson, Hirayama and Rogers19 and now maintained by the National Research Council of Canada, and distributed for free (https://github.com/nrc-cnrc/EGSnrc). The code runs on Linux, macOS and Windows operating systems and is able to model the transport of electrons, positrons and gammas with kinetic energies in the range 1 keV to 10 GeV. The code uses an implementation of the condensed history technique for charged particle propagation. Reference Kawrakow20 The code includes user codes with user-friendly graphical user interfaces that simplify the process of modelling the treatment heads of medical linear accelerators and performing patient dose calculations. BEAMnrc comprises component modules that facilitate the modelling of the geometry of components (e.g. target, primary collimator, flattening filter, jaws, MLC, etc.) found in the linear accelerator treatment head. Reference Rogers, Faddegon, Ding, Ma, We and Mackie21 DOSXYZnrc enables the calculation of dose deposited in voxelised rectilinear geometries including patient models derived from CT data. Reference Kawrakow and Walters22 The EGSnrc toolkit maintains flexibility through the inclusion of a wide range of C++ classes, known as egspp, that facilitate the modelling of more complex geometries,. Reference Kawrakow23
GEANT4
The GEANT4 toolkit was developed at CERN for modelling the passage of particles through matter Reference Agostinelli, Allison and Amako24–Reference Allison, Amako and Apostolakis26 (https://geant4.web.cern.ch/). Despite its original purpose being the simulation of high energy physics experiments and detectors, it has been extensively used for radiotherapy applications that include X-ray and particle beam therapy, micro and nano-dosimetry and radiation protection. Reference Arce, Bolst and Cutajar27 Electromagnetic physics is extended down to energies below 1 keV and up to the TeV range. Reference Chauvie, Guatelli and Ivanchenko28 The code is a C++ toolkit that makes use of contemporary object-oriented software engineering principles including the implementation of multi-threading on multi-core computer architectures. The power and flexibility that this design and implementation methodology gives GEANT4 come at the expense of the user being required to have significant prior knowledge and skills of developing applications using a modern C++ toolkit. This has limited its use in the radiotherapy community and has motivated the development of a number of more user friendly software tools that act as an interface or wrapper that makes GEANT4 more accessible to those without object-oriented programming expertise. These include GAMOS, Reference Arce, Rato, Canadas and Lagares29,Reference Arce, Lagares and Harkness30 GATE, Reference Jan, Benoit and Becheva31 PTSIM Reference Aso, Kimura, Kameoka, Murakami, Sasaki and Yamashita32,Reference Aso, Kimura, Yamashita and Sasaki33 and TOPAS. Reference Faddegon, Ramos-Méndez and Schuemann34,Reference Perl, Shin, Schumann, Faddegon and Paganetti35 A further advantage of the GEANT4 code is the GEANT4-DNA extension (http://geant4-dna.in2p3.fr/index.html) that enables modelling of the step-by-step discrete interactions of ionising particles in water at the cellular length scale. Reference Incerti, Baldacchino and Bernal36–Reference Incerti, Douglass, Penfold, Guatelli and Bezak38 Physical, chemical and biological effects of ionising radiation interactions in water can be modelled Reference Peukert, Incerti and Kempson39,Reference Sakata, Belov and Bordage40 using GEANT4-DNA.
GATE
The GEANT4 Application for Tomographic Emission (GATE) (http://www.opengatecollaboration.org/) is based on the GEANT4 toolkit and currently enables the simulation of Emission Tomography (PET and SPECT), computed tomography (CT), Bioluminescence and Fluorescence Imaging and Radiotherapy geometries. Reference Jan, Benoit and Becheva31,Reference Jan, Santin and Strul41–Reference Cuplov, Buvat, Pain and Jan44 GATE was originally developed for the nuclear medicine community with a primary aim of enabling the end user to model nuclear medicine systems without any requirement for prior knowledge of C++. Instead users use a more intuitive scripting language for creating geometries and setting simulation parameters that can then be run interactively or in batch mode. A feature of GATE is its capability for simulating dynamic or time-dependent aspects of an imaging experiment, for example, a decaying source, source and/or detector movement or breathing motion of a patient. Reference Santin, Strul and Lazaro45 In more recent times, GATE has been extended to include bioluminescence and optimal imaging Reference Cuplov, Buvat, Pain and Jan44,Reference Kang, Song, Han, Kim and Hong46 and radiotherapy, including particle therapy. Reference Jan, Benoit and Becheva31,Reference Sarrut, Bardiès and Boussion43,Reference Zarifi, Ahangari, Jia and Tajik-Mansoury47,Reference Zarifi, Ahangari, Jia, Tajik-Mansoury, Najafzadeh and Firouzjaei48
TOPAS
The TOPAS (Tool for PArticle Simulation) Monte-Carlo code was also developed to make it easier for the medical physicist to perform simulations of radiation transport using the GEANT4 code Reference Faddegon, Ramos-Méndez and Schuemann34,Reference Perl, Shin, Schumann, Faddegon and Paganetti35 (www.topasmc.org). Little or no knowledge of the GEANT4 code or the C++ programming language is required by the user. TOPAS simulations are controlled through a user-friendly TOPAS Parameter Control System that wraps the GEANT4 code while maintaining the full functionality of the underlying code including 4D time-dependent simulations. For more advanced users with C++ experience, there is the opportunity to develop their own extensions for integration into the TOPAS code. TOPAS was originally developed to facilitate the simulation of proton and carbon-ion therapy systems and has been used to develop models of passive and pencil beam scanning systems. Reference Liu, Zhang and Chen49–Reference Huang, Kang and Souris51 Despite its particle therapy roots, it is also able to model the more widely used photon and electron beams to the extent that an MV linac example is now included as part of the more recent releases. Examples of uses for modelling Brachytherapy sources are also provided. TOPAS has also been extended to include the Geant4-DNA radiobiology capabilities, Reference Schuemann, McNamara and Ramos-Méndez52–Reference McNamara, Geng and Turner55 through TOPAS-nBio (github.com/topas-nbio/TOPAS-nBio).
PENELOPE
PENELOPE is able to simulate electron and photon transport Reference Baró, Sempau, Fernández-Varea and Salvat56 utilising a mixed technique for modelling electron and positron collisions. The latest version, PENELOPE2018, is distributed by the OECD-NEA (https://www.oecd-nea.org/tools/abstract/detail/nea-1525/). The tools PENGEOM and penGUIn are available to simplify the definition of geometries and running simulations, respectively. PENELOPE has been successfully used to model medical linear accelerators through the extension PENLINAC Reference Rodríguez57 as well as the more specialised Tomotherapy Reference Sterpin, Salvat, Cravens, Ruchala, Olivera and Vynckier58,Reference Sterpin, Chen, Chen, Lu, Mackie and Vynckier59 and Leksell Gamma Knife systems. Reference Moskvin, DesRosiers, Papiez, Timmerman, Randall and DesRosiers60,Reference Moskvin, Timmerman and DesRosiers61
PRIMO
The PRIMO software (https://www.primoproject.net/primo/) enables the simulation of medical linear accelerators and patient dose calculations Reference Rodriguez, Sempau and Brualla62 using a user-friendly graphical interface. It is slightly different from other radiotherapy Monte-Carlo software in that it is available for free, but is not open source, instead being distributed as a compiled executable that runs in a 64-bit Windows environment. Parallel processing is supported on systems with multi-core processors. Reference Rodriguez and Brualla63 Beneath the intuitive graphical user interface, radiation transport is performed using the PENELOPE and DPM Monte-Carlo codes. Reference Sempau, Wilderman and Bielajew17,Reference Baró, Sempau, Fernández-Varea and Salvat56,Reference Rodriguez, Sempau, Bäumer, Timmermann and Brualla64
The ‘dose planning method’ (DPM) is a code for simulating the transport of electrons and photons in the context of radiotherapy. Reference Sempau, Wilderman and Bielajew17 The DPM code is designed to offer accurate 3D dose calculations in a fraction of the computational time of some of the other widely used codes. A mixed method for simulating electron and positron interactions is employed, with the choice of charged particle method (interaction-by-interaction or continual energy loss) chosen depending on the magnitude of the energy loss of the charged particle in the interaction. Photon interactions are modelled in an analogue manner.
Earlier versions (still available for download) of the PRIMO software supported both Elekta and Varian accelerator models, but recent versions only contain the Varian models that include a reverse-engineered TrueBeam model known as FakeBeam. Reference Rodriguez, Sempau, Fogliata, Cozzi, Sauerwein and Brualla65 Fakebeam models both flattening filter and flattening filter-free modes. Reference Belosi, Rodriguez and Fogliata66 PRIMO supports the import of DICOM RT Structure and Plan files enabling the simulation and evaluation of clinical IMRT and VMAT treatments. Reference Esposito, Silva, Oliveira, Lencart and Santos67,Reference Paganini, Reggiori and Stravato68 A further study also demonstrated the possibility of calculating patient dosimetry using Varian dynalog files. Reference Rodriguez and Brualla69 Comparisons of PRIMO with other Monte-Carlo codes have shown good agreement. Reference Aamri, Fielding and Aamry70,Reference Lloyd, Gagne, Bazalova-Carter and Zavgorodni71
MCNP
The Monte-Carlo N-Particle (MCNP) code (https://mcnp.lanl.gov/) is a general purpose radiation transport simulation code that is able to track 37 different particle types over a broad range of energies and up to 1 TeV/nucleon. Reference Kulesza, Adams and Armstrong72 The MCNP code is available to the worldwide user community, subject to USA national security restrictions and distributed through the Radiation Safety Information Computational Center, part of Oak Ridge National Laboratory. MCNP has been applied to a wide range of medical physics problems Reference Solberg, DeMarco and Chetty73 including the modelling of medical linear accelerators, Reference Mesbahi, Reilly and Thwaites74–Reference Kim, Siebers, Keall, Arnfield and Mohan78 and imaging systems such as radiotherapy electronic portal imaging devices (EPIDs), Reference Juste, Miro, Diez, Campayo and Verdu79–Reference Juste, Miró, Morera, Díez, Campayo and Verdú81 CT Reference Ding, Gu, Trofimov and Xu82 and PET/CT. Reference Waeleh, Saripan, Musarudin, Mashohor and Ahmad Saad83 MCNP has also been used to develop models for kV intraoperative radiotherapy, Reference Aghdam, Baghani, Mahdavi, Aghamiri and Akbari84,Reference Aghdam, Siavashpour, Aghamiri, Mahdavi and Nafisi85 Brachytherapy, Reference Furstoss, Reniers and Poon86–Reference Reniers, Verhaegen and Vynckier88 proton therapy Reference Jette and Chen89,Reference Verburg, Shih and Seco90 and gamma irradiator devices Reference Rodrigues, Grynberg and Ferreira91,Reference Moradi, Khandaker, Abdul Sani, Uguru, Sulieman and Bradley92 including the Leksell Gamma Knife Reference Trnka, Novotny and Kluson93,Reference Banaee, Asgari and Nedaie94
Conclusions
This paper is the first of two to give an overview of the use of Monte-Carlo simulation techniques and their application to radiotherapy. This first part has introduced the different codes that are available and currently supported with the aim of assisting the reader who wishes to develop their own use of Monte-Carlo for clinical or research applications. It is worth noting that many of the commercially available treatment planning systems, including Raystation (RaySearch Laboratories), Eclipse (Varian Medical Systems) and Monaco (Elekta) also now have Monte-Carlo algorithms as part of a suite of algorithms for electron, photon and more recently proton beams. Reference Paudel, Kim and Sarfehnia11,Reference Ali, Alsbou and Ahmad95–Reference Richmond, Angerud, Tamm and Allen97 Acceptable computational processing times are in some cases achieved through the implementation of these algorithms using GPU technology. The result as we move forward is that the radiotherapy practitioner will increasingly find themselves using Monte-Carlo techniques as part of treatment planning and treatment verification.