1.

Introduction

The diffuse photon density wave (DPDW) methodology, developed in the mid-nineties,^1–5 consists of probing tissue with frequency modulated near-infrared (IR) light and continues to be of great interest due to its biomedical applications. The ability to determine the depth and degree of cutaneous and subcutaneous tissue damage is critical for medical applications such as diagnosis and classification of pressure ulcers, deep tissue injury, and burns. In practice, it is difficult to differentiate the superficial stage of pressure ulcers from deep tissue injury that originates at the bone-soft tissue interface, especially in patients with darker skin tones, because both conditions manifest themselves as intact skin with redness or discoloration.⁶ When diagnosing burn injury, the depth of tissue damage determines the course of treatment, yet clinical assessment of burn depth is currently based on qualitative evaluation of the tissue surface appearance.⁷

Previously, we have demonstrated the ability to determine the degree of tissue damage at depths of 3 to 6 mm in acute and diabetic wounds.^8,9 We have reported preliminary results¹⁰ obtained within the frequency-domain DPDW methodology for noninvasive measurements of reduced scattering and absorption coefficients at the depths of several millimeters. A multifrequency (50 to 400 MHz) noncontact DPDW system with one light source and one detector was constructed so that incident light is focused onto the tissue surface. The ability to digitally control source–detector separation permits obtaining the precise intensity and phase data for a virtually unlimited number of source–detector separations, including small separations, beyond the verified applicability of the diffusion approximation of the radiative transfer equation generally used for interpretation of measurements.

Presently, we apply a novel Monte Carlo procedure reported previously¹¹ and based on the stochastic simulation of iterative terms of the Bethe–Salpeter equation in scattering orders for calculating DPDW signals.

Monte Carlo modeling has been widely used to simulate photon migration in tissue and tissue phantoms,^12–16 and is most often performed using the Monte Carlo modeling of light (MCML) method described in Ref. 17. This method requires quite a large sampling size due to the fact that a very small share of the incident photons is detected. In the approach presented here, every photon contributes to the detected signal, thus greatly diminishing the number of modeled photons and reducing the run-time by order, enabling the simulations to be performed using commonly available personal computers.

Let a trajectory be simulated using the conventional MCML algorithm.¹⁷ Contributing to the detected signal, the photon packet is to escape the tissue. Thus, in the conventional approach, including implementations that use parallel algorithms,^18,19 only those trajectories containing escaped photons contribute to the signal. We use the same procedure for the trajectory evolution. However, within our approach, a trajectory containing $n$ vertexes produces $n$ inputs into the signal compared with the sole input in the conventional approach. Thus, we obtain a greater efficiency as compared with the conventional method, since in our approach, we do not wait until a photon crosses the boundary. It permits to reduce the sampling volume and computation time by at least an order, while obtaining plots of the same quality. We generalized our scheme, alternatively calculating a DPDW signal as a number of photons escaping the medium, which is the main feature of a conventional approach. We have shown that the developed algorithm exhibits a signal-to-noise ratio much higher than the algorithm based on the count of escaped photons for the same sampling volume.

Here, we present results of the Monte Carlo calculations, based upon the Bethe–Salpeter equation for the scattering intensity, and compare them with the measurement data for intralipid aqueous solutions, as well with the diffusion approximation. We perform measurements of phase and amplitude of DPDW for intralipid solutions of different concentrations. Studying parameters of the DPDW as functions of the source–detector distance, we obtain an excellent quantitative agreement between experimental data, and numerical simulations, accounting for the Fresnel reflection at the boundary in form of the known extrapolated boundary conditions,²⁰ and the diffusion approximation for a half-space geometry.

The measured DPDW amplitude and phase delay are the frequency-domain parameters. There are two methods for their determination, as is discussed in Ref. 16: the Fourier transformation of the time histogram and the short-cut method of obtaining the frequency-domain parameters,²¹ based on a solution of the radiative transfer equation (RTE) in the frequency domain. Our approach also presents a direct method of calculation in the frequency domain based on the field-theoretical approach, applied to the actual experimental setup.²²

