1. Introduction

A growing area of research within the field biomedical optics involves the use of radiation-induced light emission in biological media, and can be attributed to the Čerenkov effect, a phenomena which results in broadband optical emission when a charged particle travels through a dielectric medium at a speed greater than the local phase velocity of light [1,2]. To date, the primary application of this light emission has been in Čerenkov luminescence imaging (CLI), a novel imaging modality capable of non-invasively monitoring the kinetic distribution of radionuclides in vivo [3–6]. Several studies have also proposed the extension of CLI to three-dimensional optical tomography [7–9]. A number of recent of studies have also demonstrated the feasibility of collecting Čerenkov emission from a turbid medium that is excited by an external radiation beam [10–13]. The Čerenkov radiation is predominately emitted in the ultraviolet (UV) and blue wavebands due to an inverse square wavelength dependence, and therefore several groups have investigated Čerenkov-induced fluorescence, an approach which allows researchers to shift the inherent emission to the near-infrared (NIR) wavelengths which are more favorable to light transport through biological tissue [10,14]. Conceptually, these measurements have the potential to impact clinical practice, as they could be used to optically monitor and optimize radiation therapy treatments by providing non-invasive molecular information to the clinician. As these research directives are in their relative infancy, it is necessary to develop and validate a simulation software package that can be used to investigate the radiation-induced light transport in biological media.

Modeling of the Čerenkov light emission process in scattering media cannot be easily described by analytical models of radiation transport (such as the diffusion approximation to the RTE), and instead requires the use of stochastic Monte Carlo (MC) methods. Multiple versions of publicly available MC codes exist that are often used to describe photon propagation in turbid media for biomedical applications. Most publicly available MC packages extend from the original work by Wang and Jacques on Monte Carlo for Multi-Layered media (MCML), which simulates to photon propagation through layered tissue. Later-generation packages evolved to consider time-resolution (tMCing), GPU-based acceleration (MCX), and mesh-based grid construction (MMC) [15–18]. However, these packages are not well suited to accurately consider all aspects of the generation and transport of radiation-induced photons, which requires description of the of high-energy radiation transport (i.e., either from an external beam or exogenous radionuclide), generation of the subsequent Čerenkov photons, and propagation of the photons through the medium while tracking photon interactions with scatterers, absorbers, and fluorophores. Customization of the existing software codes to describe these processes, while possible, would require substantial editing and subsequent re-validation.

The present paper presents a simulation package that is capable of accurately describing all aspects of the Čerenkov process in turbid media. At the core of this package is the Geometry and Tracking (GEANT4) software, which is an object-oriented toolkit for the simulation of particle propagation through matter [19]. It utilizes a plethora of physics models to simulate radiation transport of various particles through matter, has been extensively validated, and has been utilized to describe transport in the areas of high-energy, nuclear, space, and medical physics, among others [20–23]. Despite the power and vast simulation capabilities of the GEANT4 architecture, the package is inherently difficult to use and requires an extensive knowledge of C++ to design personalized simulations. However, the GEANT4 architecture for medically oriented simulations (GAMOS) provides a framework for users to easily interface with GEANT4 using only text-based scripts [24]. The GAMOS project has reported almost one thousand registered users and has been validated by several groups for various radiation-based applications [25,26]. In addition, while GEANT4 includes a set of optical photon physics processes, to the best of the authors’ knowledge simulation of diffuse light transport in GEANT4 has neither been fully characterized nor validated. The need for an expanded tissue optics toolbox is suggested by a recent study that used GEANT4 for modeling Čerenkov radiation for medical applications (e.g. as in [27]), which did not consider the light transport of the Čerenkov optical photons following generation; such investigations would require rigorous and accurate handling of optical photon-tissue interactions (e.g. absorption, scattering, fluorescence, etc.), and would allow studies to consider the influence of optical properties on the observed optical remission.

The present paper presents the development of a tissue optics plug-in that interfaces with the combined GEANT4/GAMOS architecture to extend new and enhanced capabilities for users to simulate optical photon transport through turbid media. The resulting package allows users to easily simulate the stochastic transport of both the high-energy particles and radiation-induced optical photons within biological media, while simultaneously allowing users to take advantage of the many pre-existing features of GEANT4/GAMOS, (i.e., the ability to specify advanced source detector distributions, generate complex heterogeneous geometries, and record a wide range of data outputs). Therefore, the aims of this study are to: (1) to describe a new tissue optics plug-in to the GEANT4/GAMOS package for radiation-induced light transport in biological media, (2) to validate the processes governing the optical photon-medium interactions within the context of its intended applications, and (3) to provide examples indicative of recently explored applications of the Čerenkov effect within the medical context. The results presented in this paper demonstrate how the resulting tissue optics modified GEANT4/GAMOS package may be of great use to, and extend beyond radiation-induced light transport to more conventional biomedical optics simulations.

