## Abstract

A heuristic method for estimating the reduced scattering coefficient (µ_{s}’) of turbid media using time-resolved reflectance is presented. The technique requires measurements of the distributions of times-of-flight (DTOF) of photons arriving at two identical detection channels placed at unique distances relative to a source. Measured temporal shifts in DTOF peak intensities at the two channels were used to estimate µ_{s}’ of the medium using Monte Carlo (MC) simulation-based lookup tables. MC simulations were used to compute temporal shifts in modeled reflectance at experimentally employed source-detector separations (SDS) for media spanning a wide range of optical properties to construct look up tables. Experiments in Intralipid (IL) phantoms demonstrated that we could retrieve µ_{s}’ with errors ranging between 6-25% of expected (literature) values, using reflectance measured across 650-800 nm and SDS of 5-15 mm. Advantages of the technique include direct processing of measured data without requiring iterative non-linear curve fitting. We also discuss applicability of this approach for media with low scattering coefficients where the commonly employed diffusion theory analysis could be inaccurate, with practical recommendations for use.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Time-resolved reflectance spectroscopy (TRS) has the potential to directly provide bulk scattering and absorbing coefficients of a turbid medium without making assumptions about the medium’s composition or structure [1–4]. In TRS, picosecond laser pulses (fired at MHz rates) are injected into a turbid medium and the multiply scattered and attenuated diffusely reflected pulses are detected using a fast single photon counting photodiodes that are usually coupled to time-correlated single photon counting (TCSPC) boards to measure the photon distribution of times-of-flight (DTOF) [4]. These measurements are then quantified using theoretical (or numerical) approaches to extract the absorption (µ_{a}) and reduced scattering (µ_{s}’) coefficients of the medium [2–5]. In biological media, accurately and efficiently recovering the optical transport coefficients can help parametrize a variety of functional and structural properties of biomedical and clinical interest [6,7].

