Diffuse correlation spectroscopy (DCS) is a well-established method that measures rapid changes in scattered coherent light to identify blood flow and functional dynamics within a tissue. While its sensitivity to minute scatterer displacements leads to a number of unique advantages, conventional DCS systems become photon-limited when attempting to probe deep into the tissue, which leads to long measurement windows (∽1 sec). Here, we present a high-sensitivity DCS system with 1024 parallel detection channels integrated within a single-photon avalanche diode array and demonstrate the ability to detect mm-scale perturbations up to 1 cm deep within a tissue-like phantom at up to a 33 Hz sampling rate. We also show that this highly parallelized strategy can measure the human pulse at high fidelity and detect behaviorally induced physiological variations from above the human prefrontal cortex. By greatly improving the detection sensitivity and speed, highly parallelized DCS opens up new experiments for high-speed biological signal measurement.

Various optical technologies have been developed to noninvasively detect dynamic biological events in deep tissues, such as changes in blood flow and other related fluctuations. For instance, methods such as diffuse optical tomography (DOT),1 diffuse optical spectroscopy (DOS),2 and functional near-infrared spectroscopy (fNIRs)3 all illuminate a tissue’s surface and measure the resulting scattered light at different locations to reconstruct images of absorption and scattering properties. Recent extensions now utilize time-gated measurements to probe even deeper by isolating scattering from non-superficial layers.4 By examining scattered light’s spectral interference fringe pattern,5 the related technique of interferometric near-infrared spectroscopy (iNIRS) can also extract additional path-length-dependent information to accurately measure deep blood flow.6 The above techniques are distinct from other measurement strategies such as Doppler OCT,7 OCT angiography,8 and laser speckle contrast imaging;9 in that their primary aim is to penetrate multiple millimeters into the tissue to extract useful signals.

Another well-established method to measure such deep tissue dynamics, termed diffuse correlation spectroscopy (DCS), utilizes the relatively simple strategy of illuminating the tissue with highly coherent light and collecting back-scattered interference using a photon-counting detector. Dynamic events, such as blood flow or other small movements, are then estimated by calculating and processing the scattered light’s temporal autocorrelation curve.10–12 DCS offers a non-invasive phase-sensitive measurement strategy13,14 to detect blood flow and dynamic changes within the tissue without requiring a reference beam and thus, offers a promising means to achieve a simple and low-cost functional imaging technique with high sensitivity. DCS has been actively studied and translated into a useful clinical technique for analysis of the brain, breast, and skin.15–19 

DCS systems are typically implemented in a reflection geometry, where a source and a detector are placed at a finite distance d apart and detected light is assumed to travel through a “banana-shaped” point-spread function (PSF) that describes the light’s most probable path through the scattering tissue10 (Fig. 1). As the mean penetration depth of this PSF typically falls between one-third and one-half of the source–detector separation d11 for typical measurement arrangements, large separation values d are generally required to probe deep within the tissue. However, due to the forward-scattering nature of a biological tissue, most photons do not reach a detector in standard DCS measurement geometries, especially when the source and detector are placed at a large distance apart to probe deep events. For example, when illuminating the strongly scattering tissue with 200 mW/cm2 visible light (at the ANSI safety limit20), a straightforward simulation21 specifies that the average number of photons obtained from a speckle grain in each exposure time is much less than 1 at a 2 cm source–detector separation and at a sampling rate of 1 MHz (see supplementary material 1). Probing variations that arise beneath the human skull requires large, multi-cm source–detector separations,10,11,22,23 indicating that current DCS systems are severely photon-limited when outfitted to measure brain dynamics. To partially overcome this low photon budget, current DCS systems must integrate measured signals over a long time period—often more than a second—reducing overall system temporal resolution.14,24,25 While DCS can, in principle, capture the effects of very small changes within the brain, such as cerebral blood flow (CBF),26 single-detector implementations of DCS quickly run into an issue of not having enough photons when aiming to detect deep tissue events.

FIG. 1.

Schematic diagram of parallelized DCS detection. The light from a 670 nm long-coherence laser is delivered by a multimode fiber to the tissue surface. Photons that penetrated deep enough will be reflected by the moving red blood cells and other scatterers. A small fraction will scatter back to the detection multimode fiber and are directed to the SPAD array. The distance between the center of the illumination fiber tip and detection fiber tip is d. The SPAD array records the temporal intensity fluctuations of the speckles at a sampling rate of 333 kHz across all 1024 channels. The autocorrelation curves g2 of each SPAD are computed and combined. Perturbations induced by blood flow are detected by comparing the decorrelation rate of the final averaged result.

FIG. 1.

Schematic diagram of parallelized DCS detection. The light from a 670 nm long-coherence laser is delivered by a multimode fiber to the tissue surface. Photons that penetrated deep enough will be reflected by the moving red blood cells and other scatterers. A small fraction will scatter back to the detection multimode fiber and are directed to the SPAD array. The distance between the center of the illumination fiber tip and detection fiber tip is d. The SPAD array records the temporal intensity fluctuations of the speckles at a sampling rate of 333 kHz across all 1024 channels. The autocorrelation curves g2 of each SPAD are computed and combined. Perturbations induced by blood flow are detected by comparing the decorrelation rate of the final averaged result.

Close modal

A clear solution to this problem is to simply detect light from more optical modes; measuring light from N independent speckle grains simultaneously will yield a proportional increase in the detected signal and a N increase in noise (under the assumption of shot-noise), which implies a N gain in the overall signal-to-noise ratio (SNR). Since the optical modes emerging from the tissue are nearly mutually incoherent (i.e., decorrelating independently), sampling multiple speckle grains on a single detector reduces measurement contrast and is not a helpful strategy (beyond integrating over a few modes per channel27,28). DCS, thus, inherently benefits from a parallelized detection strategy, wherein multiple independent detectors could be used to measure N fluctuating speckle grains simultaneously, resulting in an increase in the signal-to-noise ratio (SNR) by a factor of N28,29 (see supplementary material 2). Several multi-channel DCS systems have previously been developed to simultaneously detect more signals. For example, a number of four-channel DCS systems have been reported in prior work to detect blood flow both in humans30,31 and mice;32 Staples et al. constructed an eight-channel DCS system using one PMT and seven single photon counting modules,33 Johansson et al. developed a 5 SPAD × 5 SPAD configuration to demonstrate an improved SNR on a milk phantom and an in vivo blood occlusion test,27 and Dietsche et al. implemented a 28-channel DCS setup with multiple single-mode fibers to measure blood flow in the deep tissue.28 The DCS system in the last work has the most parallelized detection channels that we are aware of to date.

The recent development of integrated SPAD array technology, wherein thousands or more individual SPAD detectors can be integrated on a single CMOS chip, opens up a new regime for massively parallel single-photon detection, without dramatically increasing the complexity of the DCS measurement setup. SPAD arrays are solid-state detectors now fabricated in standard CMOS technology that, when manufactured at scale, offer a low-cost photon counting solution with unrivaled performance via chip or pixel-level processing capabilities.34 SPAD arrays have recently been implemented in many biological imaging experiments to achieve a wide variety of interesting tasks. These include depth ranging and imaging,35 fluorescent lifetime imaging microscopy,36–39 endomicroscopy,40 spectroscopy,41 Raman spectroscopy,42–44 ultrafast optical imaging,45 and localization-based super-resolution imaging,46 to name a few.

In this work, we use an integrated SPAD array that contains 1024 independent detectors to demonstrate a massively parallelized DCS system with ∼40 times more channels than any prior work (to the best of our knowledge). We discuss a unique software pipeline to integrate these parallelized measurements and quantitatively demonstrate the sensitivity and temporal resolution enhancements of parallelized detection with a novel phantom that consists of a digital micromirror device (DMD) placed directly behind a rapidly decorrelating tissue-mimicking phantom (Fig. 2). The DMD is employed to generate high-speed temporal perturbations with known parameters, which allows us to carefully measure and assess system performance. This phantom allowed us to determine our parallelized DCS system’s sensitivity (i.e., its smallest detectable perturbation size) through thick tissue phantoms under various conditions (e.g., average scattering mean-free path, source-to-detector separation, perturbation speed, and signal integration time). For example, we demonstrate the ability to detect a 2.7 × 2.7 mm2 perturbation that lasts 0.1 s under 1 cm of a tissue-like phantom (μa = 0.07 cm−1 and μs′ = 3.6 cm−1) with an 80% accuracy. We also verify that our system can capture a high fidelity in vivo blood flow signal from the adult forehead at a 30 ms temporal resolution and successfully detect the change in decorrelation time from above the human prefrontal cortex at a similar rate when participants switched between the behavioral tasks of resting vs reading text.