2. Tissue optics plug-in for GEANT4/GAMOS

Modeling of radiation transport for biomedical research applications requires an advanced simulation architecture due to the vast number of relevant physics processes. In addition, analytical methods are not well suited for modeling such complex particle interactions, and therefore stochastic MC methods must be used. To date, several radiation transport software packages have been developed solely for the determination of dose in radiation therapy (e.g., EGSnrc, BEAMnrc, MCNPX). On the other hand, GEANT4 represents an architecture originally developed for high-energy particle physics simulations, with a secondary application to medical physics-oriented simulations.

As such, of the available MC codes, GEANT4 is the only framework in which both the radiation transport and optical Čerenkov process have been implemented, as historically the Čerenkov effect has been utilized for a number of particle physics applications [28–33]. Therefore GEANT4 provides a unique setting for stochastically modeling and investigating the emerging biomedical applications of the Čerenkov effect. The complexity of the simulation scenario is depicted in Fig. 1, which highlights the relevant physics processes necessary for the accurate simulation of the Čerenkov effect in biological media.

 

Fig. 1 In (a) the relevant radiation processes used in radiation induced light propagation modeling. Primary photons undergo interactions with the medium, resulting in charged particles that emit optical photons via the Čerenkov process. In (b) the implemented optical photon processes in the tissue optics plug-in for GEANT4/GAMOS are shown.

Download Full Size | PDF

In Fig. 1(a), the relevant physics processes ultimately resulting in the generation of optical photons via the Čerenkov effect for an incident x-ray photon within the medically relevant energy range, (i.e., keV – MeV) is presented. Incident photons may undergo three primary radiation interactions, Compton scattering (moderate energy), pair production (high energy), and the photoelectric effect (low energy). Note, in Compton scattering the photon is only scattered and may continue to propagate and undergo more radiation events, whereas in the latter two processes the photon is terminated. All processes lead to the generation of a secondary charged particle (in the case of pair production, both a positron and electron are generated), which subsequently undergo a separate group of radiation processes. During propagation (under the condition that the electron’s energy is greater than the minimum threshold energy for Čerenkov radiation in the given medium, see [21]), optical photons are continuously generated along the electrons incremental steps. Electrons may also undergo electron scattering, liberating more secondary electrons, or the Bremsstrahlung process in which additional x-ray photons are generated. In the case of positrons, the same physics processes are encountered, in addition to annihilation in which two 0.511 MeV photons are generated upon the positron annihilating with an electron.

The light transport physics utilized in the proposed tissue optics plug-in is shown in Fig. 1(b). The implemented physics processes represent modifications to pre-existing GEANT4 optical physics processes, as well as new processes, the details of which are summarized in Table 1. Optical photons may undergo optical absorption (and subsequent remission via fluorescence), scattering, and refraction or reflection. In the case of scattering, four scattering models have been implemented, each with a different sampled phase function, P(θ,λ). These include Rayleigh scattering, Mie scattering (in which the scattering phase function is based upon the Henyey-Greenstein approximation and the spectral dependence is modulated by the anisotropy, g(λ), a modification over the pre-existing GEANT4 Mie process which accepts only a constant anisotropy for all wavelengths), Modified Henyey-Greenstein (MHG) scattering (in which the scattering represents a proportional combination of Rayleigh and Mie scattering components through α(λ), implemented as described in [34]), and a user-defined scattering process (in which the wavelength-dependent scattering phase function is explicitly defined by the user). This latter scattering process may be of use in cases where the diffusion approximation is not valid [35].

Table 1. Implemented Tissue Optics Plug-in Classes

View Table | View all tables in this article

An important feature of any optical MC code is the sampling of step sizes for optical photon propagation [36]. The implementation in the GEANT4 architecture is as described with detail in [37], the photon step size is sampled from an exponential distribution calculated from the defined scattering coefficient, μ_s(λ), and optical absorption along the propagation path is independently sampled from a distribution defined using the absorption coefficient, μ_a(λ), with the absorption treated as a binary event (i.e., no photon packeting or partial photon weight is implemented). The fluorescence process acts in a similar manner to the optical absorption process, with a separate fluorescence absorption coefficient, μ_af(λ), and fluorescence emission spectrum, I_f(λ). In addition, the probability for fluorescent photon generation following absorption can be controlled through definition of the quantum yield, ϕ, a feature not available in the pre-existing GEANT4 class, and the time profile of fluorescence emission can be designated as a delta function or exponential distribution with a specified lifetime, t.