Uncoupling µ_{a} and µ_{s}’ using experimentally measured DTOF is known to be a difficult problem that needs careful measurements and calibrations [1–5,8–10]. The measured DTOF is a convolved response of the theoretical temporal point spread function (TPSF) and the instrument response function (IRF) of the experimental system [11]. The TPSF represents the (Green's function) response of the tissue medium to an incident Dirac delta pulse, while the IRF represents the finite temporal profile of the incident laser pulse measured as it propagates through the detection optics [11]. Extraction of the optical properties is typically done by iteratively convolving a theoretical TPSF with the experimentally measured IRF and fitting the convolved response via non-linear least squares to the measured DTOF [8,12,13]. Since the theoretical TPSF is computed from known optical properties, the process yields optical properties for the medium at convergence.

The above approach requires both having appropriate theoretical estimates to accurately model the TPSF, and accurate estimates of the system IRF. Since the IRF is a function of the laser source, fiber optics, and detection electronics used to deliver and collect signals, the IRF must be measured for all detection channels and wavelengths used experimentally [14,15]. Curve fitting of the DTOF using diffusion theory (DT) is usually the most common method for extracting optical properties of turbid media using TRS [16]. Although the process requires iterative reconvolution based reconstructions and are usually performed off-line, growth of computational resources could render these methods to operate in real-time [16]. However, applications of such models to analyze biological tissues shown large variance in extracted coefficients [1–3,15,17] particularly in scattering as reported recently [18]. Inverse methods using DT can also be computationally intensive and have been reported to show crosstalk between the derived optical coefficients [1,3,9,10]. Thus, a method such as the one we have presented could potentially facilitate improved fitting with DT.

Limiting crosstalk between recovered coefficients is critical for accurate estimation of optical properties, which in turn impact clinical and diagnostic utility of these methods. Different approaches have been described to constrain inverse fitting algorithms and increase quantitative accuracy of recovered optical parameters [19–21]. For example, measurements at multiple wavelengths were used to constrain spectral properties of the medium in order to improve the accuracy of recovered parameters [21]. Thus, methods that can extract (or even constrain) the optical coefficients of the medium can help DT based models in estimation of optical properties.

Here, we present a simple-to-use, heuristic technique that facilitates estimation of the µ_{s}’ of a medium directly from experimentally measured DTOFs. This method uses lookup tables constructed using MC simulations and translates measured peak-to-peak time differences (Δt) in DTOF obtained at two distinct SDS into µ_{s}’ of the medium. We describe implementation of the method and demonstrate its application to experimental data obtained from well-described tissue simulating phantoms [13,22,23].

## 2. Materials and methods

#### 2.1 Instrumentation

Figure 1(a) shows the schematic of the experimental system, where the input pulse from a super-continuum laser (SC400, NKT Photonics, DK) was spectrally filtered via a band-pass filter (SuperK VARIA, NKT Photonics, Denmark) and coupled into one end of an optical (source) fiber (diameter=400 µm; NA=0.22; length=1 m). The distal end of the source fiber illuminated a medium of interest and formed the sensing head of a custom-fabricated optical probe (Gulf Photonics, FL). The sensing head had three detection channels (made from three optical fibers identical to the source fiber) which were epoxied at distances of 5, 10, and 15 mm from the source fiber – thus, our sensing probe had a fixed geometry. Reflectance measured from a selected channel (identified uniquely using the SDS) and was coupled into a single-photon avalanche diode (SPAD) detector (PMD-050, MPD, Italy), that was electronically coupled to a time-correlated single photon counting (TCSPC) board (SPC-130, Becker & Hickl, Germany). The laser repetition rate was 40 MHz and an electronic sync signal from the laser was used to trigger the TCSPC for signal acquisition.

Figure 1(b) shows representative time-resolved reflectance measurements at the three experimental SDS used. Figure 1(c) shows a magnified view of the DTOF peaks. The IRF was measured to ensure equality for each channel across all wavelengths by reflecting the source pulse into the detecting fiber using a mirror with a white piece of paper placed over the detecting fiber. Figure 1(d) shows the IRF measured at 650 nm for all three detector channels used. The average root mean square error (RMSE) calculated over three orders of magnitude for all IRFs measured across all channels and wavelengths used was 0.008. For comparison, repeated scans from a fixed target surface had RMSE of nearly 0.001, while DTOFs at differing absorption and scattering coefficients produce RMSE values larger than 0.1. Temporal stability of the peak position was also monitored over several hours of continuous acquisition and was shown to vary less than 0.7% per hour.

#### 2.2 Monte Carlo simulations

Although DT can be used to model TPSF measurements at multiple SDS, it is known to be inaccurate for modeling reflectance at early times and small SDS [1,4,12,13,17]. Therefore, a previously developed time-resolved Monte Carlo (MC) model [24] was used to calculate the TPSF in a semi-infinite, homogenous geometry for 48 different tissue models. Each TPSF had distinct µ_{s}’ and µ_{a} values spanning 3–18 cm^{-1} and 1.0 × 10^{−3−}0.5 cm^{-1}, respectively. Tissue models were generated by permuting six different µ_{s}’ values (3, 6, 9, 12, 15, 18 cm^{-1}) with eight different µ_{a} values (1.0×10^{−3}, 0.026, 0.050, 0.075, 0.18, 0.29, 0.39, 0.50 cm^{-1}). Absorption values were chosen on a log-scale to more accurately sample changes in absorption over several orders of magnitude.

Simulations were run to generate the TPSF (with photons being launched into and detected from the medium using optical fibers matching the experimental probe) at SDS of 5, 10 and 15 mm for all 48 tissue models. A temporal resolution of 2 ps was used for simulations at SDS of 5 mm, while the temporal resolution was 10 ps at SDS of 10 and 15 mm. Each MC run used 3×10^{8} photons, an anisotropy equal to 0.7 (using the Henyey-Greenstein phase function), an index of refraction equal to 1.35 (to match the IL phantom) [22,25,26], and fiber diameters and numerical apertures of 400 µm and 0.22, respectively, to match experimental measurements. Simulated TPSF data were resampled (using a cubic spline) to have a resolution of 0.1 ps. The peak-times calculated from the cubic spline using two different temporal resolutions in generating the TPSF were compared to ensure no errors were introduced by resampling.

Representative TPSFs from MC simulations are shown at two SDS in Fig. 2(a) (SDS = 5 mm) and Fig. 2(b) (SDS = 15 mm) for three different media (circles: µ_{s}’ = 6, µ_{a} = 0.18 cm^{-1}; squares: µ_{s}’ = 9, µ_{a} = 0.18 cm^{-1}; triangles: µ_{s}’ = 6, µ_{a} = 0.39 cm^{-1}). Vertical lines show t_{max} for each simulated reflectance (derived from the resampled data). These data show that in media with identical µ_{s}’ an increase in µ_{a} causes a decrease in t_{max} (by 7% at SDS of 5 mm and 13% for SDS of 15 mm for media shown by circles vs. triangles), while in media with identical µ_{a} an increase in µ_{s}’ causes an increase in t_{max} (by 11% at SDS of 5 mm and by 35% at SDS of 15 mm for media shown by circles vs. squares). These figures also show that at short SDS and low scattering, the absorption coefficient only weakly influenced t_{max}.

The difference (Δt) between t_{max} at SDS pairs of 5-10 mm, 10-15 mm and 5-15 mm, were obtained for each of the 48 simulations. For each SDS pair, Δt values across the 48 simulated tissue models were linearly interpolated to create 2D lookup tables of Δt values for µ_{a} spanning 1.0 × 10^{−3} - 0.5 cm^{-1} (in 1.0 × 10^{−3} cm^{-1} increments) and µ_{s}’ spanning 3-18 cm^{-1} (in 0.01 cm^{-1} increments). Several interpolations methods were tested (e.g. linear, cubic, spline) to create lookup tables and showed lesser than 0.1% variation in generated tables. These computed Δt distributions are shown for two different SDS pairs – in Fig. 2(c) (SDS pair of 5-10 mm) and Fig. 2(d) (SDS pair of 5-15 mm). Symbols for Δt_{5,15} values computed for data shown in Fig. 2(a) and 2(b) are marked on the Fig. 2(d) (shape and color of the marker identifies the optical properties). The horizonal grey line in Fig. 2(c) illustrates how a given Δt value (calculated using SDS of 5 and 10 mm) confines the range of µ_{s}’ (shown by the shaded rectangle). In other words, Δt_{5,10} = 50 ps could represent media with µ_{s}’ ranging between 7.1–10.7 cm^{-1} and µ_{a} varying from 1.0×10^{−3} and 0.5 cm^{-1}.

#### 2.3 Lookup method

Data shown in Fig. 2(c) and 2(d) form the crux of the lookup method. In the forward direction, specific values for µ_{s}’ (x-axis) and µ_{a} (color) identifies a single point in the shaded region – the ordinate of this point is Δt. In the inverse sense, given a Δt value (using a given SDS pair) it may be associated with several pairs of µ_{s}’ and µ_{a} values (identified by the horizontal line within the shaded figure in Fig. 2(c)). Thus, translation of a measured Δt for a pair of SDS into µ_{s}’ would require knowledge of µ_{a}, or vice-versa. As is evident from these figures, for fixed Δt, the limits on µ_{s}’ are much more restrictive than the limits on µ_{a}. For the range of optical properties considered here, the required Δt for a SDS pair could be produced by media with µ_{a} that span the full range of simulated values (1.0 × 10^{−3−}0.5 cm^{-1}) but had a much more confined range for possible µ_{s}’ values. In other words, the uncertainty in estimation of µ_{s}’ given an incorrect guess of µ_{a} would be much lower than the uncertainty in estimating µ_{a} given an incorrect guess of µ_{s}’. Thus, our lookup method translates the measured Δt (from a given SDS pair) directly into µ_{s}’ via the interpolated 2D MC lookup table for the corresponding SDS, but needs an input value for the µ_{a} of the medium to function.

#### 2.4 Estimating the absorption coefficient

Two different approaches were used to estimate the µ_{a} required for using the lookup method. The first approach, referred to here as the ‘tail method’, uses a limit derived from DT where the slope of the natural logarithm of the DTOF at long time scales is assumed to be directly proportional to µ_{a} (Eq. (1)) [27]. $R({r,t} )$ shown in Eq. (1) represents the measured DTOF and for longer times, the value of the estimated absorption from Eq. (1) would improve. Using this technique to accurately estimate µ_{a} would require deconvolving the IRF from the measured DTOF and having the reflectance signal span many nanoseconds. In practice, a limited dynamic range in the detection system and the drawbacks in deconvolving the IRF from the DTOF reduces the accuracy of this approach in estimating µ_{a}.

However, as we expect Δt to only weakly depend on µ_{a}, the tail method provides a simple and direct way to extract (even if crudely) a value for use with the lookup method to estimate µ_{s}’. Further, due to a limited dynamic range of our instrumentation, the DTOF time scale was shifted such that the peak occurred at $t$=0, and Eq. (1) was applied to analyze the measured reflectance signal from its peak value till it fell to 0.1% of the peak value. Derivatives were numerically calculated between successive times, smoothed and translated to calculate µ_{a} for each measured DTOF using Eq. (1) [27]. Because calculated numerical values were noisy and showed variations for different SDS (in the same media), µ_{a} was estimated as the averaged tail-method estimate for every scan, at all three SDS used.

In the second approach, referred to here as the ‘transmittance method’, the intrinsic absorption coefficient of hemoglobin was calculated by measuring the collimated transmittance using a spectrophotometer at the varying concentrations of bovine hemoglobin used experimentally (described in Sec. 2.5). µ_{a} was then calculated from the percent transmittance using Beer’s Law for each concentration, and the intrinsic absorption coefficient (${\epsilon _a}$) was determined from the slope of the linear fit of the calculated µ_{a} at each concentration. Expected (true) absorption coefficients of the medium distinct from values obtained using Eq. (1) were then calculated using ${\epsilon _a}$ and the mass concentration of absorber used. Dissolving large amounts of solid hemoglobin (∼255 mg for each concentration change) as well as uncertainties in cuvette calibrations in transmittance measurements could impact the accuracy of the estimated µ_{a} by this method.

#### 2.5 Phantom preparation and measurements

Liquid phantoms were prepared to experimentally validate the performance of the developed method. Phantoms were prepared by mixing 40 mL of 20% Intralipid (IL) (Sigma-Aldrich; MO, USA) with 750 mL of de-ionized water while absorption was independently introduced by serial additions of 250 mg of dry bovine hemoglobin (Hb) (H3760; Sigma-Aldrich; MO, USA) to the above solution. Liquid phantoms were prepared in a cylindrical container (radius = 6 cm, height = 12 cm) and preliminary experiments showed finite-boundary effects at low scattering solutions (< ∼1cm^{-1}) and shallow depths (< ∼4 cm). Therefore, the minimum µ_{s}’ used for the experiments was 5 cm^{-1} and the minimum distance between optical channel and container wall was ensured to be greater than 5 cm to best approximate a semi-infinite medium as modeled by the MC lookup tables.

Phantoms were created to have constant scattering coefficient with ∼5% IL by volume [28,29]. Since the optical properties of IL have been well described to be stable and reproducible [22,23], we could also use reference values from previous reports [13,22,23,25] to directly compare the reduced scattering coefficient we estimated to those reported previously. Hb was added five times to create a set of 6 different phantom media having fixed scattering increasing absorption coefficients (corresponding to hemoglobin concentrations ranging between 0-250 µM). This produced phantoms with expected µ_{a} values varying between 0.003-0.4 cm^{-1} for 650 nm, and between 0.022-0.14 cm^{-1} for 800 nm.

Each phantom solution was prepared in the container and mixed gently with a magnetic stirrer during measurements. The sensing head of the probe (mounted on a custom probe holder) was lowered until the surface tension of the solution was broken. The probe was clamped in place and the following sequence of measurements were obtained for each phantom. First, the laser wavelength was set to 650 ± 5 nm, the optical fiber for the channel with SDS of 5 mm was manually coupled to the SPAD (via a standard SMA fiber socket) and three repeated scans were acquired. Next, the bandpass was adjusted to illuminate the sample at 700 ± 5 nm, 750 ± 5 nm and 800 ± 5 nm and three repeated TRS measurements were acquired. The same sequence of steps was repeated by manually coupling fibers for the 10 and 15 mm SDS channels into the SPAD, across the same four source illumination wavelengths. Each TRS acquisition was obtained by allowing the TCSPC signal to acquire the signal for 30 s.

Thus, a total of 36 TRS measurements were obtained (3 scans per wavelength, for 4 wavelengths across 3 channels), for each phantom solution. TRS scans collected at 650 nm are shown for all detection channels used (with SDS of 5, 10 and 15 mm) in Fig. 1(b). The difference between t_{max} of measured DTOF at 5 and 10 mm SDS is also shown and labelled as Δt_{5, 10} in Fig. 1(c). A similar naming convention was adopted when peak differences were calculated using SDS of 5 and 15 mm (Δt_{5, 15}) and SDS of 10 and 15 mm (Δt_{10, 15}), as denoted in Fig. 1(c). Experimental DTOF measurements were obtained with time-resolutions of 12 ps.

## 3. Results

As described above, the MC lookup table was used to convert the measured Δt (for a given SDS pair) into µ_{s}’, given an estimate of µ_{a}. We first tested the developed method to extract µ_{s}’ for phantom that was mostly scattering with very low absorption – i.e. a solution of 5% IL in water. Figure 3 shows extracted µ_{s}’ as a function of illumination wavelength for the phantom with 0 µM hemoglobin (5% IL in water) using the lookup method at all three experimental SDS pairs possible (circles: Δt_{5, 10}; diamonds: Δt_{5, 15}; squares: Δt_{10, 15}; stars: literature average from Refs. [13,22,23]). Input values for µ_{a} were estimated using the tail method from measured data. Derived µ_{s}’ values in Fig. 3 show expected wavelength dependent decreases in scattering for all the three SDS pairs and tracked expected values from literature. Although Δt_{5, 10} produced µ_{s}’ closest to reference values, both Δt_{5, 15} and Δt_{10, 15} showed agreement (to better than 25%) with the reported intrinsic scattering parameters of IL [22]. Larger errors are observed at smaller SDS as the uncertainty in derived µ_{s}’ was high due to the uncertainty in correctly estimating the peak time with the experimental temporal resolution of 12 ps (TCSPC bin-width). Error bars shown represent the standard error from three repeated measurements for each concentration at each wavelength and SDS.

The values of µ_{a} used to derive µ_{s}’ shown in Fig. 3 were estimated using the tail method and could yield inaccurate results. We therefore examined the impact of changing the input µ_{a} on the extracted value of µ_{s}’, for any given DTOF. This was done by first, using an initial µ_{a0} (obtained for example using the tail method) to derive the scattering coefficient µ_{s0}’ using the lookup method. This value of µ_{a0} was then decreased by 20 - 100% in five equal intervals (where a decrease of 100% was the smallest simulated µ_{a} = 1.0×10^{−3} cm^{-1}) to determine an updated absorption coefficient µ_{a1} and then using µ_{a1} as the input to the lookup method to determine the corresponding updated scattering coefficient µ_{s1}’.

Percent errors in µ_{s1}’ (relative to µ_{s0}’) are shown in Fig. 4, as a function of percent changes in input µ_{a0}, for the three experimental SDS pairs (bars) used here. Data in Fig. 4 represent DTOFs obtained from four different (representative) media – Medium 1 (Fig. 4(a)): µ_{s0}’ = 10, µ_{a0} = 0.05 cm^{-1}; Medium 2 (Fig. 4(b)): µ_{s0}’ = 10, µ_{a0} = 0.16 cm^{-1}, Medium 3 (Fig. 4(c)): µ_{s0}’ = 13, µ_{a0} = 0.05 cm^{-1}; Medium 4 (Fig. 4(d)): µ_{s0}’ = 13, µ_{a0} = 0.16 cm^{-1}). DTOFs to estimate µ_{s0}’ in Fig. 4(a) and 4(c) were obtained from the pure scattering phantom (no hemoglobin) for illumination with 800 and 650 nm, respectively. DTOF used in Fig. 4(b) corresponds to the phantom with 256 µM Hb for 800 nm illumination, while DTOF for Medium 4 (Fig. 4(d)) was from the phantom with 53 µM Hb with illumination at 650 nm (these phantoms were selected because their absorption coefficients for the wavelengths and concentrations used were comparable to each other).