Further, we model intralipid solutions as suspensions of monodisperse hard spheres with a given diameter, and we successfully fit experimental data obtained at two wavelengths, $\lambda =685$ and 830 nm, with the same sphere diameter value. This has the potential to provide additional information on the size of inhomogeneities in tissue when performing measurements with multiple wave lengths. Our simulated numerical results for multiple concentrations, modulation frequencies, and wave lengths, supplemented with calculation of scattering coefficients within the Born approximation, are in agreement with the experimental measurements, and do not require the use of a fitting parameter.

The paper is organized as follows. In Sec. 2, the system used for the frequency-domain DPDW measurements is described. Section 3 contains the formalism of the radiative transfer theory in terms of the Bethe–Salpeter equation; in Sec. 4, the MC algorithm based upon the simulation of the Bethe–Salpeter equation is outlined. Section 5 contains the results of simulations and in Sec. 6, we evaluate the spatial dimensions of scattering inhomogeneities compared to the measurements obtained for different wavelengths of radiation. Section 7 contains our conclusions.

2.

DPDW Instrumentation

Our device is capable of measuring a virtually unlimited number of source–detector separations by using a digital actuator ($\text{resolution}=0.025\mathrm{mm}$) to move the incident light across the surface of the probed medium and a fixed detector fiber (2.5-mm ferrule diameter). The ability to select the frequency of modulation up to 400 MHz allows for better determination of a phase shift at small source–detector separations, where the shift of the phase becomes difficult to detect using smaller frequencies. The ability to digitally control the source–detector separation enables precise selection of the volume and depth of tissue that will be characterized. We can analyze the experimental data to extract information from more precise depth ranges within the overall penetration range of 2 to 10 mm. As mentioned above, the determination of optical properties at very small source–detector separations is not possible using a diffusion-based model. Monte Carlo simulations would enable the extraction of tissue optical properties at depths where the diffusion approximation is no longer valid. The setup schematic is shown in Fig. 1.

Fig. 1

Schematic of the probe assembly. The digitally controlled actuator allows multiple source–detector separations for better depth of measurement resolution.

Conventional diffusive near-infrared spectroscopy (DNIRS) utilizes relatively large-sized probes to secure multiple optical fibers at known distances. One assumption used in conventional DNIRS when solving the diffusion equation is that the measurement geometry is semi-infinite, i.e., there are no scattering or reflection events outside of the medium. This assumption is violated when light is reflected from the surface of the probe back into the medium, allowing for additional light to reach the detector and making measurements of optical properties less accurate. Our modified instrument utilizes lenses for light delivery and a single detector fiber for the collection of backscattered light, creating a practically ideal semi-infinite geometry that is less susceptible to surface reflections. Variations in the amount of pressure applied during probe placement are also minimized with our set up. Lastly, correction coefficients are needed during conventional DNIRS measurements to compensate for discrepancies in optical transmission of different fibers as well as differences in phase shift within different light source channels. Our single source, single detector interface eliminates the need for correction coefficients to adjust measured parameters, reducing system variability and minimizing error.

Let the light modulated with frequency $\omega$ be incident upon a medium; the intensity of scattered radiation turns out to be the sum of DC and AC terms. Let the AC component which is dependent on time $t$ in the form $I_{\mathrm{AC}}\sim \cos (\omega t+\theta )$, where $\theta$ is a displacement of the phase due to the optical path traveled by light in a medium, be mixed with the reference signal, $I_{\mathrm{ref}}\sim \cos (\omega t)$. Thus, assuming that the high frequency part is cut off, there appear two frequency-domain DPDW signals:

