Experimental validation of a spectroscopic Monte Carlo light transport simulation technique and Raman scattering depth sensing analysis in biological tissue.

SIGNIFICANCE: Raman spectroscopy (RS) applied to surgical guidance is attracting attention among scientists in biomedical optics. Offering a computational platform for studying depth-resolved RS and probing molecular specificity of different tissue layers is of crucial importance to increase the precision of these techniques and facilitate their clinical adoption. AIM: The aim of this work was to present a rigorous analysis of inelastic scattering depth sampling and elucidate the relationship between sensing depth of the Raman effect and optical properties of the tissue under interrogation. APPROACH: A new Monte Carlo (MC) package was developed to simulate absorption, fluorescence, elastic, and inelastic scattering of light in tissue. The validity of the MC algorithm was demonstrated by comparison with experimental Raman spectra in phantoms of known optical properties using nylon and polydimethylsiloxane as Raman-active compounds. A series of MC simulations were performed to study the effects of optical properties on Raman sensing depth for an imaging geometry consistent with single-point detection using a handheld fiber optics probe system.
RESULTS: The MC code was used to estimate the Raman sensing depth of a handheld fiber optics system. For absorption and reduced scattering coefficients of 0.001 and 1  mm  -  1, the sensing depth varied from 105 to 225  μm for a range of Raman probabilities from 10  -  6 to 10  -  3. Further, for a realistic Raman probability of 10  -  6, the sensing depth ranged between 10 and 600  μm for the range of absorption coefficients 0.001 to 1.4  mm  -  1 and reduced scattering coefficients of 0.5 to 30  mm  -  1.
CONCLUSIONS: A spectroscopic MC light transport simulation platform was developed and validated against experimental measurements in tissue phantoms and used to predict depth sensing in tissue. It is hoped that the current package and reported results provide the research community with an effective simulating tool to improve the development of clinical applications of RS.

Introduction

Raman spectroscopy (RS) has been under considerable investigation in biomedical optics during the last several decades. Compared to more established biomedical probing and imaging techniques [e.g., fluorescence imaging, magnetic resonance imaging (MRI), computed tomography (CT), nuclear imaging], RS offers several advantages. It relies on nonionizing radiation; it can be integrated into standard medical workflows because of its potential integration into compact fiber optics systems; it can be either label-free or applied to image-specific targeted ligands; and furthermore, it offers easy sample preparation and compatibility with aqueous solvents. However, there are two major drawbacks to the application of label-free RS in biomedical optics. Because of its low conversion rate, the Raman scattering cross section is weak compared to other light–tissue interaction mechanisms, which makes Raman imaging a slow process. Moreover, Raman scattering is more likely to occur in the Stokes regime where photons are redshifted, which makes the Raman photon detection process prone to contamination by potentially large fluorescence background from native biomolecules in tissue and cells.

Practical clinical applications of RS for tissue imaging (e.g., surgical guidance, targeted biopsy collection, and treatment monitoring) were mostly developed using fiber optics systems without regard to controlling tissue sensing depth. As a result, these approaches can be inherently imprecise when used for surgical guidance or targeted biopsy collection. Moreover, ambiguities in tissue sampling depth make the whole process prone to modeling errors when using pathology analysis to train predictive models based on supervised machine learning. Therefore, deriving a rigorous relation between the Raman sensing depth for specific imaging systems and tissue optical properties (absorption and elastic scattering) will be essential for the development of real-world RS applications in medicine. More recently, advances were made in depth-resolved techniques and developing spatial-offset Raman spectroscopy (SORS) with various types of probes and optical setups to harvest Raman signals at depth. These innovative techniques would also greatly benefit from robust simulation tools for modeling the depth sensitivity of RS.

Attaining a rigorous relation between Raman depth and optical properties in an RS system requires solving the radiative transfer equation (RTE) and calculating photon diffusion with predetermined boundary conditions. Light transport simulations in tissue optics traditionally involve modeling three competing interaction mechanisms: absorption and photon energy conversion into heat (resonant phenomenon), fluorescence (resonant, but redshifted re-emission), and elastic scattering (nonresonant)., These phenomena can be modeled by analytically or numerically solving the RTE; this can be simplified to the diffusion equation ( approximation) or to the simplified spherical harmonics equation ( approximation) in highly scattering media. However, these analytical/numerical approaches are limited when modeling complex imaging domains with realistic heterogeneities or curved geometries.

Biological tissue can be approximated as a random medium in which light is mainly diffused by going through many scattering incidents. This problem is either analytically impossible to solve for realistic imaging configurations or it can be computationally intensive in numerical simulations to account for absorption, elastic scattering, and fluorescence. An alternative to the analytical/numerical studies is the Monte Carlo (MC) method, which was proposed for the first time by Wilson and Adam and improved by several groups to solve the problem of light transport in tissue. Due to its stochastic nature, MC is essentially capable of dealing with any level of complexity concerning optical or geometrical properties in the RTE, as long as appropriate computational platform (hardware and software) is provided. The advent of strong microprocessors boosted with modern graphical processing units (GPU) accompanied by robust parallelization algorithms in recent years have paved the road for MC to produce rigorous and reliable solutions to light transport in tissue. In MC simulations, photons undergo sequential random walks and, during each step, the corresponding optical parameters (absorption, elastic scattering, fluorescence, and Raman scattering) can be randomly sampled to decide which of the four competing events happens along the walk.

