2.1.

SFDI Hardware

The hardware design and assembly of the SFDI arm of the system were based on the openSFDI guide²⁵ by Applegate et al.,²⁶ with adjustments to enable continuous acquisition. Since the typical mouse heart rate ranges from 5 to 14 Hz, we designed QHI to achieve raw acquisition frame rate of 50 Hz to capture pulsatile signals using LSI. This frame rate allowed for a slower effective SFDI frame rate of 2.8 Hz for three wavelengths and two spatial frequencies. On the detection arm, a fast CMOS camera (BFS-U3-32S4M-C, FLIR, Oregon, United States) was mounted so that the imaging axis was orthogonal to the sample plane [Fig. 1(a)]. We used a 50 mm fixed focal length lens (67-717, Edmund Optics, New Jersey, United States) at maximum zoom together with a 10 mm spacer (CML10, Thorlabs, New Jersey, United States) to achieve a $\sim 10\times 10\mathrm{mm}$ FOV (magnification = 0.42) at $1200\times 1200\text{pixel}$ resolution. The FOV was chosen to match the dimensions of a mouse cranial window. The working distance was $\sim 130\mathrm{mm}$. A notch filter at 632.8 nm (86-126, Edmund Optics, New Jersey, United States) was used to minimize possible crosstalk with the LSI wavelength [Fig. 1(a)].

Fig. 1

Overall schematic of the QHI system. (a) Hardware setup of QHI with LSI and SFDI components [AL, aspheric lens; AcL, achromatic lens; CL, condenser lens; D, diffuser; DMD, digital micromirror device; F, appropriate notch (for SFDI) and line filters (for LSI) at LSI wavelength, IL, imaging lens with spacers; L, plano-convex lens; LED, light emitting diode; LPDM, long-pass dichroic mirror; P, linear polarizer (cross-polarized between illumination and detection arms); and SPDM, short-pass dichroic mirror]. (b) Triggering scheme of the system where the computer sends a serial trigger to activate an automatic Arduino triggering scheme for the LED and the DMD. The electrical output of the DMD was then used to trigger the camera. The exposure time was set to 10 ms with a 20 ms period per full cycle, allowing for an effective 50 Hz raw frame rate. The Arduino triggering scheme can be deactivated with another serial trigger from the computer.

The SFDI illumination arm was configured with an illumination axis at $\sim 24\mathrm{deg}$ with respect to the imaging arm [Fig. 1(a)] to minimize distortion artifacts (e.g., keystone artifacts) on projected SFDI patterns within the FOV. All optical components were set up at 190 mm away from the mounting breadboard to provide ample space for animal imaging. The projection system included three LEDs at 660, 780, and 850 nm (M660L4, M780LP1, and M850L3, respectively; Thorlabs, New Jersey, United States) coupled into a digital micro-mirror device (DMD) projector (DLP4500 NIR, Keynote Photonics, Texas, United States) using appropriate dichroic mirrors (FF697-SDi01-25x36, Semrock, New York, United States, and DMSP805R Thorlabs, New Jersey, United States) and condenser lenses (ACL2520U-B, Thorlabs, New Jersey, United States) [Fig. 1(a)]. The three wavelengths were chosen to enhance sensitivity to the chromophores of interest (i.e., HbR and ${\mathrm{HbO}}_2$) based on their absorption spectrum and to approximately match the probing depth of the LSI wavelength (633 nm). To achieve sufficient light intensity for fast acquisition (i.e., 50 Hz), we added a plano-convex lens (LA1708-AB, Thorlabs, New Jersey, United States) in front of the combined light beam before the DMD [Fig. 1(a)]. The plano-convex lens position was carefully adjusted to avoid creating a focused image of the LED on the DMD. A 400 to 1100 nm achromatic lens (AC254-075-AB-ML, Thorlabs, New Jersey, United States) was positioned $\sim 100\mathrm{mm}$ from the DMD so that a focused image of roughly 3× magnification of the DMD’s display was formed on the focal plane. We then made fine adjustments to ensure that the illumination and detection arms were confocal using projection alignment steps listed in the openSFDI guide.²⁷ A barrier was then added to block unwanted diffraction patterns from the DMD from the imaging axis. We controlled the LEDs using dedicated drivers (LEDD1B, Thorlabs, New Jersey, United States) with triggering signals from an Arduino Uno (1050-1024-ND, Digi-Key Electronics, Minnesota, United States). Details regarding the electrical configuration are discussed in the Acquisition Software section below.

Finally, linear polarizers (LPNIRE100-B, Thorlabs, New Jersey, United States) were added to achieve crossed polarization between the illumination and detection arms while maintaining appropriate intensity counts (i.e., 40000 to 45000 counts for 16-bit images) from a static silicone phantom of known optical properties (${\mu }_a=0.019$, 0.019, and $0.019{\mathrm{mm}}^{-1}$, ${{\mu }_s}^{\prime }=1.07$, 0.83, and $0.72{\mathrm{mm}}^{-1}$ at 660, 780, and 850 nm, respectively).