FIG. 2.

Detailed diagram of our parallelized DCS phantom setup for sensitivity and temporal resolution characterization. (a) Images of the tissue phantom consisting of polystyrene microsphere-in-water suspension in a customized cuvette; (b) the DMD phantom setup with DMD placed immediately behind the tissue phantom. The distance between the centers of the illumination fiber and the detection fiber d was set to be approximately twice the cuvette thickness; (c) the DMD panel. Square colored boxes represent different perturbation areas, varying from 1.37 mm2 × 1.37 mm2 to 9.58 mm2 × 9.58 mm2; (d) example dynamic variation patterns, where pixels in the middle square area switch between “on” and “off” states at a rate of multiple kHz [corresponding to the square color box in (c)], while the peripheral area is a constant random pattern; (e) the schematic layout of the 32 SPAD × 32 SPAD array illustrated with representative detected speckle patterns. The average speckle size at the SPAD array was tuned by adjusting the detection fiber distance to approximately match the SPAD pixel active area (black circle). (f) Representative raw intensity measurements, Ij(t), from each SPAD pixel with a temporal resolution of 3 µs. (g) Corresponding intensity autocorrelation curves of each SPAD pixel, g2(j)(τ), integrated over 30 ms and (h) the final averaged autocorrelation curve g2(τ) using all 1024 SPADs.

FIG. 2.

Detailed diagram of our parallelized DCS phantom setup for sensitivity and temporal resolution characterization. (a) Images of the tissue phantom consisting of polystyrene microsphere-in-water suspension in a customized cuvette; (b) the DMD phantom setup with DMD placed immediately behind the tissue phantom. The distance between the centers of the illumination fiber and the detection fiber d was set to be approximately twice the cuvette thickness; (c) the DMD panel. Square colored boxes represent different perturbation areas, varying from 1.37 mm2 × 1.37 mm2 to 9.58 mm2 × 9.58 mm2; (d) example dynamic variation patterns, where pixels in the middle square area switch between “on” and “off” states at a rate of multiple kHz [corresponding to the square color box in (c)], while the peripheral area is a constant random pattern; (e) the schematic layout of the 32 SPAD × 32 SPAD array illustrated with representative detected speckle patterns. The average speckle size at the SPAD array was tuned by adjusting the detection fiber distance to approximately match the SPAD pixel active area (black circle). (f) Representative raw intensity measurements, Ij(t), from each SPAD pixel with a temporal resolution of 3 µs. (g) Corresponding intensity autocorrelation curves of each SPAD pixel, g2(j)(τ), integrated over 30 ms and (h) the final averaged autocorrelation curve g2(τ) using all 1024 SPADs.

Close modal

To provide a preliminary assessment of highly parallelized DCS based upon an integrated SPAD array; we constructed an easily reconfigurable phantom setup to create realistic kHz-rate perturbations deep beneath a decorrelating scattering slab. To achieve this goal, we turned to a non-biological phantom arrangement that offered a useful testing platform to quantitatively evaluate our new system’s performance gains. Multiple tissue-mimicking phantoms have previously been developed to simulate the dynamic scattering and absorption properties of human tissue such as that from the brain and skin47–52 (made of intralipid,49 glycerol with fat emulsion,48 or silicon rubber with small embedded scatterers,50,52 for example). The dynamic optical properties of these prior phantoms have been typically tuned by varying the phantom material’s temperature, viscosity, or with pump-induced motion.48,51,53 Since it is challenging to precisely vary the size and decorrelation rate of a hidden perturbation within these prior phantom setups, we instead created a digitally addressable phantom model that utilizes a DMD as the experimental perturbation source,54 which was placed immediately behind a thick and rapid decorrelating tissue phantom (5 mm–10 mm thickness) with scattering properties that closely match the living tissue [Fig. 2(b)]. We varied the relative strength and speed of perturbations by changing both the size and frequency of projected patterns within a small DMD area [Fig. 2(c)], which provided a flexible means to mimic the effects of additional decorrelation induced by blood flow (for example). Measured decorrelation lifetimes arising from our tissue phantom generally fell between 10−4 s and 10−5 s, with DMD-induced perturbations possible up to 22 kHz. We note that the DMD does not produce perturbations equivalent to any particular biological event. We instead used it here to facilitate quantitative analysis of the relative performance variations of parallelized DCS detection. In the supplementary material, note 3, we also show that with standard simulation tools,55 the generated response of our DMD perturbations can be approximately related to the expected size of a biologically realistic inclusion embedded within a semi-infinite medium.56–59 

The schematic diagram of the DCS system for our DMD phantom study is shown in Figs. 1 and 2. The decorrelating tissue phantom was created by a 1 μm polystyrene microsphere-in-water suspension (Polysciences, Inc., Warrington, PA) in custom made cuvettes. The front and back walls of the cuvette perpendicular to the fiber tips were made of microscope coverslips with a thickness of 0.17 mm. The back wall coverslip was placed directly against the surface of a DMD (Digital Light Innovations; Austin, TX; 1024 micromirrors × 768 micromirrors with a micromirror pitch size of 13.68 μm) and completely covered its surface [see Fig. 2(c)]. We aligned the center of the active region of the DMD to lie at the midpoint between the collection and illumination fibers. In this study, we tested two different phantom tissue thicknesses using cuvettes with thicknesses of 5 mm and 10 mm [Fig. 2(b)]. We also tested scattering tissue phantoms with three different reduced scattering coefficients, μs′ = 12.0 cm−1, 7.2 cm−1, 3.6 cm−1, created by varying the microsphere-in-water concentration (please refer to the supplementary material, note 4, for phantom optical properties). These coefficient values are close to previously reported reduced scattering coefficients of human tissues such as skin, brain, and other soft tissues.60 

The DMD could generate binary patterns. The middle square area (“perturbation area”) switches between the “on” and “off” states at a rate of multiple kHz to simulate changes caused by deep tissue motion, while the peripheral area is a random pattern that remains constant. We varied the size of the central perturbation area from 1.37 mm2 × 1.37 mm2 to 9.58 mm2 × 9.58 mm2, as shown in Fig. 2(c). The measurement without the DMD-induced perturbation is achieved by setting the entire DMD area to be a random pattern that remained constant.

Light from a 670 nm long coherence laser (Opto Engine LLC, Midvale, UT) with 150 mW input power was delivered by a multimode fiber (core diameter = 50 μm and NA = 0.22, Thorlabs, Newton, NJ) to the surface of the tissue phantom. To match power conditions required during in vivo experiments, we also followed the American National Standard for Safe Use of Lasers (ANSI) power limit for our phantom studies, which requires an irradiance less than 200 mW/cm2.20 We placed the illumination fiber tip 22 mm away from the cuvette surface to form a 10 mm diameter spot on the phantom surface, which provided our system with a peak irradiance of about 191 mW/cm2. As the mean penetration depth of collected photons is roughly one-third to one-half of the of the source–detector separation d,11 we set the distance between the center of the illumination spot and the detection fiber (core diameter = 1000 μm and NA = 0.5, Thorlabs, Newton, NJ) to be within this range. Specifically, for the 5 mm thickness phantom, the source–detector separation was set to 14 mm, and for the 10 mm thickness phantom, the source–detector separation was set to 21 mm. However, we note that variations in the precise definition of the penetration depth certainly exist and care should be taken in optimizing the value d.61 The SPAD array (Photon Force PF-32) we used consists of 32 × 32 single SPADs, and each single SPAD is 50 μm × 50 μm in size, with an active area of ∼7 μm in diameter [Fig. 2(e)]. The sampling rate of each SPAD is 333 kHz with a bit depth of 8.

To maximize the measurement SNR, it is beneficial to detect one to several optical modes (i.e., speckle grains) per SPAD.62 To enable this criterion across the SPAD array, we ensured that the average speckle size contained within the optical field reaching the array is equal to, or a large fraction of, the size of each SPAD active area. This was achieved by adjusting the distance (R) between the detection fiber and the photosensitive surface of the SPAD array. Please note that the distal end of the detection fiber is held in place via an optical mount such that it is positioned at desired distance R from the surface of the SPAD array active area with free space between. Free-space propagation enlarges the speckle pattern emerging from the detection fiber to ensure effective sampling by the SPAD array. Selecting a distance (R) that is too small will cause each SPAD to contain many independent speckle grains, which leads to a reduced dynamic speckle contrast, while selecting a distance that is too large will yield a low SNR. To match the average speckle size to the SPAD pixel active area, which yields maximum contrast, we calculated the speckle size ds generated from a fiber as a function of the distance and fiber core diameter using the following equation:63 

