2.1.

Monte Carlo Models

All Monte Carlo codes used in this study simulated a single layer of tissue with homogeneous optical properties and a refractive index of 1.4. The same procedures were used for selecting photon step size, absorption/scattering, and boundary interactions as those described by Wang et al.²³ The Henyey–Greenstein phase function³³ was used to calculate the scattering angle of photons^23,34,35 and all simulations assumed an anisotropy factor ($g$) of 0.9. This was considered appropriate as the anisotropy factor of brain tissue is close to 0.9 across the spectral range of interest.³⁶

The use of a single anisotropy coefficient is a common approach taken when modeling reflectance data.^14,15,27 Although decreasing $|g|$ means that scattering is more isotropic and therefore shifts spatial reflectance curves toward shorter source–detector separations, it has been demonstrated that when $g>0.8$, changes to the anisotropy factor have a minor effect on the reflectance.^24,25

The absorption and reduced scattering coefficients (${\mu }_a$ and ${\mu }_s^{\prime }$) used in these simulations were chosen to be similar to those of human brain tissue (both gray and white matter) at various wavelengths. The wavelengths investigated in this study were from 360 to 950 nm, as these are the limits of the spectrometer to be used in the final DRS instrument. The optical properties of brain tissue at these wavelengths were taken from results published in a study by Yaroslavsky et al.,³⁶ and representative values are shown in Table 1. These values were used to define the range over which a DRS instrument used in the brain would be expected to operate, and we therefore ensured that the optical property algorithms were effective within this range.

Table 1

Optical properties of brain tissue (Ref. 36).

A number of Monte Carlo–based techniques are described in this section. The first is the Mono Monte Carlo method, a technique for determining optical properties from SR-DRS data using the results of a single Monte Carlo simulation. This technique was used to simulate how removing detector fibers from a DRS fiber-optic probe design would affect the recovery of optical properties. As will be shown in Sec. 3, simulations indicated that the close proximity of detector fibers to the edge of the probe invalidated the assumption of a uniform refractive index at the tissue surface. This is a common assumption when modeling reflectance in tissue^23,37,38 and is necessary for the spatial invariance required by the Mono Monte Carlo method. This motivated the development of a modified algorithm for the recovery of optical properties from DRS data, one that considers the true geometry of the probe and, most importantly, the different refractive indices encountered by photons leaving the tissue surface (described in Sec. 2.5).

To test the validity of these Monte Carlo methods, results of the basic Monte Carlo code were compared with those from the MCML code published by Wang et al.²³ and the results of diffusion theory calculations. The extended source diffusion theory model described by Farrell et al.³⁷ was used.

2.2.

Optical Property Recovery from the Results of a Single Monte Carlo Simulation

An inverse Monte Carlo algorithm was used to recover optical properties from SR-DRS spectra. This algorithm uses the iterative Marquardt–Levenberg curve-fitting routine to fit Monte Carlo–generated spatial reflectance curves to experimental data (in this work, simulated data were used in place of experimental data). Since running repeated Monte Carlo simulations is prohibitively slow, a Mono Monte Carlo technique was used, based on the condensed simulation approach introduced by Graaff et al.²⁴ and adapted by others.^25,38

This approach used a single reference simulation with optical properties ${\mu }_{a\mathrm{ref}}$ and ${\mu }_{s\mathrm{ref}}^{\prime }$ to determine the number of interactions of each photon within the medium. In such a simulation, changing the scattering coefficient only affects the path lengths followed by photons between interactions. For this reason, the results of a single Monte Carlo simulation can be applied for a sample with new optical properties so long as the distances in the original simulation can be scaled and the new albedo $c={\mu }_s/({\mu }_a+{\mu }_s)$ is less than or equal to $c_{\mathrm{ref}}$, the albedo of the reference simulation.²⁴

The reference simulation was run for ${2\times 10}^9$ photons with optical properties ${\mu }_{a\mathrm{ref}}=0$ and ${\mu }_{s\mathrm{ref}}^{\prime }$. Note that any value of ${\mu }_{a\mathrm{ref}}$ could be used, but setting the absorption to 0 guaranteed that $c_{\text{ref}}$ was always higher than or equal to $c$. Upon exiting the tissue, each photon was binned in a grid location defined by ${\rho }_{\mathrm{ref}}$, the photon’s radial distance from the source, and $n_{\mathrm{scat}}$, the number of scattering events it experienced while travelling through the tissue. This provided a table of the reflectance $R_{\mathrm{ref}}$ at discrete values of ${\rho }_{\mathrm{ref}}$ and $n_{\mathrm{scat}}$ as illustrated in Fig. 1.

