Probing depth in diffuse optical spectroscopy and structured illumination imaging : a Monte Carlo study

Diffuse optical spectroscopy (DOS) and its modification employing structured illumination are widely used in monitoring biotissue oxygenation. In such measurements it is important to know the probing volume for definite source-detector configuration; however, it cannot be measured directly. Monte Carlo simulations allow to trace the probing depth of each individual photon contributing to the signal, which provides a numerical solution for this problem. In this study we investigate distributions of photons over maximal depth reached in turbid media (probing depth) with optical parameters typical for cutaneous tissues at the wavelength of 600 nm. Different configurations of probing illumination are considered, such as collimated point source, one-dimension sinusoidal and rectangular patterns. For collimated point source and zero source-detector separation the number of collected photons monotonously decreases with the probing depth while a pronounced maximum in the distribution is manifested with the increase of source-detector separation. The position of this maximum shifts to higher depths with the decrease of μa. For one-direction sinusoidal and rectangular illumination patterns it is shown that when the photons are collected near the center of a bright stripe, the peak of the distribution remains close to the surface. When the photons are collected near the center of a dark stripe the peak shifts towards higher depths with the decrease in spatial duty cycle and spatial frequency of the illumination pattern. Employment of rectangular illumination pattern seems more efficient for DOS applications due to wider abilities for controlling probing depth. © 2017 Journal of Biomedical Photonics & Engineering.


Introduction
Diffuse optical spectroscopy (DOS) is a modern biomedical diagnostics technique based on registration of scattered spectra for a single or several pairs of source and detectors.Reconstruction of absorption spectra from these measurements allows to evaluate chromophore content within the studied tissue such as oxy-and deoxyhemoglobin, lipid and water as well as to retrieve hemoglobin oxygen saturation index [1,2].DOS is a perspective tool for a broad range of medical diagnostics tasks including breast cancer diagnostics [4][5][6], human brain functional diagnostics [7], monitoring of tumor in laboratory animals in experimental oncology [8].
Numerical models of DOS signals are usually based on diffusion approximation of transport theory [9].It is well known that in reflectance mode typical photon trajectories from source to detector form so-called banana-shaped zones which characterize the most sensitive volume to locate inclusions in the medium ("measurement volume").Maximum depth of the most probable photon trajectory can be considered as the probing depth of the particular source-detector configuration.It is clear that the probing depth and the measurement volume depend on tissue optical properties and measurement conditions.Although the probing depth and the measurement volume are the characteristics of great importance for clinical applications, they hardly can be determined experimentally.To our knowledge, analytical evaluation of probing depth has been performed only for the simplest source-detector configuration and is unavailable in the case of complex illumination modalities.In this respect Monte Carlo technique is more appropriate for estimation of the measurement volume and the probing depth.Based on modeling of random photon walk this approach provides information about exact trajectory of an individual photon in the studied medium and allows to analyze the total set of photons trajectories contributing to DOS signal for particular source-detector separation.For instance, Monte Carlo simulations of photon trajectories from source to detector are carried out in studies [10] for semi-infinite medium with a plane boundary and [11] for complex multilayered medium model and demonstrate banana-shaped zone of typical photon trajectories.
The modification of DOS employing structured illumination can be efficiently used for mapping of biotissue absorption coefficient [2,12,13].Variation of illumination pattern, in particular, changing its spatial frequency, potentially allows to control measurement volume.At the same time, in this modality knowing probing depth is also of high importance.
In this study we report on results of Monte Carlo study of probing depth for DOS and its modality with structured illumination.

Simulated measurement configurations
Monte Carlo modeling was carried out for a uniform slab sized 10x10 mm in transversal directions and 5 mm in-depth.DOS simulation was performed for two arrangements of the medium illumination: illumination by a collimated point source and illumination by a planar source with one-dimension striped pattern.In the first case the point source is located at the medium boundary and is directed normal to it.Detection of scattered light is performed by a wide-angle detector sized 0.05x0.05mm 2 located at the medium boundary at the distance in the range 0-10 mm from the source which covers the practical range for DOS applications.
Probing with structured illumination is simulated for two types of the incident light spatial modulation: i) one-dimension sinusoidal pattern with spatial frequency varying in the range of 0.3-2.5 mm -1 and ii) one-dimension rectangular pattern with fixed width of stripes (0.1 mm) and duty cycles ranging between 0.03 and 0.25.Detection of photons was performed by the detector described above.