It is clear to see from Fig. 4 that the derived µ_{s}’ using Δt_{5, 10} varied less than 8% as the input µ_{a} values were varied by 100% in media with lower absorption (Fig. 4(a) and 4(c)), as was expected from the MC simulations. Additionally, even for the longest SDS pair (10-15 mm) used, the change in extracted µ_{s}’ was lower than 20% as the input absorption coefficients changed by 100%. However, when µ_{s0}’ was obtained with higher absorption coefficients (µ_{a0} of 0.16 cm^{-1}, Fig. 4(b) and 4(d)) the changes in the medium’s absorption translated to larger variations in extracted µ_{s}’. In general, the larger the error (or uncertainty) in knowledge of a medium’s true µ_{a} the larger the error in the recovered µ_{s}’.

Next, we tested the performance of the lookup technique to estimate µ_{s}’ in phantoms as the absorption coefficient was varied. Acquired DTOF’s using all available SDS and illumination wavelengths (with three SDS pairs and four wavelengths) were analyzed using the MC lookup method to translate the measured values of Δt into µ_{s}’, in each of the 6 experimental phantoms prepared. Figure 5(a) shows the µ_{s}’ extracted by the lookup method (for each SDS pair) in phantoms with varying hemoglobin concentrations for illumination with 750 nm and obtaining the input µ_{a} coefficients from the DTOFs using the tail method. Figure 5(b) shows these data for the same phantoms shown in Fig. 5(a) but used the transmittance method (described in Section 2.4) to determine input µ_{a}. Both Fig. 5(a) and 5(b) indicate that that the µ_{s}’ values obtained from Δt_{5,10} with the lookup method matched the expected (literature) values (shown as dashed horizontal lines [13,22,23]). Further, all three SDS pairs extracted consistent µ_{s}’ values in all the six phantom media that had different absorption properties.