Fig. 1

Photon scoring grid used to set up the Mono Monte Carlo method. ${\rho }_{\mathrm{ref}}$ is the radial distance from the source at which a photon exited the tissue and ${\mathrm{n}}_{\mathrm{scat}}$ is the total number of scattering events the photon experienced.

This grid could then be used to determine the reflectance for a sample with optical properties ${\mu }_a$ and ${\mu }_s^{\prime }$ by applying Eq. (1). First, each value of $R_{\mathrm{ref}}$ was multiplied by the new albedo raised to the power of $n_{\mathrm{scat}}$. This could be done since the albedo provides a measure of the fraction of photons that are scattered out of a collimated beam in one path length. The grid was then collapsed down to one dimension by taking the sum along each column.

The reflectance at a given source–detector distance (SDD) $R(\rho )$ could be found after first scaling the distances in the collapsed grid by the ratio of the new mean free path ($1/{\mu }_t$) to the mean free path used in the reference simulation ($1/{\mu }_{t\mathrm{ref}}$). Because scoring was performed radially, the boundaries of each bin were concentric circles around a central (source) point. The scaling procedure moved the boundaries of the bins inward or outward and therefore changed the bin positions and widths (areas). To compensate for this, each extracted reflectance value was normalized by the total area of the scaled bin in which it was scored. This area was found by subtracting the area of the circle defined by the inner boundary from the area of the circle defined by the outer boundary.

Values of $\rho$ were used that corresponded to the positions of optical fibers in a DRS probe, and extracted spatial reflectance curves were fitted to SR-DRS data using the Marquardt–Levenberg curve-fitting algorithm. When used with experimental reflectance data, this inversion procedure allows determination of the optical properties of a measured sample.

2.3.

Reducing the Probe Size

Practically, to allow ease of use in small surgical cavities, the diameter of the fiber-optic probe will be kept quite small. This meant that any data being fitted could contain only a few short SDDs near the front of a radially resolved reflectance curve. A simulation-based procedure was used to observe the precision in the recovered optical properties as the outermost SDD was reduced.

To do so, a low noise reflectance curve was generated (using Monte Carlo methods), and random Gaussian noise was added with standard deviations of 5, 10, and 15% commensurate with experimental conditions observed during preliminary investigations. These data were fitted using the Mono Monte Carlo method to recover ${\mu }_a$ and ${\mu }_s^{\prime }$. This was repeated 1000 times, and the standard deviation in the recovered optical properties was determined. Previous investigations suggested that one thousand repetitions were sufficient to ensure that the uncertainties in the results approached their limiting values. This entire process was repeated as the maximum radial distance was reduced, for various input optical properties. In these simulations, reflectance data were generated at evenly spaced distances of 0.3 mm, corresponding to typical fiber diameters.

Reducing the maximum SDD in the probe had the additional effect of changing the depth of tissue probed by reflectance measurements. This was investigated by simulating a probe containing a single source and detector and placing the detector at varying distances from the source. For any photon that was eventually collected by the detector, its depth was recorded each time it passed below a target radial distance. These results were then used to find the average depth probed below a target radial distance versus detector position.

2.4.

True Probe Forward Monte Carlo Model

Because of its ability to recover the reflectance at arbitrary SDDs, the Mono Monte Carlo method was effective for investigating how the number of fibers and their positioning affects optical property recovery in a custom probe. However, certain simplifying assumptions were required during the Monte Carlo simulation. These included

To keep within the small diameter constraint while maximizing the outermost SDD, detector fibers were placed very close to the outer edge of the probe. It was therefore expected that inhomogeneous boundary conditions (the difference in the refractive index of the probe and surrounding air) would have an effect on the reflectance. To investigate this effect, a forward Monte Carlo model was developed, which explicitly simulates the true probe design as closely as possible.

Each fiber was given a circular profile with a refractive index of 1.46, a representative value provided by an optical fiber manufacturer (CeramOptec Industries Inc., East Longmeadow, MA). Each fiber used for diffuse reflectance measurements was given a core diameter of 200 microns. This size was chosen as a tradeoff between maximizing the SDD of the outermost fiber and minimizing the required exposure time. Larger fibers would reduce the maximum achievable SDD and affect the performance of the optical property recovery algorithm, while smaller fibers would require longer measurement times to achieve an acceptable signal-to-noise ratio.