2.2.

LSI Hardware

For the LSI detection arm, an identical CMOS camera and lens-spacer system to that of the SFDI arm were used to achieve a matching FOV. We combined the two detection arms using a 638 nm long pass dichroic mirror (DMLP638R, Thorlabs, New Jersey, United States) [Fig. 1(a)]. For LSI, we adjusted the imaging lens aperture to a f/# of 4 to achieve a proper speckle-to-camera pixel size ratio (approximately $9\mu \mathrm{m}$ to $3.45\mu \mathrm{m}$, respectively).²⁸ A 632.8 nm laser line filter (FL632.8-3, Thorlabs, New Jersey, United States) was used to minimize crosstalk with illumination from the SFDI arm [Fig. 1(a)].

For the LSI illumination arm, a long-coherence 633 nm laser (115-81059-554, Coherent, California, United States) was coupled with an aspheric lens (C330TMD-B, Thorlabs, New Jersey, United States) to quickly expand the beam. A 1500-grit diffuser (DG-10-1500-MD, Thorlabs, New Jersey, United States) was used to homogenize the beam. The laser was then mounted to achieve homogeneous illumination over the camera’s FOV with the beam concentric with the FOV [Fig. 1(a)]. We adjusted the laser’s illumination power using a polarizer (LPNIRE100-B, Thorlabs, New Jersey, United States) to achieve the appropriate intensity count (i.e., 60-70 for 8-bit images) on the same static silicone phantom used for the SFDI setup. Similar to the SFDI arm, cross-polarization with linear polarizers was used to minimize specular reflectance [Fig. 1(a)].

2.3.

Computer Specification and Electrical Wiring

The performance and effective frame rate of the system also depend on the specifications of the acquisition computer. Here, we used a laptop with Intel^® Core™ i7-9750H CPU (12 CPUs) processor, 16 GB RAM, and a 1TB built-in M.2. solid state drive (SSD). The LED drivers were connected to the Arduino triggering digital pins using BNC cables with the ground leads connected to the Arduino GND pin. On the DMD, the “TRIG IN” connector J7 was connected to the Arduino DMD trigger pin and the “TRIG OUT” connector J11 to the cameras’ trigger-in pin through GPIO cables. All ground outputs from the DMD and the cameras were connected to the Arduino GND pin. Micro-B-to-A-male USB-3.0 cables were used to connect the cameras and the computer. One B-to-A male USB-2.0 cable was used to connect the Arduino to the computer.

4.1.

SFDI Processing Pipeline

We used MATLAB (R2019b, MathWorks, Massachusetts, United States) to process SFDI data following a previously published protocol.¹⁹ The first step of the processing scheme was pixel-wise spatial frequency demodulation using the three phases collected to obtain the AC signal of the reflectance

The next step of processing included an inverse lookup table (LUT) scheme to match the measured pairs of $R_d$ to a pair of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ values. Here, the LUT for $R_d(f_x)$ at $f_x=0$ and $0.3{\mathrm{mm}}^{-1}$ was generated using the same WMC code mentioned above with ${\mu }_a$ range of 0.0001 to $0.2{\mathrm{mm}}^{-1}$ and ${{\mu }_s}^{\prime }$ range of 0.001 to $3{\mathrm{mm}}^{-1}$ for 512 linearly spaced values. We employed a rapid LUT method using an interpolation scheme for scattered data from a built-in MATLAB function (scatteredInterpolant) to generate an inverse LUT for a pair of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ at each pixel.

After the extraction of ${\mu }_a$ and ${{\mu }_s}^{\prime }$, we calculated chromophore concentrations using ${\mu }_a$ at all three wavelengths. Here, we calculated an inverse matrix of extinction coefficients at each wavelength for ${\mathrm{HbO}}_2$ and HbR, which was then multiplied with ${\mu }_a$ at each pixel. In Matlab, this method is analogous to least square fitting when the matrices have mismatched dimensions. Finally, we used HbR and ${\mathrm{HbO}}_2$ concentration to calculate physiologically relevant parameters, such as HbT and ${\mathrm{StO}}_2$ using their definitions

4.2.

LSI Processing Pipeline

We employed a typical spatial speckle contrast algorithm with a 5x5 sliding window for analysis.⁹ We used MATLAB to perform the computations using the following equation:

For each LSI measurement, $\beta$ is a normalization factor that takes into account the mismatch between detector and speckle sizes and polarization. Typically, $\beta$ can be estimated using speckle contrast $K$ (where ${\tau }_c\gg T$) measured from a static phantom as $\beta \approx K_{\text{phantom}}^2$.⁹ We then use $\beta$ to correct for SFI

4.3.

Obtaining Absolute Flow Values by Correcting for Optical Properties

In this section, we describe in detail a protocol used to obtain an absolute measurement of flow by correcting speckle contrast $K$ with optical properties. We first started with an equation relating $K$ to electric field autocorrelation $G_1$ as described by Mazhar et al.¹⁵

The associated solution for $G_1$ in the spatial frequency domain is