Figure 2 presents the measurement scheme of signals in the near-IR range in the frequency domain. In our setup, the optical radiation from a laser diode, which can be modulated with a frequency from 50 to 400 MHz, is focused onto the surface of a tissue or a tissue phantom. The light scattered backward is multiplied in a demodulator by a reference radio signal (LO signal). As a result, after passage through a low frequency filter, two DC signals are measured: the $I$-signal, produced in $D_1$ demodulator mixing the useful signal with the reference one with the phase coinciding with the incident modulation phase, and the $Q$-signal, for which the phase of the reference signal is shifted by a modulation wave quarter. The signals measured give the amplitude $A=\sqrt{I^2+Q^2}$ and phase $\theta =\arctan (Q/I)$ of the DPDW.

Fig. 2

Schematic setup for measuring signals of photon density waves in the frequency domain. Amplitude-modulated light from a laser diode is delivered to the system under test through a lens that can be moved using a digitally controlled actuator. Backscattered light is delivered to an APD using an optical fiber. The amplified APD signal enters an $IQ$ demodulator (DM) where it is compared to the reference signal. The $I$ and $Q$ signals are send to a data acquisition board and personal computer.

Our device consists of two separate parts: (1) a “module assembly” containing the opto-electronic components and (2) an optical probe assembly containing a lens-coupled source fiber mounted on a digital actuator to carry light from the module assembly to the optical system being assessed, and a fixed detector fiber to carry light from the medium to the module assembly. The module assembly is enclosed within a standard 12-in. rack with six shielded NIM box modules. The first module contains the power supply, which provides four specific DC voltages (3.3, 5, 12, and 15 V) to adequately power the different components of our device. The second module contains a locking programmable sine wave generator (Novatech Instruments, LPO400A) that can produce sinusoidal radio frequency (RF) signals ranging from 200 kHz to 400 MHz with 1-Hz steps. During measurements, frequencies from 50 to 400 MHz were used. As the modulation frequency changes the impedance of our electrical components varies, causing inconsistent modulation depths (VPk-Pk/VRMS) across different frequencies. To compensate for this difference, a digitally controlled variable gain amplifier (Analog Devices, AD8375) is used to maintain a constant modulation depth of laser radiation.

The RF signal coming from the variable amplifier modulates the intensity of the semiconductor laser diodes (685 or 830 nm) inside the third module. Custom laser drivers are used to mix the DC bias of our laser diodes with the AC signal from our RF generator in order to obtain intensity modulated light; using an integrated circuit (SHARP Microelectronics, IR3CO7) and bias-tee (Mini-Circuits, TCBT-2R5G+). In the fourth module, a $4\times 1$ MEMs optical switch (DiCon Fiberoptics, Inc.) allows us to cycle between the two laser diodes and an offset (no light) position for each individual distance/frequency.

The RF signal coming from the variable amplifier modulates the intensity of two semiconductor laser diodes (685 or 830 nm) inside the third module. Custom laser drivers are used to mix the DC bias of our laser diodes with the AC signal from our RF generator in order to obtain intensity modulated light, using an integrated circuit (SHARP Microelectronics, IR3CO7) and bias-tee (Mini-Circuits, TCBT-2R5G+). In the fourth module, a $4\times 1$ MEMs optical switch (DiCon Fiberoptics, Inc.) allows us to cycle between the two laser diodes and an offset (no light) position for each individual distance/frequency. Variable optical attenuators (Oz Optics, Ltd.) are connected between each laser diode output and the corresponding optical switch input to provide control over the intensity of light delivered to the tissue. A $62.5/125\text{micron}$ multimode optical fiber is used to transport light from the optical switch to a lens assembly which first collimates and then focuses the laser light onto the medium. The lens assembly is attached to a digital actuator, which allows the position of the focused light to be digitally controlled by the user with source–detector separations ranging from 2 to 20 mm. Light that is backscattered to the surface is collected by a detector fiber (1-mm core) directly touching the medium and is transported to an avalanche photodiode (APD) module (Hamamatsu C5658) inside the fifth module. Figure 2 shows the optical probe assembly and set up.