The source fiber was simulated by emitting photons randomly at any point within its core, at any angle within its numerical aperture. Similarly, only photons incident on a detector fiber core at a (refracted) angle within the numerical aperture were scored by the model. Photons incident on the epoxy region of the probe were removed from the simulation (absorbed by the epoxy), while those reaching the tissue surface outside of the probe diameter were reflected according to the Fresnel equations at a boundary with mismatched refractive indices of 1.4 and 1.

The results of this true probe simulation were fitted using the previously described Mono Monte Carlo method to see if it was capable of recovering the input optical properties. This was done over a range of optical properties.

2.5.

Hybrid Reflectance Recovery Algorithm

To accurately model the reflectance that will be measured with our probe design, a true probe Monte Carlo simulation is required. Because it is not spatially invariant, however, such a simulation cannot be used to generate the scoring grid in the Mono Monte Carlo model, and therefore a new algorithm for rapid generation of reflectance curves was required. One possible solution would be to generate a two-dimensional look-up table using the forward Monte Carlo model of Sec. 2.4 over a range of ${\mu }_s^{\prime }$ and ${\mu }_a$ values. Achieving sufficient resolution over the wide range of optical properties listed in Table 1, however, would require a very large number of simulations, and the time required to implement this was considered impractical. Instead, a hybrid look-up table/Mono Monte Carlo method was developed for rapid generation of reflectance curves. This allowed the use of true probe simulations but required fewer simulations than a pure look-up table approach.

The Monte Carlo simulations described in Sec. 2.4 were run for selected values of ${\mu }_s^{\prime }$, but with ${\mu }_a$ set to 0. Detected photons were scored according to the detector in which they were collected and $x_t$, the total three-dimensional path length traversed while in the tissue. Combining the results of these simulations produced the three-dimensional grid shown in Fig. 2: the reflectance in each detector for each different value of ${\mu }_{s\mathrm{ref}}^{\prime }$ and $x_t$.

Fig. 2

The photon scoring grid used in the hybrid reflectance recovery algorithm. A number of different simulations were run for different values of ${\mu }_s^{\prime }$. Photons in each simulation were scored according to the detector bundle in which they were collected and $x_t$, the total three-dimensional path length traveled in the tissue.

In a process similar to the one in Sec. 2.2, this grid could later be used to recover the reflectance corresponding to a measurement on a sample with arbitrary optical properties. The absorption coefficient of the sample was incorporated by applying Eq. (2), thus collapsing the grid in one dimension. This created a look-up table of the reflectance collected by each detector versus ${\mu }_{s\mathrm{ref}}^{\prime }$. Finally, reflectance values were extracted from this table by interpolating between ${\mu }_{s\mathrm{ref}}^{\prime }$ values by cubic splines. This procedure is shown graphically in the Appendix.

In an experimental situation, the optical properties of a tissue sample can be recovered by generating successive spatially resolved reflectance curves with this algorithm and using Marquardt–Levenberg curve-fitting to find the best fit to measured data. By changing the geometry used in the initial Monte Carlo simulations, this algorithm could be applied to any desired probe design.

To ensure accuracy in reflectance curves generated with this algorithm, a number of parameters were investigated and optimized. Because photons were scored in discrete bins based on the total path length, both the maximum path length scored and the width of the bins (path length resolution) were parameters of the model. Simulations were run in which photons were scored at very high resolutions of ${1\times 10}^4\mathrm{bins}/\mathrm{mm}$ and maximum path lengths of 150 mm. While these values would be impractical for use in a final algorithm due to time and file size constraints, they were used in initial investigations to allow for independent evaluation of each parameter. The optical properties used in these investigations were ${\mu }_a=\{0.01,0.5\}$ and ${\mu }_s^{\prime }=\{0.5,13\}$. The four combinations of these values covered the range of optical properties of brain tissue listed in Table 1.

The influence of the maximum path length ($x_m$) was investigated by finding the total reflectance recovered for a series of $x_m$ values and comparing this to the reflectance when the path length was not limited. The path length resolution was investigated by progressively summing the contents of adjacent bins to simulate coarser binning. Finally, since each simulation used to generate the grid was run with a different value of ${\mu }_s^{\prime }$; the resolution in ${\mu }_s^{\prime }$ was another parameter of the model (the range was previously defined as 0.5 to ${13\mathrm{mm}}^{-1}$). This resolution was progressively improved by running additional simulations and observing plots of the total reflectance versus ${\mu }_s^{\prime }$. This procedure was continued until intermediate data points could be accurately interpolated by cubic splines. The error in interpolation was estimated by consecutively removing each point from the data (except the endpoints) and interpolating for the value of the missing data point.

2.6.

Optical Property Recovery Performance