Monte Carlo technique
A custom developed MATLAB-based code adopted from our previous Monte Carlo solutions [14] was employed in the study.Heyney-Greenstein scattering phase function of the medium was applied.The code is optimized for MATLAB environment.Since MATLAB is efficient in arrays processing, the MATLAB realization of the Monte Carlo algorithm is based on the simultaneous processing of a definite number of photons chosen as 10 6 to ensure optimal PC performance and the use of RAM.The photon parameters are stored as arrays and are transformed in accordance with MC algorithm rules, i.e. 10 6 photons steps are processed at each iteration of the algorithm.At the stage of simulation initialization the initial parameters of each photon are generated (such as photon position, direction and weight) and stored in corresponding arrays.The photon trajectories are considered to be completed when the photon leaves the sample or the critically low weight is reached.This conditions are checked at each algorithm iteration and the parameters of photons satisfying these conditions are stored.The completed trajectories in the array are replaced with the new ones initialized at the source position.The iterative process is performed until the number of completed photon trajectories reaches a predefined number equal to 10 7 in this study.
To model the effect of structured illumination the photon distribution over probing depths obtained for a point source was convolved with the considered illumination pattern resulting in the distribution corresponding to this pattern.In this calculations it was assumed that the overall probing power delivered to the object is the same for all the considered patterns.

Optical properties of studied tissue
Optical properties of the medium used in our MC simulations correspond to human skin at the wavelength of 600 nm [1,15] (see Table 1).In order to study the effect of optical properties on the distribution of photons over probing depth each of parameters was varied in the range of ±60% with respect to the basic values.The entire range of parameters is summarized in Table I.The refractive index of the medium amounted 1.4.
Table 1 Typical optical properties of human skin at the wavelength of 600 nm (basic value) and their variations employed in simulations (basic value -60% / basic value / basic value + 60%).µa, mm -1  µs , mm -1 g 0.088 / 0.22 / 0.35 8 / 20 / 32 0.8 expected, stronger attenuation of intensity with the distance from the source corresponds to the case of higher absorption coefficient.Typical distributions of photon trajectories contributing to reflectance signal over maximal depth reached in the medium (probing depth) are demonstrated in Fig. 2 as color-encoded maps where abscissa corresponds to probing depth, ordinate corresponds to source-detector separation and the logarithm of the respective fraction of photons is encoded with color.The results are demonstrated for the set of basic optical properties and for cases when the values of scattering and absorption coefficient are varied by ±60%.The maps demonstrate that at small source-detector separation (R) the contribution of photons from the smaller depth (up to 0.2 mm) prevails while when the source-detector separation increases the most probable probing depth also increases reaching about 1 mm for R = 5 mm.

