The use of tetrahedral mesh geometries in Monte Carlo simulation of applicator based brachytherapy dose distributions

Accounting for brachytherapy applicator attenuation is part of the recommendations from the recent report of AAPM Task Group 186. To do so, model based dose calculation algorithms require accurate modelling of the applicator geometry. This can be non-trivial in the case of irregularly shaped applicators such as the Fletcher Williamson gynaecological applicator or balloon applicators with possibly irregular shapes employed in accelerated partial breast irradiation (APBI) performed using electronic brachytherapy sources (EBS). While many of these applicators can be modelled using constructive solid geometry (CSG), the latter may be difficult and time-consuming. Alternatively, these complex geometries can be modelled using tessellated geometries such as tetrahedral meshes (mesh geometries (MG)). Recent versions of Monte Carlo (MC) codes Geant4 and MCNP6 allow for the use of MG. The goal of this work was to model a series of applicators relevant to brachytherapy using MG. Applicators designed for 192Ir sources and 50 kV EBS were studied; a shielded vaginal applicator, a shielded Fletcher Williamson applicator and an APBI balloon applicator. All applicators were modelled in Geant4 and MCNP6 using MG and CSG for dose calculations. CSG derived dose distributions were considered as reference and used to validate MG models by comparing dose distribution ratios. In general agreement within 1% for the dose calculations was observed for all applicators between MG and CSG and between codes when considering volumes inside the 25% isodose surface. When compared to CSG, MG required longer computation times by a factor of at least 2 for MC simulations using the same code. MCNP6 calculation times were more than ten times shorter than Geant4 in some cases. In conclusion we presented methods allowing for high fidelity modelling with results equivalent to CSG. To the best of our knowledge MG offers the most accurate representation of an irregular APBI balloon applicator.

Accounting for brachytherapy applicator attenuation is part of the recommendations from the recent report of AAPM Task Group 186.To do so, model based dose calculation algorithms require accurate modelling of the applicator geometry.This can be non-trivial in the case of irregularly shaped applicators such as the Fletcher Williamson gynaecological applicator or balloon applicators with possibly irregular shapes employed in accelerated partial breast irradiation (APBI) performed using electronic brachytherapy sources (EBS).While many of these applicators can be modelled using constructive solid geometry (CSG), the latter may be difficult and timeconsuming.Alternatively, these complex geometries can be modelled using tessellated geometries such as tetrahedral meshes (mesh geometries (MG)).
Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The recently published report of AAPM TG186 on brachytherapy dosimetry recommends accounting for the non-water equivalence of treatment geometries (Beaulieu et al 2012).The report recommends performing model based dose calculations such as the ones based on Monte Carlo (MC) simulations, finite element modelling or collapsed cone convolution.This departure from the water kernel based dose calculation approach of AAPM TG43 (Rivard et al 2004) entails adequate modelling of applicators employed in source delivery for brachytherapy for both low energy (<50 keV) and high energy (>50 keV) photon sources.The modelling of such applicators can be suboptimal when using a voxel based geometry due to the subvoxel dimensions of specific components.This may lead to volume averaging of the details of the geometry in coarse voxels, and may therefore lead to dose calculation errors propagating in the whole geometry.The combination of a voxelized Cartesian grid (representing the CT-derived patient geometry) and constructive solid geometry (CSG) describing the applicator allows applicator modelling.However applicator modelling using CSG can be tedious, may not allow complete fidelity or may be highly impractical, as in the case of deformable balloon applicators employed in accelerated partial breast irradiation (APBI).
The use of tessellated surfaces, defined by a collection of 2D tiles (e.g.triangular) of varying dimensions, or tessellated volumes defined by a collection of varying 3D elements (polyhedrons) can be used to describe complex geometrical shapes and offers an alternative to CSG modelling.This is especially attractive when manufacturer CAD designs are available.This methodology has been employed by commercial deterministic particle transport software capable of handling mesh geometries (MG) (Gifford et al 2006, Lloyd andAnsbacher 2013).Recent versions of general purpose MC codes have the ability to simulate radiation transport in tessellated or mesh geometries (MG), thus potentially facilitating the modelling of complex brachytherapy applicators.The aim of this work was to evaluate MG modelling by comparison to CSG modelling of a selection of brachytherapy applicators: the Fletcher Williamson GYN 192 Ir HDR brachytherapy applicator, successfully modelled using CSG techniques by several groups (Williamson 1990, Weeks 1998, Markman et al 2001, Gifford et al 2005, Price et al 2006, Mikell et al 2013), a shielded vaginal HDR applicator and an accelerated partial breast irradiation (APBI) balloon applicator used with a 50 kV electronic brachytherapy source (EBS).Dose distributions were obtained using the Geant4 (Agostinelli et al 2003) and MCNP6 (Goorley et al 2012) general purpose MC codes.

