Time domain diffuse correlation spectroscopy with a high coherence pulsed source: in vivo and phantom results

Diffuse correlation spectroscopy (DCS), combined with time-resolved reflectance spectroscopy (TRS) or frequency domain spectroscopy, aims at path length (i.e. depth) resolved, non-invasive and simultaneous assessment of tissue composition and blood flow. However, while TRS provides a path length resolved data, the standard DCS does not. Recently, a time domain DCS experiment showed path length resolved measurements for improved quantification with respect to classical DCS, but was limited to phantoms and small animal studies. Here, we demonstrate time domain DCS for in vivo studies on the adult forehead and the arm. We achieve path length resolved DCS by means of an actively modelocked Ti:Sapphire laser that allows high coherence pulses, thus enabling adequate signal-tonoise ratio in relatively fast (~1 s) temporal resolution. This work paves the way to the translation of this approach to practical in vivo use. © 2017 Optical Society of America OCIS codes:(170.5280) Photon migration; (170.6920) Time-resolved imaging; (170.3890) Medical optics


Introduction
Time domain near infrared spectroscopy, also known as time-resolved reflectance spectroscopy (TRS), measures the distribution of photon path lengths (or photon time-offlights) in tissue to characterize the probed volume in terms of its wavelength dependent reduced scattering ( ′) and absorption ( ) coefficients [1,2]. By collecting this information at various wavelengths, one can estimate the concentration of chromophores like water, oxyand deoxy-hemoglobin, lipids and collagen, and give information about tissue structure at the cellular level [2]. Furthermore, TRS allows the separation of early and late photons. Since the probability distribution for the maximum penetration depth moves towards deeper layers when photon time-of-flight increases [3], the contributions from deeper tissue volumes can be better separated from superficial ones if the longer path length (longer time-of-flight) photons are analysed. This ability to path length-resolve is crucial for the clinical success of this technique for many scenarios such as neuro-monitoring and neuro-imaging due to overlaying extracerebral tissue contamination of the signals.
Another emerging near-infrared technique is diffuse correlation spectroscopy (DCS). DCS has been proven to be a reliable tool for the assessment of blood flow in the microvasculature in a non-invasive way [4][5][6]. DCS exploits the speckle fluctuations of coherent light that underwent multiple scattering and quantifies its statistics by the measurement of the intensity auto-correlation function. This can then be used to estimate the electric field auto-correlation function [4][5][6]. In a diffusive medium, the correlation diffusion equation describes the propagation of the electric field temporal auto-correlation function through tissues and allows us to quantitatively relate its temporal decay to the motion of red blood cells and, therefore, to the blood flow. Traditional DCS uses continuous wave (CW) lasers and it does not readily contain a way to separate or gate short and long path lengths in a similar manner as TRS.
The extension of DCS to path length-resolved DCS (i.e. time domain DCS) is non-trivial. This is mainly due to the physics of pulsed light, which implies a limited coherence length for a short (about 100 ps) pulse. In other words, when the coherence time is much shorter than the distribution of photon time-of-flight in tissue, the bulk of the speckles decorrelate due to this long path length and random scattering, which washes out the contributions due to the motion of the scatterers.
Early on, Yodh et al. [7] demonstrated that by splitting the input pulse and using one path as a reference with a controlled delay, it is possible to carry out time-or path length-resolved DCS (or diffuse wave spectroscopy, as it is known in other fields). The crucial difference of this method in comparison to other path length-resolved laser speckle methods is that it is applicable for the highly multiple scattering regime. Their approach used non-linear effects and, therefore, it is unsuitable for in vivo measurements due to the need for a large amount of detected light.
Recently, Sutin et al. [8] were able to extend the concept to a relevant domain. They achieved this by using a pulsed light on a tissue phantom and, during the analysis, by considering only the photons that have travelled paths of similar length. This allowed them to relax the requirement for high coherence length in comparison to the case in which all the photon paths are considered. They have demonstrated that they are able to separate the contribution of photons coming from shallow and deep layers, that is, those arriving at early or late time gates. These experiments were carried out on phantoms and on small animals. One drawback of this approach is that the use of narrow time gates to calculate the autocorrelation limits the signal-to-noise ratio due to the scarcity of photons within a certain gate; this in turn limits its applicability to in vivo experiments on humans.
In this wo allow the use to-noise barri arm in a cont approach as c domain appro photons with separate cont approach the heterogeneity conceptually position while [11]. The use used with the Besides, time distribution o scattering from this approach and DCS, thu technique in superficial blo temporal reso dynamic in vi use.    [3,9]. Theref layers, and as d in a way tha icroscopic leve h resolution: knowing the ti separations is c to provide a f s an elegant tons in tissue terers, for bette measurement ture of blood the adult brai he changes du with a good e way to the tra up for the time do