The electrical signal from the APD is passed through two fixed amplifiers (Mini-Circuits, ZFL-500LN and Mini-Circuits, ZFL-500HLN). Then the signal is fed to an in-phase/quadrature ($I/Q$) demodulator (MERRIMAC IQM-9B-500), which compares the detected backscattered signal to a reference signal from the RF generator. The outputs of the $I/Q$ demodulator are the cosine and sine low frequency components of amplitude and phase shift relative to the reference signal. In the sixth and final module, a data acquisition board (National Instruments, USB-6251) is used to digitize the output from the $IQ$ demodulator and to control all digital components (RF generator, variable amplifier, optical switch, and actuator driver). A personal computer equipped with LabVIEW and MATLAB software provides¹⁰ the control interface for the overall system and analysis/fitting of the experimental data to calculate optical properties and hemoglobin concentrations within tissue.

3.

Bethe–Salpeter Equation

The Bethe–Salpeter equation for the product of a pair of spectral components of complex-conjugated scalar fields, responsible for the AC part of scattering, can be presented as

The optical theorem relates the reciprocal of the scattering length, or the scattering coefficient ${\mu }_s$, and the scattering cross section; for the scalar field, it takes the form

Defining the phase function as the normalized cross section of scattering

For $p(t)$, dependent on the angle of scattering, with $t=\cos \theta$, we use the Henyeye–Greenstein phase function.

Let $z$ be the Cartesian coordinate, $\mathbf{r}=({\mathbf{r}}_{\perp },z)$, normal to the boundary of a semi-infinite medium, $z>0$. For the normal incidence of a pencil-like beam in the vicinity of ${\mathbf{r}}_1=0$ upon the boundary $z=0$ and pencil-like normal backscattering in a small area $\sigma$ at source–detector separation $\rho =|{\mathbf{r}}_{1\perp }-{\mathbf{r}}_{2\perp }|$, the DPDW intensity can be defined as

4.

Monte Carlo Procedure

Within the stochastic method, the scattering intensity is presented as a statistical average over a sampling of $N_{\mathrm{ph}}$ incident photons

The weight $w_n^{(i)}$ arises from a random value of the integrand of the multifold spatial integral which is the $n$’th-order iteration of the Bethe–Salpeter equation [Eq. (10)]. Calculating it one simulates a stochastic sequence, or trajectory, of scattering points ${\mathbf{r}}_1,\dots {\mathbf{r}}_n$. The weight $w_n^{(i)}$ also accounts for the boundedness of a medium and guarantees the observation point to be in the vicinity of the detection.

Note that studying the photon migration in the time-domain technique, one must calculate the Fourier transform of a time histogram^12,16 for determining the frequency-domain DPDW signal. Equation (12), used previously in our work,²² yields the DPDW signals immediately in the frequency domain.

Using the well-known algorithm of radiative transfer simulation, the relative distance $r^{\prime }=|{\mathbf{r}}_j-{\mathbf{r}}_{j-1}|$ is changed to random variable $\xi =\exp (-\mu r^{\prime })$ and cosine $t=\cos \theta$ of the scattering angle is changed to $\chi =2\pi {\int }_{-1}^{t}p(t^{\prime })\mathrm{d}{t}^{\prime }$; thus, one transforms the three-dimensional spatial integral over relative coordinate ${\mathbf{r}}^{\prime }={\mathbf{r}}_j-{\mathbf{r}}_{j-1}$ in its own coordinate frame for vector ${\mathbf{r}}_{j-1}$ as

Accounting for the boundedness of the medium, one returns a photon into the medium if it reaches the boundary due to the reflection law multiplying the weight factor by the Fresnel reflection coefficient.²³ Thus, the weight $w_n^{(i)}$ is the product of the reflection coefficients and factor ${({\mu }_s/\mu )}^n$; it turns out to be unit for ${\mu }_a=0$ and no reflection.