Geant4 simulation setup
The MC toolkit Geant4.9.5.p02 was used in this work, using the G4EmLivermorePhysics class for low energy electromagnetic physics and the layered mass geometry method (Enger et al 2012).The code used was based on (Enger et al 2012) where a voxel grid is used to represent the patient geometry derived from CT imaging and a parallel world contains the model of a brachytherapy source.In our work we substituted the brachytherapy sources by an applicator model and used a voxelized water box of dimensions (20 cm) 3 with (2 mm) 3 voxels to represent the patient geometry in the main world.In locations where voxels and applicator model overlapped, photons were transported in the latter.Energy depositions occurring in the applicator model were discarded in this study.The dose in voxels partially covered by the applicator model was not corrected for the fractional voxel mass which received energy depositions, which affects only the response of these voxels.
To perform MG modelling we installed the Geant4 library CADMesh version 0.6.2(Poole et al 2012a) which enables the import of tessellated surfaces or volumes.It has been reported (Poole et al 2012b) that the use of tessellated surfaces based on the G4TessellatedSolid class requires longer simulation times than tessellated volumes based on an assembly of G4Tet (Geant4 class for modelling tetrahedra).This was confirmed in our work; therefore the latter approach was used.Several file formats are supported by CADMesh; we opted for an .eleand .nodedescription of volume meshes composed of a collection of tetrahedrons8 .The .node file contains a list of vertices which are grouped in the .elefile to describe volume elements.We will refer to volumes represented by a collection of tetrahedrons as MG.Upon import via CADMesh in Geant4 the .ele/.node tetrahedral mesh yields a G4Assembly volume (a class allowing to group a number of volumes such as G4Tet) to which a position in the parallel world and a material can be assigned (Poole et al 2012b).
The kerma per event from photon sources (electron transport was not simulated) was scored in the 2 mm 3 voxels of the water geometry using track length scoring (Williamson 1987) using the Geant4 class of (Afsharpour et al 2012).We will refer to kerma as dose in this paper.Energy deposition inside the applicator was ignored.

MCNP6 simulation setup
MCNP6 was used in this work, using the MCPLIB84 photon cross-section library in Mode P which means secondary electrons were not transported (therefore, kerma was scored).The patient geometry was approximated by a water box with the same dimensions as adopted for the Geant4 simulations.Kerma was scored per event using track length estimation with mass energy absorption coefficients from NIST (Berger et al 2005), which are also used in the Geant4 simulations.A virtual grid, FMESH, with the same resolution as the Geant4 voxels was used to define the scoring geometry.This was done because MCNP6 does not allow a voxel geometry built with CSG to overlap with mesh geometries (Martz 2012).
The capability of handling mesh geometries was recently included in the MCNP6 beta 2 release, which can handle first and second order tetrahedral, pentahedral and hexahedral elements defined though text files directly generated by two commercial programs, Abaqus (Dessault Systèmes, France) and ATTILA (Transpire Inc. Gig Harbor, WA) or by converting the volume elements generated by other programs, such as ENGRID or GMSH (Remacle and Geuzaine 2009).We opted for tetrahedral meshes defined using the .ele/.node files used for Geant4 simulations converted to the MCNP6 format (Martz 2012).

Geometries of interest
Three types of applicators were considered in this work; two GYN applicators used with an 192 Ir HDR source and a balloon applicator used with an EBS operated at 50 kVp.The 192 Ir spectrum was taken from National Nuclear Data Center (Baglin 2012) and the 50 kV spectrum from Liu et al (2008) using a Geant4 model based on the work of Rivard et al (2006).Additionally, simpler geometries such as a water cube and a spherical APBI applicator were employed for estimation of calculation efficiency and validation.All geometries were simulated with both MC codes using the same MG.Although there are differences between the input formats, both codes use the same information, which are nodes (defined using cartesian coordinates) and faces (defined through node connections).