A change implemented in the tissue optics plug-in with respect to the pre-existing GEANT4 physics processes is the handling of units. As a simulation architecture initially designed for high-energy physics research, the global GEANT4 unit system is based on the class library for high-energy physics (CLHEP), which operates on a predefined set of unit classes. The implications of this are that for biomedical optics simulations, the attenuation coefficients must be defined as attenuation lengths, and optical photon wavelengths must be defined in terms of optical photon energy (i.e., eV). To improve the ease of use and functionality of the GEANT4/GAMOS framework for biomedical optics simulations, the tissue optics plug-in has enabled the use of inverse length and wavelength units for the definition of attenuation coefficients and optical photon energies respectively. The flexibility of the designed new unit system also facilities a multitude of unit prefixes, (i.e., a user can define the inverse lengths as /mm, /cm, /um, and the optical photon wavelength as μm, nm, pm, etc.). To compliment these unit changes, a new recordable data type has also been created which allows users to record an optical photon’s wavelength, rather than energy.

In addition, one of the novelties of performing optical simulations using the GEANT4/GAMOS tissue optics plug-in is the ability to easily perform broadband, white-light simulations. Optical properties (e.g., μ_s, g, μ_a, n, etc.) are defined at discrete wavelengths, λ, and initiated photons can be sampled from a user-defined source spectrum or uniform wavelength distribution, S(λ), within a defined waveband. However, when performing a white light simulation, the number of initiated photons necessary to converge to a statistically relevant result increases proportionally. Therefore, in the tissue optics plug-in a modification has been made to allow for two sampling methods, SM, in all photon generation methods (i.e., source, fluorescence, and Čerenkov processes). These methods include continuous and fixed sampling, where in the former, linear interpolation is used to continuously sample the generation distributions, and in the latter, optical photons are only generated at discrete wavelengths corresponding to the wavelengths where the optical properties have been defined. The latter method serves as a variance reduction technique, allowing for white light simulations to be performed within a waveband and user-chosen spectral resolution.

Finally, although a number of pre-existing scoring mechanisms are currently available in the GAMOS framework, none are equipped to record the spatial fluence distribution, a parameter of interest to the biomedical optics community. Therefore an additional scoring method has been created which effectively scores the fluence by summing and recording the track length of each propagating photon’s step by the local absorption coefficient within a given voxel. In voxels at the boundary of two differing materials, the scoring is performed such that the step length in each material is treated separately, and each is divided by the appropriate absorption coefficient before being scored in the given voxel.

3. Results

3.1 Validation simulations

To validate the optical physics of the tissue optics plug-in, we start by conducting a simulation analogous to the validation of MCML by Wang et al. in which the total diffuse reflectance, R_d, and total transmittance T_d (including unscattered transmittance) are calculated for a semi-infinite slab with the following optical properties: relative refractive index n_rel = 1.0 (matched boundary condition), absorption coefficient μ_a = 10 cm⁻¹, scattering coefficient μ_s = 90 cm⁻¹, anisotropy factor g = 0.75, and slab thickness d = 0.02 cm. To approximate a semi-infinite medium the slab cross section was set to 10 x 10 cm².

Ten Monte Carlo simulations were performed using 10⁷ photons. The difference between the calculated mean and the results from van de Hulst’s table, as well as results reported by Prahl et al. and MCML are shown in Table 2 [18,38,39].

Table 2. Validation of the Total Diffuse Reflectance and Transmittance for a Matched Boundary Condition Semi-infinite Slab to Results of van de Hulst’s Table, MCML, and Prahl et al.

View Table | View all tables in this article

For a mismatched boundary condition, we perform a similar simulation for a semi-infinite slab with the following optical properties: n_rel = 1.5 (mismatched boundary condition), absorption coefficient μ_a = 10 cm⁻¹, scattering coefficient μ_s = 90 cm⁻¹, anisotropy factor g = 0 (isotropic scattering), and slab thickness d = 0.02 cm. To approximate a semi-infinite medium the slab cross-section was set to 10 x 10 cm². Ten Monte Carlo simulations were performed using 10⁷ photons, the results of which are presented in Table 3 in comparison to results reported by MCML, Prahl et al., and Giovanelli [18,39,40].

Table 3. Validation of the Total Diffuse Reflectance and Transmittance for a Mismatched Boundary Condition Semi-infinite Slab to Results of Giovanelli, MCML, and Prahl et al.