Time dom
The source w by means of ngitudinal mo 4]. It is pum se repetition ra ence pulsed so w ns) time gat urements from advantages of tivated us in th e -which can fore, a single opposed to a at is less sensit el [10,11] We have operated the laser at a wavelength ( ) of 785 nm. We have used a stack of three metallic variable attenuators (Thorlabs, Germany) to ensure that the illumination is below the maximum-permissible skin exposure limits. We have achieved the synchronization of the photon detection with the laser source pulses by splitting off a small fraction (<5%) of its light to be sent to a photodiode (Becker & Hickl, Germany). The remaining laser light was then launched into a 200 μm core step index fiber that delivered light perpendicularly to the sample. At a variable distance, depending on the experiment, a 5 μm core single mode fiber collected the diffuse light and delivered it to a single-photon avalanche diode (SPAD) detector (PDM, Micro Photon Devices, Italy) with a photon detection efficiency of about 20% at the wavelength of interest. This detector has two outputs. A high temporal resolution nuclear instrumentation module (NIM) output (nominally, with about 35 ps precision) was connected to a standard time correlated single photon counting (TCSPC) module (Pico Harp 300, PicoQuant, Berlin, Germany). The TCSPC module records the time stamps (time of arrival) and the delay with respect to the pulse emission (pulse time), with an 8 ps resolution, of each detected photon in a file. The coarser temporal resolution transistor-transistor logic (TTL) output (with few hundred ps precision) of the SPAD detector was connected to a standard hardware correlator (Correlator.com Flex-05, 8 channels, USA). The hardware correlator evaluates the intensity auto-correlation function assuming that the light is continuous wave and was used in a subset of experiments (see below).
We note that the TCSPC board can accept synchronization signals up to a maximum 84 MHz but our laser repetition rate was higher. To deal with this, we have developed an inhouse fast circuit (SYNC divider) that reduces the number of synchronization pulses by accepting just one in three, thus recording three distribution of times of flight (DTOF) curves in the same window. To avoid any loss of photon counts, these curves are combined in data analysis. Furthermore, in some experiments, by switching the pulsed source to a continuous wave high coherence diode laser (Crystal Laser, Reno, NV, US), we have tested the scatterers dynamics of the phantoms in an equivalent manner to a standard DCS setup.
The instrument response function (IRF) was acquired by facing the source and the detector fibers directly in front of each other and by introducing a sub-millimeter thick polytetrafluoroethylene sheet in between. The full width at half maximum (FWHM) of the IRF was about 100 ps.
The FoCuS-point software correlator [15], that is available on Github [16], was used for the post-processing of the data. It relies on the algorithm for fast photon auto-correlation from Wahl et al. [17]. We have adapted it to correlate the time stamps of the detected photons that were recorded by the TCSPC within a time window or a gate. This allowed us to gather the photons that arrive in an adjustable time gate and calculate the gated auto-correlation function by using time stamps of just a portion -gate-of the DTOF at a certain delay with respect to the laser emission, offline, in a similar fashion to the work of Sutin et al. [8].

