2.1.

Hybrid DSCT and LSCI System

The DSCT acquires boundary data using the same hardware as the scDCT [Fig. 1(a)], which has been reported extensively in our previous publications.^{19–25,32,33} Briefly, a high-speed scanning galvanometer mirror positioning system (maximum scanning angle: $\pm 12.5\mathrm{deg}$; switching time: $<1\mathrm{ms}$; GVS002, Thorlabs) remotely and sequentially delivered coherent focused-point NIR light (780 nm, CrystaLaser) to multiple source positions on the region of interest (ROI) [Figs. 1(b) and 1(c)]. The point light was focused on the tissue surface via an achromatic lens (AC127-019-B, Thorlabs), and the light spot size was adjusted and minimized using a lever-actuated iris (SM05D5, Thorlabs). A sCMOS camera (pixels: $2048\times 2048$; frame rate: 30/s; quantum-efficiency: 50% @800 nm; ORCA-Flash4.0, Hamamatsu Photonics) was used as a high-density 2D detector array to collect re-emitted diffuse light from the tissue for boundary measurement of spatial diffuse speckle contrasts in the selected ROI. A zoom lens (Zoom 7000, Navitar) was connected to the camera for adjusting the camera’s focus on the ROI. A high-performance long-pass filter ($>750\mathrm{nm}$; #84-761, EdmundOptics) was used to minimize the influence of ambient light on DSCT measurements.²⁶ No obvious influence of ambient light was observed on the obtained raw intensity images or reconstructed BFI maps from the experiments. A pair of polarizers (LPNIRE050-B and LPNIRE200-B, Thorlabs) crossing the source and detection paths were added to reduce the specular reflection directly from the scanning light sources on the tissue surface.

Fig. 1

Hybrid DSCT and LSCI system for 2D mapping of CBF in rats. (a) The rat was anesthetized with 1% to 2% isoflurane and placed on a heating blanket with its head fixed on a stereotaxic frame. The hybrid instrument was set up overhead for 2D mapping of CBF distributions. (b) Optical design of DSCT. (c) NIR point source was focused on the exposed skull for DSCT measurements. (d) Optical design of LSCI. Widefield illumination was generated by manually placing a diffuser in the source path for LSCI measurements.

For comparisons with LSCI, a hybrid DSCT and LSCI system was assembled [Fig. 1(a)]. LSCI and DSCT measurements were performed sequentially in one rat at the baseline. For LSCI measurement, an engineering diffuser (ED1-S20-MD, 20° Square Engineered Diffuser, Thorlabs) was manually placed in the source path of DSCT to generate a widefield illumination without the need for point source scanning [Fig. 1(d)]. After LSCI measurement, the diffuser was manually removed to facilitate DSCT measurement. Both DSCT and LSCI measurements utilized the same sCMOS camera with an exposure time of 2 ms. The $F$-number of the detection zoom lens was set as 8 to meet the Nyquist sampling rule.¹⁴

2.2.

LSCI Method for 2D Mapping of Superficial CBF

The movement of light scattering particles (i.e., red blood cells) in the measured tissue volume results in spatial fluctuations of laser speckles on the tissue surface, which appear as blurred speckle patterns when recorded by the camera.³⁴ Using the LSCI technique with widefield illumination [Fig. 1(d)], the raw intensity image taken from the selected ROI is converted to a spatial speckle contrast image. The spatial speckle contrast is quantified by calculating the ratio of the standard deviation ($\sigma$) of the intensity ($I$) to its mean intensity ($⟨I⟩$) in a pixel window of $N_{\text{pixels}}\times N_{\text{pixels}}$:^14,15,35 $K_s=\sigma /⟨I⟩=\sqrt{(⟨I^2⟩-⟨I{⟩}^2)}/⟨I⟩$. The pixel window size of $3\times 3$, $5\times 5$, or $7\times 7$ is usually selected to balance the detection sensitivity and spatial resolution.³⁵ Larger values of $N_{\text{pixels}}$ assess spatial speckle contrasts better in terms of sensitivity at the expense of lower spatial resolutions.¹⁴

Although the exact relationship between $K_s$ and BFI is nonlinear, BFI at each pixel widow can be approximated as the inverse square of the speckle contrast: $\mathrm{BFI}\sim 1/K_s^2$.^34,36 Due to its numerical simplicity, this relationship enables fast data processing of $K_s$ values over all pixel windows. Traditionally, two nested for-loops are utilized to calculate $K_s$ on each pixel window of $N_{\text{pixels}}\times N_{\text{pixels}}$. Then, a third for-loop iterates the pixel window through all pixels on the ROI to generate a 2D map of BFI.

We implemented a fast computation method for analyzing the LSCI and DSCT data, which takes advantage of functions from Image Processing Toolbox™ and Parallel Computing Toolbox™ in MATLAB to dynamically calculate the 2D map of BFI. Instead of using the inefficient for-loops to calculate $K_s$ values within the raw image, we applied MATLAB’s built-in 2D convolution function (conv2) from the Image Processing Toolbox™ to implement an implicit multi-threading approach for the fast computation.³⁷ Figure 2(a) shows the steps of the new fast parallel computation method to obtain a speckle contrast image from the raw intensity image, in which the kernel matrix was defined as a square matrix of $N_{\text{pixels}}\times N_{\text{pixels}}$ with all elements equal to $1/N_{\text{pixels}}^2$.

Fig. 2

Flowchart for data processing in LSCI and DSCT methods. (a) New fast parallel computation method for LSCI data analysis. The procedures within the green dashed box were executed on each raw intensity image to calculate $K_s$ values. Parallel Computing Toolbox™ in MATLAB was then used to process multiple raw intensity images simultaneously and generate multiple flow maps. (b) Determination of the S-D mask in DSCT for extracting a $K_s$ map at the depth of $d$. The detectors within the specified S-D distances of [$d1$: $d2$] are selected for each source to calculate $K_s$ values at the depth of $d$. (c) New fast parallel computation method for DSCT data analysis. The procedures within the blue dashed box were executed on raw intensity images of each source ($S_i$) to calculate corresponding $K_s(i)$ values. Parallel Computing Toolbox™ in MATLAB was then used to process raw intensity images of all sources simultaneously. The steps depicted in purple boxes were implemented to reduce the artifacts arising from the summation of $K_s(i)$ in the overlap regions of the masks and ultimately generate one $K_s$ map.

Because calculations of multiple speckle contrast images are independent, MATLAB’s parfor control flow statement in the Parallel Computing Toolbox™ was used to execute loop iterations in parallel to compute speckle contrast images simultaneously. Specifically, in this study, MATLAB’s codes were executed on a desktop with Intel(R) Core (TM) i9-10980XE CPU @ 3.00 GHz (18-Core), NVIDIA Quadro RTX 5000 GPU. MATLAB (version: R2021a) automatically utilized all available CPU cores by default, which in this case was 18 workers, when executing parallel commands. As a result, total computation time with the new parallel computation method was reduced dramatically compared with the use of inefficient for-loops to calculate $K_s$ values from the raw intensity images sequentially, e.g., reducing from 6.5 h to 30 s to generate 400 flow maps with the total pixels of $2048\times 2048$ in each map.

2.3.

DSCT Method for 2D Mapping of CBF at Different Depths

Following our established method,^20,26 raw intensity images were prescreened to ensure that the detected light intensities (counts) were within the linear range of our camera and above the minimal SNR requirement (i.e., $\mathrm{SNR}>3$). The DSCT method first converted each of the raw intensity images taken at different source locations within the ROI into a speckle contrast image based on the fast LSCI analysis utilizing the Image Processing Toolbox™ and Parallel Computing Toolbox™ (Sec. 2.2). Because the penetration depth ($d$) of DSCT with point source illumination is approximately one half of the S-D distance, the Euclidean distances between each source and detectors (i.e., pixel windows of $N_{\text{pixels}}\times N_{\text{pixels}}$) were calculated. To extract a $K_s$ map at the depth of $d$, the pixels/detectors within the S-D distances of [d1:d2] from the source were averaged to improve the SNR [Fig. 2(b)]. Here, $d1=2d-\mathrm{\Delta }d$, $d2=2d+\mathrm{\Delta }d$, and $\mathrm{\Delta }d$ was selected empirically as 1 mm to balance the SNR and spatial resolution. More generally for each source $(S_i)$, a binary belt shape mask denoted as $M_{S-D}(i)$ was defined such that the pixels corresponding to the detectors located within S-D distances of [$d1$: $d2$] were assigned a value of “1,” and all other pixels were assigned a value of zero. This mask identified the paired detectors for each source to extract corresponding speckle contrast values of $K_s(i)$. These calculations were implemented through all sources and the resulting speckle contrast images were summed to generate one $K_s$ map in the selected ROI [Fig. 2(c)]. To mitigate the artifacts from the summation of $K_s(i)$ in the overlap regions of masks, the inverse of mask summation was multiplied by $K_s(i)$ summation. Note that DSCT data analyses are based on our fast parallel LSCI analysis (Sec. 2.2), thus also decreasing the total computation time dramatically compared with using inefficient for-loops methods.

The DSCT data analysis requires the input of source locations. Following our established method,²⁶ the intensity image captured at each source location was analyzed to precisely determine its location. This method enabled accounting for the potential shift of each source location caused by tissue geometric factor. Briefly, the intensity image taken at each source location was first converted to a binary image. This was done by replacing all intensity values above a globally determined threshold with “1” and setting all other values to “0.” The threshold value was determined by Otsu’s method to minimize the intraclass variance of the thresholded black and white pixels.³⁸ After thresholding, those pixels with a value of “1” that connected to 8 pixels in the local neighborhoods (so called eight-connected) were found and grouped together. The group with the maximum number of pixels was identified, and the rest of the groups were removed from the 2D binary image. Morphological operators, such as opening, closing, and filling, were then applied to the binary image for removing thin protrusions and isolated areas to reduce noise. Finally, the centroid and radius of the point source were determined by estimating the center of the remaining group. This automatic process was repeated for determining all source locations.²⁶

2.4.

Depth Sensitivity Tests Using New Head-Simulating Layered Phantoms

The use of tissue-simulating phantoms with known optical properties is a commonly accepted approach for the validation of new methods and technologies.^31,39,40 In this study, head-simulating layered phantoms with known optical/geometry properties were fabricated using 3D printing techniques for testing the sensitivity of the DSCT method in mapping flow distributions at different depths. Two layers of head tissues (skull and brain) were simulated using the solid phantom with zero flow and the Intralipid liquid phantom with particle flow, respectively. Figures 3(a) and 3(b) show the 3D solid phantom (no flow) with empty channels bearing the University of Kentucky (UK) logo, fabricated by a 3D printer (SL1, Prusa). Three solid phantoms with UK logo channels were printed with the top layer thicknesses of 1, 2, and 3 mm, respectively, to mimic skulls with varied thicknesses. The solid phantoms were made of titanium dioxide (${\mathrm{TiO}}_2$), India ink (Black India, Massachusetts), and clear resin (eSUN Hard-Tough).⁴¹ The empty UK logo channels were then filled with liquid phantom solutions composed of Intralipid solution (Fresenius Kabi, Sweden), India ink, and water [Figs. 3(c) and 3(d)].²⁰ India ink concentration regulates the tissue absorption coefficient (${\mu }_a$), whereas ${\mathrm{TiO}}_2$ and Intralipid concentrations regulate the reduced scattering coefficient (${\mu }_s^{\prime }$). The optical properties of both the solid and liquid phantoms were set as ${\mu }_a=0.03{\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }=9{\mathrm{cm}}^{-1}$.^19,40 Brownian motion of Intralipid particles inside the UK logo channels provides the particle flow to mimic the motion of red blood cells inside vessels (i.e., blood flow).^20,42

Fig. 3

New head-simulating layered phantoms with known optical/geometry properties. (a) 3D design of solid phantom (no flow) with empty channels bearing the UK logo. The holes connected to the UK logo channels were used to fill in liquid phantoms. (b) 3D printed solid phantom with empty UK logo channels, fabricated by a 3D printer with the top layer thicknesses of 1 mm. The dashed square indicates the scanning ROI of $40\times 40{\mathrm{mm}}^2$ on the phantom surface. (c) Bottom view of solid phantom with the empty UK logo channels. (d) The empty UK logo channels were filled with liquid phantom solutions and then sealed with plastic film and a hot glue gun.

To evaluate the impact of sampling density on spatial resolution when using DSCT, $10\times 10$, $15\times 15$, $20\times 20$, and $30\times 30$ source locations were scanned sequentially inside the ROI of $40\times 40{\mathrm{mm}}^2$ on the UK logo phantom with the top layer thickness of 1 mm. Then, the UK logo phantoms with the top layer thicknesses of 1, 2, and 3 mm were imaged by the DSCT (with $30\times 30$ source locations) and LSCI sequentially [Fig. 3(b)]. For both DSCT and LSCI, $K_s$ was computed over a pixel window of 49 pixels (i.e., $N_{\text{pixels}}=7$) to balance the imaging sensitivity and spatial resolution. Five sequential flow maps over time were averaged to improve SNRs. Relative flow values were calculated by normalizing flow indices to their mean within the ROI. Finally, the results from the DSCT and LSCI were compared to demonstrate the depth sensitivity of DSCT in contrast to LSCI.

2.5.

In-Vivo Experiments in Rats

To evaluate the ability of DSCT for continuous monitoring of CBF variations, adult rats with intact skulls (scalps were retracted) underwent ${\mathrm{CO}}_2$ inhalation, and transient carotid arterial ligations were imaged. It is well known that ${\mathrm{CO}}_2$ inhalation leads to global CBF increases, whereas transient arterial ligations lead to regional CBF decreases. All experimental procedures involving animals were approved by the University of Kentucky Institutional Animal Care and Use Committee (IACUC).

Six male adult Sprague–Dawley rats (age: 10 to 11 months) were imaged by the DSCT. In addition, one animal was imaged by the LSCI for comparisons. For both DSCT and LSCI, $K_s$ was computed over a pixel window of nine pixels (i.e., $N_{\text{pixels}}=3$) to balance the imaging sensitivity and spatial resolution. The DSCT scanned over 441 source positions ($21\times 21$) in the ROI of $20\times 20{\mathrm{mm}}^2$ on the rat head to balance temporal and spatial resolutions. 2D maps of CBF at the depths of 1 to 3 mm were generated based on the 1 to 2 mm skull thickness of adult rats.²⁰ The relative time-course changes in CBF (rCBF) were calculated by normalizing BFI data to the baseline value before physiological changes.

Unlike the phantom studies in which temporal resolution is not critical for imaging constant particle flow contrasts (Sec. 2.4), a fast sampling rate is essential to capturing dynamic CBF variations caused by pathophysiological changes. Using a conventional 2D mapping method with 441 sources ($21\times 21$), the total sampling time to generate one 2D map of CBF was 60 s. In this study, a new method based on a moving window was adapted to improve the sampling rate, in which the topographic reconstruction of each CBF map was obtained when one new source position data point became available.^4,26 As a result, the total sampling time was reduced $\sim 440$ times (i.e., from 60 to 0.136 s).

2.6.

Statistical Analysis

SPSS software (version 29) was used for statistical analysis in animal studies. The Kolmogorov–Smirnov test was performed to assess the normality of rCBF distributions. Also, repeated measures analysis of variance (ANOVA) was used to evaluate differences in rCBF variations between the different phases of stimuli (${\mathrm{CO}}_2$ inhalations and transient arterial occlusions) and the different hemispheres. A $p\text{-}\text{value}<0.05$ is considered significant for all statistical analysis.

3.1.

DSCT Performance Depends on Sampling Density

Figure 4 shows the impact of sampling density on spatial resolution when using DSCT to image the UK logo phantom with the top layer thickness of 1 mm. Increasing scanning source numbers improved spatial resolution at the costs of a prolonged data acquisition time, an increased number of images for processing, and increased computation and storage demands (Table 1). More specifically, the maximum data acquisition time for obtaining one raw intensity image at one source position was $\sim 150\mathrm{ms}$, including 2 ms for exposure, 33 ms for capturing one image by the camera, and 115 ms for saving data to the hard drive.

Fig. 4

Effect of sampling density on spatial resolution of DSCT. (a)–(d) Numbers of source locations on the UK logo phantom with a top layer thickness of 1 mm: 100 ($10\times 10$), 225 ($15\times 15$), 400 ($20\times 20$), and 900 ($30\times 30$). (e)–(h) Resulting 2D flow maps at the depth of 2 mm with the source numbers of 100, 225, 400, and 900, respectively. Five sequential flow maps were averaged to improve the SNR.

Table 1

Impact of sampling density on DSCT performance.

3.2.

DSCT Enables 2D Mapping of Phantom Flow Contrasts with Depth Sensitivity

Figure 5 shows the results using the DSCT and LSCI to image UK logo phantoms with top layer thicknesses of 1, 2, and 3 mm, respectively. For DSCT measurements, $30\times 30$ source locations were scanned over the selected ROI of $40\times 40{\mathrm{mm}}^2$ on the phantoms. Higher flow contrasts were observed at deeper depths with larger S-D separations, indicating the depth sensitivity of the DSCT method. SNRs decreased with the increases of imaging depth and top layer thickness (i.e., from 1 to 3 mm). These results are expected as the deeper penetration and thicker top layer resulted in fewer diffused photons being detected, thus leading to lower SNRs. By contrast, SNRs and spatial resolutions of superficial LSCI measurements were much worse than those of DSCT measurements, especially on the phantoms with thicker top layers.

Fig. 5

DSCT and LSCI measurement results for UK logo phantoms. (a), (b) Resulting 2D maps of Intralipid particle flow contrasts on three phantoms with the top layer thicknesses of 1, 2, and 3 mm, respectively. Flow indices were normalized to their mean values to generate relative flow maps for comparisons.

3.3.

DSCT Enables 2D Mapping of rCBF at Different Depths in Rats

Figures 6(a) and 6(b) show white light images of a rat head (scalp was removed) and 441 source positions ($21\times 21$) on the rat skull for DSCT measurements. Figure 6(c) shows rCBF map of a representative rat (rat #4), imaged by the LSCI. Figures 6(d)–6(h) show rCBF maps of rat #4 at depths varying from 1 to 3 mm, imaged by the DSCT. As observed in phantom measurements (Sec. 3.2), DSCT enables 2D mapping of rCBF distributions at various depths, with higher spatial resolutions and SNRs than LSCI.

Fig. 6

2D mapping of rCBF distributions with LSCI and DSCT. (a) A white light image of a rat head with nose up (rat #4, scalp was retracted). The yellow dotted line indicates the central line that separates the left and right hemispheres (RH). (b) Sources distribution (red dots) on the ROI of $20\times 20{\mathrm{mm}}^2$. (c) 2D map of rCBF, imaged by LSCI. (d)–(h) 2D maps of rCBF at different depths from 1.5 to 3.0 mm, imaged by DSCT.

3.4.

DSCT Captures Dynamic rCBF Responses During ${\mathrm{CO}}_2$ Inhalations in Rats

Figure 7(a) shows typical blood flow maps in a representative rat (rat #2) at the depth of 3.0 mm before, during, and after $8\%{\mathrm{CO}}_2$ inhalation, measured by the DSCT with 441 source positions ($21\times 21$). Figure 7(b) shows time-course changes in rCBF in one rat (rat #2), calculated by the moving window method and normalizing BFI data to their baseline values (assigned as “100%”) before ${\mathrm{CO}}_2$ inhalation. Figure 7(c) compares the results of time-course rCBF changes between the conventional reconstruction (red curve) and moving window (blue curve) methods. As expected, both methods generated similar trends in rCBF, although the sampling rate of the moving window method was much higher than the standard reconstruction method. Figure 7(d) compares time-course changes in rCBF at different depths ranging from 1.0 to 3.0 mm. As expected, rCBF changes in deeper brain tissues ($\text{depth}\ge 2.0\mathrm{mm}$) were higher than those in the skull ($\text{depth}\le 1.5\mathrm{mm}$). Figure 7(e) shows average time-course changes in rCBF (moving window method) over six rats. Table 2 and Table S1 in the Supplemental Material summarize the average changes in rCBF and corresponding significance levels ($p$-values), respectively, during different phases of ${\mathrm{CO}}_2$ inhalation. The $8\%{\mathrm{CO}}_2$ led to a significant increase in rCBF at the endpoint of inhalation compared with the baseline (mean ± standard error: $117\%\pm 5\%$; repeated measures ANOVA: $p\text{-}\text{value}=0.024$), which agrees with previously reported changes in rCBF.²⁰ There was also a significant difference (repeated measures ANOVA: $p\text{-}\text{value}=0.046$) between the $8\%{\mathrm{CO}}_2$ inhalation and the endpoint of the recovery phase.

Fig. 7

DSCT measurements of rCBF responses to ${\mathrm{CO}}_2$ inhalation in rats. (a) 2D maps of rCBF at the depth of 3 mm before, during, and after $8\%{\mathrm{CO}}_2$ inhalation in an illustrative rat (rat #2). The dashed green box indicates the selected brain region for data analysis. (b) Time-course changes in rCBF at a depth of 3 mm, averaged over the brain region (dashed green box) and reported as mean value ± standard deviation in one rat ($n=1$, rat #2). The shaded area represents standard deviations of spatial changes. (c) Comparison of time-course changes in rCBF (mean values) between the conventional reconstruction (red curve) and moving window (blue curve) methods (rat #2). (d) Comparison of time-course changes in rCBF (mean values) in one rat ($n=1$, rat #2) at different depths ranging from 1.0 to 3.0 mm (moving window method). (e) Average time-course rCBF changes (moving window method) at a depth of 3 mm over six rats (mean value ± standard error). The shaded area represents standard errors.

Table 2

Average rCBF variations (mean ± standard error) from their baselines during CO2 inhalation over six rats (moving window method).

3.5.

DSCT Captures Dynamic rCBF Responses to Transient Arterial Ischemia and Hypoxia in Rats

Figure 8(a) shows typical blood flow maps in a representative rat (rat #3) at the depth of 3.0 mm before, during, and after sequential CCA ligations, measured by the DSCT with 441 source positions ($21\times 21$). Figure 8(b) shows time-course changes in rCBF at the right hemisphere (RH) and left hemisphere (LH) of one rat (rat #3), which were calculated by the moving window method and normalizing BFI data to the baseline values prior to CCA ligations. Figure 8(c) shows average time-course changes in rCBF over five rats. Table 3 and Table S2 in the Supplemental Material summarize the average changes in rCBF and corresponding significance levels ($p$-values), respectively, during different phases of CCA ligations. As expected, the sequential right and left CCA ligations as well as subsequent releasing resulted in significant alterations in rCBF in both hemispheres compared with their baselines (repeated measures ANOVA: $p\text{-}\text{value}<0.05$). In particular, the bilateral ligation led to significant reductions in rCBF in both hemispheres, compared with all other phases of transient CCA ligations (repeated measures ANOVA: $p\text{-}\text{value}<0.05$). These results agree with our previous study in adult rats utilizing scDCT with 3D reconstruction.²⁰

Fig. 8

DSCT measurements of rCBF responses to sequential unilateral and bilateral CCA ligations in rats. (a) 2D maps of rCBF at the depth of 3 mm before, during, and after unilateral and bilateral CCA ligations in an illustrative rat (rat #3). The dashed green boxes indicate selected brain regions at RH and LH for data analyses. (b) Time-course changes in rCBF (moving window method), averaged over the brain regions (dashed green boxes) and reported as mean value ± standard deviation (rat #3). The shaded area represents standard deviations. (c) Average time-course rCBF changes over five rats (mean value ± standard error). The shaded area represents standard errors.

Table 3

Average rCBF variations (mean ± standard error) from their baselines during unilateral and bilateral CCA ligations over five rats (moving window method).

Figure 9(a) shows typical blood flow maps in a representative rat (rat #4) at the depth of 3.0 mm before and during $100\%{\mathrm{CO}}_2$ euthanasia, measured by the DSCT with 441 source positions ($21\times 21$). Figure 9(b) shows time-course changes in rCBF, which were calculated by the moving window method and normalizing BFI data to the baseline values prior to asphyxia. Figure 9(c) shows average time-course changes in rCBF over five rats. At the end, rCBF values were close to “zero,” meeting the physiological expectation.

Fig. 9

DSCT measurements of rCBF variations during asphyxia ($100\%{\mathrm{CO}}_2$) in rats. (a) 2D maps of rCBF at the depth of 3 mm before and at the end of asphyxia challenge in an illustrative rat (rat #4). The dashed green box indicates the selected brain region for data analysis. (b) Time-course changes in rCBF, averaged over the brain region (dashed green boxes) and reported as mean value ± standard deviation. The shaded area represents standard deviations. (c) Average time-course rCBF changes over five rats (mean value ± standard error). The shaded area represents standard errors.