Diffuse correlation spectroscopy (DCS) is an optical technique that measures the scatterer’s motion in diffusive media, traditionally, by injecting a coherent continuous wave (CW) laser beam in a turbid sample and by characterizing the speckle fluctuations by their intensity autocorrelation function. Afterward, using a correlation diffusion model, the motion of red blood cells in the microvasculature is estimated, which has been shown to be proportional to local blood flow (BF).1 The technique is an advancement over diffuse wave spectroscopy,2 which is widely used in the study of nonbiological media. Over the last two decades, DCS has emerged as a promising technique for many clinical applications since it provides a direct measure of BF in deep () tissues.
In standard DCS, which uses a CW light source, photons that travel all the path lengths in the tissue are detected without discrimination. On the contrary, using a pulsed light source, it is possible to measure the photon time-of-flight (ToF) and, thus, to select a specific range of path lengths. By exploiting the relationship between photon ToF and mean penetration depth,3,4 depth-discrimination in the BF can be enabled, which is the first advantage of the newly emerging time-domain (TD) DCS with respect to standard (i.e., CW) DCS.5,6 A second advantage is the parallel estimation, using the same probe, of the tissue optical properties. This has been shown to increase the accuracy of the BF estimation7 and allows, by using more than one wavelength, the determination of oxy- and deoxyhemoglobin concentrations alongside BF.
The first proof-of-concept of TD-DCS was given by Yodh et al.8 In this work, by using a gating scheme based on an optical delay line and a nonlinear crystal, they showed that it is possible to measure path length resolved autocorrelation functions. However, their scheme was quite complex and needed a high optical power, unsafe for in vivo applications. Recently, a demonstration on phantoms and small animals was given by Sutin et al.,5 where a pulsed laser and a time-correlated single-photon counting (TCSPC) module were used to measure the BF with narrow (tens of picoseconds) time gates. This allowed the researchers to overcome an important barrier, namely the limited temporal coherence of a pulsed light source. Later, Pagliazzi et al.6 overcame the signal-to-noise ratio (SNR) barrier for human in vivo measurements, by using a custom-made mode-locked pulsed light source and broader (few nanoseconds) temporal gates. A hardware gating acquisition scheme, which uses fast-gated single-photon detectors enabling very short source–detector (SD) separation measurements, was also recently demonstrated.9 Recently, a theoretical model for TD-DCS in multilayer turbid media has been proposed for improving data analysis in nonhomogeneous biological tissues.10
Another method to achieve path length resolved DCS, initiated by the works of Tualle et al.,11,12 has been recently shown to be available in vivo on small animals.13 While in TD-DCS, path lengths are resolved based on their ToF, and interferometric measurements achieve the same goal by sweeping source wavelength in time. This promising technique, which can directly measure the electric field autocorrelation function, has not yet been translated to humans.
The current challenges of the technique are both on the technological side, regarding the development of effective source and detection schemes, and on the theoretical side, for the improvement of the physical understanding, necessary for data analysis. In this context, our work focuses on the study of finite temporal resolution effects in TD-DCS. These effects were very recently studied using a model that is applicable only to narrow time gates.14
Here, we argue, as in our earlier work,6 that broad gates are crucial to achieve a good SNR and, motivated by this argument, build a more complete mode. This approach would also permit movement toward an accurate depth-resolved BF measurement, and can be useful for driving hardware design. As we show, these effects have an important impact on the accuracy of the BF estimates. We also provide and validate a model for BF estimation in the case of a realistic (i.e., nonideal) TD-DCS experiment.
In the following, we derive the model, based on the instrument response function (IRF), for describing gated autocorrelation functions (Sec. 2); then, we describe the Monte Carlo (MC) simulator, the experimental system, and the phantom model exploited (Sec. 3); then the validation of the model with a numerical simulation (Sec. 4) and a phantom experiment (Sec. 5) are presented. Finally, some discussion and conclusions are reported (Sec. 6).
In classical CW-DCS, for ergodic systems and infinitely coherent sources (i.e., sources with coherence length much greater than the mean path length of photons in the sample), the normalized intensity autocorrelation function , where is the correlation delay time, can readily be related to the normalized electric field autocorrelation function , by means of the so-called Siegert relation:15,16
Here, is a constant, equal to one for an ideal experiment, that is proportional to the inverse of the number of detected modes and decreases for decreasing source coherence length. We note that halved by the use of the typical unpolarized source–detection pairs in DCS experiments.
In a TD-DCS experiment, Eq. (2) is modified to account for the selection of the path lengths in a specific range , where is the gate opening path length and is the gate width. The parameter is chosen to probe the BF at different depths: with earlier or later gate positions, superficial or deeper BF is measured, respectively. Thus, we can write as follows:
In a realistic TD-DCS experiment, it is worth noting that the injected light pulse has a finite temporal duration and the detector has a finite temporal response/resolution. When the photon ToFs are recorded (e.g., using a TCPSC module) and are subsequently time-gated via software, the system can be assumed to be linearly time-invariant and, thus, is fully characterized by measuring its IRF. IRF is typically obtained by facing the injection and the collection fibers against each other while taking care to fill all the modes. Notably, the assumption of temporal invariance may not be valid when a fast-gated detector is used for gating the photons via hardware, which does so by switching on the detector only in the desired temporal window.18,19
We note that Eq. (3) holds only in the case of an ideal system, for which , where is the Dirac delta function. In the case of an ideal system, it is assumed that it is possible to measure the photon path length with a perfect accuracy and thus select, by temporal gating of the distribution of time-of-flight (DTOF), a specific range of path lengths. On the other hand, in a real system, the path length selection is affected by inaccuracies due to the system’s nonideal temporal response. This is because the path length that a photon ideally that undergoes in the medium (i.e., the ideal path length, ) can differ from the path length that we estimate using the measured ToF (i.e., the measured path length, ) in a stochastic manner. In time-resolved spectroscopy,19 the spreading of path length is well known and affects the expression of the normalized distribution of path lengths: , where is the convolution by the IRF of the theoretical time-resolved diffuse reflectance .
In our approach, the expression of , the normalized electric field autocorrelation for a single path length, in the realistic case is calculated as a weighted average of its value in the ideal case. Assuming a fixed measured path length , the number of detected photons within an infinitesimal interval of path lengths between and can be expressed as . We, therefore, have7), assuming for normalization, we obtain our main result as follows:
The function can be interpreted as the probability, for a photon with ideal path length , to be mapped into the range or equivalently in the corresponding time gate. This function models the fact that, due to finite temporal resolution, also path lengths outside the gate can be assigned to it, contributing to the measured autocorrelation decay rate.
The measured IRF in Fig. 1 is a composition of the bell-shaped excitation pulse and the detector temporal response. The IRF width [full width at half-maximum (FWHM)] is determined mainly by the input pulse, having a duration of since the time jitter of the detector and the temporal broadening of the graded-index and single-mode fibers we used are much smaller than the width of the laser pulse. The long diffusion tail is a characteristic of the single-photon avalanche diode that was used and depends on the detector architecture.18 From Fig. 1, it is also possible to see how the ideal gating probability (equal to 1 inside the gate and 0 outside) is modified by the IRF. In fact, the EGF describes, for each time point, how much this ideal gating probability is distorted by the finite temporal resolution of the system. Due to the characteristic tail of the IRF toward positive times, many short path length photons are mapped into the late gate and contribute to its autocorrelation decay rate (see the leading edge of EGF in Fig. 1). This effect is non-negligible due to the higher number of early photons compared to late ones.3,4
We note that Eqs. (7) and (8) are applicable for any gate width and predict a distortion of the autocorrelation function also in the narrow gate limit used in Ref. 5. In fact, in this limit, by means of the well-known mean value theorem for integrals, we can write [compared to Eq. (8)] as follows:14 except for the coherence length term, which we are ignoring. Thus, our model extends the previous results to a broader range of applicability (i.e., broad gates).
Materials and Methods
Monte Carlo Simulator
To study the signals involved in a TD-DCS experiment, we have performed MC simulations using software developed at ICFO based on the MCML algorithm.20 In particular, we have simulated a semi-infinite homogeneous medium surrounded by air in the reflectance geometry by assigning the reduced scattering coefficient (), the absorption coefficient (), and the refractive index (). We have considered a collimated beam perpendicular to the surface and placed, at a fixed SD separation , a 1 mm/side cubic detector. For every detected photon packet , the path length in the medium and the arrival weight are recorded in a history file. In postprocessing, the ToF is calculated starting from the path length. To introduce the timing inaccuracy of the system, for every detected packet, we have added to the before mentioned ToF a time shift obtained by performing an MC sampling20 of a realistic IRF. First, the IRF cumulative distribution function (CDF) is computed by numerical integration, then the , such that the CDF is equal to a given random number , is selected. Finally, the gated autocorrelation function is computed, in diffusive approximation, in the following way:15
For experimental validation, we have used the custom-made Ti:Sapphire laser source, operated in active mode-locking, that was deployed in our previous works.6,9 The acousto-optic modulator permits the phase-locking of the longitudinal modes, enabling a pulse repetition frequency of 100 MHz. The wavelength was tuned to 785 nm. The light was sent to the surface of the sample using a core step-index fiber and the diffused light, recollected with a core single-mode fiber (780 HP, Nufern, East Granby, Connecticut), was delivered to a single-photon avalanche diode detector (PDM, Micro Photon Devices, Bolzano, Italy). We have implemented optical coupling of the fibers with laser and the detector, respectively, by using a combination of lenses and aperture irises, which were adjusted manually. A synchronization signal was generated using a beam-splitter and photodiode (OCF-401, Becker & Hickl, Berlin, Germany) placed at the source’s output. A TCSPC module (PicoHarp 300, PicoQuant, Berlin, Germany) recorded, for every detected photon, two quantities: the photon time stamp, so-called arrival time, and, with a resolution of 8 ps, the delay with respect to the synchronization signal, so-called pulse time. To stay within maximum sync rate of the TCSPC module, which is 84 MHz, a SYNC divider circuit (SDC, developed in our laboratory) was used to reduce the sync rate by two times to 50 MHz. The SDC operates on the photodiode electronic signal by removing the signal generated by exactly one pulse from every two pulses. The resulting TCSPC histogram is composed of two identical replicas, separated by a laser pulse period (10 ns), which are summed together in postprocessing in order not to lose useful photons. More details on the setup can also be found in Ref. 6. The IRF was acquired by facing the input and output fibers, separated by a thin polytetrafluoroethylene sheet. Its FWHM was , see Fig. 1.
Phantom Studies and Data Analysis
We have mixed a tissue-mimicking homogenous liquid phantom by adding 6.5% in mass of Lipofundin (B. Braun Melsungen AG, Germany) to water.21 The phantom optical properties were estimated with our system to be and by fitting the measured DTOF curve with the theoretical time-resolved diffuse reflectance for a semi-infinite homogeneous medium, convolved with the IRF.4 The optical coefficients are then used to compute the normalized path length distribution by normalizing the theoretical reflectance to the unitary area. We note that can also be estimated directly by deconvolving the IRF from the measured DTOF.22 This method could be an alternative to estimate the path length distribution, but a high number of counts than the used approach is necessary for a correct inversion to avoid noise in the results. After the experiment, the gated autocorrelation function is computed using the software correlator, as was used in our previous work.6 Its working principle can be divided conceptually into two parts. First, the photons with pulse time belonging to the desired temporal window are selected. Then, the arrival times of the selected photons are correlated with the fast photon autocorrelation algorithm described in Ref. 23.
To validate the proposed model, we have performed an MC simulation, as described in Sec. 3.1. We have simulated a system with the same optical properties of the phantom described in Sec. 3.3. The simulated values of the other parameters were , , and . We have launched photons and randomly divided the detected ones into 26 groups (to match the number of groups used in the experiments) to obtain a collection of autocorrelation functions for evaluating their variability with respect to the sampling noise and used for random sampling the IRF of our experimental system. The high number of photons was chosen to obtain good convergence of the simulated reflectance curves. We note that in order to simulate realistic experimental noise in the autocorrelation functions, a proper TD-DCS noise model is needed that takes the correlator structure into account in a similar manner done for CW-DCS.15
We have considered two broad temporal gates: an early gate, starting well before the rise of the DTOF to where it falls below 2% of its peak value (i.e., from to 0.9 ns, where we set at the peak of the IRF), and a late gate, from that point to where the DTOF is smaller than 0.1% of its peak value (i.e., from 0.9 to 5.9 ns). After computing the autocorrelation functions using Eq. (1), we have performed a fit for using Eqs. (3) and (7) for an uncorrected and IRF-corrected fit, respectively. The fit minimizes the mean squared error between the simulated and the model, using the MATLAB (MathWorks) function “fminsearch,” which exploits the Nelder–Mead simplex method. We have chosen a fitting region from to the point were , for the late gate, and , for the early gate. The starting and ending points of the fit are values typically used in TD-DCS experiments, to reject after-pulsing effects and regions that are effected by the probe contact, superficial heterogeneities and high noise, respectively.6,9 The different ending limit for the late gate was chosen to reject contaminations from small path length photons, which are characterized by slower decay rates and thus decorrelate at a larger .
Figure 2 shows a comparison between the theoretical autocorrelations, using the IRF-corrected and uncorrected models, and the simulated ones in the presence of the IRF, averaged over the 26 groups. In addition, a comparison between the retrieved with an uncorrected and IRF-corrected fit of the simulated autocorrelations is shown, also for the case of ungated acquisition (from to 5.9 ns).
From Fig. 2, it is possible to see that the uncorrected model does not describe the behavior of the gated autocorrelation functions correctly, in particular in the late gate case. On the other hand, the corrected model better matches the simulated autocorrelations for both early and late gates. By making use of the corrected model for a fit of the simulated autocorrelations, the estimated becomes very close to the expected value () for both the gates and also ungated acquisition. In the late gate case, the one more affected by the IRF, the error is reduced from to .
Simulations for Different Optical Properties and Source-Detector Separations
To further validate the model, the simulation has been repeated for different values of optical properties and source–detector separations. We simulated six media with different combinations of absorption and reduced scattering coefficients, as specified in Table 1. In each simulation, photons were launched and recollected at three different SD separations: , 2, and 3 cm. The simulated value of was kept to in all cases. After the simulations were launched, the gated autocorrelations (early and late) were computed using the same temporal limits of the previous section. Figure 3 shows the fitted values of the BFI for the two temporal gates and ungated acquisition, using the uncorrected and IRF model.
Simulated values of reduced scattering coefficient, absorptions coefficient, and Brownian diffusion coefficient for the six simulated media.
Using the uncorrected model, the discrepancy between early and late gate BFI is high, especially for small SD separations and high absorption coefficients, probably due to the smallest width of the path length distribution. The dependence of the error on the reduced scattering coefficient is lower. Conversely, by using the IRF-corrected model, it is possible to recover a good agreement between the gated BFIs for all the considered values of optical properties and SD separation, with an error always smaller than 3%.
We have then performed a homogenous phantom experiment, using the liquid mixture with optical properties specified in Sec. 3.3 and an SD separation . We have considered five temporal gates, in which the first one is centered at the peak of the IRF (i.e., ). Each gate had a width of 1 ns and overlapped 500 ps to the previous one. Using the measured arrival times and ToFs, for each gate, we have computed 20 autocorrelations with 1-s integration time for 20 s and averaged them, repeating this 26 times, for a total of 520 s of acquisition. The group duration was chosen to increase the SNR of the autocorrelation curves while maintaining a sufficient number of groups to obtain reliable statistics from the fitted BFI. For all the autocorrelations, the fitting interval was chosen from to the point where .
Figure 4 shows the measured autocorrelation functions averaged over all the temporal windows and the corresponding fits, uncorrected and IRF-corrected, for an early and a late gate, centered at and , respectively. We also show a comparison between the retrieved BFI using the uncorrected and IRF-corrected models, together with error-bars indicating their standard deviations along the 20 s blocks, and the average ungated value retrieved fitting the corresponding autocorrelations with the CW solution of the correlation diffusion equation.1
As it is possible to see from Fig. 4, the uncorrected and IRF-corrected models give similar fit quality in the considered regions. The obtained value of the parameter decreases gradually from 0.35 for the first gate to 0.22 for the last. The parameter changes depending on the gate position, as already reported in Refs. 5 and 6, due to the different degree of interference arising from the gating of the DTOF and the finite laser coherence length.14 The tail in the measured late gate at large may be due to larger noise induced by the lower number of photons collected at later ToFs, further study is required to clarify this point. However, the impact of the tails on the results is negligible since they were excluded by the fitting area. After the IRF correction, the discrepancy between the retrieved from the gated autocorrelations and the ungated ones is strongly reduced. For the later gates (), which show a larger distortion, the error is reduced from 37% to 3% on average. Since the medium is homogeneous, we do not expect any difference between earlier and later gates. Therefore, the IRF correction largely compensates the discrepancy in the early-late result caused by the nonideal acquisition system.
Discussion and Conclusion
In this work, we have proposed a model to describe a realistic TD-DCS experiment, characterized by a finite temporal resolution, based on measurable or known quantities, such as the IRF and the time/path length gate limits. We have introduced a function, called EGF, that quantifies the path length selection capabilities of the experimental system. In this framework, the distortions in the autocorrelation functions are explained as contaminations from path lengths in the proximity of the gate. In particular, the use of this model for IRF correction permitted to reduce the error in the retrieved BFI for the late, broad gate from 34% to 3% when MC simulations are considered. Moreover, in a homogeneous phantom experiment, the use of the proposed model reduced the BFI difference between gated and ungated autocorrelations from 37% to 3%.
Regarding the limits of validity of our method, it is worth noting that the model is applicable to any gate width since it takes into account the time interval considered for the computation of the autocorrelation function. A further improvement for our method could be to consider the finite coherence length effects as in Ref. 14 by adapting the theory described in Ref. 16. While the coherence effects are not taken care of, the results in Fig. 4 and in Sec. 7 (Fig. 5) demonstrate that as long as the laser is sufficiently coherent, as indicated by relatively high value, the model is able to use wide gates and account for the IRF effects. In Sec. 7, we show that the model is able to correctly retrieve the dynamic properties down to a pulse-width of 200 ps. Further modeling of coherence effects should improve this but it is beyond the scope of this work due to the complexity of accurate measurements of the coherence length of this versatile, tunable laser. We note that the method was checked for homogeneous media, but it may be generalized to heterogeneous samples, for instance, to a layered medium.
To conclude, the method described in this work, which can also be applied to in vivo data, paves the way to an accurate depth-resolved BF measurement. This is a critical step for many biomedical applications, where the discrimination of deeper BF from superficial BF is necessary. In addition, this work can provide useful guidance for hardware design since it permits one to understand the role of different parameters, such as detector diffusion tail, on the overall system performance.
Derivation of the Expression for the Nonideal Electric Field Autocorrelation Function
In this appendix, we report the calculations to derive Eq. (6) from Eq. (5). Starting from Eq. (5) and substituting the expressions of and , we obtain
In the above passages, we have used the relations and ; this last is valid if . Then, Eq. (12) can be rewritten as follows:
In the above passages, we have used the relations and the commutative property of convolution. Note that in the last step, we started the integration from since the ideal path length distribution is equal to zero for negative path lengths.
Coherence Length Effects on Amplitude and Decay-Rate of the Autocorrelations
To study the effects of limited coherence length in the autocorrelation functions, we have performed a phantom experiment decreasing gradually the width of the laser pulse. Through increasing the radio-frequency power of the acousto-optic modulator, we obtained five different values of IRF FWHM, from down to 150 ps. We have prepared a calibrated liquid phantom with the same receipt used in Sec. 5 and used an SD separation , measuring 30 ungated autocorrelations with 10-s integration time each, for a total acquisition time of 300 s. Figure 5 shows the resulting intensity and, using the Siegert relation, electric field autocorrelation functions. Also, we show the dependence of the measured and the estimated BFI, using the CW solution of the correlation diffusion equation, on the IRF FWHM.
Decreasing the pulse width degrades the value of as expected by the lowering of the coherence length. However, down to IRF FWHM of the order of 200 ps, the shape of the curve and the BFI estimation are not significantly distorted compared to the initial configuration (maximum IRF FWHM). Thus, for pulse widths above that threshold, neglecting the coherence length term should be a good approximation, especially for gate widths smaller than the width of the path length distribution, where is larger compared to the ungated case.
The research was funded by Fundacio CELLEX Barcelona; Ministerio de Economia y Competitividad/FEDER (PHOTODEMENTIA, DPI2015-64358-C2-1-R), Institutes de Salud Carlos III/FEDER (MEDPHOTAGE, DTS16/00087), the “Severo Ochoa” Program for Centers of Excellence in R&D (SEV-2015-0522), the Obra Social “la Caixa” Foundation (LlumMedBcn), Institució CERCA, AGAUR-Generalitat (2017 SGR 1380), LASERLAB-EUROPE IV, the European Union’s Horizon program (LUCA, 688303H2020-ICT-2015), and the EU Marie Curie Initial Training Network “OILTEBIA” (317526).
Lorenzo Colombo received his MS degree in engineering physics from Politecnico di Milano, Italy, in 2017. Currently, he is enrolled in the PhD program in physics at Politecnico di Milano. His research interests focus on photon migration in diffusive media, and the development of diffuse optical systems and methods for noninvasive monitoring of biological tissues.
Marco Pagliazzi is a PhD candidate in the medical optics group at ICFO, Spain. He graduated in physics at the University of Pisa, Italy. The focus of his research is developing diffuse optical instrumentation, algorithms, and protocols, and applying them to the study of human physiology.
Sanathana Konugolu Venkata Sekar obtained his PhD (cum laude) from Politecnico di Milano, Milan, Italy, in 2017. Currently, he is a postdoc at Biophotonics@Tyndall National Institute (IPIC), Cork, Ireland. His current area of research and interest includes upconverting nanoparticles for deep tissue imaging and optogenetics applications, GAs in scattering media absorption spectroscopy (GASMAS) for infant lung vitals monitoring, time-domain diffuse optics, diffuse correlation spectroscopy, Raman spectroscopy, imaging and nonimaging probes development, and entrepreneurial activities in biophotonics.
Davide Contini received his MS degree in electronic engineering and his PhD in physics from the Politecnico di Milano, Italy, in 2004 and 2007, respectively. He is an associate professor in the Department of Physics, Politecnico di Milano. He is the author of more than 100 papers in international peer-reviewed journals and conference proceedings. His research activity is focused on time-resolved spectroscopy of highly diffusive media for applications in biology and medicine.
Alberto Dalla Mora graduated summa cum laude in electronics engineering from the Politecnico di Milano, Italy, in 2006 and received his PhD summa cum laude in information and communication technology from the same university in 2010. Since 2011, he has been an assistant professor in the Department of Physics, Politecnico di Milano. He has authored more than 30 papers in international peer-reviewed journals. Currently, his research interests include time-resolved diffuse spectroscopy instrumentation and applications for biomedical diagnosis.
Lorenzo Spinelli graduated in physics from the University of Milan, Italy, in 1994. In 1999, he received his PhD in physics from the University of Milan. Since November 2001, he has been a researcher at the National Research Council in service at the Institute of Photonics and Nanotechnologies—Milan section. His research field has moved toward the study of photon propagation in turbid media, in order to develop optical systems for medical diagnostics.
Alessandro Torricelli received his MS degree in electronic engineering from Politecnico di Milano, Italy, in 1994 and his PhD in physics from Politecnico di Torino in 1999. He is a full professor in the Department of Physics, Politecnico di Milano. He is the author of more than 100 papers in international peer-reviewed journals. His current research interests include photon migration in diffusive media, fNIRS, and noninvasive diffuse spectroscopy with time-domain systems.
Turgut Durduran was trained at the University of Pennsylvania. In 2009, he moved to ICFO—The Institute of Photonic Sciences, Spain, where he leads the medical optics group. His research interests revolve around the use of diffuse light to noninvasively probe tissue function. The group develops new technologies and algorithms and routinely translates them to preclinical, clinical, and industrial applications.
Antonio Pifferi received his MS degree in nuclear engineering from Politecnico di Milano, Italy, in 1991 and his PhD in physics from Politecnico di Torino in 1995. He is a full professor in the Department of Physics, Politecnico di Milano. His research is directed toward the development of laser techniques and instrumentation for diagnosis and the study of light propagation in diffusive media, with applications to optical biopsy, optical mammography, and functional brain imaging.