Table 1 shows the mean percent errors between lookup table derived µ_{s}’ and literature values averaged across all wavelengths, for each SDS pair while using the tail method to determine the µ_{a} used in the lookup table. Percent errors reported in Table 1 varied lesser than 7% of the values obtained when the transmittance method was used to estimate µ_{a}. Largest deviations between the two approaches were observed in media with highest absorption where the influence of the IRF impacts the derived values most significantly.

At wavelengths of 750 and 800 nm, estimating µ_{a} via the tail method provided improved accuracy in estimating µ_{s}’. However, at wavelengths of 650 and 700 nm, calculating µ_{a} using the tail method decreased accuracy in estimating µ_{s}’ (µ_{a} obtained via transmittance data yielded lower errors in extracted µ_{s}’ in these cases). Best accuracy was obtained for media with low µ_{a} and by using SDS pairs closest to the source (these data cells are highlighted in Table 1).

## 4. Discussion

We have described a heuristic technique capable of assessing the reduced scattering coefficient of a homogenous medium directly from experimentally measured TRS signals. We have shown that an experimentally measured peak-time difference between DTOFs at two different SDS can be translated to obtain a robust estimate of the medium’s scattering coefficient using preconstructed MC lookup tables. Although our method requires an input (estimated) value for the absorption coefficient, the derived µ_{s}’ using our approach is only weakly dependent on the input absorption coefficient. Application to experiments with tissue phantoms showed that we could estimate the expected (true) scattering coefficients in media with errors ranging from 5% -25%.