In Eq. (9), ${\mu }_{tr}={\mu }_a+{{\mu }_s}^{\prime }$ is the transport coefficient, and $P_0$ is the illumination power. Here, $A$ was defined as

In the spatial frequency domain, ${{\mu }_{\mathrm{eff}}}^{\prime }$ is defined as

In this equation, $f_x$ is the spatial frequency of the illumination pattern; $n$ is the refractive index; $\lambda$ is the LSI wavelength; and $D_b$ is the “effective” Brownian diffusion coefficient describing the motion of scatterers.³³ This equation is often used in dynamic light scattering theory to describe the mean-squared displacement of moving particles. Furthermore, the Brownian model yielded a better fit for the observed correlation decay curves in various tissue types.³⁴ In this work, we assumed that LSI illumination is planar; thus, $f_x$ is assumed to be 0.

Our goal was to retrieve $D_b$ using ${\mu }_a$, ${{\mu }_s}^{\prime }$, and speckle contrast $K$. The overall schematic of deriving $D_b$ from measured parameters is described in Fig. 2. First, we used concentrations of ${\mathrm{HbO}}_2$ and HbR determined with SFDI to calculate ${\mu }_a$ at 633 nm. Then, to determine ${{\mu }_s}^{\prime }$ at 633 nm, we fit an exponential function in the form of the Mie-scattering exponential decay function³⁵

Fig. 2

Flow chart of the method of obtaining absolute flow value using the Brownian diffusion coefficient $D_b$ from measured $K$, ${\mu }_a$, and ${{\mu }_s}^{\prime }$. Here, a $\beta$ factor was approximated using a static silicone phantom.

To perform the fit more efficiently, we rewrote the decay function, as shown in Eq. (12). We then solved for $\log (A)$ and $b$ using MATLAB “\” operator using the ${{\mu }_s}^{\prime }$ measured at 660, 780, and 850 nm. From $\log (A)$ and $b$, we obtained a value of ${{\mu }_s}^{\prime }$ at 633 nm with a forward solution. We used a $\beta$ factor to create a LUT of $K^2$ values as a function of $D_b$, ${\mu }_a$, and ${{\mu }_s}^{\prime }$ by numerically integrating Eq. (8). The values of ${\mu }_a$ and ${\mu }_s\mathrm{n}$ in the LUT were chosen to cover the range that included the values at 633 nm (i.e., ${\mu }_a$ from 0.0001 to $0.3{\mathrm{mm}}^{-1}$ and ${{\mu }_s}^{\prime }$ from 0.001 to $3{\mathrm{mm}}^{-1}$ with 128 equally spaced values), and the range of $D_b$ was set to 0 to ${20}^{-6}{\mathrm{mm}}^2/\mathrm{s}$ with 256 equally spaced values. Using the MATLAB function scatteredInterpolant, we then obtained the value of $D_b$ from the set of measured $K$, ${\mu }_a$, and ${{\mu }_s}^{\prime }$ at 633 nm.

5.1.

Flow Sensitivity Using LSI with a Flow Phantom Setup

We used a flow phantom setup to characterize the flow sensitivity of the LSI arm of QHI as previously described.^36,37 The static part of the phantom was made using silicone with no added absorbers and ${\mathrm{TiO}}_2$ as scatterers. The optical properties of the static portion were measured using SFDI, with ${\mu }_a=0.0049$, 0.0059, and $0.0062{\mathrm{mm}}^{-1}$ and ${{\mu }_s}^{\prime }=0.87$, 0.67, and $0.60{\mathrm{mm}}^{-1}$ at 660, 780, and 850 nm, respectively. A $25\mu \mathrm{L}$ glass capillary pipet tube (Kimble 71900, inner radius $\sim 330\mu \mathrm{m}$) was embedded in the silicone phantom so that the tube was flush with the surface of the phantom [Fig. 3(a)]. A 1% intralipid solution was prepared and infused into the capillary tube using a precision pump (70-4504, Harvard Apparatus, Massachusetts, United States) at incremental speeds of 0, 0.5, 1, 2, 3, 4, and $5\mathrm{mm}/\mathrm{s}$ for 60 s at each speed. The imaging session lasted 7 min in total. We repeated the experiment for 5 trials, with 15 min of wait time between each trial. This wait time is implemented to ensure that the flowing fluid is restored to a baseline state.

Fig. 3

Flow phantom experiment to characterize the sensitivity of the LSI arm. (a) A schematic of the silicone flow phantom with an embedded glass capillary pipet tube where a solution of 1% intralipid was flowed at different speeds using a syringe pump. A representative flow map of the captured FOV was also shown with the chosen ROI (black rectangles) for static and flow regions. (b) Time-course data of the experiment where the flow speed was increased at each 1-min interval (top table). The time course of ${\mathrm{SFI}}_{\beta \text{-corrected}}$ within the flow region showed an increasing pattern with increasing speed between 0 and $\sim 4\mathrm{mm}/\mathrm{s}$. A plateau pattern was observed at higher speeds ($>4\mathrm{mm}/\mathrm{s}$). For the static region, ${\mathrm{SFI}}_{\beta \text{-corrected}}$ showed no change during the experiment.