Water cube.
A water cube of (20 cm) 3 was modelled as a MG with a number of tetrahedrons varying between 12 and 191 514.In Geant4 this cube was treated as an applicator and placed in the parallel world.In Geant4 no voxels were used in the main world to compare calculation times fairly against MCNP6, where it is not currently possible to combine voxels and MG.Again a point source with an 192 Ir spectrum was positioned at the origin and the time to transport 10 7 photons was recorded for Geant4 and MCNP6 using a single core on the same computer.The simplistic model was used to measure the speed of both codes while dose distributions were compared to verify the consistency between the codes.For the latter 1 bn photons were simulated.

Idealized APBI applicator.
A spherical shell of thickness 0.4 mm and outer radius 2.5 cm was modelled using two MeshLab (Visual Computing Lab, ISTI, CNR) generated tessellated spherical surfaces.A python script (tetnest.py9) was then used to convert the two surfaces into a MG (.ele and .nodefiles).Three levels of detail were investigated with MG containing 4472, 17 824 and 69 408 tetrahedrons.Figure 1 shows an example MG.It was necessary for the G4Tet class of Geant4 to disable a G4Exception caused by degenerate polyhedrons having volumes smaller than a set threshold.These MG were used for direct comparison with simulations performed with CSG representation of the same geometry.The material used is a polymer loaded with barium for x-ray imaging purposes.The composition and percentage by weight are as follows: H -4.640%; C -30.970%; N -7.230%;O -1.096%; Si -16.540%,Cl -36.570%; S -0.549%; Ba -2.350% and the mass density is 1.2 g cm −3 .For this case, a photon point source was located at the origin of the spherical shell.In addition to the 50 kV and 192 Ir spectra, monoenergetic photons with energies 20, 30, 40, 50 and 100 keV were used to verify the agreement for different energies.
Using the 50 kV spectrum, the time to transport 10 7 photons was recorded for Geant4 using the MG of different elements as well as the CSG.

APBI balloon applicator.
An APBI breast balloon applicator (Xoft Inc., an iCad subsidiary, San Jose, CA) consists of a polymer balloon with the same composition and mass density as in the previous section.Figure 2(a) shows the enhanced visibility on CT images due to the presence of barium in the balloon wall as well as the shape of the applicator.The inserted balloon's inner contour was manually outlined on each CT slice of an APBI patient (White et al 2014), forming a cloud of points.A second set of points was generated by expanding the initial set by 0.4 mm about the centre of mass of the initial set.Each set was imported in MeshLab where a tessellated surface was generated using the convex hull function for each cloud of points (inner and outer surfaces of balloons with normal facing inwards and outwards respectively) and exported in .plyformat.A convex hull is the minimum geometrical volume that contains all of the contour points such that the vector between any two of the contour points does not lie outside of the volume.An example MG is presented in figure 2(b).A photon point source with a 50 kV photon spectrum was positioned at the centre of mass of the cloud of points used to generate the balloon.The contour points were also used to create a water phantom with the balloon wall represented by voxels with the same resolution as the reference patient CT images 0.82 × 0.82 × 2.00 mm 3 , evaluated using the same simulation parameters as the MG simulations.As can be seen in figure 2 the APBI applicator also includes a high density source channel end cap which is used to compensate for the lack of attenuation resulting from the air channel.This end cap and the catheter were evaluated in a separate study (White et al 2014) since they can be modelled using CSG and the goal of this study is MG modelling of the balloon.

Shielded HDR vaginal applicator.
A shielded cylindrical vaginal applicator consisting of a PMMA (H 8 C 5 O 2 , density 1.19 g cm −3 ) cylinder containing a tungsten half cylinder shield and a central stainless steel channel (see table 1 for the composition of the channel and the shield) for an 192 Ir source and a steel drive cable was modelled (figure 3).The applicator was also modelled by CSG for validation.The geometry for this applicator was based on a published model (Petrokokkinos et al 2011).The mesh geometry was created using Abaqus to first create a surface mesh.The surface mesh was then exported in the STereoLithography (STL) file format to Engrid to create the tetrahedral elements forming a volume mesh.Two levels of refinement were investigated (MG with 16 530 and 129 860 elements).Additionally a treatment plan obtained from a clinical case with 7 dwell positions spaced 0.5 cm apart was also considered.
The shield of the applicator was modelled in a separate simulation to estimate the time required to simulate 10 7 primary photons from the 192 Ir source in Geant4 and MCNP6.The number of elements composing the shield was varied from 101-59 365.In this simulation, as for the MG water cube, the water box present in the main world was not voxelized.