A principal advantage of the presented technique is its inherent speed and simplicity of use (i.e. directly being able to translate measurements into the reduced scattering coefficient of the medium). However, it is worth noting that there are important instrumentation calibration and prerequisite conditions that must be satisfied for it to function accurately. Although our method does not use the IRF for translating measurements into the medium's scattering coefficient, it is not independent of it. It is critical for the IRF (i.e. the intrinsic temporal shape of the incident laser pulse) to remain stable across the duration of experiments and also be identical for each measurement channel used. In our case, each channel used the same length and type of fiber and were all detected by a common detector and as shown in Fig. 1(c) were identical across channels for all wavelengths. Systems that use different fiber lengths and/or employ multiple photodetectors that seek to use this technique might need to consider such issues carefully during development of lookup tables.

Phantom results showed that errors between literature (expected) values and derived values were (on average) lower than 6% when data from SDS of 5 and 10 mm were used and were as high as 25% when data from SDS of 10 and 15 mm were used. Thus, utilizing smaller SDS would yield the best estimates of scattering, consistent with previous reports [5,30–32]. We note that accurately resolving peak time-differences for small SDS requires higher temporal resolution of measured data. Further, non-linear distortions in t_{max} (from the convolved IRF response) could impact longer SDS measurements more and may explain our findings.