Two regions of interest (ROIs) were chosen over the dynamic and static region of the phantoms and median values of ${\mathrm{SFI}}_{\beta \text{-}\text{corrected}}$ were calculated at each imaging time point and filtered using a 50-frame ($\sim 1\mathrm{s}$) moving average filter. The resulting time course data sets are shown below in Fig. 3(b). Here, the dynamic ROI showed the expected increase in ${\mathrm{SFI}}_{\beta \text{-}\text{corrected}}$ values with increasing speed, with the sensitivity to flow reaching a plateau starting at $\sim 4\mathrm{mm}/\mathrm{s}$, similar to what has been reported previously with $T=10\mathrm{ms}$.³¹ The static ROI showed no change in SFI despite the increasing speeds.

5.2.

Optical Property Accuracy Assessment of SFDI

We fabricated four static silicone phantoms and compared optical property measurements between the SFDI arm of QHI with a commercial SFDI system. The phantoms were made using a base of silicone, ${\mathrm{TiO}}_2$ powder as scatterers, and India ink as absorbers. One phantom had a specific set of ${\mu }_a$ and ${{\mu }_s}^{\prime }$, and the other three phantoms were fabricated to have the following optical property values relative to the first phantom: (1) same ${\mu }_a$, $2\times {{\mu }_s}^{\prime }$; (2) same ${{\mu }_s}^{\prime }$, $2\times {\mu }_a$; and (3) $2\times {\mu }_a$ and $2\times {{\mu }_s}^{\prime }$. Here, the values of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ measured using QHI were compared with measurements from a commercial SFDI system (ReflectRS, Modulim, California, United States). A common circular ROI with a diameter of 3.3 mm was selected and median values of the optical properties were extracted. We performed a linear regression analysis using the median ${\mu }_a$ and ${{\mu }_s}^{\prime }$ values from the four phantoms measured using QHI and the ReflectRS for 660 and 850 nm. This wavelength choice resulted from the fact that the ReflectRS does not have a 780 nm LED. Here, we observed a high linear correlation between QHI and the ReflectRS ($R^2>0.98$) for both ${\mu }_a$ and ${{\mu }_s}^{\prime }$ (Fig. 4). Using a Bland–Altman plot, we observed a mean difference of $-0.0059$ and $0.0316{\mathrm{mm}}^{-1}$ for ${\mu }_a$ and ${{\mu }_s}^{\prime }$, respectively (Fig. 4). Collectively, these data suggest that QHI slightly underestimated ${\mu }_a$ and slightly overestimated ${{\mu }_s}^{\prime }$; however, this discrepancy was linear for various combinations of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ (as shown through high $R^2$). One reason for this discrepancy could be the mismatched FOV of the two systems, as QHI had a much smaller total FOV (i.e., $1\times 1\mathrm{cm}$ versus $20\times 15\mathrm{cm}$, respectively, $\sim 300\times$ smaller) due to its higher magnification, with an increased sensitivity to any heterogeneities within the phantom (e.g., locally dense regions of ${\mathrm{TiO}}_2$), which could skew the median value.

Fig. 4

Correlation (top) and Bland-Altman plots (bottom) between QHI and a commercial ReflectRS system from Modulim of four static silicone phantoms. Here, we compared the median values obtained from the chosen ROI. Only 660 nm and 850 nm wavelengths were compared, as 780 nm is not available with the ReflectRS. For the correlation plot, $R^2$ values were also listed for the linear fit.

5.3.

Intra-Day Inter-Session Precision Measurements of Static Silicone Phantom

To investigate the stability of the system during a session and among multiple sessions within one day of imaging, we imaged phantom 1 in six separate trials. Here, each trial lasted 5 min, and the time between the start of each trial was 30 min. The total length of the experiment was 3 h. A circular diameter with a radius of 200 pixels ($\sim 1.6\mathrm{mm}$) was chosen at the center of the FOV and median values of ${\mu }_a$, ${{\mu }_s}^{\prime }$, and $K$ were recorded at each respective wavelength at each time point. A whisker bar graph of all median values for the five minutes of imaging for each measured optical parameter at each respective wavelength is shown in Fig. 5. We also calculated the coefficient of variation across all trials using the mean value of each trial’s time course. For SFDI measurements of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ of all wavelengths, the standard deviations for each of the trials were $<0.001$ while the coefficient of variation across all trials was $<0.5\%$. For LSI measurements of $K$, we observed standard deviations ranging from 0.0014 to 0.0013 for each of the trials with a coefficient of variation across all trials of 0.16%. These values were listed in Table S1 in the Supplementary Material and indicated good intraday stability of the system.

Fig. 5