Shielded HDR Fletcher Williamson applicator.
The shielded HDR Fletcher Williamson applicator consists of three stainless steel channels and two polysulfone (PSU) ovoids with tungsten shields included to reduce the radiation dose to the bladder and the rectum (figure 4).This figure shows how modelling the applicator using CSG may prove challenging if all components including screws and holders need to be modelled.The compositions were obtained from the vendor under a non-disclosure agreement.The applicator was modelled using a CAD (Computer-Aided Design) file provided by Nucletron (Elekta AB, Stockholm, Sweden) to create the mesh surfaces with Abaqus.
An 192 Ir point source was located at 20 dwell positions inside each of the three channels with inter-dwell distance of 0.5 cm and loading 5.5 cm (from the tip) of the central channel and 1.5 cm of each ovoid channel.The source positions and the dwell times were obtained from a treatment plan to create a more realistic dose distribution.
CSG modelling of this applicator proved challenging since MCNP6 requires a torus to be rotationally symmetric around an axis parallel to one axis of the main system (x, y or z).To allow a straightforward CSG representation the lateral channels were aligned with the central channel by a 4.5° rotation.This rotation was also applied to the ovoids and shields.An equivalent MG was generated.

Simulation details
One billion primary photons were simulated for each case studied, except for the simulations performed for the computation time investigation, which used 10 7 primary photons instead.The Geant4 calculations were divided in 600 jobs and distributed using the HTCondor queuing system (Thain et al 2005) on a heterogeneous cluster of 9 computers comprising a total of 194 cpus and running Ubuntu and OSX operating systems, each having the same  version of Geant4 and CADMesh installed.Calculation time estimations were performed on a single core of an Intel Xeon X5650 processor running at 2.67 GHz.
The MCNP6 calculations were performed using a single Intel i7 (2860QM) with four cpus of 2.5 GHz running Windows 7.For calculation time estimations the same machine and setup was used as for the Geant4 calculations.
The use of MG was validated as follows: for cases with the idealized APBI spheres, dose distributions from Geant4/MCNP6 obtained with the mesh representations with different numbers of tetrahedrons and with the CSG representation of Geant4/MCNP6 were compared using dose ratios.The same was done for the gynaecological applicators (MG versus CSG).In this case only an MCNP6 CSG representation was employed.
The balloon case cannot be easily simulated since the balloon has a deformable irregular shape which is patient-dependent.Thus only simulations with mesh geometries in Geant4/ MCNP6 were considered.

Water cube
The water cube represented by a MG model of 12-191 514 elements yielded dose distributions which agreed with those obtained from CSG representation within 1%. Figure 5 shows  calculation times for 10 7 primary photons from Geant4 and MCNP6 as a function of the number of elements in the MG.We observed that for MG with less than 10 4 elements Geant4 took approximately 1.5 times longer than MCNP6 with differences up to 3.5 times for the highest number of elements.Calculation times increased with the number of mesh elements for both codes.Above a threshold between 10 4 and 10 5 elements the Geant4 calculation time increased at a greater rate than MCNP6.It is currently not clear to us why Geant4 behaves this way, although it may be due to different tracking algorithms since MCNP6 creates a neighbour list at the beginning of the simulation, which may result in a more efficient tracking algorithm at a cost of a longer time to process the input files at the beginning of the simulation.