View Table | View all tables in this article

We reproduce the validation used to verify the MCML package and use the GEANT4/GAMOS tissue optics plug-in to compute the angularly resoled diffuse reflectance and transmittance of a semi-infinite slab with the following optical properties: n_rel = 1.0, absorption coefficient μ_a = 10 cm⁻¹, scattering coefficient μ_s = 90 cm⁻¹, anisotropy factor g = 0.75, and slab thickness d = 0.02 cm. To approximate a semi-infinite medium the slab cross section was set to 10 x 10 cm². For the simulation, 10⁸ photons were used and the global coordinate system was set such that photons were launched in the –z direction, with the x-y plane parallel to the semi-infinite medium surface. The pre-existing data recording capabilities of GEANT4/GAMOS were used to record the directional cosines of each photon emerging from the slab as either reflectance or transmittance (i.e., for a photon emerging with a directional vector, $v=v_1\widehat{x}+v_2\widehat{y}+v_3\widehat{z}$ the values $v_1$, $v_2$ and $v_3$ were recorded, where the exit angle is defined as $cos(\theta )=v_3$, where $\sqrt{v^2{}_1+v^2{}_2+v^2{}_3}=1$). The results were used to generate a histogram using 30 bins between 0 and 90 deg. and are compared to data from van de Hulst’s table in Fig. 2. In addition to recording the diffuse transmittance, the unscattered transmittance was also calculated and found to be 0.13534, in agreement with the theoretically expected unscattered transmittance which is given by exp[-(μ_s + μ_a)d] = 0.13534 for the given optical properties.

 

Fig. 2 In (a) the angularly resolved diffuse reflectance, R_d(θ) is shown as a function of exit angle in comparison to van de Hulst’s table. In (b) the diffuse transmittance, T_d(θ) is plotted.

Download Full Size | PDF

For comparison to MCML, we performed simulations and recorded the internal fluence in the depth direction for a semi-infinite medium with a matched and mismatched refractive index boundary condition for the following optical properties: n_rel = 1.0 or 1.37, absorption coefficient μ_a = 0.01 cm⁻¹, scattering coefficient μ_s = 100 cm⁻¹, and anisotropy factor g = 0.9. The fluence was scored using 200 bins between 0 and 1 cm, and calculated by scoring the average track length of photons traversing each voxel divided by the absorption coefficient of the corresponding voxel. To avoid boundary effects at a depth of 1 cm, the depth of the actual medium was extended to 10 cm. In addition, to approximate a semi-infinite medium the medium cross section was set to 10 x 10 cm²_. The results are shown in Fig. 3 normalized to those reported by MCML.

 

Fig. 3 The depth resolved fluence for a semi-infinite medium with a refractive index matched (n_rel = 1.0) and mismatched (n_rel = 1.37) boundary are shown in comparison to results from MCML.

Download Full Size | PDF

To validate the time-resolving capabilities of the GEANT4/GAMOS tissue optics plug-in, a simulation similar to Fang et al. was performed [16]. A 60 x 60 x 60 mm³ homogenous domain with an origin located at the center of the domain was generated with a point source located at (0,0,60) mm incident with an initial direction of (0,0,-1) and the following optical properties: n_rel = 1.0, absorption coefficient μ_a = 0.005 mm⁻¹, scattering coefficient μ_s = 1 mm⁻¹, and anisotropy factor g = 0.01. The entire cubic domain was split into 2 x 2 x 2 mm³ voxels, and the temporal fluence was recorded for the voxel with a centroid located at (0,-16,20) mm by scoring the average track length of photons traversing the voxel by the absorption coefficient of the voxel. Due to the pre-existing GEANT4/GAMOS data recording capabilities, the temporal fluence was calculated by considering the initial global time that a photon entered the voxel of interest.

The resulting hits were used to generate a histogram of 50 bins between 0 and 5 ns. For the simulation, 10⁸ photons were used. The results are normalized and plotted in Fig. 4 in comparison to the analytical expectation from diffusion theory.

 

Fig. 4 The temporal fluence for a homogenous cubic medium in comparison to the time-domain response function (TPSF) predicted by diffusion theory.

Download Full Size | PDF