Measurements of QHI precision were performed for six trials (5 min each) within the same day on the same static silicone phantom (i.e., phantom 1). The mean and standard deviation are plotted here. The recorded values for ${\mu }_a$, ${{\mu }_s}^{\prime }$, and $K$ and the calculated coefficient of variation across all trials were shown in Table S1 in the Supplementary Material.

5.4.

Longitudinal Precision Measurement of Static Silicone Phantom

Using the four static silicone phantoms listed above, we conducted longitudinal precision measurements over 5 days. A single full set of LSI and SFDI frames (a total of 18 frames) was taken and analyzed. A circular ROI with a radius of 200 pixels ($\sim 1.6\mathrm{mm}$) was chosen at the center of the FOV and median values of ${\mu }_a$, ${{\mu }_s}^{\prime }$, and $K$ were determined at each respective wavelength. These median values across the 5 days were then used to calculate a coefficient of variation for each measured parameter. Here, we noted that all coefficients of variations were $<4\%$ (Fig. 6, Tables S2–S4 in the Supplementary Material), indicating good precision of the system across multiple days. Here, $K$ will theoretically be 1 for a completely static phantom. However, the measured $K$ from the four phantoms ranged between 0.5 and 0.7 due to potential environmental, experimental, or instrumental sources. For phantom 2, $K$ increased with increasing ${\mu }_a$ values compared to those of phantom 1 due to the decreasing probability of detection of a photon hitting a dynamic scattering event (i.e., a dynamic photon). On the contrary, phantom 3 showed a decrease in contrast $K$ with increasing ${{\mu }_s}^{\prime }$ due to the increasing probability of detection of a dynamic photon. Finally, phantom 4 showed similar $K$ values to those of phantom 1 when both increasing ${\mu }_a$ and ${{\mu }_s}^{\prime }$. The measured $K$ values differed across the four phantoms due to ${\mu }_a$ and ${{\mu }_s}^{\prime }$ changes as expected from Eqs. (8)–(11) (Fig. 6, red dotted line for $K$ graphs).

Fig. 6

Measurements of QHI precision over five days. The median values within the selected ROI for each day are plotted and further used in the analysis. The recorded values for ${\mu }_a$, ${{\mu }_s}^{\prime }$, and $K$ and the calculated coefficient of variation across all five days are shown in Tables S2–S4 in the Supplementary Material. $K$ is higher for the phantom with higher ${\mu }_a$ (phantom 1 versus phantom 2) and lower for the phantom with higher ${{\mu }_s}^{\prime }$ (phantom 1 versus phantom 3). With higher ${\mu }_a$ and ${{\mu }_s}^{\prime }$ (phantom 4), $K$ is similar to that of phantom 1. This agrees with the findings of Mazhar et al.¹⁵

6.1.

Intra-Animal CBF Measurement During Hypercapnia Gas Challenge

To assess the impact of the processing methods on dynamic measurements, we subjected a female mouse (C57BL/6, 5 months old) to a cranial window preparation for hypercapnia. Anesthesia was initiated at 4% isoflurane using an induction chamber and the animal was maintained at 2% for the surgical period at $1\mathrm{L}/\min$ flow rate. For the gas challenge, we used 1.5% isoflurane to maintain anesthesia and the following experimental protocol: 5 min of baseline with balance room air, 4 min of 5% ${\mathrm{CO}}_2$ and balance room air, followed by 6 min of the initial isoflurane-room air condition. We extracted the time course of SFDI and LSI data using a semicircular ROI on the left hemisphere (radius of 175 pixels, $\sim 1.4\mathrm{mm}$). The ROI was chosen with its position along the midpoint between bregma and lambda and an offset of 100 pixels ($\sim 0.8\mathrm{mm}$) from the midline. The representative CBF maps of $D_b$ and SFI with the chosen ROI can be seen in Fig. 7(a). We then normalized $D_b$ and SFI measurements to time $t=0$ and obtained a percentage ΔCBF. Due to normalization, $\beta$ correction is not required for SFI measurements. The time courses were then filtered using a moving average of 200 data points, approximating a 4-s period. In this experiment, as expected we observed a hyperemic response during the hypercapnic period [Fig. 7(b)]. However, there was a slight difference between the filtered SFI and $D_b$ time courses [Fig. 7(c)]. This difference ranged from $-1.5$ to 0.5% during the experimental period, corresponding to the range of changes of $\sim 0.005{\mathrm{mm}}^{-1}$ ($\sim 7\%$) in ${\mu }_a$ and $\sim 0.025{\mathrm{mm}}^{-1}$ ($\sim 2\%$) in ${{\mu }_s}^{\prime }$ [Figs. 7(c) and 7(d)]. This decrease in ${\mu }_a$ and ${{\mu }_s}^{\prime }$ resulted in a minimal change in CBF dynamics between ${\mathrm{D}}_b$ and SFI measurements. In fact, Eq. (9) predicts that ${\mathrm{K}}^2$ remains constant with a defined $D_b$ if the changes of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ happen with the same multiplicative factor. For example, if ${\mu }_a$ and ${{\mu }_s}^{\prime }$ both increase by 50%, $K^2$ stays the same. Intuitively, the loss of correlation due to scattering cancels the increase of correlation from the loss of the traveling photon due to absorption. However, in the case where only one of the two optical properties remains constant while the other is changing, $K^2$ will change accordingly. A simulation result of this phenomenon was shown in Fig. S1 in the Supplementary Material.