Data analysis
The Siegert relation, ( ) = 1 + | ( )| , when the assumptions behind it are fulfilled [18], can be applied in order to relate the measured normalized intensity auto-correlation function ( ( )) to the normalized electric field auto-correlation function ( ( )), in order to quantify the motion of scatterers. Here is the correlation delay time and is a constant, determined by the collection optics. is equal to 1 for ideal experimental conditions but decays with an approximately inverse proportionality to the number of modes detected at the same time. It also decreases by decreasing the source coherence length. As first proposed by Pine et al. [18] for diffusing wave spectroscopy, the theoretical model for ( ) can be computed by evaluating the integral Where ( )is the probability distribution of the photon path lengths ( ) in tissue and is a quantity that describes the decay rate of the auto-correlation function. The integration over the path lengths is infinite for an ungated experiment, but, in gated experiments, just the path lengths between two extremes of the gate, and , are used. The other term in the integrand, (2) is the contribution to the correlation decay from a single path (of a multiply scattered photon) of length . The decay rate is given by the formula is the square of the wavenumber of the light in the medium and is the fraction of moving to static scatterers [4]. The decay rate, therefore, represents the rate of loss of correlation per unit path length in the medium and is proportional to the scatterer Brownian diffusion coefficient ( ) of the sample. By applying the above-mentioned Siegert relation we are left with two parameters to be estimated, namely the optical parameter and either (in phantoms, = 1) or the blood flow index (BFI = ). The BFI is assumed to be related to the shear-induced Brownian motion of red blood cells [20] in case of in vivo measurements, and was shown to be proportional to both absolute and relative microvascular blood flow [4][5][6]21]. Eq.
(1) tells us that, if the gate is sufficiently narrow compared to the decay of ( ), then ( ) could be treated as a constant and simplified. Afterwards, the average D or the BFI of a collection of narrow adjacent gates can be recovered from the slope of the decay of the autocorrelation function at different gate positions or path lengths (s) [7,8]. However, for broader gates, this approximation is not anymore valid and further steps need to be taken. In particular, in order to solve Eq. (1) numerically, ( ) must be estimated. As we assume an uniform refractive index in the sample, the path length is the time-of-flight of the photons times the speed of light in the sample and ( ) can be rewritten, in the limits of validity of the photon diffusion approximation, as the time-resolved reflectance [22,23] for the relevant boundary conditions. For all our experiments, we have assumed a semi-infinite medium. In order to correctly recover the path lengths ( ) from the photon time-of-flight, we have set the zero time as the time at which the IRF has a peak. In principle, the reflectance could also be estimated from the measured DTOFs and IRF by the deconvolution of the latter from the former [21] and by using a proper regularization to reduce the deconvolution noise. In our case, we have chosen to compute ( ) from theory using the absorption and reduced scattering coefficients that have been measured by knowing the DTOF and IRF using a semiinfinite medium model with extrapolated boundary by convolving the theoretical model with the IRF [22].
Here, we have tested both the narrow (160 ps) and broad gate approaches on tissue simulating phantoms. The broad gates are used for fast, light-starved in vivo measurements of blood flow in humans. We chose the width of the broad gates depending on the shape of the DTOF. We have considered the photons with times of flight up to a certain fraction of the DTOF in an early gate and the rest in a late gate, for a total of two broad gates. The width of the narrow gates was chosen as a compromise between the detected photon counts within the gates and β value. For the broad gates, the early gate extended from the time at which the IRF peaked (t ) to the time at which the measured DTOF dropped to 90% on the falling edge of the curve, with respect to its maximum, for a total width of about 800 ps. The late gate extends from this point to the time at which the DTOF is no longer distinguishable from the noise floor which results in a width of between 2 and 3 ns. Then, we numerically solve the full expression in Eq.
(1) to model the auto-correlation.  We have recovered the optical properties ( ′and ) of the liquid phantoms by fitting the convolution of the IRF and the solution of the radiative transport equation in the diffusion approximation with the corresponding DTOF. The DTOF was acquired for each phantom separately in a homogeneous one-litre container.

2.3Phantom
As discussed in section 2.2, we have first computed the electric field auto-correlation of narrow, 160 ps wide gates, overlapping for 80 ps from one to the next. The normalized intensity auto-correlation function for each gate was then fitted against the delay time with the single path length exponential decaying correlation (Eq. 2) to obtain • . For each gate, is the average path length. We have plotted the fitted • against and calculated the slope, which is , in two regions, separately: early, between 80% of maximum counts on the rising edge of the DTOF and 90% on the falling edge, and late, between 60% and 10% on the falling edge of the DTOF.