Idealized APBI applicator
Our initial results showed a dose discrepancy in the region inside the idealized APBI applicator when using MCNP6 which was traced to the loss of photons scattered towards  the interior of the balloon (figure 6) which caused an underestimation of the dose there.This happens because MCNP6 considers the inner part of the spherical shell as a void gap and eliminates particles going towards this region.Although MCNP6 only accepts MG within a background cell, such as a water cube, and tracks particles in this cell once they leave the MG, the background cell is not considered if a mesh gap is identified.Hollow surfaces or gaps between different mesh parts may be treated as a void cell leading to missing particles.In this study this issue was corrected by filling the spherical shell with a mesh water sphere instead of the background water phantom, ignored due to the hollow volume.
Figure 7 presents the results for the 50 kV point source located at the origin of the idealized APBI balloon applicator.We observe agreement within 1% for the majority of the voxels between MG (using 4472 elements) and CSG representations of the geometry in Geant4 as well as good agreement between Geant4 and MCNP6 when using MG in both codes.Similar agreement was observed when using photon point sources of energies 20, 30, 40, 50 and 100 keV.Changing the number of elements in the MG from 4472-17 824 and 69 408 did not alter the results but yielded calculation times which were ~2, ~3 and ~8 times longer when using Geant4 compared to the CSG representation.

APBI balloon applicator
Figure 8 presents similar results as above for the irregular APBI balloon applicator defined from TPS contour points representing the balloon geometry during a clinical case irradiation.Agreement within 1% is observed within the 25% isodose surface between Geant4 and MCNP6.The MG contained 5195 tetrahedrons.This irregular balloon cannot easily be modelled using CSG geometries.
Figure 9 shows a voxelized APBI balloon applicator slice and the dose ratio between the results obtained with voxels and with MG.Differences of up to 30% were observed due to the low geometric accuracy of the voxel model whose resolution exceeds the balloon thickness by a factor of 2 and 5 for the voxel width/height and for the slice thickness, respectively.Moreover, the balloon wall may be represented by more than one voxel depending on the region.A more accurate model could be obtained with higher resolution voxels, which was not evaluated in this work since irregularly shaped mesh elements provided an optimal geometry representation validated against CGS models.

Shielded HDR vaginal applicator
Figure 10 shows the results of comparing dose distributions obtained from Geant4 and MCNP6 using both MG (two levels of refinement) and CSG representations of the shielded vaginal applicator.In general the agreement between Geant4 and MCNP6 and between MG and CSG are within 1% as shown.When using the MG representation with more elements some artefacts were observed in the dose distribution from MCNP6 (figure 10(c), aft of applicator) which were traced to lost particles in the simulation.With mesh geometries there is a possibility of overlapping tetrahedrons and/or void gaps between the tetrahedrons.A certain degree of overlap is accepted by both codes, however MCNP6 may kill some particles due to complete or substantial overlaps or gaps leading to results as shown in figure 6. Geant4 prints a track stuck warning and shifts the particle by a small displacement and keeps tracking.We verified that only a low fraction of simulated particles caused track stuck warnings to be printed.No killed track messages were observed in Geant4.In addition, Geant4 warning messages can be used to track the geometry problems as they provide information on the tetrahedron responsible for the stuck track.Stuck tracks were not observed in the CSG model.Increasing the number of elements again caused longer simulation times and may not result in a higher accuracy since higher number of elements does not assure the absence of overlaps or gaps.Figure 11 shows the simulation time for 10 7 primary photons using a  different number of elements in the shield for both codes.Again Geant4 requires longer calculation times than MCNP6.

Shielded HDR Fletcher Williamson applicator
Figure 12 presents the result of Geant4 and MCNP6 MG simulations of the shielded Fletcher Williamson applicator.We observed agreement within 1% between the dose ratios of both codes, except for differences up to 8% observed at some points close to the bottom of the applicator and far from the dwell positions.
The agreement observed between MCNP6 and Geant4 MG models was also obtained between a CSG and a MG model obtained for the same applicator with lateral channels rotated 4.5 degrees to align with the central channel as shown in figure 13.Most of points are within 1% with maximum difference of 4.1%.