Fig. 7

CBF dynamics slightly differentiate between $D_b$ and SFI measurements for a mouse undergoing a hypercapnia gas challenge. (a) Representative images of $D_b$ and SFI during hypercapnia with the chosen semicircular ROI (black) on the left hemisphere. Scale bars are 1 mm. b) CBF measurements normalized to $t=0$ where both SFI and $D_b$ analysis was performed. Median data are chosen from the left hemisphere to extract the time course of SFI and ${\mathrm{D}}_b$ values. The data is further filtered using a moving filter of 200 frames ($\sim 4\mathrm{s}$). (c) A difference of $-1.5$ to 0.5% between the normalized $D_b$ and SFI curves is observed throughout the gas challenge experiment. (d) Optical properties derived at 633 nm showing a slight decrease in both ${\mu }_a$ and ${{\mu }_s}^{\prime }$ across the time course.

6.2.

Inter-Animal CBF Measurement at Resting State

We measured the resting-state CBF of a cohort of mice ($n=13$, B6SJL background, 15 months old). Anesthesia was initiated using 4% isoflurane (balance oxygen) within an enclosed chamber and the animals were maintained at 2% isoflurane (balance room air) at $1\mathrm{L}/\min$ flow rate. We then performed an acute intact-skull cranial window surgical procedure to gain access to the cortical layer of the brain.³⁸ The scalp was shaved, depilated, and sterilized with alcohol and povidone-iodine wipes. We then performed a midline incision, retracted the scalp, and removed the fascial layer to expose the relatively transparent skull. The window was then filled with saline and a coverslip was applied to create a translucent skull for imaging the cortical surface.

Resting-state data were taken for one minute at 2% isoflurane and a common semi-circular ROI (radius of 175 pixels, $\sim 1.4\mathrm{mm}$) on the left hemisphere was chosen [Fig. 8(a)]. The ROI was chosen with its position along the midpoint between bregma and lambda and an offset of 100 pixels ($\sim 0.8\mathrm{mm}$) from the midline. To calculate SFI and $D_b$, we chose the median value within the ROI to generate time course data for all measured parameters (e.g., $K$, ${\mu }_a$, and ${{\mu }_s}^{\prime }$). We applied the processing methods discussed above to generate optical properties at 633 nm and convert SFI to $D_b$ [Fig. 8(a)]. Here, we only applied the processing methods on the 1-min time-course data and not on every pixel within the FOV.

Fig. 8

Baseline CBF measurements from 13 mice (2% isoflurane for 1 min) highlighted the effect of optical properties on the interpretation of CBF measurements using LSI. (a) Median values of data obtained from the left hemisphere were used to generate time course data for the 1-min imaging period. Both the pulsatile CBF waveform and breathing-related motion artifacts could be resolved. On the representative image on the right, the scale bar is 1 mm. (b) Mean values of the SFI time course during each imaging period were then used to rank the animals in ascending order (blue to red circles). Animal IDs (1 to 13) were then assigned based on this ranking. Here, the $D_b$ values differed from SFI in ranking for nine out of 13 animals.

To assess the impact of optical property correction on CBF, we calculated the mean of the time courses of SFI and $D_b$. An initial ranking of resting-state CBF was generated using SFI values [Fig. 8(b), blue to red for lowest to highest SFI measurements]. The animal IDs were also assigned using this ordering (i.e., animal 1 was the one with the lowest SFI). To facilitate a more direct comparison of SFI and $D_b$, we normalized the individual SFI and $D_b$ values for all animals to those of animal 1. Due to normalization, $\beta$ correction is not required for SFI measurements. Using this comparison method, we observed that the ranking of CBF using $D_b$ differs from that using SFI. Specifically, 9 out of 13 animals showed a change in the ranking [Fig. 8(b)], thus highlighting the potential effect of optical properties on speckle-based measurements of CBF.

Figure 9(a) showed representative normalized CBF images that were generated using a pixel-wise correction algorithm. All pixels in the $D_b$ and SFI images were normalized to the respective values of the center pixel in each image. We also calculated a percentage difference between normalized $D_b$ and SFI maps [Fig. 9(a)] and noted higher CBF values in larger vessels (red arrows) in the $D_b$ image. Currently, this pixel-by-pixel method of correction is currently not applied to time-course analysis due to the long processing time for images with such high pixel resolution (1.44 MP). We further investigated if these spatial differences would also yield a change in temporal dynamics. On the same representative animal, we used the percentage difference map to define three groups of pixels within the chosen ROIs mentioned in Fig. 8. These three groups correspond to the pixels that showed a decrease [$<-5\%$, Fig. 9(b) blue], minimal change [$-5$ to 5%, Fig. 9(b) green], or increase [$>5\%$, Fig. 9(b) red] in the spatial percentage difference map. We then calculated rCBF as relative CBF from $t=0$ for $D_b$ and SFI time courses using the median value of all the pixels within each respective region. As shown in Fig. 9(b) below, the $D_b$ and SFI time courses showed similar dynamics at the resting state for all regions. The power spectra, not shown here, did not yield notable differences in the pulsatile signal.