The calculation time depends strongly on the number of terms $n_{\mathrm{sc}}$ in sum [Eq. (12)], as well as on the sampling volume $N_{\mathrm{ph}}$; in Ref. 14, the number of scattering events was traced up to $n_{\mathrm{sc}}={10}^5$. At fixed transport length $l_{\mathrm{tr}}=1/{\mu }_s^{\prime }$, the number of scattering events to be taken into account increases for smaller scattering length $l_s=l_{\mathrm{tr}}(1-g)$, i.e., for larger scattering anisotropy $g$. Evaluating the argument of the oscillating factor in Eq. (12) as $c^{-1}\omega R_n^{(i)}\sim c^{-1}\omega n{l}_s\sim 1$ we obtain an estimate $n\approx {\omega }^{-1}{\mu }_s^{\prime }c{(1-g)}^{-1}$.

To compare the present approach with the conventional approach, we generalize the code and alternatively calculate the signal of the scattering intensity, as the number of photons escaping the medium within the lines of the MCML method.¹⁷ Namely, we change Eq. (12) to

In Fig. 3, we present the simulation flow chart, describing the two ways of signal accumulation, either with Eq. (12), or with Eq. (14), for steady-state radiation, $\omega =0$. For brevity, we neglect the inner reflection and adsorption. The calculation yields two two-dimensional arrays of the $n$’th scattering order inputs for a series of discrete values of the source–detector distances.

Fig. 3

Flow chart for Monte Carlo simulation for a normal backscattering from a half-space; BS denotes the calculation based on the present approach and MCML-like denotes the calculation based on the count of escaped photons.

Note that the scattering intensity calculated with Eq. (12) can be interpreted as an average of exponentials $\exp [-\mu z_n^{(i)}]$, which describe the extinction of the radiation returning from the medium to the boundary after the $n$’th act of scattering, and with Eq. (14) we calculate the mean number of photons escaping the medium also after $n$ scattering events.

To evaluate statistical noise inherent in the two approaches, we calculated the coefficient of variation $c_v$, or the relative standard deviation, for distributions of two random variables, either $J_0^{(i)}(\rho )$ or $J_{\mathrm{MCML}}^{(i)}(\rho )$. We found that the $c_v$ value for variables [Eq. (14)] significantly exceeds the value for variables [Eq. (12)]. Thus, we conclude that the commonly used approach based upon the calculation of the number of escaping photons requires much larger sampling than that for the approach based on the Bethe–Salpeter equation.

For comparison, we calculated the DPDW parameters using alternative definitions [Eqs. (12) and (14)]. We found that plots obtained with Eq. (14) for sampling $N_{\mathrm{ph}}={10}^6$ exhibit the same level of noise as those obtained within the present approach, Eq. (12), for sampling $N_{\mathrm{ph}}={10}^5$. Thus, we conclude that for the same noise, the MCML-like approach will require about a 10 times longer run-time than the considered one.

5.

Simulation Results

We performed measurements of the phase and magnitude of DPDW in dependence on the source–detector separation for intralipid–water solutions, varying the concentration from 0.5 to 2 vol. %, for two modulation frequencies, ${(2\pi )}^{-1}\omega =100$ and 200 MHz, and for two wavelengths of incident radiation, $\lambda =685$ and 830 nm.

Simulating, we varied the water absorption coefficient from ${\mu }_a=0.0045{\mathrm{cm}}^{-1}$ to ${\mu }_a=0.01{\mathrm{cm}}^{-1}$ for $\lambda =685\mathrm{nm}$ with no noticeable effect, and took ${\mu }_a=0.04{\mathrm{cm}}^{-1}$ for $\lambda =830\mathrm{mn}$, due to known data.²⁴

First, we discuss the data for $\lambda =685\mathrm{nm}$. Adjusting diffusion approximation curves, we have varied ${\mu }_s^{\prime }$ values in a wide range; fitting data for the lower concentration, $c=0.5\%$, with ${\mu }_s^{\prime }=5.5{\mathrm{cm}}^{-1}$ we consider further scattering coefficients to be proportional to the concentration within the studied range. The same values of ${\mu }_s^{\prime }$ were used for the Monte Carlo simulation.