To validate the accuracy of the GEANT4/GAMOS tissue optics plug-in in resolving 2D fluence in a homogenous medium, a simulation was performed with an identical geometry to that used for the temporal resolved fluence with optical properties of: n_rel = 1.0, absorption coefficient μ_a = 0.005 mm⁻¹, scattering coefficient μ_s = 1 mm⁻¹, and anisotropy factor g = 0.01. However, in this case the fluence was recorded in a 30 x 30 grid of 2 x 2 x 2 mm³ voxels centered at y = 0 mm in the 60 x 60 x 60 mm³ homogenous medium domain. For the simulation, 10⁸ photons were simulated, and the fluence was calculated by scoring the average track length of photons traversing the voxel by the absorption coefficient of the voxel. The results are shown in Fig. 5(a) for contours produced at 10 dB spacing in comparison to that of MMC [16].

 

Fig. 5 In (a) a comparison between the 2D fluence calculated by GEANT4/GAMOS and MMC for a homogenous medium. In (b) the comparison between the 2D fluence calculated by GEANT4/GAMOS and the fluence results returned by MMC for an identical inhomogeneous medium. In both figures the contour lines are shown at 10 dB spacing.

Download Full Size | PDF

To validate the accuracy of the GEANT4/GAMOS tissue optics plug-in in resolving 2D fluence in an inhomogeneous medium, a simulation was performed with a similar geometry to that used for the homogenous 2D fluence validation: n_rel = 1.0, absorption coefficient μ_a = 0.002 mm⁻¹, scattering coefficient μ_s = 1 mm⁻¹, and anisotropy factor g = 0.01. However, a spherical inclusion with a radius of 10 mm was centered inside the cubic medium with optical properties: n_rel = 1.0, absorption coefficient μ_a = 0.05 mm⁻¹, scattering coefficient μ_s = 5 mm⁻¹, and anisotropy factor g = 0.9. The fluence was recorded in a 30 x 30 grid of 2 x 2 x 2 mm³ voxels centered at y = 0 mm in the 60 x 60 x 60 mm³ homogenous medium domain. For the simulation, 10⁸ photons were simulated, and the fluence was calculated by scoring the average track length of photons traversing the voxel by the absorption coefficient of the voxel. The results are shown in Fig. 5(b) for contours produced at 10 dB spacing in comparison to that of MMC [16].

3.2 Example radiation induced optical simulations

To demonstrate the flexibility of the GEANT4/GAMOS tissue optics plug-in in simulating radiation-induced light transport in biological media, two example simulations were performed. Each was designed to mimic relevant scenarios involving the practical utilization of optical photons generated by the Čerenkov effect.

In the first scenario, optical photons are generated by a 6 MeV external x-ray photon beam, which targets a spherical inclusion, indicative of an external beam radiation therapy (EBRT) treatment. For simplicity of the example, the beam was assumed mono-energetic, although accurate simulation of a linear accelerator’s poly-energetic particle spectrum can be accomplished by generating primaries from a phase space file. The time profile of the beam was set such that all x-ray photons were generated at time t = 0 ns, analogous to an impulse function. The geometry is shown in Figs. 6(a) and 6(b). As the x-ray photon beam enters the tissue volume, secondary electrons are liberated within the medium, which subsequently result in optical photon generation via the Čerenkov effect.

 

Fig. 6 In (a) and (b) the simulation geometry for the external radiation beam scenario. In (c) and (d) the simulation geometry for the PET agent scenario.

Download Full Size | PDF

In the second scenario, shown in Figs. 6(c) and 6(d), optical photons are generated via a radionuclide [in this case a positron-emitting tomography (PET) agent], ¹⁸F, which produces Čerenkov photons during radioactive decay, primarily due to the propagation of the emitted positrons.

In both cases, the geometry is the same and consists of: a 4 × 4 × 4 cm³ tissue volume with the origin located at its center. A 1 cm radius spherical inclusion is placed with its center located at (0, 0, 0.75) cm. In addition, a 0.1 cm radius, 4.0 cm length cylindrical vessel with its center located at (−1.0, 0, 1.75) cm rotated 63.5 degrees in the x-y plane is placed in the tissue volume. The main tissue volume contains an absorption coefficient, μ_a, and scattering coefficient, μ_s. The spherical inclusion contains a separate absorption coefficient, μ_a2, and the same scattering coefficient. In addition, the spherical inclusion contains a fluorophore with absorption coefficient, μ_af, and emission spectrum I_f. The quantum yield of the fluorophore, ϕ, was set to 1.0, and the time profile set to an exponential with a lifetime, t, of 0.1 ns. Finally, the vessel absorption coefficient was set to μ_a3, and the same scattering coefficient as the background and inclusion. For all three volumes, the refractive index was set to a spectrally constant 1.41, and the anisotropy, g, was set to a spectrally constant 0.9. In all simulations only the Mie scattering process was used. The simulated waveband was constrained to 550 – 850 nm, and the spectral characteristics of the defined attenuation coefficients are shown in Fig. 7. Note two spectra for μ_a2 are shown, as the absorption coefficient of the inclusion was simulated for a fully oxygenated and deoxygenated case.

 

