Fast and sensitive diffuse correlation spectroscopy with highly parallelized single photon detection

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),¹ diffuse optical spectroscopy (DOS),² and functional near-infrared spectroscopy (fNIRs)³ 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.⁴ By examining scattered light’s spectral interference fringe pattern,⁵ the related technique of interferometric near-infrared spectroscopy (iNIRS) can also extract additional path-length-dependent information to accurately measure deep blood flow.⁶ The above techniques are distinct from other measurement strategies such as Doppler OCT,⁷ OCT angiography,⁸ and laser speckle contrast imaging;⁹ 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 strategy^13,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 tissue¹⁰ (Fig. 1). As the mean penetration depth of this PSF typically falls between one-third and one-half of the source–detector separation d¹¹ 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/cm² visible light (at the ANSI safety limit²⁰), a straightforward simulation²¹ 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),²⁶ single-detector implementations of DCS quickly run into an issue of not having enough photons when aiming to detect deep tissue events.

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 channel^27,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 $N$^28,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 humans^30,31 and mice;³² Staples et al. constructed an eight-channel DCS system using one PMT and seven single photon counting modules,³³ 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,²⁷ and Dietsche et al. implemented a 28-channel DCS setup with multiple single-mode fibers to measure blood flow in the deep tissue.²⁸ 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.³⁴ 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,³⁵ fluorescent lifetime imaging microscopy,^36–39 endomicroscopy,⁴⁰ spectroscopy,⁴¹ Raman spectroscopy,^42–44 ultrafast optical imaging,⁴⁵ and localization-based super-resolution imaging,⁴⁶ 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 mm² perturbation that lasts 0.1 s under 1 cm of a tissue-like phantom (μ_a = 0.07 cm⁻¹ and μ_s′ = 3.6 cm⁻¹) 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.

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 skin^47–52 (made of intralipid,⁴⁹ glycerol with fat emulsion,⁴⁸ 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,⁵⁴ 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⁻⁴ s and 10⁻⁵ 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,⁵⁵ 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⁻¹, 7.2 cm⁻¹, 3.6 cm⁻¹, 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.⁶⁰ 

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 mm² × 1.37 mm² to 9.58 mm² × 9.58 mm², 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/cm².²⁰ 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/cm². As the mean penetration depth of collected photons is roughly one-third to one-half of the of the source–detector separation d,¹¹ 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.⁶¹ 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.⁶² 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 d_s generated from a fiber as a function of the distance and fiber core diameter using the following equation:⁶³ 

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. d_s 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:

where I_j(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 g₂ 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,

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 g₂ computation, as we observed that curves decayed to ∼1 when the lag time τ was larger than 200 µs.

The autocorrelation function, g₂(τ), is related to the normalized electric field autocorrelation function, g₁(τ), by the Siegert relation,⁶⁴ 

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

and report t_d defined as the delay time when g₂ decays to 1/e of its initial value [g₂(3 µs)] as the “decorrelation time.” We sample t_d 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⁻¹ 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 g₂(τ) from all M = 1024 SPADs within the array. The results show that with longer autocorrelation curve rates, the g₂(τ) 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.

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 g₂(τ) curve. We compared the averaged g₂(τ) curve that resulted from using 1, 9, and 1024 SPADs to detect a 4.1 mm² × 4.1 mm² perturbation area, presented at 5 kHz, behind a 5 mm thick tissue phantom with μ_s′ = 12.0 cm⁻¹. 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 g₂(τ) 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 g₂ 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.

To quantify the detection sensitivity, we computed a simple classification accuracy metric by alternately measuring the decorrelation time t_d with and without a DMD-induced perturbation and comparing the two adjacent measurements. This alternate measurement method can minimize the disturbance of t_d caused by phantom and experimental environment changes. If the t_d determined with the perturbation present was smaller than the t_d 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 g₂ curves (5 mm thick tissue phantom with μ_s′ = 12.0 cm⁻¹) 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 mm² × 3.42 mm² perturbation via the g₂ 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 mm² × 2.74 mm² and 2.05 mm² × 2.05 mm², respectively.

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 g₂ 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⁻¹, 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⁻¹. With an autocorrelation curve rate of 30 ms, the g₂ curve clearly dips with increasing DMD frequency, and it is easier to differentiate g₂ curves corresponding to different perturbation frequencies for larger perturbation areas, as expected.

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/cm²). 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).

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 g₂(τ) [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 g₂(τ) [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.⁶⁵ 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.⁶ 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 g₂(τ) and the associated decorrelation time t_d 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 t_d 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, t_d, 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.⁶ 

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 available⁶⁶), and recent consumer-oriented developments have led to combined sensing and processing layers via 3D stacking that open up new functionalities.³⁴ 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 dynamics^14,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.⁶⁸ 

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.