In this article, a new MC package is introduced that can simulate elastic and inelastic scattering of photons in biological tissue. First, a review of the previously reported work on MC analysis of photon transport is provided and the advantages of the package presented in this article are highlighted. Then, the stochastic analysis methodology of elastic and inelastic diffusion of photons is presented. To validate the new package, simulation results in realistic tissue phantoms are presented and compared with those of experimental phantoms of known absorption and scattering coefficients. Detailed Raman sampling depth analysis is then presented to establish relationships between sensing depth, tissue optical properties, and imaging geometry parameters. This work sets the stage for the development of more controlled experimental protocols leading to the clinical translation of RS systems into medical applications.

Monte Carlo Simulation Strategy

Overview of MC Simulators in Biomedical Optics

Various MC algorithms have been implemented to stochastically solve the RTE and study light transport in turbid media, including biological tissue. Monte Carlo Multi-Layered (MCML) was the first MC package developed by Wang et al. in standard C language to model photon diffusion in layered tissue. MCML initially had two limitations: the first was its long computational time, which could take many hours for a single simulation, and the second was its restriction to model only layered geometries and its inability to simulate nonlayered structures. Concurrent with MCML, Wang et al. released convolutional MCML (CONV MCML) to simulate the interaction between a finite-width light beam and multilayered tissues. The advent of more powerful computational methods as well as the development of multicanonical MC simulators helped to overcome the first drawback of MCML. Hybrid models were also introduced to model light transport in more complex geometries, allowing inclusion of simple geometrical shapes such as cuboids, spheres, and cylinders in layered structures., Among the reported results obtained from the hybrid models were MC study of optimized light delivery for tumor laser treatment, MC modeling of light delivery and focusing in tissue with blood vessels, arteries, and capillaries, MC analysis of photoacoustic imaging of sentinel lymph nodes, MC study of near-infrared (NIR) light propagation within adult and neonatal head models, and grid-based modeling of skin under laser irradiation. In recent decades, the hybrid models were superseded by more versatile techniques such as CUDAMCML, mesh-based MC, MC eXtreme, tetrahedron-based inhomogeneous MC optical simulator, and voxel-based MC.

In addition to the traditional MC simulation techniques, where mainly elastic scattering and absorption of light in tissue are modeled, a substantial number of MC studies investigating fluorescence and the Raman effect in turbid media have been reported. Modification of the primary MC approaches incorporating the quantum yield and emission spectra of specific fluorophores were used to stochastically model fluorescence in scattering and absorbing media. In two consecutive papers,, MC simulations were applied to empirically model the fluence rate and fluorescence re-emission as a function of effective penetration depth and diffuse reflectance. Beuthan et al. showed that (NADH) concentration in turbid media can be estimated through simultaneous detection of fluorescence and light backscattering. This work was extended by Minet et al. to estimate the concentration of NIR fluorophores. In another work, a fast Fourier method and scaled fitting procedures were used to improve the speed and accuracy of the MC simulations in reconstructing the fluorescence spectra.,

A deficiency of existing MC tools is that they are unable to simultaneously account for all tissue optics competing mechanisms, including not only absorption/elastic scattering and fluorescence, but also nonelastic scattering, i.e., the Raman effect. Enejder et al. used an MC method to model the Raman-scattered light within blood samples in a quartz cuvette and quantified the analytes in whole blood. A semi-analytical study facilitated by an MC simulation on infinite and finite single-layer media to correct the turbidity-induced distortions in Raman spectra was presented by Shih et al.. Mo et al. presented an optical fiber Raman probe coupled with a ball lens and developed an MC simulation to study the depth-resolved Raman signal collected from a two-layer epithelial tissue. Hokr and Yakovlev developed an MC model to consider elastic scattering, absorption, and spontaneous Raman scattering in a turbid medium and they showed that an enhancement in elastic scattering leads to the growth of forward and backward Raman signals. However, due to the extensive computational costs of the proposed MC model, the authors chose artificially large values (on the order of ) for the effective Raman cross section and made the assumption that the laser pump depletion is negligible, which may hinder the given conclusion in realistic inelastic scattering light–tissue interaction measurements. An MC approach for calculating the depth sensitivity of single-fiber and multi-fiber Raman probes interrogating skin was presented by Reble et al.. Modeling skin as a two-layer geometry with a finite epidermis layer backed by a semi-infinite dermis layer, the authors demonstrated that skin with nonmelanoma cancer reveals higher sampling depth compared to normal skin. In 2014, two groups separately developed MC simulations with Raman modeling capabilities. Wang et al. built an eight-layer skin model and performed an MC calculation to study how different layers affect the measured in vivo Raman spectrum. Hokr et al. presented a model based on the MCML package and investigated the stimulated Raman scattering and nonlinear dynamics of light in turbid media. Both the proposed models were, however, limited to layered structures, which limits their applications for more complex geometries. Periyasamy et al., presented Raman MC simulations based on the MCML package with embedded objects with spherical, cuboidal, ellipsoidal, and cylindrical shapes in layered structures. MC calculations were also conducted to analyze SORS systems.,,