Fig. 7 In (a) the spectral characteristics of the defined attenuation coefficients. The vessel background absorption, μ_a2 is shown for the completed oxygenated and deoxygenated cases. In (b) the fluorescence absorption and emission are shown.

Download Full Size | PDF

To highlight the ability of the GEANT4/GAMOS tissue optics plug-in to simulate relevant optical systems, the recorded data outputs were chosen such that an imaging and fiber based system could be modeled and investigated. To facilitate both of these systems in a single simulation, a phase plane was designated at the top surface of the tissue volume above the vessel and inclusion. Each optical photon generated within the medium due to the Čerenkov effect, emerging from the tissue volume, and traversing the phase plane was recorded and its exiting position (x,y), direction (u,v,w), creator process (Čerenkov or fluorescence), and wavelength were saved in CSV format. The data was then post-processed to yield the desired outputs.

For the external radiation beam scenario, a single simulation of 10⁷ incident x-ray photons was performed and the SO₂ = 100% μ_a2 was used. To generate the camera-based measurements, the recorded optical photon histories were binned into 1 × 1 mm² pixels for all emerging photons regardless of exiting direction, (i.e, all were assumed to be captured by an externally placed camera lens focused on the tissue surface). The resulting remission was further normalized to the number of incident x-rays, i.e., the values displayed in the images is in [captured optical photons mm⁻² incident x-ray⁻¹]. The resulting images for the total collected light, as well as that created due to the Čerenkov effect, and that due to Čerenkov-excited fluorescence are shown in Fig. 8.

 

Fig. 8 In (a) the total remitted light captured by an externally placed camera for the external radiation beam scenario. In (b) and (c) the components of the remitted light due to the Čerenkov effect and Čerenkov-excited fluorescence are shown. The white outlines denote the locations of the vessel and inclusion.

Download Full Size | PDF

In the case of the external radiation beam, Čerenkov light is generated in the horizontal direction along the entire volume of the beam traversing the tissue volume, see Figs. 6(a) and 6(b). Less remitted light is seen within the first several millimeters due to the existence of a buildup region where the number of secondary electrons generated by the incident x-ray photons has not yet reached equilibrium due to their non-zero propagation distance. In addition, the Čerenkov image appears brightest at the left side of the tissue volume due to the forward dominance of x-ray, electron scattering and forward dominated optical photon scattering, g = 0.9 for all volumes in the simulation, (i.e., optical photons generated to the right side of the reflectance image continue to propagate to the left resulting in a buildup of intensity). As expected, the fluorescence appears localized at the inclusion, and the absorbing vessel shows up as an attenuating feature in the captured images.

Utilizing the wavelength information for each recorded photon, the data processing can be taken a step further to mimic spectral imaging. To investigate the concept, the images were further processed in wavelength using 20 nm binning between 550 – 850 nm. The resulting images at three wavebands of interest, 550 – 570 nm, 690 – 710 nm, and 770 – 790 nm are shown in Fig. 9. Between 550 – 570 nm, the remitted light is weakest due to the enhanced absorption of the hemoglobin spectra, see μ_a2 and μ_a3 in Fig. 6. Spatially, the dark portions of the images in this waveband also correspond with the locations of the vessel and inclusion. Furthermore, as expected no remitted fluorescence is observed.

 

Fig. 9 In (a-c) the captured reflectance images between 550 – 570 nm for an external radiation beam incident on the tissue volume. In (d-f) and (g-i) reflectance images are shown between 690 – 710 nm, and 770 – 790 nm respectively.

Download Full Size | PDF

In the second waveband, between 690 – 710 nm, both the inclusion and vessel appear to have little effect on the captured images. This is due to the minimal absorption of hemoglobin in this waveband. Furthermore, in this waveband there exists no additional absorption due to the fluorophore, and no fluorescence emission due to the fluorophore, see Fig. 6(b). Finally, in the third waveband between 770 – 790 nm, the remitted fluorescence within the inclusion can be clearly seen.

From the same simulation and recorded data, fiber based measurements were also investigated. In order to investigate the ability of the fiber measurements to spectrally detect changes in the oxygenation of the inclusion, an identical external radiation beam simulation with the same data outputs was performed, in this case with the SO₂ = 0% μ_a2. The three fiber locations selected for analysis are shown in Fig. 10(a). The first location was chosen to be above the vessel and inclusion, the second location centered on the inclusion, and the third location symmetric to the first, but in this case not above the vessel.

 