Our findings are consistent with previous reports showing that the peak arrival time (t_{max}) of the DTOF is strongly sensitive to µ_{s}’ and only weakly to µ_{a} [15,17], especially for shorter SDS. It also builds off separate work that seek to mitigate the influence of the IRF in TRS analysis [8,11,33]. Data shown in Fig. 2(a) and 2(b) reflect these findings, where both µ_{s}’ and µ_{a} impact the peak time to varying degrees, but as denoted in Fig. 2(c) (shaded bar), a fixed value for Δt (for a given SDS pair) tightly binds the range of allowed µ_{s}’ values, but not µ_{a}. Although the technique needs an input estimate for µ_{a}, this can be derived by applying Eq. (1) directly to the measured DTOF. Accurately using Eq. (1) to estimate µ_{a} would require a large dynamic range in the detecting system and deconvolving the IRF.

Here, we have shown (Fig. 4) that only a rough estimate of µ_{a} is needed as it weakly affects the peak position. Utilizing smaller SDS decreases the effect of absorption on the peak time difference leading to more accurate estimates of µ_{s}’ by reducing the probed volume. Therefore, small SDS would be more sensitive to shallower layers of the medium, while longer SDS could differentially be used to probe deeper tissues. The extension of this method to inhomogeneous media will be addressed in future investigations. Although not experimentally studied here, any fiber geometry with at least two detector channels with unique SDS would be amenable for use with the lookup method presented. Further, using this method to derive the scattering by setting µ_{a}=0 could be used to obtain a lower bound on µ_{s}’.

We also note that the SDS must be well defined (within 0.5 mm) to accurately estimate µ_{s}’. Peak times between DTOF’s at 10 and 11 mm SDS could vary by ∼ 50% translating into recovered error of ∼ 80% in µ_{s}’. Additionally, although MC simulations and experimental geometries must be matched (i.e. we need experimental measurements from large volumes to represent semi-infinite media) we note a specific advantage of this method in “finite” geometries as only the earliest arriving photons determine Δt – and thus, could be robust boundary effects.

The described method could be used to bootstrap (or confine the range of µ_{s}’) other analytical inverse solvers to quantify optical properties from TRS measurements, especially in weakly scattering media where DT-based approaches are typically more vulnerable to higher errors [1–3,15,17]. Uncertainty in the index of refraction of the medium also weakly affects the relative peak-time difference. Lastly, since this method works by establishing the temporal peak of measured DTOF reflectance profiles, it should be possible to sparsely sample the DTOF histogram to acquire data faster. These topics will be subject of future investigations.

## 5. Conclusion

An easy-to-use approach to extract µ_{s}’ of turbid media that overcomes current limitations in TRS analysis was presented. By using a pair of short SDS (<10 mm), we show that µ_{s}’ can be estimated within 6% of reference values with reasonable estimations of µ_{a}. The approach could be further extended to include timing differences of later arriving photons as well as aid traditional inverse solvers in further parameterizing optical properties. The method shows great promise in recovering a medium’s optical transport properties via TRS in regimes where common techniques fail, while facilitating quantitation in real-time without directly measuring the IRF.

## Acknowledgements

We acknowledge support from Jens Mueller at the Research Computing Group at Miami University for discussion and guidance regarding setting up parallel runs of MC simulations on the supercomputing cluster.

## Disclosures

The authors declare no conflicts of interest related to this article.

## References

**1. **V. G. Cubeddu R, A. Pifferi, P. Taroni, and A. Torricelli, “Experimental test of theoretical models for time-resolved reflectance,” Med. Phys. **23**(9), 1625–1633 (1996). [CrossRef]

**2. **A. Pifferi, A. Torricelli, A. Bassi, P. Taroni, R. Cubeddu, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, J. Swartling, T. Svensson, S. Andersson-Engels, R. L. P. Van Veen, H. J. C. M. Sterenborg, J.-M. Tualle, H. L. Nghiem, S. Avrillier, M. Whelan, and H. Stamm, “Performance Assessment of Photon Migration Instruments: The MEDPHOT Protocol,” Appl. Opt. **44**, 2104–2114 (2005). [CrossRef]

