2.1.

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.^26,27 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,^29–31 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^40–51 and the Raman effect^52–62 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.^40–45 In two consecutive papers,^46,47 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 ${\mathrm{C}}_{21}{\mathrm{H}}_{27}{\mathrm{N}}_7{\mathrm{O}}_{14}{\mathrm{P}}_2$ (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.^50,51

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 ${10}^{-2}$) 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.^59,60 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.^16,61,62

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.

2.2.

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.

2.3.

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 ${\mu }_s$ 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:¹⁸

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 $\mathrm{\Omega }$ being scattered into direction ${\mathrm{\Omega }}^{\prime }$, i.e., $p(\mathrm{\Omega },{\mathrm{\Omega }}^{\prime })$. For the case of unpolarized light, the probability of scattering is equally distributed for all the angles in the azimuthal plane, i.e., $p_{\varphi }(\varphi )=\varphi /2\pi$ and hence $p(\mathrm{\Omega },{\mathrm{\Omega }}^{\prime })=p(\mathrm{\Omega }·{\mathrm{\Omega }}^{\prime })=p(\cos \theta )$. To determine the azimuthal direction of the scattering, a random number ${\xi }_2$, distributed uniformly between 0 and 1, is generated and the azimuthal direction of scattering is obtained as $\varphi =2\pi {\xi }_2$. Then, by assigning the Henyey–Greenstein phase function^17,18 to $p(\cos \theta )$, the CDF of the scattering angle in the polar plane is

The absorption coefficient ${\mu }_a$ specifies the probability of absorption per unit distance traveled by a photon with corresponding CDF,

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 $[{\lambda }_a,\lambda ]$ and the wavelength of illumination is ${\lambda }_i$, the fluorescence and Raman CDFs, which are specifying the probability of fluorescence or Raman shifts from ${\lambda }_i$ to $\lambda$, are

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 $n_1$ and $n_2$, then the package generates a random number ${\xi }_7$, between 0 and 1, and compares it with the reflectance of the interface, given by the following Fresnel relation:¹⁸

2.4.

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., ${\xi }_1$ to ${\xi }_7$ 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.

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 ($\varphi =2\pi {\xi }_2$) 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 $\varphi$ and $\theta$. 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 ($\varphi$ and $\theta$) are updated. If there is no interface, then the photon travels along the path determined by $l_{\mathrm{diff}}$ and the initial $\varphi$ and $\theta$. 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 $\varphi$ and $\theta$, 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.

2.5.

Physical Environment Modeling

The optical properties (${\mu }_s$, ${\mu }_a$, $g$, $n$, ${\rho }_F$, and ${\rho }_R$) 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 ${\mu }_s$, $g$, and $n$ 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 $a_1$ and power $a_2$

The dependence of absorption, fluorescence, and Raman scattering (${\mu }_a$, ${\rho }_F$, and ${\rho }_R$) 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.

3.1.

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 (TiO₂) powder for scattering.^65,66 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 $913\text{-}\mu \mathrm{m}$ 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 $20\mu \mathrm{m}$ in air. Two OCT B-scans were acquired by rotating the sample by 90 deg around the $z$ axis normal to the cross section of the sample; one along $\widehat{x}$ direction and one along $\widehat{y}$ direction. The PDMS thickness was determined to be uniform ($\pm 10\mu \mathrm{m}$) and $913\text{-}\mu \mathrm{m}$ 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 $d=913\mu \mathrm{m}$). Measured (b) absorption and (c) reduced scattering coefficients of nylon and synthetized PDMS samples.

Characterization studies were then conducted to determine the absorption coefficient (${\mu }_a$) and reduced scattering coefficient [${\mu }_s^{\prime }=(1-g){\mu }_s$] 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, ${\mu }_a$ ranges from 0.11 to $0.075{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }$ varies from 0.65 to $0.45{\mathrm{mm}}^{-1}$ between 500 and 1000 nm. For the nylon phantom, ${\mu }_a$ changes from 0.02 to $0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }$ fluctuates from 1.05 to $0.25{\mathrm{mm}}^{-1}$ within the wavelength range 500 to 1000 nm.

3.2.

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 ($500\text{-}\mu \mathrm{m}$ 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).^6,10,69,70 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.

3.3.

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 $600\text{-}\mu \mathrm{m}$ diameter surrounding the illumination area modeled by a $400\text{-}\mu \mathrm{m}$ 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 $400\mu \mathrm{m}$ indicates the source port, while the green region with a diameter of $600\mu \mathrm{m}$ 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: $256\times 256\times 256$ for the single-layer phantoms I and II and $501\times 501\times 203$ 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 $n=1.4$ 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 (${\rho }_F$ and ${\rho }_R$) 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 $7.2\times {10}^8$ 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 $0.1\mu \mathrm{m}$. 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.

4.1.

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 $d$ 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., $z<0$). The illuminating source was modeled by a monochromatic point source ($\lambda =785\mathrm{nm}$), which had a Gaussian radial distribution with a standard deviation of $0.1\mu \mathrm{m}$ 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 $2\times 2\times 2{\mathrm{cm}}^3$, which was modeled by a 3-D grid composed of $256\times 256\times 256$ voxels. For each simulation, $6\times {10}^7\text{photons}$ 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.

4.2.

Impact of Inelastic Scattering Conversion Rate on Raman Depth Sampling

The effect of the Raman conversion rate ${\rho }_R$ 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 ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and ${\mu }_a=0.001{\mathrm{mm}}^{-1}$, and Raman conversion rate ${\rho }_R$ varying between ${10}^{-6}$ and ${10}^{-3}$. 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 ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and ${\mu }_a=0.001{\mathrm{mm}}^{-1}$. (b) 75% and 90% Raman depths versus conversion rate.

4.3.

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 ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and various absorption coefficients. In this figure, it is shown that for a highly absorptive medium (${\mu }_a=0.1{\mathrm{mm}}^{-1}$), 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 $10{\mathrm{mm}}^{-1}$ 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 ${\mu }_s^{\prime }$ value, increasing ${\mu }_a$ results in a smaller decrease in the 75% and 90% Raman sensing depth compared to the depth decrease observed by increasing ${\mu }_s^{\prime }$ for a fixed ${\mu }_a$ value. Furthermore, as ${\mu }_s^{\prime }$ grows, the 75% and 90% Raman depth curves versus ${\mu }_a$ 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 ${\mu }_s^{\prime }$, the Raman sensing depth curves versus ${\mu }_a$ behave like parallel lines, which shows that for stronger scattering coefficients, the amount of change in the sensing depth is progressively independent of ${\mu }_a$. 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 ${\rho }_R={10}^{-5}$ and (a) ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$, (b) ${\mu }_s^{\prime }=10{\mathrm{mm}}^{-1}$.

Fig. 7

(a) and (c) 75% and 90% Raman depth versus absorption coefficient for various ${\mu }_s^{\prime }$ and ${\rho }_R={10}^{-6}$. (b) and (d) 75% and 90% Raman depths versus scattering coefficient for various ${\mu }_a$ and ${\rho }_R={10}^{-6}$.