Fig. 10 In (a) the three chosen fiber locations are shown. In (b-d) the corresponding spectra from the three positions are plotted. The solid lines denote the total collected light, and the two dotted type lines the spectral components due to Čerenkov light and Čerenkov-excited fluorescence. The oxygenated (red) and deoxygenated (blue) inclusion absorption coefficient cases are also denoted.

Download Full Size | PDF

In order to extract the appropriate subset of the total number of photons recorded in the phase plane, the photons captured by any fiber were filtered based on location and direction. For example, for the center fiber position, assuming a 1 mm diameter fiber with a numerical aperture (NA), of 0.5 and refractive index, n, of 1.5, the position data was used to select only photons with emerging positions of $(x^2+y^2)\le 0.5$mm, and directions of $co{s}^{-1}(w)\le si{n}^{-1}(NA/n)$ where the photons exiting angle is calculated as $cos(\theta )=w$where the photon directions are normalized such that $u^2+v^2+w^2=1$.

The results for the three fiber position spectra are shown in Figs. 10(b)-10(d). In each plot, the total collected spectra is shown by a solid line, and the spectral components from Čerenkov and fluorescent light by the two dotted lines. In addition, all three lines are shown for the oxygenated (red) and deoxygenated (blue) inclusion absorption coefficients.

At the first fiber position, Fig. 10(b) the spectral components due to the Čerenkov light and excited fluorescence can be clearly seen. Due to the presence of the vessel, the wavelengths below 650 nm are highly attenuated. At the second fiber position, Fig. 10(c), due to its centering on the inclusion, twice the fluorescence intensity is observed. Finally, at the third fiber position, Fig. 10(d), the spectra are nearly identical to that of the first position, except the wavelengths below 650 nm are less attenuated due to the absence of the vessel. More captured light is also observed in this waveband relative to the second position due to its increased distance from the inclusion, which contains the additional fluorophore absorption in this waveband.

Finally, to compare the Čerenkov light from an external radiation beam to that of the ¹⁸F radionuclide, a single simulation was performed using 10⁸ decay particles and the oxygenated inclusion absorption coefficient, SO₂ = 100% μ_a2. Note, an additional order of magnitude in primary particles was simulated due to the less energetic charged particles produced by ¹⁸F relative to the external radiation beam, and therefore fewer Čerenkov photons produced. For comparative purposes, only the camera based white light images are shown in Fig. 11.

 

Fig. 11 In (a) the total remitted light captured by an externally placed camera for the PET tracer scenario. In (b) and (c) the components of the remitted light due to the Čerenkov effect and Čerenkov-excited fluorescence are shown. The white outlines denote the locations of the vessel and inclusion.

Download Full Size | PDF

The difference in spatial characteristics can be clearly seen. In contrast to the external radiation beam, the localization of the incident particles within the inclusion leads to Čerenkov light and Čerenkov-excited fluorescence which are only localized within the inclusion (i.e., no remitted light is visible to the left side of the images where the external radiation beam Čerenkov light was peaked). Furthermore, due to the less energetic charged particles produced by the radionuclide relative to the external radiation beam, the remitted light, in units of [captured optical photons mm⁻² decay particle⁻¹] is approximately an order of magnitude weaker.

4. Discussion

In this study we have presented and validated a tissue-optics software plug-in for the GEANT4/GAMOS architecture that facilitates the simulation of radiation-induced light transport in biological media. While the GEANT4/GAMOS package is a well-validated simulation tool for high-energy particle transport through matter, to the best of the authors’ knowledge, this is the first study to rigorously validate the light transport in comparison to accepted standards within the biomedical optics community. The results of validation simulations performed in this study confirm the ability of the combined simulation package to accurately calculate the total diffuse reflectance and transmittance for a refractive index matched and mismatched boundary, as well as angularly, temporally, and spatially resolve the characteristics of optical photon propagation in tissue.