Illumination with a collimated point source
With the increase of scattering coefficient (from Fig. 2a to Fig. 2c) the distribution shrinks to the region of smaller probing depths and smaller source-detector separations.Higher scattering coefficient with all other parameters preserved results in longer photon pathways required to reach particular depth (due to multiple scattering) and, hence, stronger absorption at lower depth.Similar tendency can be observed for increase of absorption coefficient (from Fig. 2d to Fig. 2f): higher absorption results in shrinkage to lower probing depths and smaller source-detector separations.Similar response of the studied distribution to the variation of scattering and absorption coefficient can be explained by prevailing of diffusive regime of light propagation in the considered media.In diffusive regime the light propagation is governed by the effective attenuation coefficient µ  = √µ  (µ  + µ  (1 − )) ; given µ  ≪ µ  (1 − ) the effect of variation of µs and µa becomes equivalent which is confirmed by the obtained results.
In order to study quantitatively the effect of optical properties on the photon distribution over probing depth, the cross-sections of the corresponding distributions were plotted separately for different source-detector separations.Fig. 3 demonstrates the distributions of photon fractions over probing depth for several source-detector separations and for varying scattering coefficient.Distributions for R = 0 demonstrate monotonous decrease of the dependence on Z while for R = 1 mm and R = 2 mm the dependence demonstrates pronounced maximum.This maximum can be associated with the centers of the banana patterns [3] indicating DOS measurement volume.With the increase of R this maximum is shifted to the region of larger probing depths with simultaneous decrease of the fraction of contributing photons.For R = 0 (Fig. 3a) the fraction of photons with small probing depth (below 0.8 mm) attenuates with µs decrease while the fraction of photons with higher probing depth (above 0.8 mm) increases.This fact is explained by the overall weaker attenuation of probing radiation with depth for smaller values of µs.For R = 1 mm the fraction of photons in the vicinity of the distribution peak is the highest for the basic parameters compared to higher and lower µs values (Fig. 3b).This fact is explained by two controversial effects: on the one hand, higher µs.ensures higher scattering and, hence, contribution to backscattering signal; on the other, longer photons travel paths in medium results in higher absorption and, hence, lower photon fractions.As for the most probable probing depth, it is about 0.35 mm and.For R = 2 mm the effect of intensity attenuation due to absorption at longer pathlengths start  to prevail and the amount of collected photons decreases with the increase in µs (Fig. 3c).The most probable probing depth for this R is about 0.6 mm.
On the contrary to variation in µs, variation in µa results in monotonous dependencies: the fraction of collected photons rises with the decrease in µa for all three considered R values (Fig. 4a).For R = 1 and 2 mm the most probable probing depth shifts to higher values with the decrease in µa (Figs.4b and 4c).Overall, the most probable probing depth demonstrates stronger dependence on µa compared to µs.

Illumination with spatially modulated light
One of the modifications of DOS employs illumination with structured patterns of light in order to control measurement volume.Different parameters of the pattern provide different distribution of probing radiation within the object.In this study we consider one-dimension sinusoidal and rectangular spatial patterns which are usually employed in measurements DOS systems.Sinusoidal patterns are more common [12,13] since they allow for analytical solution of diffusion equation with further ability to reconstruct a map of absorption coefficient [12].It is assumed that probing volume significantly depends on the spatial frequency of the sinusoidal pattern employed.
Typical distributions of the detected photons over probing depth (abscissa) and the position of the detector (ordinate) for sinusoidal patterns with different spatial frequency are shown in Fig. 5. From this figure one can observe that for smaller probing depths (below 1 mm) the distribution dependence on R resembles the illumination pattern while for higher depths the  Similar trend was obtained for the case of rectangular patterns with the decrease in spatial duty cycle (Fig. 6), however, in this case the change in the structure of the probing depth map is more significant.This indicates higher potential of employing rectangular patterns for controlling measurement volume compared to sinusoidal patterns.Cross-sections of photon distributions over probing depth were plotted for transversal coordinates of the detector corresponding to centers of bright and dark stripes for both types of illumination pattern and are shown in Figs.7 and 8, respectively.For sinusoidal patterns when the detector is placed at the center of bright stripe (Fig. 7a) the distribution maximum defining the most probable probing depth is located at the depth of about Z = 0.1 mm for spatial frequency f = 0.3 mm -1 and shifts to the surface of biotissue sample with the increase of frequency.When the detector is placed at the center of dark stripe (Fig. 7b) the maximum of distribution is located at Z = 0.25 mm for f = 0.3 mm -1 and, similarly, it shifts to the surface of biotissue sample with the frequency value increase.Note that the fraction of photons detected at the center of bright stripe is highest for f = 0.3 mm -1 and decreases with the frequency increase.On the contrary the fraction of photons detected at the center of the dark stripe increases with the frequency increase.This is explained by the fact that in the case of low spatial frequency local probing intensity is higher in the vicinity of the bright stripe center and is lower in the vicinity of the dark stripe center as compared to the case of higher frequency.
Similar results were obtained for rectangular illumination patterns (Fig. 8).The fraction of photons detected at the center of bright stripe is highest for the spatial duty cycle D = 0.03 and decreases with the duty cycle increase.On the contrary, at the center of dark stripe the fraction of collected photons increases with the duty cycle increase.As in case of the sinusoidal patterns, this fact originates from higher local probing intensity in the vicinity of bright stripe center compared to the dark stripe.In the case of rectangular pattern the difference is stronger which yields the shift of the distribution peak to the sample surface (Z = 0 mm) (Fig. 8a).When the detector is placed in the center of dark stripe, the position of the distribution peak demonstrates strong dependence on the duty cycle: the peak shifts from Z = 0.6 mm to Z = 0.1 mm when the duty cycle alters from 0.03 to 0.25.