The model based upon the Bethe–Salpeter equation contains the scattering length $l_s$ as the starting parameter. Thus, with ${\mu }_s^{\prime }$ being fixed, we need to preset the mean cosine $g$ for the determination of the scattering coefficient, using the definition ${\mu }_s^{\prime }={\mu }_s(1-g)$. We present simulation plots for two mean cosine values. As is seen, the plots for a fixed value of the reduced scattering coefficient do not differ noticeably for different scattering anisotropy parameters, in line with the diffusion approximation, within the simulation error. Calculating the amplitude and phase of the DPDW within the diffusion approximation, we used the closed equation derived in Ref. 23 accounting for the boundedness of the medium with the extrapolated boundary conditions.

In Fig. 4, we present the DPDW phase shift as a function of separation $\rho$, in comparison with simulations and with diffusion approximation results, for modulation frequency ${(2\pi )}^{-1}\omega =100\mathrm{MHz}$. The measurement data are shown for two solutions, with concentrations $c=0.5$ or 2%, which are simulated as media with ${\mu }_s^{\prime }=5.5{\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }=22{\mathrm{cm}}^{-1}$, correspondingly. The simulation plots are shown for two anisotropy parameters, $g=0.5$ or 0.8.

Fig. 4

Phase via source–detector separation with experimental data for two concentrations: open circle—concentration $c=0.5\%$ and open triangle—$c=2\%$. Monte Carlo simulations: lower plots, ${\mu }_s^{\prime }=5.5{\mathrm{cm}}^{-1}$: closed circle—anisotropy parameter $g=0.8$, and times—$g=0.5$; upper plots, ${\mu }_s^{\prime }=22{\mathrm{cm}}^{-1}$: open triangle—$g=0.8$ and plus—$g=0.5$; lines present diffusion approximation: solid line—${\mu }_s^{\prime }=5.5{\mathrm{cm}}^{-1}$ and dashed line—${\mu }_s^{\prime }=22{\mathrm{cm}}^{-1}$. The modulation frequencies ${(2\pi )}^{-1}\omega =100\mathrm{MHz}$, wavelength $\lambda =685\mathrm{nm}$.

As a whole, we admit an excellent agreement of measurements, diffusion approximation, and numerical simulations. As is seen, the Monte Carlo curves as well as the diffusion approximation quite perfectly follow the experimental data, up to the range of small separations of order of several transport lengths. We also admit that when using the transport length as the spatial scale, one does not observe a noticeable dependence of the studied DPDW parameter on the scattering anisotropy.

Comparing statistical errors, defined as the standard deviations related to the diffusion approximation values, we calculated, as an example, two plots of phase shift for source–detector separations from 2 to 10 mm and ${\mu }_s^{\prime }=1.1{\mathrm{mm}}^{-1}$, using two algorithms [Eqs. (12) or (14)]. We found relative errors 2.2% and 8.4% with the sampling volume $N_{\mathrm{ph}}={10}^4$, and relative errors 0.25% and 0.65% with $N_{\mathrm{ph}}={10}^6$, respectively. As is seen, the errors are significantly smaller for our algorithm.

In Fig. 5, the measurement data, diffusion approximation plot, and Monte Carlo simulation results for function $F(\rho )=\lg ({\rho }^2\sqrt{I^2+Q^2})$ for a modulation frequency of ${(2\pi )}^{-1}\omega =200\mathrm{MHz}$, which are usually exploited for presentation of an intensity decrease with the source–detector separation for a half-space geometry, are also presented. Experimental data are given for concentration $c=0.5\%$, in comparison with the numerical result for a model with ${\mu }_s^{\prime }=5.5{\mathrm{cm}}^{-1}$. Plots with anisotropy parameters $g=0.5$ and $g=0.8$ are practically identical; it means that for the fixed transport length the light scattering does not depend on the scattering anisotropy, beyond the relationship ${\mu }_s^{\prime }={\mu }_s(1-g)$.