**3. **S. A.-E. and T and S. Erik Alerstam, “Improved accuracy in time-resolved diffuse reflectance spectroscopy,” Opt. Express **16**(14), 10440–10448 (2008). [CrossRef]

**4. **A. Torricelli, D. Contini, A. Pifferi, M. Caffini, R. Re, L. Zucchelli, and L. Spinelli, “Time domain functional NIRS imaging for human brain mapping,” NeuroImage **85**, 28–50 (2014). [CrossRef]

**5. **A. Pifferi, D. Contini, A. D. Mora, A. Farina, L. Spinelli, and A. Torricelli, “New frontiers in time-domain diffuse optics, a review,” J. Biomed. Opt. **21**(9), 091310 (2016). [CrossRef]

**6. **R. H. Wilson, K. Vishwanath, and M.-A. Mycek, “Optical methods for quantitative and label-free sensing in living human tissues: principles, techniques, and applications,” Adv. Phys.: X **1**(4), 523–543 (2016). [CrossRef]

**7. **M. Chandra, K. Vishwanath, G. D. Fichter, E. Liao, S. J. Hollister, and M.-A. Mycek, “Quantitative molecular sensing in biological tissues: an approach to non-invasive optical characterization,” Opt. Express **14**(13), 6157 (2006). [CrossRef]

**8. **D. Milej, A. Abdalmalak, D. Janusek, M. Diop, A. Liebert, and K. St. Lawrence, “Time-resolved subtraction method for measuring optical properties of turbid media,” Appl. Opt. **55**(7), 1507 (2016). [CrossRef]

**9. **A. Pifferi, A. Torricelli, P. Taroni, D. Comelli, A. Bassi, and R. Cubeddu, “Fully automated time domain spectrometer for the absorption and scattering characterization of diffusive media,” Rev. Sci. Instrum. **78**(5), 053103 (2007). [CrossRef]

**10. **V. Ntziachristos, X. ma, A. G. Yodh, and B. Chance, “Multichannel photon counting instrument for spatially resolved near infrared spectroscopy,” Rev. Sci. Instrum. **70**(1), 193–201 (1999). [CrossRef]

**11. **S. Wojtkiewicz, A. Gerega, M. Zanoletti, A. Sudakou, D. Contini, A. Liebert, T. Durduran, and H. Dehghani, “Self-calibrating time-resolved near infrared spectroscopy,” Biomed. Opt. Express **10**(5), 2657 (2019). [CrossRef]

**12. **H. Wabnitz, D. R. Taubert, M. Mazurenka, O. Steinkellner, A. Jelzow, R. Macdonald, D. Milej, P. Sawosz, M. Kacprzak, A. Liebert, R. Cooper, J. Hebden, A. Pifferi, A. Farina, I. Bargigia, D. Contini, M. Caffini, L. Zucchelli, L. Spinelli, R. Cubeddu, and A. Torricelli, “Performance assessment of time-domain optical brain imagers, part 1: basic instrumental performance protocol,” J. Biomed. Opt. **19**(8), 086010 (2014). [CrossRef]

**13. **L. Spinelli, F. Martelli, A. Farina, A. Pifferi, A. Torricelli, R. Cubeddu, and G. Zaccanti, “Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. Time-resolved method,” Opt. Express **15**(11), 6589 (2007). [CrossRef]

**14. **S. Konugolu Venkata Sekar, A. Dalla Mora, I. Bargigia, E. Martinenghi, C. Lindner, P. Farzam, M. Pagliazzi, T. Durduran, P. Taroni, A. Pifferi, and A. Farina, “Broadband (600-1350 nm) Time-Resolved Diffuse Optical Spectrometer for Clinical Use,” IEEE J. Sel. Top. Quantum Electron. **22**(3), 406–414 (2016). [CrossRef]

**15. **V. Ntziachristos and B. Chance, “Accuracy limits in the determination of absolute optical properties using time-resolved NIR spectroscopy,” Med. Phys. **28**(6), 1115–1124 (2001). [CrossRef]

**16. **M. Giovannella, D. Contini, M. Pagliazzi, A. Pifferia, L. Spinelli, R. Erdmann, R. Donat, I. Rocchetti, M. Rehberger, N. Konig, R. Schmitt, A. Torricelli, T. Durduran, and U. M. Weigel, “BabyLux device: a diffuse optical system integrating diffuse correlation spectroscopy and time-resolved near-infrared spectroscopy for the neuromonitoring of the premature newborn brain,” Neurophotonics **6**(02), 1 (2019). [CrossRef]

**17. **R. Cubeddu, M. Musolino, A. Pifferi, P. Taroni, and G. Valentini, “Time-Resolved Reflectance: A Systematic Study for Application to The Optical Characterization of Tissues,” IEEE J. Quantum Electron. **30**(10), 2421–2430 (1994). [CrossRef]