All the cited studies on the use of MC methods to analyze Raman scattering are either limited to layered structures with simple embedded shapes or suffer from a high computational cost in simulating Raman photons. The work presented here overcomes those difficulties by enabling rapid GPU-boosted simulation of all four competing mechanisms, including the Raman phenomenon, in turbid tissues with any arbitrary geometry or optical properties over any desired spectral range. The availability of this tool is important as it will highlight the roadmap toward the fabrication of more realistic optical phantoms, determine the Raman depth sensitivity across all possible tissue parameters, and make a quantitative assessment of expected levels of Raman signal-to-noise ratio (SNR), which is dramatically impacted by the relatively high levels of fluorescence and heterogeneity of biological tissues.

New MC Simulator and Its Advantages Over Currently Available Methods

In this section, a newly developed MC package is presented for simulating elastic and inelastic tissue light scattering. The current MC package was developed based on parallelization of the algorithm on graphics cards, which is obtained using the Open Graphics Library (OpenGL). Compared to other published Raman simulators, the developed package offers several important features that are emphasized here.

First, a marching cube algorithm facilitated by several shaders from the OpenGL shading language was used in the implementation of the medium that allows robust rendering of three-dimensional (3-D) curved geometries with locally refined structures that may contain multiple inclusions of different optical properties. In generating the geometries, the voxelated 3-D objects are saved in stack of two-dimensional (2-D) images and the marching cube algorithm performs the divide-and-conquer approach to extract isosurfaces (3-D regions with identical optical properties). More details on the applied marching cube algorithm can be found in Refs. 63 and 64. Second, the package was developed essentially by incorporating several shaders of OpenGL, including vertex shaders, fragment shaders, and geometry shaders. By means of these shaders, all the MC calculations regarding light transport, processing the initial and final positions, directions, and wavelengths of photons, reflection and refraction of photons, geometry rendering, as well as camera/sensor implementation are performed on the GPU. Consequently, the running time of each simulation is shorter compared to CPU-based packages. Furthermore, due to the flexibility of OpenGL and its compatibility with many graphic cards, the developed package does not require a very specific architecture and can be run on any computer with a graphic processor supporting OpenGL 4.1 or above. Thus, significant benefit in portability, implementation, and maintenance of the package is brought forth. Third, the package parallel implementation is robust and dynamic such that the wavelength shift of photons, due to inelastic scattering and fluorescence events, and their diffusion at the shifted wavelengths can be simulated concurrently with the propagation of photons at the launching wavelength. Therefore, the MC package is capable of simulating both fluorescence and Raman spectra, i.e., phenomena at multiple wavelengths. Fourth, the demonstrated package offers easy access to information regarding the final location and direction of photons, their final wavelengths, as well as the locations where inelastic scattering, fluorescence, or absorption events occurred. Having access to this information provides users with the ability to perform various types of post-processing analysis, such as spatial or temporal filtering. Fifth, the package allows for inputs from structural imaging modalities such as MR and CT, which allows simulations in realistic situations through different organs and within contrast-enhancing tumors. Sixth, the MC simulator provides the user with several light source and sensor options. For illumination, different types of sources such as point source, wide field source with the Gaussian profile, and patterned illumination for spatial frequency-domain imaging can be modeled. The detection features include single-point fiber optics detection and wide-field camera-based detection.

Stochastic Modeling of Light–Tissue Interaction and Random Variables

In a turbid medium, photons are traveling random paths composed of a sequence of straight segments. Each path can be interrupted by one of the following four phenomena: elastic scattering (Rayleigh or Mie), absorption, fluorescence, or inelastic (Raman) scattering. To include each of these events in the MC analysis, their respective phenomenological physical constants are used to build a probability density function (PDF) and hence, sample random numbers. For elastic scattering, the scattering coefficient represents the scattering probability per unit distance along each photon trajectory. Based on a modified version of the Beer–Lambert law, the PDF and cumulative distribution function (CDF) for elastic scattering are defined as follows: The function is the CDF, which represents the probability that a photon has gone through a scattering event after having traveled a distance . Then for the length of the photon path between two diffusion events (), a random number , between 0 and 1, is generated and following the elastic scattering CDF, we have which gives

In addition to the diffusion length, the direction of elastic scattering is determined from the phase function, which is essentially the angular probability density of a photon coming from solid angle direction being scattered into direction , i.e., . For the case of unpolarized light, the probability of scattering is equally distributed for all the angles in the azimuthal plane, i.e., and hence . To determine the azimuthal direction of the scattering, a random number , distributed uniformly between 0 and 1, is generated and the azimuthal direction of scattering is obtained as . Then, by assigning the Henyey–Greenstein phase function, to , the CDF of the scattering angle in the polar plane is By letting , where is distributed between 0 and 1, the scattering orientation along the polar plane is, thus, where is the anisotropy coefficient.