In vivo experiments
For all the in vivo studies, the laser power was limited by using an optical attenuator to meet the maximum permissible exposure limit for human skin. The protocol was approved by the Ethical Committee of Politecnico di Milano and it was conducted in agreement with the Declaration of Helsinki. All subjects gave written consent before their participation.
First, we have taken measurements on the forehead of a healthy subject (male, 30 years old). He was instructed to lay down supine on a bed and remain alert for the duration of the experiment. The probe was a black plastic fiber holder (5 cm by 5 cm by 1 cm) with the fibers inserted 1 cm apart and pointing perpendicularly to the skin. A cotton mesh wrap encircling the head was used to secure the probe to the right side of the forehead far from the midline. We have repeated the experiment two times. The second time we have manually applied a constant increased pressure (the maximum the subject was comfortable with for all the duration of the data acquisition). We have acquired, for each pressure scenario (no pressure, pressure), 2.7•10 7 photons over 300 s, and we have used the broad gate approach for the analysis.
Second, we have used the same source-detector separation (1cm) for measuring the BFI in the arm of another healthy subject (male, 30 years old) before, during and after an arterial cuff occlusion. The black plastic fiber holder was applied, with the fibers perpendicular to the skin, above the brachioradialis muscle in the right arm and was secured with a loose bandage. The tourniquet for compressing the affluent arteries was applied around the arm right below the shoulder joint. The cuff occlusion started at 245 s and lasted for 180 s. During the cuff occlusion the tourniquet pressure was raised to about 180 mmHg which is well above the arterial blood pressure of the subject. In order to gain temporal resolution, we have divided the 600 s measurement period in 1 s epochs and analysed them separately. The BFI was then estimated by using the numerical integration method in the early and late gates separately for each epoch. At the same time, we have acquired correlation curves independently on the time of arrival of the photons with 1 s time resolution by using the hardware correlator. These were analysed with the standard continuous wave DCS theory and the results were compared with the gated analysis.

Results
The first step was to fit the measured DTOF (Fig. 3, top) from both liquid mixtures (see section 2.1) to determine the optical properties. The reduced scattering coefficients were 10.6 cm -1 for the phantom with 30% glycerol and 10.5 cm -1 for the phantom without glycerol. Absorption coefficients were 0.026 cm -1 and 0.025 cm -1 , respectively. This confirms that the contrast between the two compartments would be due to the differences in viscosity and not due to differences in the reduced scattering coefficient, which is tightly coupled to in the correlation diffusion model. Fig. 3, bottom, shows the results for the dynamically heterogeneous and dynamically homogeneous two-layer phantom. As expected, both regions (early and late) reveal the same decay constant, k, and varying the thickness of the superficial layer (S) does not change k in the dynamically homogeneous case (DH, triangles). All the measured auto-correlation curves are reported, for completeness, in Fig. A1 in the Appendix, for the DH, δ=5 mm case. As expected (Section 2.2), by increasing the average gate path length, the decay constant, k·s, of the auto-correlation curves increases. The small distortions of the early gates with respect to the late in the DH case are due to the contributions of instrument response function, the contributions of the photons that underwent a small number of scattering events and the finite gate width. The effects of these contributions are different depending on the position of the center of the gate with respect to the DTOF curve, which is a topic for future study. The value of becomes significantly different between the early and late regions when we have filled the superficial (S) layer with the reduced phantom with 30% glycerol (DA, circles). Notably, in this case the retrieved is lower (see Table 1) in the early compared to the late region, and decays faster when the thickness of layer S is increased (52%, 36% and 28% of the dynamically homogeneous case in early region; 100%, 71% and 50% in late region).  with the nsity autoble 1, just ing much relative to %). As we increase the superficial layer thickness δ, the estimated relative decrease of is practically independent on the evaluation methods (numeric integration or k·s-fitting).
The results from the adult-human head are shown in Fig. 4. In Fig. 4 (a), we report the measured intensity auto-correlation ( ) functions for the early and late broad gates and for the ungated case. It is evident how, even in the pulsed case, we have achieved a good signalto-noise ratio for the estimated correlation curves, with =0.41 and 0.25, respectively for the early and late gates, and =0.29 for the ungated case. We note that the values of observed in these experiments are higher compared with the ones reported in Sutin et al. [8], as a result of the tunable active pulse locking that our particular laser allowed. Fig. 4 (b) shows the measured electric field auto-correlation function ( ) at early gates for the baseline (control) and pressure scenarios, demonstrating that the probe pressure causes a significant change in the measurements from the early gate. The numeric fit was carried out considering the experimental points up to g ~ 0.3. We report the auto-correlations of the other gates in Fig. A2, in the Appendix. Fig. 4(c) shows the measured BFI normalized relative to the ungated value, which reflects a mixed contribution of superficial and deep photons, of the baseline (no pressure) condition. As expected [27], when pressure was applied, the ungated analysis showed a decrease of 61% of BFI. Since the blood flow in the scalp is reduced, the contrast between the early and late gates increases when pressure is applied to the probe. The results of the cuff occlusion experiment are summarized in Fig. 5. The measured average was 0.36 and 0.32 for the early and late broad gates, respectively. This set-up is able to achieve a good signal-to-noise ratio (count-rate and total number of photons) in a 1 s window, and we are able to discern the appropriate decay during the cuff occlusion, followed by the rapid and strong hyperemia. The features of the three time-traces are very similar. Trends in the estimated relative BFI (rBFI) with the broad gates analysis are in good agreement with the ungated estimation carried out with the standard hardware correlator.