Conclusions
In this work the effect of optical properties and probing pattern configuration on the probing depth of diffuse optical spectroscopy is analyzed.Distributions of photons over maximal depth reached in the medium associated with probing depth are derived from numerical Monte Carlo simulations.For collimated point source it is shown that the fraction of collected photons demonstrates monotonous decrease with probing depth for source-detector separation equal to zero while with increase in source-detector separation a maximum at the definite depth is manifested.The position of this maximum shifts to higher probing depths with decrease of μa.
For the problem of structured illuminations the maps of photons distribution over probing depth are derived by the convolution of Monte Carlo results for a point source with the given illumination spatial patterns.Sinusoidal and rectangular patterns are considered.It is shown that the decrease in duty cycle and spatial frequency results in significant redistribution of the maps for rectangular and sinusoidal patterns, respectively.In particular, the distribution resembles illumination patterns for probing depths below 1 mm while for higher probing depths the distribution remains almost uniform and unchanged for different pattern parameters.Distribution of photons over probing depth demonstrates opposite trends when the detector is located in the center of bright versus dark stripe of the illumination pattern.For the center of bright stripe the peak of the discussed distribution remains close to the surface while for the center of dark stripe it shifts towards higher depths with decrease in spatial frequency and duty cycle for sinusoidal and rectangular patterns, respectively.Stronger dependence of this distribution on the duty cycle for rectangular patterns indicates high potential of this type of illumination pattern for controlling probing depth in DOS applications.

Fig. 1
Fig.1demonstrates typical simulated two-dimension maps of intensity scattered within tissue from a collimated point source for a definite scattering coefficient and different absorption coefficients.As

Fig. 5
Fig. 5 Distribution maps of photons over probing depth for sinusoidal illumination pattern with different spatial frequencies f: (a) f = 2.5 mm -1 ; (b) f = 1 mm -1 ; (c) f = 0.3 mm -1 .The profile of the illumination pattern is shown at the left.Color encoding corresponds to the logarithm of the respective fraction of collected photons.

Fig. 6
Fig. 6 Distribution maps of photons over probing depth for rectangular patterns with fixed width of stripes (0.1 mm) and various spatial duty cycles D: (a) D = 0.25; (b) D = 0.1; (c) D = 0.03.The profile of the illumination pattern is shown at the left.Color encoding corresponds to the logarithm of the respective fraction of collected photons.

Fig. 7
Fig. 7 Distribution of photons over probing depth for sinusoidal illumination patterns with different spatial frequency f for the detector position (a) at the center of bright stripe and (b) at the center of dark stripe.

Fig. 8
Fig. 8 Distribution of photons over probing depth for rectangular illumination patterns with different duty cycles D for the detector position (a) at the center of bright stripe and (b) at the center of dark stripe.distributionremains almost uniform and unchanged for different spatial frequencies.The major changes are manifested when the detector is located at the minima of illumination patterns where the region of most probable probing depths shifts to higher values.Similar trend was obtained for the case of rectangular patterns with the decrease in spatial duty cycle (Fig.6), however, in this case the change in the structure of the probing depth map is more significant.This indicates higher potential of employing rectangular patterns for controlling measurement volume compared to sinusoidal patterns.Cross-sections of photon distributions over probing depth were plotted for transversal coordinates of the detector corresponding to centers of bright and dark stripes for both types of illumination pattern and are shown in Figs.7 and 8, respectively.For sinusoidal patterns when the detector is placed at the center of bright stripe (Fig.7a) the distribution maximum defining the most probable probing depth is located at the depth of about Z = 0.1 mm for spatial frequency f = 0.3 mm -1 and shifts to the surface of biotissue sample with the increase of frequency.When the detector is placed at the center of dark stripe (Fig.7b) the maximum of distribution is located at Z = 0.25 mm for f = 0.3 mm -1 and, similarly, it shifts to the surface of biotissue sample with the frequency value increase.Note that the fraction of photons detected at the center of bright stripe is highest for f = 0.3 mm -1 and decreases with the frequency