ds=1.4λRD,
(1)

where λ is the wavelength of the illumination laser, R is the distance between the detection fiber and SPAD array, and D is the core diameter of the detection fiber. ds should be similar to the size of the SPAD pixel active area of 7 μm. In our setup, we used a multimode fiber with D = 1 mm, making the ideal fiber-to-SPAD distance R = 7.5 mm. We note that the relatively low fill factor of current SPAD array technology [7 μm SPAD active area with 50 μm pixel pitch in our setup, Fig. 2(e)] does not significantly impact parallelized DCS given the above measurement SNR constraint.

The intensity fluctuation sequence of all pixels in the SPAD array is first stored in the internal memory of the camera. After a measurement window is finished, the data are transferred to the computer memory for post-processing (see the data flow chart in the supplementary material, note 5). We first computed the autocorrelation curve at each SPAD pixel separately and then averaged the computed autocorrelation curves into a final result [Figs. 2(f)2(h)]. We note that there are many other possible methods to process the SPAD array’s raw 3D video frame data that we plan to examine in the future. To quantify the scattered field’s fluctuations, we computed the normalized temporal intensity autocorrelation function at the jth SPAD pixel, g2(j)(τ), using the following equation:

g2(j)(τ)=t=0PΔtIj(t)Ij(t+τ)t=0PΔtIj(t)2,
(2)

where Ij(t) is the intensity of the jth SPAD in the array at a given time t, τ is the “lag” or correlation time, and we perform a summation over P temporal samples until reaching a desired degree of statistical averaging for one autocorrelation curve. In our phantom studies, each autocorrelation function g2 was typically calculated across P = 10 000 or 35 000 temporal points, which yields one new autocorrelation curve every PΔt = 30 ms or 105 ms, where Δt = 3 µs is the SPAD’s temporal sampling rate. We will refer to PΔt as the “autocorrelation curve rate,” which is different than the SPAD temporal sampling rate and defines the effective measurement rate for our system. After calculating g2(j)(τ) for each single SPAD [Fig. 2(g)], we then average these autocorrelations curves across j = 1 to M SPAD pixels to compute the final system autocorrelation function,

g2(τ)=j=1Mg2(j)(τ),
(3)

which yields a curve with a higher SNR, as shown in Figs. 2(g) and 2(h). In most of our tests, we used a lag time that ranged from τmin = 0 to τmax = 300 µs for g2 computation, as we observed that curves decayed to ∼1 when the lag time τ was larger than 200 µs.

The autocorrelation function, g2(τ), is related to the normalized electric field autocorrelation function, g1(τ), by the Siegert relation,64 

g2(τ)=1+β|g1(τ)|2,
(4)

where β is a coefficient factor. To attain a decorrelation measure, we fit g2(τ) using the following equation:

f(τ)=1+a×ebτ
(5)

and report td defined as the delay time when g2 decays to 1/e of its initial value [g2(3 µs)] as the “decorrelation time.” We sample td at the system autocorrelation curve rate PΔt and report the coefficient b of the index term as “decorrelation speed.”

Figure 3 demonstrates the accuracy gained by using 1024 SPADs for parallelized DCS measurement. During this acquisition, we used a 5 mm thick tissue phantom with μs′ = 12.0 cm−1 and did not induce any additional perturbation with the DMD, showing a tissue phantom decorrelation time of ∼40 µs. Figure 3(a) shows representative raw data captured by a single SPAD over 10 ms. The average number of photons per pixel per sample is ∼0.58. Figures 3(b)3(c) show the resulting autocorrelation curves generated at 10 ms and 100 ms autocorrelation curve rates. The blue curves represent the results using measurements from one SPAD, g2(1)(τ), and approximately match previous experimental measurements from common single-detector DCS setups. The red curves are the average of the autocorrelations g2(τ) from all M = 1024 SPADs within the array. The results show that with longer autocorrelation curve rates, the g2(τ) curve of a single SPAD smoothens to approximately match the curve produced by all 1024 SPADs, but when the autocorrelation curve rate is short, the results of a single SPAD not only have a low SNR but also deviate from the averaged curve created by our parallelized approach.

FIG. 3.

The average decorrelation curve generated by multiple SPAD array pixels increases the DCS accuracy. (a) Raw data of temporal light intensity fluctuation from one SPAD pixel. Autocorrelation curves calculated using a single SPAD (blue) and all 1024 SPAD pixels (red) with a measurement window (i.e., the autocorrelation curve rate) of (b) 10 ms and (c) 100 ms, respectively.

FIG. 3.

The average decorrelation curve generated by multiple SPAD array pixels increases the DCS accuracy. (a) Raw data of temporal light intensity fluctuation from one SPAD pixel. Autocorrelation curves calculated using a single SPAD (blue) and all 1024 SPAD pixels (red) with a measurement window (i.e., the autocorrelation curve rate) of (b) 10 ms and (c) 100 ms, respectively.

Close modal

To demonstrate the sensitivity gains of highly parallelized DCS, we tested the system’s ability to detect an embedded perturbation using a variable number of SPADs to compute a final g2(τ) curve. We compared the averaged g2(τ) curve that resulted from using 1, 9, and 1024 SPADs to detect a 4.1 mm2 × 4.1 mm2 perturbation area, presented at 5 kHz, behind a 5 mm thick tissue phantom with μs′ = 12.0 cm−1. We first calculated the autocorrelation curve of each SPAD pixel using an autocorrelation curve rate of 30 ms (P = 10 000 temporal measurements detected at 333 kHz) and then averaged the resulting g2(τ) curves from M = 1, 9, and 1024 SPADs, as shown in Fig. 4(a). After repeating this experiment over 420 independent trials, we computed the autocorrelation curve mean ± standard deviation (SD), displayed as solid center lines with a band of finite width in Fig. 4(b). The red and blue bands represent measurements performed with and without the perturbation present, respectively, where the measurement without a perturbation had the entire DMD area set to a random constant pattern (see Sec. II A). Here, we see that the averaged autocorrelation g2 curves are either completely indistinguishable or only partially distinguished when using either 1 or 9 SPADs but become completely distinguishable (i.e., fully separated) when using 1024 SPADs, clearly reflecting an improved system sensitivity.

FIG. 4.

Demonstration of increased sensitivity to embedded perturbations with an increased number of SPADs. (a) Autocorrelation curves calculated using 1, 9 (3 × 3), and all 1024 (32 × 32) SPADs, respectively, with an autocorrelation curve rate (PΔt) of 30 ms and a 4.1 mm2 × 4.1 mm2 DMD perturbation area; (b) autocorrelation curve separability and the corresponding perturbation detection accuracy using 1, 9, and 1024 SPADs. Banded curves represent the mean ± SD of 420 averaged autocorrelations, where the red curve is measurement with perturbation and the blue curve is measurement without perturbation; (c) the perturbation detection accuracy increases with a larger number of SPADs for three example setup configurations: 5 mm phantom thickness and 2.7 mm2 × 2.7 mm2 perturbation area (blue), 5 mm phantom thickness and 4.1 mm2 × 4.1 mm2 perturbation area (red), and 10 mm phantom thickness and 9.6 mm2 × 9.6 mm2 perturbation area (yellow).

FIG. 4.

Demonstration of increased sensitivity to embedded perturbations with an increased number of SPADs. (a) Autocorrelation curves calculated using 1, 9 (3 × 3), and all 1024 (32 × 32) SPADs, respectively, with an autocorrelation curve rate (PΔt) of 30 ms and a 4.1 mm2 × 4.1 mm2 DMD perturbation area; (b) autocorrelation curve separability and the corresponding perturbation detection accuracy using 1, 9, and 1024 SPADs. Banded curves represent the mean ± SD of 420 averaged autocorrelations, where the red curve is measurement with perturbation and the blue curve is measurement without perturbation; (c) the perturbation detection accuracy increases with a larger number of SPADs for three example setup configurations: 5 mm phantom thickness and 2.7 mm2 × 2.7 mm2 perturbation area (blue), 5 mm phantom thickness and 4.1 mm2 × 4.1 mm2 perturbation area (red), and 10 mm phantom thickness and 9.6 mm2 × 9.6 mm2 perturbation area (yellow).

Close modal

To quantify the detection sensitivity, we computed a simple classification accuracy metric by alternately measuring the decorrelation time td with and without a DMD-induced perturbation and comparing the two adjacent measurements. This alternate measurement method can minimize the disturbance of td caused by phantom and experimental environment changes. If the td determined with the perturbation present was smaller than the td without the perturbation present, we assumed successful perturbation detection, and if not, we assumed that the detection was failed. This simple classification accuracy metric, defined as the proportion of successful detections over 420 trials, is plotted as a function of the number of SPADs (M) used to compute the average decorrelation curve [Fig. 4(c)]. The classification accuracy generally increases with an increase in the number of SPADs. Different color curves represent different tissue phantom properties and DMD perturbation areas, and the accuracy increases when using a thinner tissue phantom or larger perturbation area.