The absorption coefficient specifies the probability of absorption per unit distance traveled by a photon with corresponding CDF, which indicates the probability that an absorption event occurs within a traveled distance . Some MC approaches maintain a packet weight index for each photon, which is decreased by absorption during the photon travel. Here, the “intensity” of the photons is unaltered along their paths and they are considered as quantized (albeit nonpolarized) particles. Absorption is considered on the same footing as the other interaction mechanisms with the probability of photon survival based on the Russian roulette mechanism. The random number of the roulette wheel is compared with to determine whether the photon survives or annihilates at a given iteration. The survival condition of a photon is, therefore, determined by the following expression: where is a random number between 0 and 1.

The probabilities of Raman scattering or fluorescence re-emission are related to the molecular structure of the biological tissue and can be modeled by their emission spectra. If the re-emission bandwidth of interest is and the wavelength of illumination is , the fluorescence and Raman CDFs, which are specifying the probability of fluorescence or Raman shifts from to , are where and are the conversion rates of fluorescence and Raman events, which are essentially proportional to their re-emission spectra. This is because if we consider fluorescence and Raman scattering as absorption/re-emission and scattering phenomena, respectively, then the re-emission intensities in fluorescence and scattering intensity in Raman scattering are directly proportional to the corresponding fluorescence quantum yield and Raman cross section, respectively. These two intensities are, in fact, the inherent fluorescence and Raman spectra of medium under interrogation. It should be noted that in the current package, the CDFs defined in Eqs. (9) and (10) are indicated by and , respectively. To decide whether, on each diffusion path, the inelastic scattering or fluorescence re-emission occurs or not, the MC package uses Eqs. (9) and (10) to calculate and for all the wavelengths of interest. Then it produces two random numbers, and , and compares them with all the calculated and . If the produced random numbers are smaller than the CDF value of a particular wavelength, then that wavelength will be chosen as the inelastic wavelength. If multiple wavelengths have their CDF values larger than the generated random numbers, then the shortest wavelength among them will be chosen. If the generated random numbers exceed all the CDF values of the considered wavelengths, then the diffusion would be elastic.

Finally, the reflection and refraction of photons at the interface of two layers are computed using the Fresnel coefficients. If the photon encounters an interface between two layers with refractive indices and , then the package generates a random number , between 0 and 1, and compares it with the reflectance of the interface, given by the following Fresnel relation: where is the angle of incidence and is the angle of transmittance. In Eq. (11), is the average reflectance for the two orthogonal polarizations of light, for in the current context of interest, we have assumed that light is unpolarized. Therefore, if is smaller than , then the photon is reflected, otherwise it is transmitted.

Algorithm Managing Light–Tissue Interaction Events

Implementation details of the MC algorithm are shown in Fig. 1. To simulate elastic scattering, absorption, fluorescence, or Raman scattering, their respective physical constants are used to build pertinent PDFs and hence, sample random numbers. At the start of each step, the devised algorithm samples seven different random numbers, i.e., to as summarized in Table 1, corresponding to the length of diffusion step, azimuthal and polar directions of the diffusing photon, probabilities of absorption, florescence, and Raman scattering, and specular reflection if the next diffusion increment crosses a boundary between two media. Then, the algorithm takes the generated random numbers into account and follows survival tests for each of the above-mentioned events to decide which event dominates at that step.

Fig. 1

Proposed MC algorithm including all the competing events in light–tissue interaction.

Table 1

Generated random numbers.

Name	Variable	Corresponding equation
Diffusion free path	l_diff	l_diff=−ln[(1−ξ_1)/μ_s]/μ_s
Azimuthal direction of diffusion	ϕ	ϕ=2πξ_2
Polar direction of diffusion	θ	θ=arccos{12g[1+g^2−(1−g^21−g+2gξ_3)^2]}
Russian Roulette	ξ_4	ξ_4>1−e^−μ_al_diff
Probability of fluorescence shift	ξ_5	ξ_5>ρ_F(λ)
Probability of Raman shift	ξ_6	ξ_6>ρ_R(λ)
Probability of reflection at a surface	R	ξ_7<R