Using the Marquardt–Levenberg curve-fitting routine with absolute reflectance data from the hybrid reflectance recovery algorithm will allow recovery of ${\mu }_a$ and ${\mu }_s^{\prime }$ from measured SR-DRS data. A simulation-based testing procedure was used to characterize the precision of the optical property algorithm. The relative uncertainty (defined here as the standard deviation divided by the mean of 1000 runs) in the recovered optical properties was observed when simulated reflectance curves corrupted by random noise were input to the algorithm.

3.1.

Monte Carlo Model

Results of the comparison between our basic Monte Carlo model, the MCML code,²³ and diffusion theory are shown in Fig. 3 for selected optical properties. It can be observed that at short SDDs (Fig. 3 inset) the results are consistent with MCML, while at longer distances, the results are consistent with both MCML and diffusion theory.

Fig. 3

Validation of Monte Carlo codes. Solid lines are spatially resolved diffuse reflectance spectroscopy (SR-DRS) reflectance curves generated using the Mono Monte Carlo method. Dots indicate the results of the MCML code by Wang et al.²³ Crosses show the results of diffusion theory calculations using the extended source model described by Farrell et al.³⁷

3.2.

Reduced Probe Size

The results of the investigation described in Sec. 2.3 (reducing the probe diameter) are shown in Fig. 4, which plots the percent uncertainty in the recovered optical properties versus the maximum SDD in the simulated probe design. These results indicated that the presence of random noise in the reflectance curves had a strong effect on the precision of the optical properties recovered by the Mono Monte Carlo algorithm. The relative uncertainty in the recovered optical properties was found to increase approximately linearly with the standard deviation of the added noise and to vary approximately inversely with the SDD of the outermost detector. Examination of Fig. 4 suggests that, for the optical properties shown, an outermost distance $\ge 1.5\mathrm{mm}$ would demonstrate acceptable precision. This was consistent with other simulations using different optical properties (data not shown). It can also be observed that at short SDDs, the uncertainty in ${\mu }_a$ is significantly greater than ${\mu }_s^{\prime }$. At short SDDs, collected photons have experienced short path lengths in tissue, and therefore the influence of small changes in ${\mu }_a$ is insignificant.

Fig. 4

The percent uncertainty (standard deviation divided by the mean) in the value of the optical properties recovered from 1000 reflectance curves corrupted by normally distributed random noise with a standard deviation of 5, 10, or 15% as indicated. Optical properties of ${\mu }_a=0.31{\mathrm{mm}}^{-1}$ and ${\mu }_a=2.34{\mathrm{mm}}^{-1}$ were used to generate all initial reflectance curves.

The effect that reducing the maximum SDD had on the depth probed by reflectance measurements is demonstrated by the simulation results in Fig. 5. This shows the average depth probed by photons as they passed below a radial distance of 0.89 mm when the detector was placed at various SDDs.

Fig. 5

The influence of the source to detector distance on the depth of tissue probed by DRS. The average depth of tissue probed below a radial distance of 0.89 mm is shown for optical properties of ${\mu }_a=0.31{\mathrm{mm}}^{-1}$ and ${\mu }_a=2.34{\mathrm{mm}}^{-1}$.

Using these data alongside other practical constraints, a design with three detectors placed up to a maximum SDD of 1.67 mm was chosen for use in brain tissue. The furthest detector was placed near to the outer edge of the probe to maximize its SDD, and the diameter of the probe was limited to $<2\mathrm{mm}$ to allow easy manipulation in a small surgical cavity.

Multiple detector fibers were placed at each SDD, and their outputs were bundled together to improve the photon collection efficiency. More fibers were placed in bundles at further SDDs to help compensate for the decreased signal at these distances. The design chosen for the probe is shown in Fig. 6, and the relative signal strength collected by each detector bundle is shown in Table 2.

Fig. 6

The proposed design for the fiber optic probe. DRS detector fibers are bundled together into three groups, each at the indicated distance from the DRS source fiber (measured from the centers of each fiber). The central fiber is used for both illumination and detection in time-resolved fluorescence spectroscopy measurements.

Table 2

Relative signal strength in each of the outer two diffuse reflectance spectroscopy detector bundles of the probe shown in Fig. 6, obtained through Monte Carlo simulations. This was defined as the total reflectance (photons per unit area) collected in all fibers of the indicated bundle divided by the total reflectance in the detector bundle 0.23 mm from the source.