For a specific tissue phantom, the ability to detect a perturbation depends on both the perturbation area size and autocorrelation curve rate. We, thus, varied both of these parameters while repeating our classification accuracy tests. The results of this exercise are summarized in Figs. 5 and 6. For each setting, we repeated 120 independent experiments to calculate the classification accuracy. The mean ± SD of several example g2 curves (5 mm thick tissue phantom with μs′ = 12.0 cm−1) is shown in Fig. 5, where once again red and blue curves represent the case of perturbation present vs absent, respectively. The curves within each row use the same autocorrelation curve rate, while the curves within each column have the same DMD perturbation area. The results show that when using an autocorrelation curve rate of 10 ms, it is possible to accurately identify a 3.42 mm2 × 3.42 mm2 perturbation via the g2 curve separation. When the autocorrelation curve rate increases to 30 ms or 100 ms, it is possible to detect smaller perturbations, down to sizes of 2.74 mm2 × 2.74 mm2 and 2.05 mm2 × 2.05 mm2, respectively.

FIG. 5.

Parallelized DCS detection of mm-scale size perturbations at different autocorrelation curve rates. (a)–(c) Detection results with 10 ms autocorrelation curve rate for perturbation areas of size 2.05 mm2 × 2.05 mm2, 2.74 mm2 × 2.74 mm2 and 3.42 mm2 × 3.42 mm2, respectively. Banded curves represent the mean ± SD of 120 averaged autocorrelation curves, where the red curve indicates with perturbation, blue curve band indicates without, and red box is a selected zoom-in. (d)–(f) Similar detection results using a 30 ms and (g)–(i) 100 ms autocorrelation curve rate.

FIG. 5.

Parallelized DCS detection of mm-scale size perturbations at different autocorrelation curve rates. (a)–(c) Detection results with 10 ms autocorrelation curve rate for perturbation areas of size 2.05 mm2 × 2.05 mm2, 2.74 mm2 × 2.74 mm2 and 3.42 mm2 × 3.42 mm2, respectively. Banded curves represent the mean ± SD of 120 averaged autocorrelation curves, where the red curve indicates with perturbation, blue curve band indicates without, and red box is a selected zoom-in. (d)–(f) Similar detection results using a 30 ms and (g)–(i) 100 ms autocorrelation curve rate.

Close modal
FIG. 6.

Parallelized DCS classification accuracy of different perturbation sizes, tissue phantom properties, and autocorrelation curve rates. (a) The classification accuracy increases for larger perturbation areas and lower reduced scattering coefficients. Here, μs′ = 3.6, 7.2, and 12.0 cm−1 for red, blue, and green curves, respectively, for a 5 mm tissue phantom thickness and 30 ms autocorrelation curve rate. (b) Matching results with the 10 mm thick tissue phantom. (c) Comparison of the classification accuracy of different autocorrelation curve rates with the 10 mm tissue phantom. The dotted line is for 30 ms and the solid line is for 100 ms. (d) The minimum detectable perturbation area as a function of the number of SPADs. Yellow, red, and blue curves represent results corresponding to classification accuracy thresholds at 0.8, 0.85, and 0.9 for a 5 mm tissue phantom thickness and μs′ = 12.0 cm−1.

FIG. 6.

Parallelized DCS classification accuracy of different perturbation sizes, tissue phantom properties, and autocorrelation curve rates. (a) The classification accuracy increases for larger perturbation areas and lower reduced scattering coefficients. Here, μs′ = 3.6, 7.2, and 12.0 cm−1 for red, blue, and green curves, respectively, for a 5 mm tissue phantom thickness and 30 ms autocorrelation curve rate. (b) Matching results with the 10 mm thick tissue phantom. (c) Comparison of the classification accuracy of different autocorrelation curve rates with the 10 mm tissue phantom. The dotted line is for 30 ms and the solid line is for 100 ms. (d) The minimum detectable perturbation area as a function of the number of SPADs. Yellow, red, and blue curves represent results corresponding to classification accuracy thresholds at 0.8, 0.85, and 0.9 for a 5 mm tissue phantom thickness and μs′ = 12.0 cm−1.

Close modal

Figures 6(a) and 6(b) summarizes the classification accuracy of our parallelized DCS system under different phantom scattering coefficients μs′ and for different perturbation sizes for 5 mm and 10 mm thick tissue phantoms, respectively. A 30 ms autocorrelation curve rate was used to calculate g2 in (a) and (b). It is clear that the classification accuracy increases with a larger perturbation area and a lower tissue phantom reduced scattering coefficient and thickness. The classification accuracy also improves with a longer autocorrelation curve rate [100 ms vs 30 ms, Fig. 6(c)]. Figure 6(d) plots the minimum detectable perturbation area at different detection accuracy thresholds as a function of the number of SPADs M used for parallelized detection (5 mm thickness, μs′ = 12.0 cm−1, and the autocorrelation curve rate of 30 ms). These results show the significant improvement in the system sensitivity with increased measurement parallelization.

To test the effect of the temporal perturbation rate on the average autocorrelation function of our parallelized DCS system, we varied the refresh frequency of the DMD perturbation area between 1 kHz and 20 kHz for 3 different perturbation area sizes with the results shown in Fig. 7. The tissue phantom thickness in this experiment was 5 mm with μs′ = 12.0 cm−1. With an autocorrelation curve rate of 30 ms, the g2 curve clearly dips with increasing DMD frequency, and it is easier to differentiate g2 curves corresponding to different perturbation frequencies for larger perturbation areas, as expected.

FIG. 7.

The effect of perturbation’s temporal rate on the average autocorrelation function of the parallelized DCS system. Different colored autocorrelation curves represent different perturbation frequencies. The perturbation areas corresponding to (a), (b), and (c) are 1.37 mm2 × 1.37 mm2, 4.10 mm2 × 4.10 mm2, and 6.84 mm2 × 6.84 mm2, respectively. The autocorrelation curve rate is 30 ms.

FIG. 7.

The effect of perturbation’s temporal rate on the average autocorrelation function of the parallelized DCS system. Different colored autocorrelation curves represent different perturbation frequencies. The perturbation areas corresponding to (a), (b), and (c) are 1.37 mm2 × 1.37 mm2, 4.10 mm2 × 4.10 mm2, and 6.84 mm2 × 6.84 mm2, respectively. The autocorrelation curve rate is 30 ms.

Close modal

To demonstrate that our highly parallelized DCS system can also improve the SNR of in vivo measurements, we modified our setup to perform in vivo forehead blood flow experiments [Fig. 9(a)]. To reduce motion effects, the participant rested their head on a chinrest, which was further secured by a head strap. Similar to the phantom study setup, the light was delivered by the same multimode fiber to the participant’s forehead. We fixed both the illumination and the collection fibers via fiber clamps onto two manual translation stages, which were attached to the mounted chinrest assembly. We attached a 3D printed holder to the tip of the illumination fiber, and before each experiment, we adjusted the positions of the illumination fiber holder and the collection fiber tip to touch the top right of each subject’s forehead by translating the two stages. With the distance from the fiber tip to the subject’s forehead fixed at 11 mm (via the fiber holder length), an illumination spot with a diameter of ∼5 mm was created on the forehead, and we further reduced the illumination power to 40 mW to comply with the ANSI standard (less than 200 mW/cm2). The distance between the center of the illumination spot and the detection fiber was set to 1.5 cm, and the SPAD array sampling rate was kept at 333 kHz. The study was approved by the Duke campus Institutional Review Board, and six healthy adult volunteers were enrolled after obtaining a written informed consent.

In this experiment, we got 21 s of continuous frame data for each subject and verified our measurements with a commercial pulse oximeter (EMAY, Ltd.; Hong Kong, China). We processed the data to generate decorrelation speed plots sampled at different autocorrelation curve rates, PΔt, as well as with using measurements from a variable number of SPADs, M, to demonstrate the temporal resolution advantages of our parallelized approach. Example results from these tests are shown in Fig. 8. Similar to our phantom studies, we fit the autocorrelation curves using an exponential function and extract and plot the decay coefficient b [Eq. (5)]. A larger b value represents faster decorrelation speed (i.e., faster blood flow).

FIG. 8.