More specifically, after launching photons, the package first checks if there is any Raman shift. If there is a Raman shift, it updates the wavelength and optical properties of the medium and then computes the diffusion length [Eq. (4)] as well as azimuthal () and polar [Eq. (6)] directions, respectively. If there is no Raman shift, the package computes the diffusion length, azimuthal, and polar directions based on the optical properties at the wavelength of illumination. Then, the simulator checks if there is any change of interface along the diffusion segment and in the direction of the initial and . If there is an interface on the way, then the photon is either transmitted or reflected based on the Fresnel laws, and accordingly, the azimuthal and polar directions of subsequent diffusion ( and ) are updated. If there is no interface, then the photon travels along the path determined by and the initial and . At the next step, the photon undergoes the Russian roulette process and competes with absorption. If there is no absorption, the algorithm goes back to the step of Raman shift checking. If there is an absorption event, then the algorithm checks whether the photon is annihilated or undergoes fluorescence re-emission. In the case of annihilation, the package launches a new photon. In the case of fluorescence re-emission, the algorithm updates the angles and , as fluorescence is assumed to be an isotropic process. Furthermore, the algorithm checks if the fluorescence re-emission is at the initial wavelength or not (second last box in Fig. 1). If yes, then the package returns to the Raman shift inspection step. If no, the package updates the wavelength and optical properties for the new wavelength, and then jumps back to the Raman shift inspection step. It should be mentioned that the launched photons can repeatedly go through the described algorithm in a number of iterations set initially by the user. If before reaching the maximum number of iterations, the photon is absorbed and not re-emitted, then that whole sequence is terminated, and a new photon is generated. Furthermore, the simulation domain is considered as a void (i.e., vacuum) space where all the modeled structures are contained. Hence, the exiting photons propagate along straight lines without experiencing any scattering or absorption in the surrounding void medium. More details on the technical implementation of the MC algorithm based on OpenGL can be found in Ref. 63.

Physical Environment Modeling

The optical properties (, , , , , and ) of the simulated medium need to be defined at any point in space for all the wavelengths. As is a convention in the diffuse optics literature, empirical relations for , , and are used, although the code can be adapted to be compatible with the tabulated values for these constants. In the spectral therapeutic and imaging window for biological media, the dependence of the scattering coefficient on wavelength can be empirically described by a decreasing power law. Thus, the scattering is represented by an exponential equation with coefficient and power where and are positive real numbers. The anisotropy coefficient and refractive index exhibit more linear behaviors

The dependence of absorption, fluorescence, and Raman scattering (, , and ) on wavelength cannot be modeled empirically. For this reason, the package should be provided with the values of these properties at each wavelength, adapted to any desired biological application.

The MC simulation package allows flexibility in defining the illumination and detection geometries. The light sources are represented according to their position, direction, intensity (in terms of number of photons), radius, distribution (standard deviation), and emission wavelength. Thus, a monochromatic source having a Gaussian emission profile is represented only by 10 numbers. In addition, more than one source can be modeled simultaneously by adding a supplemental definition in an input file. Thus, it is possible to model a source containing several wavelengths by superimposing several sources in the file. The source radius defines the boundary beyond which no photon is emitted, which makes the simulator suitable for optical fiber sources. The cameras and sensors used in biophysics are diverse. For this reason, it is important to leave enough flexibility at the detection port. In the current version of the MC package, the pinhole camera is supported with the generic parameters tabulated in Table 2. Allowing the user to control these parameters, the simulator offers a solution adaptable to a wide variety of applications. Furthermore, in order to implement new sensors and cameras, the package allows the user to insert the corresponding camera matrix of the new sensing system and apply it to the captured photons.

Table 2

Generic parameters of the camera/sensor implemented in the package.

Parameter	Representation space
Position	(X,Y,Z)
Orientation	(X,Y,Z)
Resolution	(X,Y)
Field of view	(Near plane, far plane, angle of view)

Experimental Validation of the Simulation Technique

In this section, experimental measurements made in tissue phantoms using a single-point RS probe are compared with simulation results for the same geometry in synthetic phantoms in order to validate the practical utility of the new MC simulation technique.

Experimental Tissue Phantoms

Solid tissue phantoms were fabricated from polydimethylsiloxane (PDMS) and nylon matrices. The absorption and scattering properties of the phantoms were controlled by adding India ink for absorption and titanium dioxide (TiO2) powder for scattering., PDMS and nylon were chosen in part because they have Raman signatures allowing them to be clearly distinguishable. Three phantoms were fabricated [Fig. 2(a)]: (I) a single layer of PDMS (thickness: 10 mm) with absorption and scattering agents, (II) a single layer of nylon (thickness: 10 mm) with absorption and scattering agents, and (III) a two-layer phantom made of a substrate of nylon layer (thickness: 10 mm) over which a layer of PDMS was deposited. The thickness of the PDMS layer in phantom III was measured using a commercial optical coherence tomography (OCT) system (Thorlabs SL1310V1). The axial resolution of the system was estimated to be in air. Two OCT B-scans were acquired by rotating the sample by 90 deg around the axis normal to the cross section of the sample; one along direction and one along direction. The PDMS thickness was determined to be uniform () and thick over the 3.71-mm width of the B-scans.

Fig. 2

(a) Structural geometry of phantoms I (PDMS), II (nylon), and III (PDMS on the top nylon with ). Measured (b) absorption and (c) reduced scattering coefficients of nylon and synthetized PDMS samples.

Characterization studies were then conducted to determine the absorption coefficient () and reduced scattering coefficient [] of the single-layer phantoms. Specifically, reflection and transmission spectra were acquired using a commercial PerkinElmer Lambda 1050 spectrometer equipped with a single integration sphere. The optical properties were computed using the Inverse Adding Doubling technique. Figure 2(b) shows the absorption coefficients and Fig. 2(c) shows the reduced scattering coefficients of PMDS and nylon phantoms, respectively. The optical properties are consistent with biological tissue in the NIR window. For the PDMS phantom, ranges from 0.11 to and varies from 0.65 to between 500 and 1000 nm. For the nylon phantom, changes from 0.02 to and fluctuates from 1.05 to within the wavelength range 500 to 1000 nm.