**18. **S. Mosca, P. Lanka, N. Stone, S. Konugolu Venkata Sekar, P. Matousek, G. Valentini, and A. Pifferi, “Optical characterization of porcine tissues from various organs in the 1697–1100 nm range using time-domain diffuse spectroscopy,” Biomed. Opt. Express **11**(3), 1697 (2020). [CrossRef]

**19. **A. Corlu, R. Choe, T. Durduran, K. Lee, M. Schweiger, S. R. Arridge, E. M. C. Hillman, and A. G. Yodh, “Diffuse optical tomography with spectral constraints and wavelength optimization,” Appl. Opt. **44**(11), 2082–2093 (2005). [CrossRef]

**20. **C. Li, S. R. Grobmyer, L. Chen, Q. Zhang, L. L. Fajardo, and H. Jiang, “Multispectral diffuse optical tomography with absorption and scattering spectral constraints,” Appl. Opt. **46**(34), 8229–8236 (2007). [CrossRef]

**21. **C. D’Andrea, L. Spinelli, A. Bassi, A. Giusto, D. Contini, J. Swartling, A. Torricelli, and R. Cubeddu, “Time-resolved spectrally constrained method for the quantification of chromophore concentrations and scattering parameters in diffusing media,” Opt. Express **14**(5), 1888 (2006). [CrossRef]

**22. **R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express **16**(8), 5907 (2008). [CrossRef]

**23. **L. Spinelli, M. Botwicz, N. Zolek, M. Kacprzak, D. Milej, P. Sawosz, A. Liebert, U. Weigel, T. Durduran, F. Foschum, A. Kienle, F. Baribeau, S. Leclair, J.-P. Bouchard, I. Noiseux, P. Gallant, O. Mermut, A. Farina, A. Pifferi, A. Torricelli, R. Cubeddu, H.-C. Ho, M. Mazurenka, H. Wabnitz, K. Klauenberg, O. Bodnar, C. Elster, M. Bénazech-Lavoué, Y. Bérubé-Lauzière, F. Lesage, D. Khoptyar, A. A. Subash, S. Andersson-Engels, P. Di Ninni, F. Martelli, and G. Zaccanti, “Determination of reference values for optical properties of liquid phantoms based on Intralipid and India ink,” Biomed. Opt. Express **5**(7), 2037 (2014). [CrossRef]

**24. **K. Vishwanath and M.-A. Mycek, “Time-resolved photon migration in bi-layered tissue models,” Opt. Express **13**(19), 7466 (2005). [CrossRef]

**25. **M. Raju and S. N. Unni, “Concentration-dependent correlated scattering properties of Intralipid 20% dilutions,” Appl. Opt. **56**(4), 1157 (2017). [CrossRef]

**26. **S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star, and M. J. C. Van Gemert, “Optical Properties of Intralipid: A Phantom Medium for Light Propagation Studies,” Lasers Surg. Med. **12**(5), 510–519 (1992). [CrossRef]

**27. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties,” Appl. Opt. **28**(12), 2331 (1989). [CrossRef]

**28. **K. Vishwanath, K. Chang, D. Klein, Y. U. Feng Deng, V. Chang, J. E. Phelps, and N. Ramanujam, “Portable, fiber-based, diffuse reflection spectroscopy (DRS) systems for estimating tissue optical properties,” Appl. Spectrosc. **65**(2), 206–215 (2011). [CrossRef]

**29. **J. E. Bender, K. Vishwanath, L. K. Moore, J. Q. Brown, V. Chang, G. M. Palmer, and N. Ramanujam, “A robust monte carlo model for the extraction of biological absorption and scattering in vivo,” IEEE Trans. Biomed. Eng. **56**(4), 960–968 (2009). [CrossRef]

**30. **A. Puszka, L. Di Sieno, A. D. Mora, A. Pifferi, D. Contini, A. Planat-Chrétien, A. Koenig, G. Boso, A. Tosi, L. Hervé, and J.-M. Dinten, “Spatial resolution in depth for time-resolved diffuse optical tomography using short source-detector separations,” Biomed. Opt. Express **6**(1), 1 (2015). [CrossRef]

**31. **A. Torricelli, A. Pifferi, L. Spinelli, R. Cubeddu, F. Martelli, S. Del Bianco, and G. Zaccanti, “Time-resolved reflectance at null source-detector separation: Improving contrast and resolution in diffuse optical imaging,” Phys. Rev. Lett. **95**(7), 078101 (2005). [CrossRef]

**32. **M. Mazurenka, A. Jelzow, H. Wabnitz, D. Contini, L. Spinelli, A. Pifferi, R. Cubeddu, A. D. Mora, A. Tosi, F. Zappa, and R. Macdonald, “Non-contact time-resolved diffuse reflectance imaging at null source-detector separation,” Opt. Express **20**(1), 283 (2012). [CrossRef]

**33. **A. Liebert, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, and H. Rinneberg, “Evaluation of optical properties of highly scattering media by moments of distributions of times of flight of photons,” Appl. Opt. **42**(28), 5785 (2003). [CrossRef]