Discussion and conclusion
We have demonstrated the feasibility of separating the Brownian diffusion coefficients of scatterers in deep and shallow layers by using a mode-locked pulsed laser source. To do so, we have used a two-layer phantom in which the two layers have a different Brownian diffusion coefficient, but same reduced scattering and absorption coefficients. When the Brownian coefficient of the superficial layer decreases, we have observed a selective drop in the typical decorrelation time corresponding to the photons arriving early to the detector. We have shown a contrast in the relative dynamics of the two layers with respect to the ungated case in both phantoms and in vivo. Even though we do not report all possible combinations of changes of dynamics and optical properties in a two layer media, this work paves the way to works of that nature by establishing the proper experimental and theoretical framework.
We have also demonstrated the feasibility of applying this technology directly in vivo on humans. This can be used to discern the blood flow of a superficial layer (e.g. scalp and skull) with respect to deeper one (e.g. the brain), in order to reduce partial volume artefacts that affect standard DCS measurements. Here, we have focused on two broad gates. However, this method can be extended to the use of an analysis that uses multiple gates and/or a full set of autocorrelation curves at different time bins. For example, Spinelli et al. [28] utilized a similar method for TRS. This could be a future extension beyond this work where the relevant physical and algorithmic aspects could be studied.
On the head of a healthy subject we have observed a reduction in the BFI corresponding to shallow tissue as we apply additional pressure. The additional pressure we have applied does not affect the deeper brain tissue, which is protected by the hard skull bone, but it causes a decrease of the blood flow in the shallow skin and scalp layers. This is in agreement also with previous studies that demonstrated that an externally applied pressure modulates the contribution of the extra cerebral layer on the BFI assessment by DCS [27,29]. In these studies, different source-detector separations were utilized to gain a pseudo-depth resolution, but here we do take advantage of the multiple gate strategy to achieve depth resolution with a single-shot measurement. We note that we have used a relatively short source detector separation (1 cm). The depth resolution is in fact achieved by the time gate only, in a way similar to the null source detector separation TRS [30]. In a similar approach, a gated detector can be used here to systematically discard the very early photons that would blind the detection electronics at very small source detector separation. This concept can be used in a multi-source, multi-detector system to increase the resolution of diffuse correlation tomography.
In the arm cuff experiment we show that the broad time gates analysis is able to achieve a good signal-to-noise ratio with 1 s resolution, and therefore it is suitable to follow fast changes of the BFI. Surprisingly, we see that the early and late gates do not differ in terms of the BFI values, and there is no statistically significant difference between the two gates, or between the ungated and gated results, during all three phases considered separately. Although we have expected a higher hyperemic peak during reperfusion in the late (deep) gate [31], due to higher oxygen consumption in the muscle, with respect to the shallower skin and fat layers, we were not able to detect it. We believe that this is due to the very thin superficial tissue (skin and fat) thickness of 2.5 mm, measured by a skin caliper, of this particular subject. Nevertheless, we were able to prove that even fast BFI-changes due to challenges like the cuff occlusion can be monitored by time domain DCS. The results are in agreement with the ones obtained with a standard hardware correlator for ungated measurements and the same pulsed light source.
In all the experiments we have used a single wavelength but, by tuning the laser wavelength sufficiently fast (i.e. every 500ms), we can gain specific sensitivity to the deoxygenated or oxygenated hemoglobin or even further, to characterize the sample multispectrally while also measuring the blood flow. We note also that multiple source-detector pairs can readily be utilized by this method by multiplexing.
In this work, we have assumed homogeneous optical properties and we have focused on demonstrating the ability to separate the dynamic properties of the superficial and deeper layers. However, in vivo measurements contain heterogeneities in both optical and dynamic properties. It is well known that the inaccuracies in the recovery of the optical properties are reflected in the recovery of dynamic properties [32]. However, since time domain DCS measurements contain the TRS information, we believe that the quality of the results will be similar as the quality of the TRS results [33] and this is a topic for future work.
While, as other works suggest [18,34], the Siegert relation may not be valid when the coherence length of the source is in the order of the distribution of photon path lengths in tissue, or even shorter, here we demonstrate that this relation still allows us to extract meaningful values of BFI from in vivo experiments. Sutin et al. [8] have justified the application of the Siegert relation to the use of very narrow gates, and we believe that the greater coherence length of our system allows us to avoid most distortions when working with broader (up to 1 ns) gates. Moreover, no evidence of distortion that may be caused by the non-validity of the Siegert relation, such as the intensity auto-correlation curves not decaying to one at infinite time, were observed.
The longer coherence length of the pulsed laser source that we have used here, which is manifest in the relatively high values of observed, allows for a time resolution that is high enough to monitor even fast hemodynamic changes like those due to the cuff occlusion in the arm. However, a dedicated setup like the one presented in this work offers several disadvantages. First, the need for a long cavity and a powerful pump laser makes the device bulky and hard to use at the bedside. Second, the fact that this is a unique system -it is a modified system adapted to operate in active mode-locking -does not allow one to easily fabricate multiple systems. The present system is a powerful laboratory workstation, suitable for exploration of basic physics of time-domain DCS and of its potential. However, the use of similar laser sources has been common in the past for human in vivo investigation [35,36] which encourages us, in order to move towards the clinics, to work on making the instrumentation more compact, by substituting the custom Ti:Sapphire laser with a modification of a commercially available laser that will provide a suitable pulse coherence and power. For the broad gate approach applied in vivo in the head (see Fig. 4) we report the early and late gate estimated auto-correlation curves (Fig. A2). This way we intend to highlight the full picture of observed changes when pressure is applied, for both gates. Here we show two gated (early, late) and ungated calculations for different applied pressure conditions: no pressure (circles) and with pressure (diamonds). The solutions to the correlation diffusion equation was used to fit the ungated and the numeric integration was used to fit the gated results (solid lines).

Disclosures
Turgut Durduran is an inventor on a relevant patent (Patent US8082015B2, "Optical measurement of tissue blood flow, hemodynamics and oxygenation"). ICFO has equity ownership in the spin-off company HemoPhotonics S.L. Potential financial conflicts of interest and objectivity of research have been monitored by ICFO's Knowledge & Technology Transfer Department. Rainer Erdmann is CEO of Picoquant GmbH. The authors declare that there are no conflicts of interest related to this article.