Experimental Raman Spectroscopy Measurements

Fluorescence and Raman scattering spectra were measured with a handheld RS probe in contact with phantoms I, II, and III. All Raman acquisitions ( spot size) consisted of 10 consecutive spectra acquired following excitation with a 785-nm light source (25 mW at the tip, 1-s exposure per spectrum).,,, The raw signals were processed to isolate the contribution associated with the Raman effect: (1) the background signal (measurement acquired with laser turned off) was subtracted and the resulting spectrum corrected for system response by normalizing it to a measurement made on a Raman standard (SRM2241; NIST, Gaithersburg, Maryland); (2) the remaining low-frequency background (mainly fluorescence) was removed using the rolling ball algorithm; and (3) the spectrum was normalized with standard normal variate (SNV) technique to have a mean of zero and standard deviation of one. Figure 3 shows the post-processed spectra for each phantom. All expected Raman peaks for PDMS and nylon were detected as seen in Figs. 3(b) and 3(d), respectively. The measured Raman spectrum for the two-layer phantom [Fig. 3(f)] shows all peaks associated with the top PDMS layer and several peaks associated with nylon in the higher wavelength range above 800 nm.

Fig. 3

Simulated versus measured fluorescence and Raman spectra of (a), (b) phantom I, (c), (d) phantom II, and (e), (f) phantom III, respectively.

Monte Carlo Simulations

An MC light transport protocol was developed [Fig. 4(d)] to run simulations consistent with the illumination/detection geometry of the handheld single-point RS probe as well as with the geometry and optical properties of phantoms I, II, and III (Fig. 2). Figures 4(a)–4(c) show the cross-sectional and side views of the simulated source and detector configurations. The detection area was approximated by a cylinder of diameter surrounding the illumination area modeled by a diameter disk. The acceptance angles (numerical apertures) for illumination and detection were 8 deg and 30 deg, respectively.

Fig. 4

(a) Cross-sectional view of the handheld interrogating probe. The red region with a diameter of indicates the source port, while the green region with a diameter of illustrates the collecting port. (b) Side-view of the detector geometry. (c) Side-view of source geometry. (d) Elastic and inelastic photon diffusion in a two-layered sample under illumination of a light beam.

Three-dimensional (3-D) grids of voxels were generated to emulate the geometry of the phantoms: for the single-layer phantoms I and II and for the two-layer phantom III. Optical properties parameters were assigned to each voxel to match conditions met in the experimental tissue phantoms. Those included a refractive index of and the measured absorption and reduced scattering coefficients reported in Fig. 2. The measured background (the residual after applying the rolling ball algorithm) and resulting Raman spectra of the single-layer PDMS and nylon were used to assign relative probabilities at each voxel associated with the fluorescence and Raman scattering processes, respectively. As noted while presenting Eq. (10), the fluorescence and Raman conversion probabilities are proportional to their corresponding spectra of re-emission and scattering, respectively. Accordingly, the acquired spectra from single layer PDMS and nylon were multiplied by a proportionality constant and used as probability conversion rates ( and ) in the CDF manner at 100 discrete wavelengths equally distributed between 810 and 932 nm. Tissue excitation was modeled by a monochromatic point source with an excitation wavelength of 785 nm. All simulations were performed on a DELL Precision 7920 Tower workstation with an 8 core Intel Xeon Silver 4110 microprocessor (2.1 GHz, 3.0 GHz Turbo), an NVIDIA Quadro P4000 graphic card, and OpenGL 4.5. As illustrated in Fig. 4(d), a total of photons in 3600 batches were launched onto the center of the top face for each phantom to simulate Raman scattering and fluorescence re-emission at 100 discrete wavelengths equally spaced between 810 and 932 nm. The required time for the MC simulating and subsequent post-processing calculations were 5.5 and 4 s per batch, respectively. The source was approximated as a Gaussian distribution with standard deviation of . A post-processing analysis was applied at the detection end to count the Raman photons collected within the aperture of the detectors.

The simulations for each phantom resulted in spectra that were post-processed using the same background removal technique as for the experimental spectra. Figures 3(a), 3(b) and 3(c), 3(d) show the SNV-normalized experimental and simulated spectra for the pure PDMS and nylon phantoms, respectively. The post-processed simulated and experimental spectra are in agreement, preliminarily demonstrating the efficacy of the method to simulate inelastic Raman scattering in a diffusive medium. Figures 3(e) and 3(f) show the simulated and experimental results for the two-layer phantom, respectively. The two spectra show overall good agreement with minor differences between the measured and simulated curves at longer wavenumber shifts, where most Raman nylon peaks of the nylon substrate are located.

Inelastic Scattering Sampling Depth

Simulation Geometry and Depth Sensing Evaluation Metrics