Human forehead pulse measurement results. Representative autocorrelation curves (blue) and their best exponential function fit (red) are shown at left [(a), (c), and (e)]. At right are associated DCS pulse measurements over 10 s [(b), (d), and (f)], which are represented by the plot of normalized decorrelation speed, extracted from autocorrelation curves after exponential fitting. (a) and (b) Using M = 10 SPADs and an autocorrelation curve rate PΔt = 100 ms, (c) and (d) M = 1024 SPADs and an autocorrelation curve rate of PΔt = 100 ms, and (e) and (f) M = 1024 SPADs and an autocorrelation curve rate, PΔt = 30 ms. (g) The Fourier transform of the DCS pulse signal (f) and the pulse signal acquired using a commercial pulse meter (h).

FIG. 8.

Human forehead pulse measurement results. Representative autocorrelation curves (blue) and their best exponential function fit (red) are shown at left [(a), (c), and (e)]. At right are associated DCS pulse measurements over 10 s [(b), (d), and (f)], which are represented by the plot of normalized decorrelation speed, extracted from autocorrelation curves after exponential fitting. (a) and (b) Using M = 10 SPADs and an autocorrelation curve rate PΔt = 100 ms, (c) and (d) M = 1024 SPADs and an autocorrelation curve rate of PΔt = 100 ms, and (e) and (f) M = 1024 SPADs and an autocorrelation curve rate, PΔt = 30 ms. (g) The Fourier transform of the DCS pulse signal (f) and the pulse signal acquired using a commercial pulse meter (h).

Close modal

The representative autocorrelation curves and the pulse signals calculated using the normalized decorrelation speed b are shown in Figs. 8(a)8(f). With an autocorrelation curve rate of PΔt = 100 ms and averaging all M = 1024 SPADs, we obtain a high-SNR autocorrelation function g2(τ) [Fig. 8(c)] with a decorrelation speed that clearly follows the periodicity of the subject pulse [Fig. 8(d)]. Two deflections caused by ventricle contractions and repolarization can be clearly observed. However, if only M = 10 SPADs were used for averaging, pulse signal retrieval remains challenging [Fig. 8(b)], as the decorrelation speed b extracted from the exponential fit is heavily affected by a low SNR autocorrelation function g2(τ) [Fig. 8(a)]. Moreover, we are still able to retrieve the pulse signal with a shorter autocorrelation curve rate, PΔt = 30 ms [Figs. 8(e) and 8(f)]. This is further validated by examining the Fourier transform of the DCS decorrelation speed signal [Fig. 8(f)] and the pulse signal simultaneously acquired by the pulse oximeter [Fig. 8(h)], shown in Fig. 8(g), which have approximately the same peak frequency at 1.1 Hz. Consistent results were obtained for six different human subjects.

As an additional demonstration of the sensitivity gains of highly parallelized DCS, we also measured dynamic variations induced by a behavioral task. The prefrontal cortex region is considered to manage in part cognitive processes such as planning, cognitive flexibility, and working memory.65 Previous studies with the complementary non-invasive optical measurement approach termed interferometric near-infrared spectroscopy (iNIRS), for example, measured and reported dynamic changes within the prefrontal cortex area for subjects who read a paragraph of unfamiliar text after a 10 min rest period.6 Here, we explored the likelihood of observing such dynamic activity using autocorrelation functions obtained from our highly parallelized DCS system. The schematic diagram of the experimental setup is shown in Fig. 9(a). The subjects first rested for 5 min and then began a 15 min test divided into three stages: reading stage 1, resting stage, and reading stage 2, where each stage lasted 5 min. Due to the data transfer limitation from the SPAD array to computer memory, we collected DCS measurements for the first 10 s of every test minute using all 1024 SPADs. The autocorrelation function g2(τ) and the associated decorrelation time td were calculated using an autocorrelation curve rate of 100 ms. A representative plot of the decorrelation time measurement of one subject is shown in Fig. 9(b). We then calculated the mean decorrelation time for each subject by averaging the 100 values of td that comprise each 10 s DCS measurement window [see the plot in Fig. 9(c)]. A total of six subjects were enrolled in this study. The mean decorrelation curve for each enrolled subject is shown in the supplementary material, Fig. 6(a). To better evaluate whether there exist unwanted motion artifacts during the experiment, we introduced the SD, σd, of all the decorrelation time, td, within each 10 s measurement window (supplementary material, Fig. 7). σd should be relatively stable, as the same behavioral task was performed within each measurement window. After obtaining the 95% confidence interval of σd from all the six subjects (red dashed lines in the supplementary material, Fig. 7, 90 total measurements), we removed experiments exhibiting a decorrelation time standard deviation σd outside this confidence interval. In this study, subject No. 3 and No. 6 included such data, suggesting that significant motion occurred, which we aim to address with a new setup in future work (see Sec. V). The mean decorrelation time of the remaining four subjects are shown in Fig. 9(d) with the mean ± SD shown as a dotted line with a banded curve. Please note that when averaging the data of different subjects, we have divided the data of each subject by the average decorrelation value of their respective measurement sequence (computed across the 20 min experiment). We performed this normalization to avoid uneven weighting of different subjects’ measurement sequences to our cross-subject analysis. The solid horizontal line represents the mean decorrelation time across each of the three test stages (the first reading stage, resting stage, and the second reading stage). The mean decorrelation time in both reading stages (blue) was generally observed to be lower than that in the resting stage (red), indicating an increase in blood flow speed within the measurement area (i.e., activation of the prefrontal cortex), matching prior results.6 

FIG. 9.

The experimental setup and results of the human prefrontal cortex activation test. (a) The schematic diagram of the experimental setup. (b) The plot of the decorrelation time (td) values over the 15 min test including two reading stages and one intermediate rest stage, where 10 s of the signal is collected every minute. (c) The plot of mean decorrelation time corresponding to (b). Each mean decorrelation time value is obtained by averaging all the decorrelation time values within the corresponding 10 s window. (d) Mean ± SD results of the mean decorrelation time of 4 subjects calculated after dividing the data of each subject by the average decorrelation value of their respective measurement sequence. Solid horizontal lines represent the average of the five normalized decorrelation times in each stage.

FIG. 9.

The experimental setup and results of the human prefrontal cortex activation test. (a) The schematic diagram of the experimental setup. (b) The plot of the decorrelation time (td) values over the 15 min test including two reading stages and one intermediate rest stage, where 10 s of the signal is collected every minute. (c) The plot of mean decorrelation time corresponding to (b). Each mean decorrelation time value is obtained by averaging all the decorrelation time values within the corresponding 10 s window. (d) Mean ± SD results of the mean decorrelation time of 4 subjects calculated after dividing the data of each subject by the average decorrelation value of their respective measurement sequence. Solid horizontal lines represent the average of the five normalized decorrelation times in each stage.

Close modal

In this study, we demonstrated the highly parallelized acquisition of the DCS signal using an integrated SPAD array that contained 32 × 32 pixels, which allowed us to average over a thousand independent speckle fluctuation measurements per frame to improve the system sensitivity and speed. We found that the SPAD array provides enhanced sensitivity as compared to a single or a small number of SPADs by testing its performance with a novel perturbation model, implemented by inserting a DMD beneath a thick tissue phantom, which demonstrated accurate detection of mm-scale perturbations under 1 cm of the tissue-like phantom. We also used this new approach to measure the blood pulse from the human forehead with a 30 ms temporal resolution and detect behavior-induced variations from above the pre-frontal cortex.

Monolithic CMOS SPAD arrays continue to rapidly increase in pixel count (e.g., an 1 megapixel array is now available66), and recent consumer-oriented developments have led to combined sensing and processing layers via 3D stacking that open up new functionalities.34 Like other integrated CMOS technology, we expect both pixel counts and functionality to increase with improvements from an ever-advancing semiconductor manufacturing industry. Prior work has shown that more than one single-photon detector can increase the DCS measurement sensitivity for the challenging task of non-invasive detection of functional brain dynamics14,67 and tomographically measure cerebral blood flow,29,51 for example. Our new approach here offers a means to implement such efforts with dozens of times more detectors and thus increase speed or sensitivity by a proportional factor, thus opening up a new suite of possible experiments.

The sampling rate of our SPAD array in this study is 333 kHz, which is slower than the typical sampling rate (1 MHz or more) of single SPADs in previously reported DCS setups. Nevertheless, as our in vivo study and phantom-based experiments demonstrate, this lower sampling rate is still fast enough to retrieve the entire autocorrelation function of the relevant biological activity such as blood flow with a source–detector distance of 1.5 cm. Our SPAD array currently operates at a bit depth of 8, which is more than sufficient for most of DCS applications. In future work, we aim to work at an increased sampling rate of 720 kHz with a decreased per-sample bit depth, which is possible with our current hardware.