Fig. 5

Diffuse photon density wave (DPDW) intensity function, $F(\rho )=lg({\rho }^2\sqrt{I^2+Q^2})$ against the source–detector separation $\rho$, in intralipid solution with concentration open circle—$c=0.5\%$ in comparison with the Monte Carlo data for ${\mu }_s^{\prime }=5.5{\mathrm{cm}}^{-1}$: closed circle—$g=0.8$ and times—$g=0.5$; solid line presents diffusion approximation for ${\mu }_s^{\prime }=5.5{\mathrm{cm}}^{-1}$. Modulation frequency ${(2\pi )}^{-1}\omega =200\mathrm{MHz}$, $\lambda =685\mathrm{nm}$.

We admit once more an illustrative agreement for data obtained with measurements, simulations, and with the diffusion equation. Some deviations are observed at quite small separations, where the determination of local intensity itself becomes uncertain, both experimentally and theoretically.

The similar plots on the source–detector separation have been obtained earlier¹⁸ using parallel computing successfully comparing the Monte Carlo simulation and the diffusion approximation, presenting results for isotropic scattering and zero modulation frequency. Here, we presented a demonstrative agreement between the Monte Carlo simulation, the diffusion approximation, and measurements, simultaneously, for different scattering anisotropies.

6.

Optical Parameters of Intralipid Suspension

Let the water–intralipid solution be described as a suspension of $N$ identical intralipid spherical particles with diameter $D$ randomly distributed in points ${\mathbf{R}}_1,\dots ,{\mathbf{R}}_N$ in a uniform medium. The local permittivity $ϵ(\mathbf{r})$ of such a suspension can be presented as

Within the Born approximation, the scattering cross section $\tilde{G}(\mathbf{q})$ is the Fourier transform of a pair correlator of random local permittivity

Thus, the optical theorem permits to present the scattering coefficient in the closed form

A similar relationship is valid for the reduced scattering coefficient

Given the refractive indices of water $n_0=1.33$ and intralipid $n_s=1.47$, we get $\delta ϵ=0.39$ for the permittivity mismatch, and calculate within the Born approximation the scattering and transport lengths for varying sizes of scatterers. The results are given in Table 1. The last line with close values of scattering coefficients presents the extrapolation of data²⁷ obtained for the 10% intralipid solution using the Mie equation.

Table 1

Scattering and reduced scattering coefficients, in cm−1, and scattering anisotropy g=cos θ¯ of intralipid-water 1% suspension for different scatterer diameters D. Permittivity mismatch Δϵ=0.39.

Performing simulations for $\lambda =685\mathrm{nm}$, we took ${\mu }_s^{\prime }=11{\mathrm{cm}}^{-1}$ and $g=0.8$; as seen from Table 1, these values correspond to the monodisperse 1% suspension of intralipid particles $0.5\mu \mathrm{m}$ in diameter, in agreement with the manufacturer’s specifications. In its turn, in the case of longer wavelength, $\lambda =830\mathrm{nm}$, for the same $0.5\text{-}\mu \mathrm{m}$ diameter intralipid particles, we get ${\mu }_s^{\prime }=8.73{\mathrm{cm}}^{-1}$ for the same 1% concentration. These results permit us to fix also the anisotropy parameter $g$.

In Figs. 6 and 7, the simulation data on the DPDW phase and intensity obtained for the IR radiation with $\lambda =830$ are presented; for the reduced scattering coefficient, we took ${\mu }_s^{\prime }=4.37{\mathrm{cm}}^{-1}$ for concentration $c=0.5\%$ and ${\mu }_s^{\prime }=8.74{\mathrm{cm}}^{-1}$ for $c=1\%$ concentration, calculated for the suspension model of the $0.5\text{-}\mu \mathrm{m}$ diameter intralipid spheres within the Rayleigh–Gans approximation. For absorption, we took ${\mu }_a=0.04{\mathrm{cm}}^{-1}$ due to Ref. 24.