A comprehensive study of the effects of absorption, elastic, and inelastic scattering on the depth of Raman conversion was conducted. Because of light diffusion in tissue, each optical measurement was associated with a distribution of photons having sampled different sensing depths. From the Beer–Lambert law, the fraction of detected photons having sampled a given depth exponentially decreases with depth. In the analysis presented here, Raman sensing depth is defined as the depth after which a specific cumulative percentage of detected Raman photons was scattered. To ensure that the detected signals represented a realistic sensing depth, two values were computed as references. Specifically, depth sensing values were computed associated with the depths after which 75% or 90% of all Raman photons were generated in the face of all other competing interaction mechanisms.

To conduct the Raman depth sampling analysis, a series of MC simulations were done. In each simulation, the in silico phantom was first assumed to be a semi-infinite space (i.e., ). The illuminating source was modeled by a monochromatic point source (), which had a Gaussian radial distribution with a standard deviation of and was located at the center of the top face of the phantom. However, to minimize computational time and required memory, the simulated phantom was approximated as a slab of , which was modeled by a 3-D grid composed of voxels. For each simulation, in 300 batches were launched, where the required time for the MC simulation and subsequent post-processing analysis were 5.5 and 4 s per batch, respectively. As for the collection part, a post-processing analysis was performed to count the Raman photons collected by an aperture identical to the detector in Fig. 4.

Impact of Inelastic Scattering Conversion Rate on Raman Depth Sampling

The effect of the Raman conversion rate on the Raman depth was first studied. Figure 5(a) shows the cumulative percentage of detected Raman photons versus the Raman penetration depth for a phantom with and , and Raman conversion rate varying between and . As observed in this figure, increasing the Raman conversion rate increases the slopes of the curves, which indicates that the Raman depth is approximately inversely proportional to the conversion rate. This relation is more clearly visible in Fig. 5(b), where the 75% and 90% Raman depths are plotted for various conversion rate values. It is observed that as the Raman conversion rate increases, most of the propagating photons undergo a Raman shift within shallower layers.

Fig. 5

(a) Cumulative percentage of detected Raman photon versus penetration depth for various Raman conversion rates within a phantom with and . (b) 75% and 90% Raman depths versus conversion rate.

Impact of Absorption and Elastic Scattering on Raman Depth Sampling

The effects of absorption and elastic scattering on Raman sensing depth are illustrated in Figs. 6 and 7. Figure 6(a) shows the photon intensity curves for and various absorption coefficients. In this figure, it is shown that for a highly absorptive medium (), the slope of the cumulative captured photon intensity curve is larger compared to those with lower absorption coefficients, indicating shallower Raman sensing depths for more absorptive media. Figure 6(b) shows the photon intensity curves for a higher scattering coefficient of and varying absorption coefficients. For the higher scattering coefficient, all the intensity curves become much sharper, indicating much shallower Raman sensing depths for more scattering media. This indicates the dominance of elastic scattering event on the depth of Raman re-emission. These two facts are more visible in Fig. 7, where the 75% and 90% Raman sensing depths are traced versus the reduced scattering and absorption coefficients, respectively. Comparing Figs. 7(a), 7(c) and 7(b), 7(d) reveals that the effect of absorption coefficient on the Raman sensing depth is less significant than that of the scattering coefficient. According to these figures, for a constant value, increasing results in a smaller decrease in the 75% and 90% Raman sensing depth compared to the depth decrease observed by increasing for a fixed value. Furthermore, as grows, the 75% and 90% Raman depth curves versus become flatter, which indicates that the Raman sensing depth is less influenced by the absorption in media with stronger scattering coefficients. In other words, for higher amounts of , the Raman sensing depth curves versus behave like parallel lines, which shows that for stronger scattering coefficients, the amount of change in the sensing depth is progressively independent of . Finally, as can be seen in Figs. 7(a) and 7(c), the blue curves are not following the general expected behavior at the last two points, where kinks appear. This is due to the fact that, at such points, the absorption coefficient is much larger than the reduced scattering coefficient, and hence, estimating the photon paths and fluence rate is out of the realm of the diffusion equation. Nevertheless, such extreme values for absorption and reduced scattering coefficients are rarely seen in biological tissues.

Fig. 6

Cumulative percentage of detected photons versus penetration depth for various absorption coefficients of a phantom with and (a) , (b) .

Fig. 7

(a) and (c) 75% and 90% Raman depth versus absorption coefficient for various and . (b) and (d) 75% and 90% Raman depths versus scattering coefficient for various and .

Discussion

Prior to being used in human subjects, optical imaging instruments should undergo an iterative design procedure facilitated by tissue phantoms with realistic geometrical and optical properties such that a reliable assessment of the imaging performance and utility can be obtained. In this regard, the optimized system-specific parameters can be estimated by optical simulation of the designed imaging instrument, e.g., by detailed ray tracing analysis of included optical components as well as their arrangements. Tissue-specific technical features are usually estimated based on tissue phantom experiments, which model a range of optical properties representing the expected inter- and intrapatient heterogeneity, potential impact of specular reflections, and possible variations in tissue surface characteristics including curvature and texture. However, in order to complement the phantom analysis and estimate a wider spectrum of optical properties of the tissue under interrogation, in diffuse optics, numerical light transport simulations have been invoked in a number of modern clinical applications. A few examples of such numerical studies are fluorescence imaging systems for surgical guidance, source–sensor geometries in diffuse optics applications, and diffuse optical tomography systems.