Fig. 9

On one representative animal, $D_b$, when compared to SFI map, showed a spatial difference in CBF measurements but not in temporal dynamics. (a) Representative normalized maps of $D_b$ and SFI highlight the spatial difference in estimated CBF. Individual pixels in each image were normalized relative to the respective center pixel (*). A percentage difference map was generated using $\%\text{Difference}=(D_{b,\text{norm}}-{\mathrm{SFI}}_{\text{norm}})/{\mathrm{SFI}}_{\text{norm}}\times 100$. $D_b$ showed increased blood flow in large vessels (red arrows) compared with SFI. Scale bars are 1 mm. (b) Within the chosen semicircular ROI, temporal dynamics (i.e., rCBF) did not show notable differences between SFI and $D_b$ measurements among regions with decreasing ($<-5\%$), minimally changing ($-5$ to 5%), and increasing ($>5\%$) spatial differences. The power spectrum of both $D_b$ and SFI (not shown here) showed negligible differences in the pulsatile waveforms.

7.1.

Current Hardware Limitations

The system has a 50 Hz raw frame rate at $1200\times 1200\text{pixel}$ resolution of a co-aligned $10\mathrm{mm}\times 10\mathrm{mm}$ FOV. This frame rate enabled an effective 50 Hz and 2.8 Hz frame rate for LSI and SFDI arms, respectively. Using in vitro phantoms, for sensitivity and accuracy measurements, we showed that QHI can detect differences in LSI flow speed from 0 to $4\mathrm{mm}/\mathrm{s}$ while also achieving a strong linear correlation of optical properties (i.e., ${\mu }_a$ and ${{\mu }_s}^{\prime }$) measurements with a commercial SFDI system ($R^2>0.95$). We observed a slight underestimation of ${\mu }_a$ and a slight overestimation of ${{\mu }_s}^{\prime }$ between QHI and the commercial system.

The openSFDI guide suggested that the exposure time of the camera should be multiples of the DMD’s exposure time to avoid demodulation artifacts.²⁶ However, despite matching our camera’s exposure time to the DMD’s exposure time (i.e., 10 ms), we did not notice significant demodulation artifacts after performing phantom calibration. We achieved this by ensuring that the trigger signal from the DMD to the camera was sent $11\mu \mathrm{s}$ before the pattern display to match the minimum $11\mu \mathrm{s}$ trigger delay time of the camera used in the current build. And as noted above, the fast frame rate at high pixel resolution may result in dropped frames during acquisition. This could be due to the bottleneck of the writing speeds of the internal SSDs. A potential solution would be saving LSI and SFDI data into separate SSDs to maximize the writing speed.

The flow phantom validation experiment showed that the system is sensitive up to $4\mathrm{mm}/\mathrm{s}$ of directed flow speed for 1% intralipid solution. In live mouse brains, researchers have used fluorescence microscopy to measure absolute flow speed on a vascular level, where the speed of moving red blood cells ranged from 0 to $1\mathrm{mm}/\mathrm{s}$ in capillaries, 1 to $10\mathrm{mm}/\mathrm{s}$ in venules, and 5 to $20\mathrm{mm}/\mathrm{s}$ in arterioles.⁴⁵ The dynamic range of our system presents a limitation in quantitative measurements of venules and arterioles. Thus, we recommend choosing the median value of a large ROI during ROI selection and processing, such as the chosen semicircles as shown above, as the quantitative factor. This median value typically represents the parenchymal regions of the brain, where capillaries dominate.

Also, the described method of deriving ${\mu }_a$ and ${{\mu }_s}^{\prime }$ at the LSI wavelength from values obtained using SFDI can be prone to errors. This would then affect the deduced $D_b$ values. One solution to this problem is to perform LSI and SFDI sequentially, where LSI frames are interlaced with SFDI frames. However, this method would significantly increase the time required for the acquisition or reduce the exposure time to accommodate for a faster acquisition speed. Decreasing exposure time can increase the detectable noise levels and also affect the spatial speckle contrast measurements when the exposure time starts to approximate the correlation time (i.e., the flow would appear more static at a faster exposure time). Additionally, the sequential triggering of SFDI frames limits the SFDI frame rate to one for every 18 LSI frames. This leads to a temporal mismatch in the correction method, where 18 LSI frames, which can span the pulsatile signal, have the same optical properties. This limitation poses a potential explanation for the lack of temporal dissimilarity between $D_b$ and SFI, as shown in Fig. 9(b) above. Furthermore, the mismatched wavelengths raised the question of whether the two modalities probe similar tissue volumes. Here, we calculated the estimated effective depth penetration of the SFDI system using the derived optical properties at 633 nm for the 13 mice measured at resting state. The depth penetration ranges from 0.5 to $\sim 2\mathrm{mm}$ depending on the spatial frequencies (Table S5 in the Supplementary Material). It is acknowledged that the penetration depth for SFDI falls within this specific range. This range overlaps with the typical penetration depth of an LSI measurement.⁴⁶ However, the exact probing volumes remain dissimilar. To address these concerns, there have recently been reports of utilizing speckle patterns as a function of spatial frequency to extract optical properties.^47–49 However, to the best of our knowledge, a method to decouple the effect of blood flow from optical properties on speckle contrast measurements has not been reported.