Fig. 6

DPDW phase as a function of source–detector separation $\rho$ for longer wavelength $\lambda =830\mathrm{nm}$. Measurement data: open circle—$c=0.5\%$ and open triangle—$c=1.0\%$; for smaller concentration ${\mu }_s^{\prime }=4.37{\mathrm{cm}}^{-1}$, closed circle—Monte Carlo data and solid line—diffusion approximation; for higher concentration ${\mu }_s^{\prime }=8.74{\mathrm{cm}}^{-1}$, closed triangle—Monte Carlo and dashed line—diffusion approximation. Modulation frequency ${(2\pi )}^{-1}\omega =100\mathrm{MHz}$.

Fig. 7

A similar plot, as in Fig. 6, for intensity decrease, $F(\rho )$.

Thus, for $\lambda =830\mathrm{nm}$, we obtain a fair agreement of measurements, Monte Carlo simulations, and diffusion approximations, without any fitting parameters; we determine the value of the reduced scattering coefficient using the knowledge of the scatterer size, obtained in the experiment with $\lambda =685\mathrm{nm}$.

7.

Conclusion

Thus, we performed the Monte Carlo simulations using the developed approach, varying the parameter of scattering anisotropy and comparing the results with measurements as well as with the diffusion approximation predictions. We observed a fair agreement across all three methods.

The Bethe–Salpeter equation, which is a microscopic justification of the radiative transfer equation, contains well-known assumptions, such as the scalar field model, ladder diagram approximation, and artificial phase function. Thus, the demonstrated agreement between experimental data and plots of the DPDW phase and amplitude calculated on the basis of the Bethe–Salpeter equation once more reaffirms these assumptions.

The diffusion approximation equations indicate that if one scales the extinction in terms of the transport length, the description of radiative transfer becomes universal and does not have any dependence on the scattering anisotropy. Comparing Monte Carlo simulation data, diffusion approximation results, and measurement data, we have shown that this is true well beyond the expected applicability of the diffusion approximation with proper accounting for the boundary conditions, at least for a half-space. We performed simulations for mean scattering cosine $g=0.5$ and $g=0.8$, and did not find a noticeable difference between corresponding phase and intensity plots, apart the well-known universal dependence on $g$ contained in the relationship between scattering and reduced scattering coefficients.

We performed measurements and simulations for two different near-IR wavelengths. We find that both measurements and calculations yield supplemental information on the size of scattering inhomogeneities. Given the size of scatterers, we achieved a fair agreement between measurements and simulations without the use of fitting parameters.

In our modeling scheme, every scattering event contributes to the detected signal. This makes our approach much more efficient as compared to the commonly used Monte Carlo algorithm wherein one simulates the migration of an incident photon until it reaches the detector. Performing comparative simulations for the two approaches, we have shown that the conventional approach based on a count of photons escaping the medium requires a sampling volume that is an order larger than the one required for our approach based upon the Bethe–Salpeter equation.

Simulating the series over scattering orders, one must satisfy rigid requirements on the number of scattering orders taken into account. Calculating one plot of the DPDW parameters dependent on the source–detector separation distance takes about 10 min for several thousand scattering events and sampling volume of ${10}^5$ photons. The procedure presented permits calculation of the DPDW parameters by means of a commonly available personal computer in a run-time which gives our method the potential for real-time monitoring of biomedical practices.

Our approach, based upon the theoretical-field Bethe–Salpeter equation, may be very appropriate for simulation of light migration while accounting for the electromagnetic nature of radiation, i.e., polarization effects and optical anisotropy of tissues. Our approach can also be readily applied for various geometries of tissues, and in the future a multilayer medium model may be implemented for studying the optical properties of multilayer tissues. We hope the considered MC algorithm may also be useful for tissue studies with the temporal domain method.