This article presented the development of a spectroscopic MC code to simulate both elastic and inelastic light interactions in biological tissue. The Raman and fluorescence spectroscopy capabilities of the MC package were validated showing agreement between simulations and experimental measurement in pure materials, here either nylon or PDMS tissue phantoms. In layered materials (a thin layer of PDMS over a substrate of nylon), differences were observed between the measured and simulated spectra in the higher wavelength regime, where most of the nylon Raman peaks are located. This could potentially be due to multiple factors including the inability of the reduced scattering approximation formula incorporated in the package (Mie scattering) to effectively model scattering at larger depths, the inaccuracy in estimating the refractive index mismatch between the two layers, and/or inefficiencies of the normalization process applied to extract the Raman and fluorescence spectra.

To illustrate the analytic capabilities of the MC package, a detailed analysis was performed to estimate the depth sensing capabilities of an existing surgical-guidance handheld probe instrument. The study was performed with the objective of quantifying the impact of tissue optical properties (absorption and elastic scattering) and the relative probability of occurrence of the Raman effect on inelastic scattering depth sampling. Although realistic values for biological tissue parameters such as absorption and elastic scattering can be found in the literature, no realistic values were available to estimate the relative probability, or physical cross section, of inelastic scattering (herein referred to as Raman conversion rate) in biological tissue. As a result, values associated with different orders of magnitude for this parameter were simulated to evaluate its impact on depth sampling (Table 3).

Table 3

Optical properties used in the sensing analysis and corresponding calculated Raman sensing depths.

Parameters				Depth of 75% sensitivity (μm)	Depth of 90% sensitivity (μm)
Absorbtion μ_a (mm^−1)	Scattering μs′ (mm^−1)	Mean free path (μm)	P_Raman
0 to 1.4	0.5	200	10^−6	270 to 600	75 to 160
	4	25		140 to 190	41 to 80
	30	3.3		10 to 80	15 to 25
0.001	1	100	10^−3 to 10^−6	105 to 225	40 to 75

It was shown that in media with higher inelastic cross-section values, the rate of Raman conversion at superficial layers is higher and the sensing depth associated with the probing instrument decreases accordingly. This is expected because when the Raman cross section increases, the likelihood of Raman conversion of photons as they diffuse across shallow layers grows, which results in a higher percentage of superficial Raman scattering events. Since the probability of two consecutive Raman shifts for a single photon is negligible, the likelihood of another Raman scattering at deeper layers for the already-Raman-shifted photons becomes vanishingly low, and hence, most of the detected Raman photons are those re-emitted from upper layers. However, for lower conversion rates, the overall intensity of the collected Raman signal is also lower. Therefore, reducing the Raman conversion rate for the sake of increasing the Raman depth is at the expense of lowering the count associated with inelastically scattered photons. As a result, there is a trade-off between Raman depth sensing and Raman SNR. It was also shown that absorption has a direct impact on Raman sensing depth: larger tissue absorption values lead to smaller sensing depths. This makes sense intuitively, since for smaller absorption coefficient values the likelihood that photons diffuse over longer distances increases, hence, a larger relative fraction of elastic scattering events occurs.

In the course of photon transport within tissue, among the four competing events, elastic scattering acts as a stimulator for other three events. In other words, the higher the elastic scattering coefficient, the shorter would be the diffusion length between two consecutive scattering events, the greater would be the number of collisions between photons and molecules, and therefore, the higher would be the chance of absorption, fluorescence, and Raman scattering. In a more quantitative description, in a medium with a large scattering coefficient, the spatial variation of electric polarizability , which mainly contributes to the amplitude of re-emitted Raman signal, is large. Therefore, a positive incremental addition to the scattering coefficient, such as , leads to a positive incremental increase in the Raman conversion rate, i.e., . Hence, as discussed in the previous section, as the scattering coefficient becomes larger, the Raman conversion rate increases, and the event of Raman shift at upper layers becomes more likely. Therefore, the Raman sensing depth in tissues or phantoms with higher scattering coefficients is less than the depth in structures with lower scattering coefficients.

As a final point, it should be mentioned that the discussed MC algorithm is a direct method of light transport calculation, which may require large number of photons, especially in situations where one of the competing processes has a very small probability to trigger this. As reviewed comprehensively in Ref. 21, there are several options to accelerate the MC modeling of light transport in turbid media, including scaling, perturbation, and convolution. In these accelerated techniques, the MC analysis is performed indirectly by modifying an existing base MC calculation. As noted earlier, the developed package gives the option of recording the final location, direction, and wavelength of the simulated photons. Hence, in principle, a user can readily use those recorded data as a base simulation and perform accelerated MC analysis in new cases. Following the given explanations and shown results, it is hoped that the developed package will bring more understanding in the area of light transport modeling and equip the contributing research community with a versatile toolset for optical detection and treatment as well as surgical guidance.