7.2.

Limitations of the Homogeneous, Semi-Infinite Tissue Assumption

A limitation of the current processing method is the employment of a semi-infinite homogeneous model for SFDI data. Furthermore, the solution from Mazhar et al. was derived for the same tissue geometry.¹⁵ Thus, the current correction scheme does not include the effects of volumetric heterogeneities such as the skull on the measurements. To improve the accuracy of QHI measurements of CBF, the skull should be modeled as a scattering layer with minimal absorption and a fixed $D_b$ value that resembles that of a static phantom.

Furthermore, it should be noted that the equations used in this publication were derived from dynamic light scattering theory with an assumption that the photons encounter multiple scattering events before being detected. This is typically true for the parenchymal region where small capillaries act as the main source of dynamics, resulting in a more diffusive behavior typically observed from multiple and unordered scattered photons.⁵⁰ However, it was also shown that larger vessels follow a more single and ordered scattered regime. Finally, the authors also noted a mixed dynamic (i.e., multiply and ordered scattering, and single and unordered scattering) in the immediate regions surrounding these vessels. Thus, our assumption of multiple and unorder scattering potentially results in the mischaracterization of the true dynamics of the larger vessels and their surrounding regions. Nevertheless, this assumption should be appropriate for the parenchymal region, where capillaries dominate.

7.3.

Understanding the Diffusion Coefficient and its Relevance in CBF Measurements

The diffusion coefficient $D_b$ is used to describe the mean squared displacement of the moving particles using the Brownian diffusion model. In diffuse correlation spectroscopy, a fiber-based analog of LSI, $D_b$ has been used to quantify CBF with assumed optical properties.^51–53 Researchers have found that this Brownian model better fits the electric field autocorrelation curves compared to a random flow model when measuring the movement of red blood cells inside tissues. A potential explanation is that, in tissue, red blood cells do not simply move from one position to another but also experience shear stresses, causing them to “roll and tumble” in the vasculature.³⁴ This is particularly true for the parenchymal region of the brain, where the imaging system is mainly probing small capillaries. From a macroscopic point of view, this highly randomized motion exhibits behaviors similar to Brownian motion. However, the typical values of $D_b$ for blood flow are magnitudes higher than that of pure thermally induced Brownian motion. As for large vessels, as discussed above, a more rigorous model should be utilized to correctly account for the more single-directional motion of the red blood cells.

The unit of $D_b$ is area over time (e.g., ${\mathrm{mm}}^2/\mathrm{s}$). This unit differs from the clinically used units of CBF such as perfusion rate per tissue volume per time (e.g., mL/g/min). However, these clinical measurements are often obtained using a volumetric scanner, such as MRI or PET. This is different from an en-face widefield optical imaging measurement, such as LSI. While LSI can probe CBF at deeper brain regions, this information is usually integrated into the superficial cortical measurements. Hence, each LSI pixel represents mainly the motion detected at the surface and partial subsurface contributions. Both contribute to the blurriness of the speckle pattern at the analyzed pixel. Thus, the unit of $D_b$ (i.e., area per time) is appropriate to measure such dynamics. In fact, it has been shown that LSI measurements of directed motions do not measure volumetric flow directly but a combination of the channel size and the speed of the particles (i.e., mostly single-directional area flux).⁵⁴ However, the work in that study specifically investigated the flow dynamics of large vessels, ranging from 20 to $100\mu \mathrm{m}$ for in-vivo and 65 to $175\mu \mathrm{m}$ for in-vitro demonstrations and is not directly applicable to measurements of the parenchymal regions as performed here. Future work is warranted to fully understand which model works best for each selected vascular region of the cortex.

Finally, we are interested in leveraging $D_b$ as a non-unitless measurement of CBF in calculating the cerebral metabolic rate of oxygen (${\mathrm{CMRO}}_2$). In studies combining LSI with spectroscopic imaging, relative measurements of ${\mathrm{rCMRO}}_2$ can be estimated with rCBF and hemoglobin concentrations using certain assumptions.^5,13,27 A new model utilizing ${\mathrm{D}}_b$ as a quantitative CBF surrogate will allow for a direct calculation of ${\mathrm{CMRO}}_2$ instead of a relative value. However, this new method requires changes to the typical ${\mathrm{CMRO}}_2$ equation (i.e., Fick’s law) due to the mismatched dimensions and will be the subject of future work.