Although our DMD-based phantom platform enables quantitative analysis of sensitivity advantages of our SPAD array DCS system, the introduced micromirror perturbation is not equivalent to any particular biological event. In addition, our phantom is also not semi-infinite. While we demonstrated both experimentally and theoretically that we can obtain decorrelation curves that are similar to those measured from real human tissues, further studies are needed to compare the dynamics of our DMD-based phantom with true biological events, such as blood flow. Another limitation of our current phantom setup is that the glass surface of the cuvette likely introduced light-guiding effects, which might lead to leakage of significant amount of photons and overwhelm the actual photons from the deep tissue. To overcome this limitation, highly scattering walls such as Mylar foils could be employed in future experiments, which will likely improve the sensitivity of the phantom-based measurement system.68 

In our in vivo experiments, one potential concern with the use of a multi-mode fiber to deliver speckle onto the SPAD array is the introduction of unwanted additional speckle decorrelation due to fiber movement (bending, twisting, etc.). To reduce this effect, we limited fiber motion by not attaching the fiber to the subject’s head using pads or stripes. Instead, we fixed the fiber to a stationary post and ensured that the collection fiber tip touched the subject’s forehead when the subject was rested on the chinrest [see Fig. 9(a)]. Because the fiber and the subject’s head were not directly coupled, we do not expect that subject motion caused additional decorrelation due to fiber movement. In addition, we expect the effects of added speckle decorrelation, induced by a moving fiber, to be quite minimal, as the decorrelation speed caused by the light’s long journey through the in vivo tissue is much greater than the decorrelation speed caused by a relatively slowly moving fiber (see the experiment in the supplementary material, Note 7). Instead, we expect that relative subject motion with respect to a fixed multi-mode fiber tip during long-term measurement introduced spurious decorrelation variations in our reading task experiment, for it will cause rapid changes in the measurement position and fiber tip pressure on the skin. We aim to address this with a modified measurement arrangement and more appropriate behavioral tasks in future work.

In the human prefrontal cortex activation test, the source–detector separation is set to 15 mm, which is shorter than the separation of 20 mm–30 mm reported in some prior DCS setups on measuring the prefrontal cortex. While we believe that a significant portion of photons we measured is from the frontal cortex, it is inconclusive whether changes in systemic parameters such as blood pressure or blood flow from shallower depths also contributed to the signal changes we observed during the reading task. In vivo studies with a larger sample number and an increased source detector separation are needed to further validate our findings in this pilot study in the future. In addition, we used a 670 nm illumination source, as the detection sensitivity of our SPAD array is higher in the visible wavelength range. However, a longer wavelength would potentially increase the penetration depth into deep tissues. The illumination wavelength still needs to be carefully optimized.

Finally, the large spatiotemporal datasets acquired per experiment suggest that a new suite of software tools is on the horizon for highly parallelized DCS setups.

W.L. and R.Q. contributed equally to this work.

See the supplementary material for the supporting content.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Research reported in this publication was supported by the National Institute of Neurological Disorders and Stroke of the National Institutes of Health under Award No. RF1NS113287 as well as the Duke-Coulter Translational Partnership. W.L. acknowledges the support from the China Scholarship Council. The authors would like to thank Kernel, Inc., for their kind assistance.

1.
D. A.
Boas
,
D. H.
Brooks
,
E. L.
Miller
,
C. A.
DiMarzio
,
M.
Kilmer
,
R. J.
Gaudette
, and
Q.
Zhang
, “
Imaging the body with diffuse optical tomography
,”
IEEE Signal Process. Mag.
18
,
57
75
(
2001
).
2.
A.
Cerussi
,
D.
Hsiang
,
N.
Shah
,
R.
Mehta
,
A.
Durkin
,
J.
Butler
, and
B. J.
Tromberg
, “
Predicting response to breast cancer neoadjuvant chemotherapy using diffuse optical spectroscopy
,”
Proc. Natl. Acad. Sci. U. S. A.
104
,
4014
4019
(
2007
).
3.
M.
Ferrari
and
V.
Quaresima
, “
A brief review on the history of human functional near-infrared spectroscopy (fNIRs) development and fields of application
,”
Neuroimage
63
,
921
935
(
2012
).
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
).
5.
D.
Borycki
,
O.
Kholiqov
,
S. P.
Chong
, and
V. J.
Srinivasan
, “
Interferometric near-infrared spectroscopy (iNIRS) for determination of optical and dynamical properties of turbid media
,”
Opt. Express
24
,
329
354
(
2016
).
6.
O.
Kholiqov
,
W.
Zhou
,
T.
Zhang
,
V.
Du Le
, and
V. J.
Srinivasan
, “
Time-of-flight resolved light field fluctuations reveal deep human tissue physiology
,”
Nat. Commun.
11
,
1
15
(
2020
).
7.
J. A.
Izatt
,
M. D.
Kulkarni
,
S.
Yazdanfar
,
J. K.
Barton
, and
A. J.
Welch
, “
In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomography
,”
Opt. Lett.
22
,
1439
1441
(
1997
).
8.
R. K.
Wang
,
S. L.
Jacques
,
Z.
Ma
,
S.
Hurst
,
S. R.
Hanson
, and
A.
Gruber
, “
Three dimensional optical angiography
,”
Opt. Express
15
,
4083
4097
(
2007
).
9.
D.
Briers
,
D. D.
Duncan
,
E.
Hirst
,
S. J.
Kirkpatrick
,
M.
Larsson
,
W.
Steenbergen
,
T.
Stromberg
, and
O. B.
Thompson
, “
Laser speckle contrast imaging: Theoretical and practical limitations
,”
J. Biomed. Opt.
18
,
066018
(
2013
).
10.
T.
Durduran
and
A. G.
Yodh
, “
Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement
,”
Neuroimage
85
,
51
63
(
2014
).
11.
E. M.
Buckley
,
A. B.
Parthasarathy
,
P. E.
Grant
,
A. G.
Yodh
, and
M. A.
Franceschini
, “
Diffuse correlation spectroscopy for measurement of cerebral blood flow: Future prospects
,”
Neurophotonics
1
,
011009
(
2014
).
12.
L.
Colombo
,
S.
Samaei
,
P.
Lanka
,
D.
Ancora
,
M.
Pagliazzi
,
T.
Durduran
,
P.
Sawosz
,
A.
Liebert
, and
A.
Pifferi
, “
Coherent fluctuations in time-domain diffuse optics
,”
APL Photonics
5
,
071301
(
2020
).
13.
M. N.
Kim
,
T.
Durduran
,
S.
Frangos
,
B. L.
Edlow
,
E. M.
Buckley
,
H. E.
Moss
,
C.
Zhou
,
G.
Yu
,
R.
Choe
,
E.
Maloney-Wilensky
 et al, “
Noninvasive measurement of cerebral blood flow and blood oxygenation using near-infrared and diffuse correlation spectroscopies in critically brain-injured adults
,”
Neurocrit. Care
12
,
173
180
(
2010
).
14.
D.
Wang
,
A. B.
Parthasarathy
,
W. B.
Baker
,
K.
Gannon
,
V.
Kavuri
,
T.
Ko
,
S.
Schenkel
,
Z.
Li
,
Z.
Li
,
M. T.
Mullen
 et al, “
Fast blood flow monitoring in deep tissues with real-time software correlators
,”
Biomed. Opt. Express
7
,
776
797
(
2016
).
15.
E. M.
Buckley
,
B. F.
Miller
,
J. M.
Golinski
,
H.
Sadeghian
,
L. M.
McAllister
,
M.
Vangel
,
C.
Ayata
,
W. P.
Meehan
 III
,
M. A.
Franceschini
, and
M. J.
Whalen
, “
Decreased microvascular cerebral blood flow assessed by diffuse correlation spectroscopy after repetitive concussions in mice
,”
J. Cerebr. Blood Flow Metabol.
35
,
1995
2000
(
2015
).
16.
A.
Rajaram
,
G.
Bale
,
M.
Kewin
,
L. B.
Morrison
,
I.
Tachtsidis
,
K. S.
Lawrence
, and
M.
Diop
, “
Simultaneous monitoring of cerebral perfusion and cytochrome c oxidase by combining broadband near-infrared spectroscopy and diffuse correlation spectroscopy
,”
Biomed. Opt. Express
9
,
2588
2603
(
2018
).
17.
L.
He
,
Y.
Lin
,
C.
Huang
,
D.
Irwin
,
M. M.
Szabunio
, and
G.
Yu
, “
Noncontact diffuse correlation tomography of human breast tumor
,”
J. Biomed. Opt.
20
,
086003
(
2015
).
18.
H. S.
Yazdi
,
T. D.
O’Sullivan
,
A.
Leproux
,
B.
Hill
,
A.
Durkin
,
S.
Telep
,
J.
Lam
,
S. S.
Yazdi
,
A. M.
Police
,
R. M.
Carroll
 et al, “