Conclusion
We have presented methods applicable to Geant4 and MCNP6 MC codes, which allow the modelling of complex brachytherapy applicator geometries in MC simulations, as an alternative to CSG representations.We have validated our results using simple cases allowing CSG representation.For complex applicators which do not allow straightforward CSG modelling such as the shielded Fletcher Williamson applicator, we validated the MG method using a modified model.For the realistic APBI balloon applicator, which has an irregular geometry which precludes easy modelling with CSG, we have observed good agreement between MCNP6 and Geant4.The use of MG generally entails a computation efficiency penalty compared to simulation times for CSG with the same code, for the cases evaluated in this work.However, simulation times depend on the Monte Carlo code and on the simulation setup with MG being faster than CSG for some cases in the literature (Yeom et al 2014).Strategies such as scoring a phase space at the surface of applicators could be employed to minimize the impact of this penalty.The methods presented here can be used in the validation process of treatment planning systems or to evaluate modifications to applicator design.The MG method should also be equally useful in modelling other complex radiotherapy devices in external beam radiotherapy with photons, electrons or light ions.

Figure 1 .
Figure 1.Example MG for the idealized APBI applicator showing the external surface and an inner section using a cutaway plane.The wall material is barium loaded polymer.

Figure 2 .
Figure 2. (a) Axial CT image of an APBI balloon applicator inserted in a post surgical breast cavity.The wall is clearly visible due to barium loading.The EBS channel is occupied by a dummy insert to identify dwell positions.The balloon is filled with a saline solution.The high intensity pixels correspond to the end cap.Neither saline solution nor end cap are modeled in this work.(b) MG for the APBI applicator showing the external surface.

Figure 3 .
Figure 3. Schematic representation of the shielded cylindrical vaginal applicator.

Figure 5 .
Figure 5. Calculation times in seconds for simulating 10 7 primary photons from an 192 Ir source in a water cube represented by a MG with varying number of volume elements.

Figure 6 .
Figure 6.(a) Inside the balloon wall there appears to be a lack of backscatter from the balloon wall and the water beyond it in the MG case.(b) However, the photon spectra directly outside the balloon agree between MG and CSG.

Figure 7 .
Figure 7. (a) Dose ratio in an axial slice intersecting the origin from dose distributions obtained with MG and CSG representations of the idealized APBI balloon applicator using Geant4 and a 50 kV photon spectrum.Isodose lines are also presented.(b) Dose ratio obtained using MG in Geant4 and MCNP6.Dose in the balloon wall was not scored in Geant4 hence the low values of the ratio.Isodoses overlap for (a) and (b) due to the small differences obtained.The first two colour maps show a histogram (black) of the distribution of values of the dose ratio over the whole phantom volume.

Figure 8 .
Figure 8.(a) The MG of the irregular balloon.(b) Axial plot of dose ratios obtained from MCNP6 and Geant4 using the MG of the balloon and a 50 kV photon point source.Isodose lines are also presented.Dose in the balloon wall was not scored in Geant4 hence the low values of the ratio.The colour map shows a histogram (extreme right) of the dose ratio distribution over the whole phantom volume.

Figure 9 .
Figure 9. (a) APBI balloon applicator wall represented by voxels.(b) Axial dose ratio obtained from MCNP6 using the voxel and MG models of the balloon.

Figure 10 .
Figure 10.(a) Dose ratio in central axial, coronal and sagittal slices from Geant4 and MCNP calculations of the dose distribution from the shielded HDR vaginal applicator represented with MG containing 16 530 elements.(b) The dose ratio when representing the applicator with a MG of 16 530 elements in MCNP and a CSG representation.Differences are not visible in these figures (b) since most of the results are within 0.5% with maximum difference around 1%. (c) Dose ratio between MG and CSG applicator models in MCNP using 129 860 elements for the MG.Isodoses inside of the applicator were not shown.

Figure 11 .
Figure 11.Calculation times in seconds for 10 7 primary photons from an 192 Ir source with the W shield from the vaginal applicator represented by a MG with varying number of volume elements.

Figure 12 .
Figure 12.Dose ratios between Geant4 and MCNP6, using MG models.(a) Axial slice.(b) Sagittal slice of the dose ratio.(c) Coronal slice.Isodose lines from both codes are also plotted.

Figure 13 .
Figure 13.(a) MG representation of the shielded Fletcher Williamson applicator used for validation purposes.The MG geometry was exported by MCNP6 as an output file (Martz 2012).Ratio between MCNP6 (CSG) and Geant4 (MG): (b) Sagittal slice of the dose ratio.(c) Axial slice.(d) Coronal slice.

Table 1 .
Material properties of the steel channel and tungsten shield of the shielded vaginal applicator (see figure3).Elemental composition expressed in percentage of weight (%w).