With validated light transport mechanics, the tissue optics plug-in provides users with access to a host of stochastic Monte Carlo (MC) modeling capabilities already included in the GEANT4/GAMOS framework, (i.e., the ability to specify advanced source detector distributions, generate complex heterogeneous geometries, and record a wide range of data outputs). The resultant GEANT4/GAMOS tissue optics package therefore provides an innovative tool to accurately simulate the measurement of Čerenkov-based optical photons that originate from high-energy radiation propagation in biological media. The Čerenkov process requires the coupled consideration of high-energy and optical photon transport (as detailed in Fig. 1), a requirement that surpasses the capabilities of all standard MC light transport packages that are currently publically available. The Čerenkov-based example simulation presented in Section 3.2 demonstrates the capability of this novel package to simulate broadband light propagation in a complex heterogeneous geometry with varying optical properties induced by a detailed source specification (i.e., external radiation beam or radionuclide of a given particle type, energy, positional, directional and temporal distribution), with flexible data outputs (e.g., data recording location, and particle position, direction, time, wavelength, and creation process). The implication of the data recording flexibility is that a given optical system of interest can be easily simulated and investigated, or multiple modalities can be simulated in a single simulation (i.e., camera and/or fiber measurements).

The presented examples serve as demonstrations of the capabilities of the simulation package, and are directly relevant to a number of already reported investigations into the medical application of Čerenkov light. For example, fiber measurements of Čerenkov light and/or fluorescence induced by an external radiation beam in tissue have been explored for oxygenation measurements and treatment monitoring during radiation therapy [10–13]. Furthermore, a number of studies have investigated the use of CLI to track radionuclides in vivo using camera-based systems [3–6]. Therefore the envisioned use of this software is to guide the design and development of optical devices used to collect Čerenkov-based optical measurements. Future work in this area is likely to address light transport related pressing questions for Čerenkov-based applications in tissue, including: optimization of detector details (e.g. type, location, orientation), improving the understanding of the location of origin of the collected photons during measurement, and characterizing influence of background tissue optical properties and heterogeneities on the collected spectral signal. Furthermore, although not detailed in the present study, radiation-induced light transport may also encompass scintillation, the luminescence process analogous to fluorescence where the absorbed energy is provided by a charged particle rather than photon. Similar to the Čerenkov effect, applications of scintillation within the field of biomedical optics are in its infancy, although a recent report has investigated the use of nano-scintillators as an internal excitation source for photodynamic therapy (PDT) [41]. Therefore, the plug-in modified software package described in the present paper may be of use in studying radiation-induced light transport of scintillation photons for an emerging number of novel applications.

While the work presented within this paper focuses on the use of the GEANT4/GAMOS tissue optics plug-in to simulate Čerenkov-based measurements, it is important to note the core aspects of the package, which makes this a unique investigative tool that can be extended to multiple areas within the field of biomedical optics. A brief list of important advantages include: (1) The specification of optical detectors at multiple (co-registered) locations that are capable of scoring photon tissue interactions from an extensive catalogue within the GAMOS framework (e.g., detectors can return interaction type, frequency, and location for collected photons); this represents an advantage over standard MC codes which may require reprogramming to perform non-standard interaction/event tracking or specification of multiple co-registered detectors. (2) Broadband simulation capabilities allow for user-specification of wavelength-dependent optical properties and model estimation of spectral response; this capability avoids the need to either initialize independent simulations at multiple wavelengths or apply perturbation scaling factors which may limit the range of scattering properties within a spectrum that can be accurately investigated [42,43]. Furthermore, the ability to sample photon generation in a fixed or continuous manner has been provided to allow users to control the spectral resolution of such simulations. (3) A geometry can be constructed using a combination of standard shapes (i.e., polygons, ellipsoids, cylinders, etc.) to yield heterogeneous volumes containing layers and/or inclusions, each of which may contain user-specified optical properties. (4) A user-defined wavelength-dependent scattering phase function can be defined, providing users with a unique framework for investigating phase functions beyond the more conventional Henyey-Greenstein approximation. Furthermore, in complex heterogeneous geometries, a separate user-defined phase function can be defined in each material, a new feature added by the tissue optics plug-in.

While these listed advantages are not necessarily unique in functionality compared with previously reported custom MC packages, the robust catalogue of user-defined inputs and outputs within the GAMOS interface allows users to adjust these simulation parameters without requiring C or C++ level programming. Furthermore, the presented plug-in and core simulation software is fully open-source and publically available unlike many modified in-house MC codes that have previously reported to achieve some of the functionality described in the present study. Despite its robust features, the GEANT4/GAMOS tissue optics package is not well-suited for all simulation purposes, as it does not offer GPU-based acceleration or mesh-based grid generation, options that limit the efficiency and accuracy in complex media. However, increased efficiency represents the goal of a future study, as the plug-in modified simulation package described herein represents the only current validated architecture capable of simulating both the high-energy radiation and subsequent diffuse light transport.

The plug-in, user guide, example files, as well as the necessary files to reproduce the validation simulations are available at http://www.dartmouth.edu/optmed/research-projects/monte-carlo-software.