Mapping breast cancer blood flow index, composition, and metabolism in a human subject using combined diffuse optical spectroscopic imaging and diffuse correlation spectroscopy
,”
J. Biomed. Opt.
22
,
045003
(
2017
).
19.
N. B.
Agochukwu
,
C.
Huang
,
M.
Zhao
,
A. A.
Bahrani
,
L.
Chen
,
P.
McGrath
,
G.
Yu
, and
L.
Wong
, “
A novel noncontact diffuse correlation spectroscopy device for assessing blood flow in mastectomy skin flaps: A prospective study in patients undergoing prosthesis-based reconstruction
,”
Plast. Reconstr. Surg.
140
,
26
31
(
2017
).
20.
A. N. S. Institute
,
American National Standard for Safe Use of Lasers
(
Laser Institute of America
,
2007
).
21.
J.
Jönsson
and
E.
Berrocal
, “
Multi-scattering software: Part I: Online accelerated Monte Carlo simulation of light transport through scattering media
,”
Opt. Express
28
,
37612
37638
(
2020
).
22.
C. G.
Favilla
,
R. C.
Mesquita
,
M.
Mullen
,
T.
Durduran
,
X.
Lu
,
M. N.
Kim
,
D. L.
Minkoff
,
S. E.
Kasner
,
J. H.
Greenberg
,
A. G.
Yodh
 et al, “
Optical bedside monitoring of cerebral blood flow in acute ischemic stroke patients during head-of-bed manipulation
,”
Stroke
45
,
1269
1274
(
2014
).
23.
R. C.
Mesquita
,
T.
Durduran
,
G.
Yu
,
E. M.
Buckley
,
M. N.
Kim
,
C.
Zhou
,
R.
Choe
,
U.
Sunar
, and
A. G.
Yodh
, “
Direct measurement of tissue blood flow and metabolism with diffuse optics
,”
Philos. Trans. Math. Phys. Eng. Sci.
369
,
4390
4406
(
2011
).
24.
E. M.
Buckley
,
N. M.
Cook
,
T.
Durduran
,
M. N.
Kim
,
C.
Zhou
,
R.
Choe
,
G.
Yu
,
S.
Shultz
,
C. M.
Sehgal
,
D. J.
Licht
 et al, “
Cerebral hemodynamics in preterm infants during positional intervention measured with diffuse correlation spectroscopy and transcranial Doppler ultrasound
,”
Opt. Express
17
,
12571
12581
(
2009
).
25.
T.
Durduran
,
C.
Zhou
,
B. L.
Edlow
,
G.
Yu
,
R.
Choe
,
M. N.
Kim
,
B. L.
Cucchiara
,
M. E.
Putt
,
Q.
Shah
,
S. E.
Kasner
 et al, “
Transcranial optical monitoring of cerebrovascular hemodynamics in acute stroke patients
,”
Opt. Express
17
,
3884
3902
(
2009
).
26.
J. J.
Selb
,
D. A.
Boas
,
S.-T.
Chan
,
K. C.
Evans
,
E. M.
Buckley
, and
S. A.
Carp
, “
Sensitivity of near-infrared spectroscopy and diffuse correlation spectroscopy to brain hemodynamics: Simulations and experimental findings during hypercapnia
,”
Neurophotonics
1
,
015005
(
2014
).
27.
J. D.
Johansson
,
D.
Portaluppi
,
M.
Buttafava
, and
F.
Villa
, “
A multipixel diffuse correlation spectroscopy system based on a single photon avalanche diode array
,”
J. Biophotonics
12
,
e201900091
(
2019
).
28.
G.
Dietsche
,
M.
Ninck
,
C.
Ortolf
,
J.
Li
,
F.
Jaillon
, and
T.
Gisler
, “
Fiber-based multispeckle detection for time-resolved diffusing-wave spectroscopy: Characterization and application to blood flow detection in deep tissue
,”
Appl. Opt.
46
,
8506
8514
(
2007
).
29.
C.
Zhou
,
G.
Yu
,
D.
Furuya
,
J. H.
Greenberg
,
A. G.
Yodh
, and
T.
Durduran
, “
Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain
,”
Opt. Express
14
,
1125
1144
(
2006
).
30.
D.
Tamborini
,
K. A.
Stephens
,
M. M.
Wu
,
P.
Farzam
,
A. M.
Siegel
,
O.
Shatrovoy
,
M.
Blackwell
,
D. A.
Boas
,
S. A.
Carp
, and
M. A.
Franceschini
, “
Portable system for time-domain diffuse correlation spectroscopy
,”
IEEE Trans. Biomed. Eng.
66
,
3014
3025
(
2019
).
31.
D.
Irwin
,
L.
Dong
,
Y.
Shang
,
R.
Cheng
,
M.
Kudrimoti
,
S. D.
Stevens
, and
G.
Yu
, “
Influences of tissue absorption and scattering on diffuse correlation spectroscopy blood flow measurements
,”
Biomed. Opt. Express
2
,
1969
1985
(
2011
).
32.
S.
Han
,
M. D.
Hoffman
,
A. R.
Proctor
,
J. B.
Vella
,
E. A.
Mannoh
,
N. E.
Barber
,
H. J.
Kim
,
K. W.
Jung
,
D. S.
Benoit
, and
R.
Choe
, “
Non-invasive monitoring of temporal and spatial blood flow during bone graft healing using diffuse correlation spectroscopy
,”
PLoS One
10
,
e0143891
(
2015
).
33.
C. J.
Stapels
,
N. J.
Kolodziejski
,
D.
McAdams
,
M. J.
Podolsky
,
D. E.
Fernandez
,
D.
Farkas
, and
J. F.
Christian
, “
A scalable correlator for multichannel diffuse correlation spectroscopy
,”
Proc. SPIE
9698
,
969816
(
2016
).
34.
C.
Bruschini
,
H.
Homulle
,
I. M.
Antolovic
,
S.
Burri
, and
E.
Charbon
, “
Single-photon avalanche diode imagers in biophotonics: Review and outlook
,”
Light: Sci. Appl.
8
,
1
28
(
2019
).
35.
S.
Chan
,
A.
Halimi
,
F.
Zhu
,
I.
Gyongy
,
R. K.
Henderson
,
R.
Bowman
,
S.
McLaughlin
,
G. S.
Buller
, and
J.
Leach
, “
Long-range depth imaging using a single-photon detector array and non-local data fusion
,”
Sci. Rep.
9
,
1
10
(
2019
).
36.
L.
Pancheri
and
D.
Stoppa
, “
A spad-based pixel linear array for high-speed time-gated fluorescence lifetime imaging
,” in
2009 Proceedings of ESSCIRC
(
IEEE
,
2009
), pp.
428
431
.
37.
L.
Parmesan
,
N.
Dutton
,
N. J.
Calder
,
N.
Krstajic
,
A. J.
Holmes
,
L. A.
Grant
, and
R. K.
Henderson
, “
A 256× 256 spad array with in-pixel time to amplitude conversion for fluorescence lifetime imaging microscopy
,” in
International Image Sensor Workshop, Vaals, The Netherlands, Memory
(
International Image Sensor Society (IISS)
,
2015
), Vol. 900.
38.
D.-U.
Li
,
J.
Arlt
,
J.
Richardson
,
R.
Walker
,
A.
Buts
,
D.
Stoppa
,
E.
Charbon
, and
R.
Henderson
, “
Real-time fluorescence lifetime imaging system with a 32× 32 0.13 μm cmos low dark-count single-photon avalanche diode array
,”
Opt. Express
18
,
10257
10269
(
2010
).
39.
G.
Giraud
,
H.
Schulze
,
D.-U.
Li
,
T. T.
Bachmann
,
J.
Crain
,
D.
Tyndall
,
J.
Richardson
,
R.
Walker
,
D.
Stoppa
,
E.
Charbon
 et al, “
Fluorescence lifetime biosensing with dna microarrays and a cmos-spad imager
,”
Biomed. Opt. Express
1
,
1302
1308
(
2010
).
40.
E.
Pedretti
,
M. G.
Tanner
,
T. R.
Choudhary
,
N.
Krstajić
,
A.
Megia-Fernandez
,
R. K.
Henderson
,
M.
Bradley
,
R. R.
Thomson
,
J. M.
Girkin
,
K.
Dhaliwal
 et al, “