A 0.4-mm diameter fiber was included in the probe design for future use in a combined diffuse reflectance/time-resolved fluorescence instrument. Although photons entering this fiber were not scored during Monte Carlo simulations, its presence was considered because it altered the reflectance profile of the probe due to a difference in the index of refraction of the fiber and the surrounding probe surface.

3.3.

Boundary Effects

The consequences of ignoring the change in refractive index at the probe’s edge (boundary effects) are demonstrated in Fig. 7. These results show that although the Mono Monte Carlo algorithm was able to generate a good fit to the (simulated) true probe reflectance, this fit corresponds to incorrect optical properties. This analysis was performed over a range of optical properties, and results are listed in Table 3. The discrepancies between the target and recovered absorption coefficients in Table 3 show the strong influence of boundary effects on this probe design. This motivated the development of a new optical property recovery algorithm, one that considers the true geometry of the probe surface. This algorithm was described in Sec. 2.5.

Fig. 7

The effects of boundary conditions on optical property recovery with the Mono Monte Carlo algorithm. The solid line is the reflectance curve predicted by the Mono Monte Carlo model for optical properties of ${\mu }_a=0.02{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.99{\mathrm{mm}}^{-1}$. Points marked with an “$x$” show the reflectance predicted by the true probe forward Monte Carlo model for the same optical properties. The dotted line is the fit generated with the Mono Monte Carlo algorithm, which corresponds to optical properties of ${\mu }_a=2.8\times {10}^{-6}{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}$.

Table 3

The effects of boundary conditions on the Mono Monte Carlo algorithm’s ability to recover optical properties.

3.4.

Parameter Optimization: Hybrid Reflectance Recovery Algorithm

Selected results from the investigation of the maximum path length parameter are shown in Fig. 8. It was found that for path lengths $<150\mathrm{mm}$, the reflectance in each detector reached $>99\%$ of its maximum value, the total when the path length was not limited. This was true for all investigated optical properties (data not shown), and therefore 150 mm was chosen as the maximum path length to record in all simulations.

Fig. 8

Investigation of the maximum path length parameter for simulations with optical properties of (a) ${\mu }_a=0.01{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=13{\mathrm{mm}}^{-1}$ and (b) ${\mu }_a=0.5{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=0.5{\mathrm{mm}}^{-1}$. Each curve corresponds to a simulated detector bundle with the indicated source-to-detector distance. $N_p(x)$ is the number of photons collected that traveled a total path length equal to $x$. $N_{p\mathrm{ref}}$ is the total number of photons collected for each source–detector distance (SDD) when the path length was not limited.

Investigation of the path length resolution parameter showed that when it was decreased to 13 bins/mm, the simulated reflectance remained within 1% of the reflectance at ${1\times 10}^4\mathrm{bins}/\mathrm{mm}$. This was true across the range of investigated optical properties, and therefore a path length resolution of $13\mathrm{bins}/\mathrm{mm}$ was selected for $x_t$. Since the maximum path length recorded was 150 mm, 1950 path length bins were used in total.

The resolution in ${\mu }_s^{\prime }$ is shown in Fig. 9 for ${\mu }_a=0.02$ and ${0.5\mathrm{mm}}^{-1}$. The error in each point was found to be $<2\%$ when the resolution in ${\mu }_s^{\prime }$ was set to ${0.5\mathrm{mm}}^{-1}$.

Fig. 9

The reflectance recovered using the method described in Sec. 2.5 for each detector bundle in the probe layout of Fig. 6. Results are shown over a range of desired ${\mu }_s^{\prime }$ values, with ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ (a) and ${0.5\mathrm{mm}}^{-1}$. (b) Data points indicate the ${\mu }_{s\mathrm{ref}}^{\prime }$ values used to generate the grid and the lines are the results of interpolation by cubic splines. The error in each data point was $<2\%$.

3.5.

Optical Property Recovery Performance

Results characterizing the precision of the optical property recovery algorithm are shown in Fig. 10. All reflectance data used in the fitting procedure corresponded to the probe design of Fig. 6. As expected, the highest uncertainty in the results occurs when recovering the absorption coefficient from data corresponding to low ${\mu }_a$ values. This is a consequence of using a probe with short SDD values; the algorithm is less sensitive to changes in ${\mu }_a$ than ${\mu }_s$. This is especially true for tissues with very low absorption because of the shape of the corresponding reflectance curve.

Fig. 10

Performance of the optical property recovery algorithm when presented with simulated reflectance curves corrupted by random noise with a standard deviation of 5%. The values on the $x$ and $y$ axes are the target optical properties, while the contour lines represent the relative uncertainty in the recovered optical properties, defined as the standard deviation of the results divided by the mean over 1000 runs.