High-speed dual color fluorescence lifetime endomicroscopy for highly-multiplexed pulmonary diagnostic applications and detection of labeled bacteria
,”
Biomed. Opt. Express
10
,
181
195
(
2019
).
41.
J. L.
Lagarto
,
C.
Credi
,
F.
Villa
,
S.
Tisa
,
F.
Zappa
,
V.
Shcheslavskiy
,
F. S.
Pavone
, and
R.
Cicchi
, “
Multispectral depth-resolved fluorescence lifetime spectroscopy using spad array detectors and fiber probes
,”
Sensors
19
,
2678
(
2019
).
42.
J.
Kostamovaara
,
J.
Tenhunen
,
M.
Kögler
,
I.
Nissinen
,
J.
Nissinen
, and
P.
Keränen
, “
Fluorescence suppression in Raman spectroscopy using a time-gated CMOS SPAD
,”
Opt. Express
21
,
31632
31645
(
2013
).
43.
J.
Blacksberg
,
Y.
Maruyama
,
E.
Charbon
, and
G. R.
Rossman
, “
Fast single-photon avalanche diode arrays for laser Raman spectroscopy
,”
Opt. Lett.
36
,
3672
3674
(
2011
).
44.
I.
Nissinen
,
J.
Nissinen
,
J.
Holma
, and
J.
Kostamovaara
, “
A 4×128 SPAD array with a 78-ps 512-channel tdc for time-gated pulsed Raman spectroscopy
,”
Analog Integr. Circuits Signal Process.
84
,
353
362
(
2015
).
45.
H. K.
Chandrasekharan
,
K.
Ehrlich
,
M. G.
Tanner
,
D. M.
Haynes
,
S.
Mukherjee
,
T. A.
Birks
, and
R. R.
Thomson
, “
Observing mode-dependent wavelength-to-time mapping in few-mode fibers using a single-photon detector array
,”
APL Photonics
5
,
061303
(
2020
).
46.
I. M.
Antolovic
,
S.
Burri
,
C.
Bruschini
,
R. A.
Hoebe
, and
E.
Charbon
, “
Spad imagers for super resolution localization microscopy enable analysis of fast fluorophore blinking
,”
Sci. Rep.
7
,
1
11
(
2017
).
47.
J. P.
O’Reilly
,
N. J.
Kolodziejski
,
D.
McAdams
,
D. E.
Fernandez
,
C. J.
Stapels
, and
J. F.
Christian
, “
A capillary-mimicking optical tissue phantom for diffuse correlation spectroscopy
,”
Proc. SPIE
10056
,
1005613
(
2017
).
48.
L.
Cortese
,
G. L.
Presti
,
M.
Pagliazzi
,
D.
Contini
,
A.
Dalla Mora
,
A.
Pifferi
,
S. K. V.
Sekar
,
L.
Spinelli
,
P.
Taroni
,
M.
Zanoletti
 et al, “
Liquid phantoms for near-infrared and diffuse correlation spectroscopies with tunable optical and dynamic properties
,”
Biomed. Opt. Express
9
,
2068
2080
(
2018
).
49.
P.
Di Ninni
,
F.
Martelli
, and
G.
Zaccanti
, “
Effect of dependent scattering on the optical properties of intralipid tissue phantoms
,”
Biomed. Opt. Express
2
,
2265
2278
(
2011
).
50.
L. A.
Dempsey
,
M.
Persad
,
S.
Powell
,
D.
Chitnis
, and
J. C.
Hebden
, “
Geometrically complex 3D-printed phantoms for diffuse optical imaging
,”
Biomed. Opt. Express
8
,
1754
1762
(
2017
).
51.
L.
Gagnon
,
M.
Desjardins
,
J.
Jehanne-Lacasse
,
L.
Bherer
, and
F.
Lesage
, “
Investigation of diffuse correlation spectroscopy in multi-layered media including the human head
,”
Opt. Express
16
,
15514
15530
(
2008
).
52.
J.
Dong
,
R.
Bi
,
J.-H.
Ho
,
K.
Lee
,
P. S.
Thong
, and
K.-C.
Soo
, “
Diffuse correlation spectroscopy with a fast Fourier transform-based software autocorrelator
,”
J. Biomed. Opt.
17
,
097004
(
2012
).
53.
J. G.
Kim
and
H.
Liu
, “
Investigation of biphasic tumor oxygen dynamics induced by hyperoxic gas intervention: The dynamic phantom approach
,”
Appl. Opt.
47
,
242
252
(
2008
).
54.
W.
Liu
,
R.
Qian
,
S.
Xu
,
P. C.
Konda
, and
R.
Horstmeyer
, “
Fast sensitive diffuse correlation spectroscopy with a spad array
,” in
Optical Tomography and Spectroscopy
(
Optical Society of America
,
2020
), p.
SM3D-3
.
55.
Q.
Fang
and
D. A.
Boas
, “
Monte Carlo simulation of photon migration in 3d turbid media accelerated by graphics processing units
,”
Opt. Express
17
,
20178
20190
(
2009
).
56.
D. A.
Boas
,
S.
Sakadžić
,
J. J.
Selb
,
P.
Farzam
,
M. A.
Franceschini
, and
S. A.
Carp
, “
Establishing the diffuse correlation spectroscopy signal relationship with blood flow
,”
Neurophotonics
3
,
031412
(
2016
).
57.
C.
Funck
,
F. B.
Laun
, and
A.
Wetscherek
, “
Characterization of the diffusion coefficient of blood
,”
Magn. Reson. Med.
79
,
2752
2758
(
2018
).
58.
G. J.
Smits
,
R. J.
Roman
, and
J. H.
Lombard
, “
Evaluation of laser-Doppler flowmetry as a measure of tissue blood flow
,”
J. Appl. Physiol.
61
,
666
672
(
1986
).
59.
D.
Arvidsson
,
H.
Svensson
, and
U.
Haglund
, “
Laser-Doppler flowmetry for estimating liver blood flow
,”
Am. J. Physiol. Gastrointest. Liver Physiol.
254
,
G471
G476
(
1988
).
60.
S. L.
Jacques
, “
Optical properties of biological tissues: A review
,”
Phys. Med. Biol.
58
,
R37
(
2013
).
61.
F.
Martelli
,
T.
Binzoni
,
A.
Pifferi
,
L.
Spinelli
,
A.
Farina
, and
A.
Torricelli
, “
There’s plenty of light at the bottom: Statistics of photon penetration depth in random media
,”
Sci. Rep.
6
,
1
14
(
2016
).
62.
T.
Dragojević
,
J. L.
Hollmann
,
D.
Tamborini
,
D.
Portaluppi
,
M.
Buttafava
,
J. P.
Culver
,
F.
Villa
, and
T.
Durduran
, “
Compact, multi-exposure speckle contrast optical spectroscopy (SCOS) device for measuring deep tissue blood flow
,”
Biomed. Opt. Express
9
,
322
334
(
2018
).
63.
J. W.
Goodman
,
Speckle Phenomena in Optics: Theory and Applications
(
Roberts and Company Publishers
,
2007
).
64.
P.-A.
Lemieux
and
D.
Durian
, “
Investigating non-Gaussian scattering processes by using nth-order intensity correlation functions
,”
J. Opt. Soc. Am. A
16
,
1651
1664
(
1999
).
65.
S. V.
Siddiqui
,
U.
Chatterjee
,
D.
Kumar
,
A.
Siddiqui
, and
N.
Goyal
, “
Neuropsychology of prefrontal cortex
,”
Indian J. Psychiatry
50
,
202
(
2008
).
66.
K.
Morimoto
,
A.
Ardelean
,
M.-L.
Wu
,
A. C.
Ulku
,
I. M.
Antolovic
,
C.
Bruschini
, and
E.
Charbon
, “
Megapixel time-gated spad image sensor for 2D and 3D imaging applications
,”
Optica
7
,
346
354
(
2020
).
67.
F.
Jaillon
,
J.
Li
,
G.
Dietsche
,
T.
Elbert
, and
T.
Gisler
, “
Activity of the human visual cortex measured non-invasively by diffusing-wave spectroscopy
,”
Opt. Express
15
,
6643
6650
(
2007
).
68.
H.
Wabnitz
,
D. R.
Taubert
,
M.
Mazurenka
,
O.
Steinkellner
,
A.
Jelzow
,
R.
Macdonald
,
D.
Milej
,
P.
Sawosz
,
M.
Kacprzak
,
A.
Liebert
 et al, “
Performance assessment of time-domain optical brain imagers, Part 1: Basic instrumental performance protocol
,”
J. Biomed. Opt.
19
,
086010
(
2014
).

Supplementary Material