Multimodal Noninvasive Functional Neurophotonic Imaging of Murine Brain‐Wide Sensory Responses

Abstract

Modern optical neuroimaging approaches are expanding the ability to elucidate complex brain function. Diverse imaging contrasts enable direct observation of neural activity with functional sensors along with the induced hemodynamic responses. To date, decoupling the complex interplay of neurovascular coupling and dynamical physiological states has remained challenging when employing single‐modality functional neuroimaging readings. A hybrid fluorescence optoacoustic tomography platform combined with a custom data processing pipeline based on statistical parametric mapping is devised, attaining the first noninvasive observation of simultaneous calcium and hemodynamic activation patterns using optical contrasts. Correlated changes in the oxy‐ and deoxygenated hemoglobin, total hemoglobin, oxygen saturation, and rapid GCaMP6f fluorescence signals are observed in response to peripheral sensory stimulation. While the concurrent epifluorescence serves to corroborate and complement the functional optoacoustic observations, the latter further aids in decoupling the rapid calcium responses from the slowly varying background in the fluorescence recordings mediated by hemodynamic changes. The hybrid imaging platform expands the capabilities of conventional neuroimaging methods to provide more comprehensive functional readings for studying neurovascular and neurometabolic coupling mechanisms and related diseases.

Introduction

Functional neuroimaging has become a primary tool in neuroscience and research into neurodegenerative diseases. Blood-oxygen-level-dependent (BOLD) functional magnetic resonance imaging (fMRI) is the mainstay for noninvasive mapping of brain function.^[1] However, interpretation of activity-related BOLD signals is challenged by their dependence on multiple factors such as blood flow, blood volume, metabolic rate of oxygen, and baseline physiological state.^[2] Functional ultrasound imaging can provide additional information on the blood flow.^[3] However, the blood flow changes similarly represent an indirect correlate of the underlying neural dynamics, which is insufficient for fully characterizing the complex interplay between neural activity and its accompanying hemodynamics.

Optical imaging approaches arguably provide the most powerful means for retrieving functional and molecular information underlying brain activity. The emergence of GCaMP-type and other genetically encoded calcium indicators has facilitated fluorescence-mediated sensing of neural activation,^[4] and has for example been widely employed to characterize resting state activity and stimulus-evoked responses in the mouse cortex.^[5] Brain hemodynamics can be further extracted by means of multispectral optical intrinsic signal imaging,^[5,6] which is however a purely planar (2D) imaging technique limited to superficial observations and commonly requires invasive procedures, such as scalp removal and skull thinning. The lack of depth resolution further compromises accuracy of the extracted individual hemodynamic components as the signals are basically unevenly integrated along the depth dimension. Since hemoglobin exhibits distinct absorption spectra in its oxygenated (HbO) and deoxygenated (HbR) states,^[7] brain hemodynamics can alternatively be visualized with functional optoacoustic (fOA) tomography across 3D space, time, and optical spectrum, thus enabling the mapping of multiple hemodynamic parameters across the entire murine brain, not accessible with other modalities.^[8,9] fOA can thus potentially massively advance our knowledge on the mechanisms of neurovascular coupling, that is, if the mapping of both the neuronal activity directly as well as its associated hemodynamic changes could be realized simultaneously. However, a detailed characterization and rigorous side-by-side validation of hemodynamic fOA readings against well-established, direct optical neuroimaging methods is currently lacking.

Here, we introduce a hybrid fluorescence and optoacoustic tomography (FLOT) platform for concurrent multiparametric noninvasive characterization of brain-wide hemodynamic and calcium responses to sensory stimulation in mice. A custom data processing pipeline inspired by statistical parametric mapping (SPM) was further developed to facilitate the functional data analysis. Simultaneous FLOT-based imaging of calcium and hemodynamic responses to electrical paw stimulation is demonstrated through the intact skull and scalp of GCaMP6f mice, providing noninvasive readouts of both direct^[5,10] and indirect correlates of neural activity. To the best of our knowledge, this represents the first successful observation of activation patterns (direct and indirect responses simultaneously) in the entire mouse cortex with optical contrast. The use of near-infrared (NIR) wavelengths in fOA imaging further enabled to reach deeper brain regions thus averting cross-talk between calcium and hemodynamic signals.

Results The Hybrid FLOT System

The dedicated hybrid FLOT imaging system employs a flexible fiberscope for capturing fluorescence images in epi-illumination mode with ≈44 µm pixel resolution. The distal end of the fiberscope was inserted into the central cavity of a spherical matrix array transducer used for acquiring 3D fOA data with nearly isotropic spatial resolution which was estimated to be 113 µm using microsphere measurements, slightly degrading at the periphery FOV (Figure 1a; and Figure S1, Supporting Information).^[11,12] Both modalities provided a 12 × 12 mm² field of view (FOV) along the lateral x-y plane covering the entire mouse cortex. Trains of unipolar rectangular electrical pulses with 5 ms duration and 1.0 mA intensity were applied to the left hindpaw of Thy1-GCaMP6f mice. The trains consisted of 32 pulses with stimulus frequency of 4 Hz (i.e., 8 s total duration) repeated every 82 s (Figure 1b). Five wavelengths (700, 730, 755, 800, and 850 nm) were optimally selected for fOA to avoid significant spectral coloring while ensuring a small condition number of the spectral unmixing matrix (Figure 1c). The unmixed hemodynamic components of HbO and HbR are shown in Figure 1d. Subsequently, both fluorescence and fOA data were analyzed with a dedicated data processing pipeline (Figure 1e). We employed regressors for analyzing the hemodynamic and calcium responses by convolving the stimulation pulse train with a hemodynamic response function (HRF) and a GCaMP calcium response function (CRF), respectively (see the Experimental Section for a detailed description).

Figure 1. Layout of the hybrid FLOT platform for multiparametric imaging of murine brain activation. a) Schematics of the imaging system incorporating the simultaneous functional optoacoustic (fOA) and fluorescence readings. b) Burst electric current stimulation paradigm applied to the left hindpaw. c) Wavelength selection for the spectroscopic fOA imaging of blood oxygenation. The molar extinction spectra for HbO and HbR are plotted. d) Spatial distribution of the spectrally unmixed HbO and HbR hemodynamic components in the brain. e) The custom FLOT data processing pipeline includes fOA image reconstruction and spectral unmixing (left column), fluorescence data reshaping, preprocessing and the general linear model (middle column), and statistical parametric mapping (right column).

GCaMP6f Calcium Response Function (CRF)

We first estimated the CRF of the rapid calcium activity using epifluorescence recordings acquired at a high frame rate of 200 Hz in Thy1-GCaMP6f mouse brains in vivo. The high frame-rate acquisitions were employed to extract well-defined response function profiles subsequently used in the fluorescence data processing pipeline. The scalp of Thy1-GCaMP6f mice was removed prior to imaging while keeping the skull intact (Figure 2a), which resulted in clearly resolvable cortical vasculature. Temporal fluorescence signal profiles from two representative 0.3 × 0.3 mm² regions of interest (ROIs) were tracked in both the contralateral (opposite side from the stimulated paw) and ipsilateral hemispheres, as indicated in Figure 2a. The recorded profiles were averaged across different stimulation bursts to reveal the activation time course (Figure 2b). The averaged contralateral and ipsilateral signal profiles manifested a high degree of correlation, consisting of both spontaneous wave activity as well as sharp contralateral activity spikes evoked by the electric pulses. The respective averaged time-lapse images further reveal the spatiotemporal patterns of stimulus-evoked neuronal activation in the somatosensory cortex (Figure 2c). In order to estimate the average stimulus-evoked fluorescence response (Figure 2d; and Figure S2, Supporting Information), we first applied a temporal notch filter with 0.2–3 Hz stop-band to the unaveraged time profiles in the selected contralateral ROI. Next, the filtered profiles were averaged for all the stimulation electric pulses with the resulting curve fitted to a gamma function y = ct^{a − 1}e^−t/b/[b^a*Γ(a)], a = 4.775, b = 0.016, c = 0.079. In the subsequent analysis, the gamma fitted curve was adopted as the characteristic GCaMP6f CRF to external electrical paw stimulation. The calculated correlation coefficient between the measured curve and the fitted gamma function is 0.991 with p-value < 10⁻⁵ for t test. Other functional parameters were further calculated across all the stimulation trials, including the activation intensity 1.3 ± 0.4%, time-to-peak (TTP) 0.063 ± 0.013 s, full-width-at-half-maximum (FWHM) of the calcium response 0.085 ± 0.024 s and decay time (T_1/2) 0.106 ± 0.025 s (Figure 2e). A linear regressor was then obtained by convolving the stimulation paradigm with the GCaMP6f CRF. Fluorescence data processing with the proposed pipeline rendered an activation map revealing the expected location of the activated region with high specificity (Figure 2f). Compared to an activation map based on the difference image between activated frames and baseline recording, the proposed method exhibits superior robustness and has better localization specificity as it is less affected by the slow changes in light intensity affected by hemoglobin absorption variations (Figure S3, Supporting Information).

Figure 2. Assessment of the GCaMP6f impulse response function to electrical paw stimulation in mice. a) Representative time-averaged epi-fluorescence image of the mouse brain. b) Temporal fluorescence signal profiles spatially averaged over two 0.3 × 0.3 mm2 sized ROIs in the contralateral (CL) and ipsilateral (IL) somatosensory cortex, as indicated in panel a). Each vertical line indicates single electric pulse. c) Time-lapse images post the pulsed electric current stimulation obtained by averaging the responses from all the stimulation trials (pulses). d) Averaged stimulus-evoked response function estimated by averaging activation curves from 160 consecutive stimulation pulses fitted to a Gamma function. e) Functional parameters extracted from all the stimulation trials. f) Activation map obtained from the proposed data processing pipeline. All scale bars: 1 mm.

FLOT Reveals Concurrent Calcium and Hemodynamic Response Maps Noninvasively

Next, we simultaneously imaged calcium and hemodynamic responses in Thy1-GCaMP6f mice (n = 6). Since the fOA imaging was performed in the near-infrared spectral window, the scalp was kept intact for 4 mice, while it was removed for another 2 mice to mitigate image artifacts due to skin pigmentation (see the Experimental Section). Activation maps corresponding to different hemodynamic components, namely HbO, HbR, total hemoglobin (HbT) and oxygen saturation (sO₂), were further rendered from the fOA image sequences. Localized responses were clearly observed in the transverse, sagittal, and coronal cross sections overlaid to the Allen P56 mouse brain atlas^[13] (Figure 3a; voxels in the activation map are statistically significant with respect to a null hypothesis of no activation: p < 0.05 one-sample t-test after false discovery rate (FDR) correction). Localized activity was detected within the primary somatosensory cortex as well as within the primary motor area on the contralateral side. No obvious activation was observed in the corresponding regions on the ipsilateral side. Since the HbR response was negatively correlated to the hemodynamic response regressor while other components were positively correlated, both positive and negative thresholds were applied to the calculated activation map (i.e., t-map) of each component. Multiple functional parameters were estimated from these fOA maps, i.e., the activation intensity, TTP, FWHM, decay time, and activation volume size (Figure 3b).

Figure 3. Concurrent measurement of calcium and hemodynamic responses in the mouse brain. a) Transverse, sagittal, and coronal views of the activation maps corresponding to HbO, HbR, HbT, and sO2 components extracted from fOA images. b) Superimposed fOA activation maps and statistics of the functional parameters from the ROI indicated by the black square situated at a depth of ≈300 µm from the surface of the cortex. c) Corresponding GCaMP and hemodynamic activation maps calculated from the fluorescence measurements overlaid with the fluorescence image of a GCaMP mouse brain with intact skull and scalp. d) Unaveraged time courses of the fluorescence and fOA signals from the ROIs indicated in b) and c). The fluorescence signal is lowpass filtered to emphasize slow trends (top row; Butterworth, cutoff frequency 0.5 Hz). e) Averaged activation time courses of fluorescence, HbO, HbR, HbT, and sO2 from the ROI indicated in b) and c). CL: contralateral, IL: ipsilateral. All scale bars: 1 mm.

We next analyzed the concurrently-acquired fluorescence data. Due to the strong light scattering from the scalp, the fluorescence images appeared blurry compared to those acquired after scalp removal. Yet, localized brain activation upon electrical paw stimulation was clearly observed through the intact scalp, corroborating the proposed technique's noninvasive functional imaging capacity. Moreover, increased light absorption by hemoglobin results in modulation of fluorescence signals by the hemodynamic responses in the activated region. A secondary hemodynamic activation map was thus obtained simultaneously from the fluorescence measurement by simply replacing the GCaMP CRF with the HRF in the construction of the regressor and overlaid onto the GCaMP activation map (Figure 3c; both activation maps are thresholded at 70% of their corresponding maximum t-value).

Both unaveraged fluorescence and fOA signal time courses are shown in Figure 3d, depicting the highly dynamic spontaneous fluctuations and stimulus-evoked responses. A significant correlation exists between the simultaneously recorded fluorescence and HbO/HbR responses, despite the fact that the signals were measured with two independent modalities, which may partially explain the slow background fluorescence signal fluctuations. These dynamics may further depend on the physiological status of the mouse, which is profoundly affected by many factors such as anesthesia and individual differences. Fractional averaged signal changes from representative ROIs in contralateral hemisphere are shown in Figure 3e. Interestingly, the lowpass filtered fluorescence signal peaked right after the paw stimulation due to the fast calcium-mediated response but then slowly declined below its baseline value gradually going back to the baseline, reflecting changes in light attenuation associated to increased local hemoglobin absorption in the activated regions.

Our results show that fOA provides diverse information on brain activation responses complementing the fluorescence readings. Further, it was observed that fractional hemodynamic signal changes significantly differed across various brain regions whereas fluorescence responses did not exhibit such spatial variability (Figure S4, Supporting Information).

Multiparametric Analysis of Stimulus-Evoked Brain Activation

Multiparametric analysis was performed to characterize the coupling between simultaneous fluorescence and fOA readings. The same data processing pipeline was repeated for all mice (n = 6) across a total of 80 stimulation trials. Activation maps were first calculated using both the fluorescence and fOA images to infer activated brain regions and thereby to select ROIs for calculating the time courses. Spontaneous neural activity largely averages out pre- and poststimulation, while the stimulation-evoked responses are preserved (Figure 4a). The fractional signal changes of the different hemodynamic components were then averaged across animals (Figure 4b–e, n = 6). Strongly correlated signal increase in HbO (Figure 4b), HbT (Figure 4d) and sO₂ (Figure 4e), and HbR decline (Figure 4c) were observed from the fOA measurements. As expected, these signal changes are closely correlated to the predicted hemodynamic response by convolving the HRF with the mean neural activity trace (green curve in Figure 4a). Based on the fractional change in HbO and HbR, the fluorescence signal contaminated by slowly varying changes in hemoglobin absorption was successfully corrected^[5] using estimates of the pathlength factors xEx = 0.56 mm⁻¹ and xEm = 0.57 mm⁻¹ (Figure 4f).

Figure 4. Multiparametric analysis of the coupling between simultaneous fluorescence and fOA readings. a) Averaged fluorescence activation curve after applying a bandpass filter (3–8 Hz) and the predicted hemodynamic response (HR) calculated via convolution with the HRF. b–e) Averaged activation curves of HbO, HbR, HbT, and sO2 across all the stimulation trials. Shaded regions show standard error of mean (SEM) across n = 6 mice in all the trials. f) Fluorescence signal compensated for hemoglobin absorption variations by considering the fractional changes of HbO and HbR from the contralateral side. g–j) Coupling between HbO, HbR, HbT, sO2, and lowpass filtered fluorescence signal from simultaneously measured single activation traces. Single trial responses were assembled into 10 equally distributed amplitude bins where mean and standard deviation values were respectively estimated. k) Statistics of the activation intensity, TTP, FWHM, and decay time (T1/2) across all the trials from 6 mice. Note that all values are statistically significant compared to their ipsilateral counterparts (p-value [less than] 0.05, two-tailed paired t-test). CL: contralateral, IL: ipsilateral.

The coupling between neural activity and the various hemodynamic response components was further studied by calculating the relationship between the peak values of the corresponding (simultaneously measured) single activation traces of different normalized response components (Figure 4g–j). To calculate these plots, the range of observed single-trial values of each respective quantity was binned into 10 equal bins, and the mean and range of the respective coupled quantity was calculated. In Figure 4g,h, the results show a strong coupling relationship between HbT and sO₂ responses versus both HbR and HbO responses, respectively. These results show that both HbT and sO₂ responses saturate for high HbO responses, while remaining linear across the full HbR response range. This is partly explained by the related observation that higher HbO responses were tied to higher (and weakly saturating) fluorescence responses than those related to respective HbR ranges (Figure 4i). An overall stronger coupling of sO₂ versus HbT responses to both HbO and HbR was observed (respective slopes are higher in Figure 4g,h). The sO₂ responses were also more tightly coupled than HbT responses to the neuronal fluorescence responses, both in terms of higher values and narrower deviation ranges (Figure 4j). The strong relationship between the different hemodynamic components extracted from the fOA measurement and the fluorescence signal corroborates the synergistic detection of neurovascular coupling phenomena in stimulus-evoked brain activation by the hybrid FLOT method.

Despite the strong correlation among fOA components, noticeable differences were observed in the calculated functional parameters (Figure 4k). sO₂ (1.7 ± 0.9%) presented stronger activation intensity as compared to HbO (1.6 ± 0.9%), HbR (−1.2 ± 0.6%), and HbT (0.7 ± 0.2%). In terms of TTP, the HbR readings (15.7 ± 5.9 s) exhibited the most delayed response followed by sO₂ (10.3 ± 1.3 s), HbO (10.3 ± 1.4 s), and HbT (9.9 ± 2.9 s). Similar pattern was observed in the FWHM values with HbT (11.5 ± 4.9 s) holding the longest response duration followed by HbR (11.3 ± 6.5 s), HbO (10.5 ± 4.0 s), and sO₂ (10.3 ± 3.8 s). In terms of decay time (T_1/2), the HbR (28.2 ± 5.9 s) had the longest value followed by HbT (27.0 ± 8.6 s), sO₂ (25.5 ± 6.7 s), and HbO (24.4 ± 7.3 s).

Discussion

Major efforts are underway to devise new neuroimaging methods enabling scalable recording of neural activity and the associated hemodynamic responses. Both calcium and hemodynamic changes can be captured with sequential (nonconcurrent) multimodal strategies combining different contrast mechanisms. However, the mammalian brain's complexity and its highly dynamic physiological status challenge the establishment of interpretable links across separately acquired multimodal readings often covering different spatiotemporal scales. In this work, calcium signals and hemodynamic responses were imaged concurrently and a tailored functional data processing pipeline was developed to facilitate the comparison of calcium and hemodynamic responses. From a methodological perspective, we present the first hybridization of ultrafast 5D (i.e., spectroscopic time-resolved volume data) optoacoustic tomography and highly sensitive macroscopic fluorescence recording for concurrent multiparametric functional neuroimaging with optical contrast. The hybrid FLOT approach thus expands the capabilities of previously reported multimodal methods to provide more comprehensive functional information on stimulus-evoked neural activation in the murine brain in an entirely noninvasive fashion.

Our results exhibit a significant correlation between the GCaMP fluorescence signal and the associated increases in the HbO, HbT, sO₂, and decrease in the HbR signals (Figure 4a–e). In previous neuroimaging studies with fMRI, it was determined that HbR changes are the chief contributors to the BOLD signal.^[14] Interestingly, this signal is affected by changes in blood flow, blood volume, and metabolic rate of oxygen, which makes the BOLD signal relatively ambiguous to interpret. In contrast, we were able here to retrieve signal changes from multiple individual hemodynamic components as well as calcium responses simultaneously. The coupling analysis (Figure 4g,h) indicates that both ΔHbO and ΔHbR contribute to the local changes in HbT and sO₂, yet a tighter coupling is manifested with the ΔHbO component that is also closely related to the corresponding low-pass-filtered calcium signal ΔF measured from the same stimulation trials (Figure 4i). These physiological findings complement previous observations on event-evoked brain activation mapping and related hemodynamic readings based on wide-field optical mapping,^[5,15] fOA,^[8,16–18] and fMRI.^[2,19,20]

When comparing to wide-field optical mapping and fMRI, the fOA technique is particularly advantageous for characterizing hemodynamic responses in the brain as it can map multiple hemodynamic parameters in 3D with high spatial resolution and ultrafast (real-time) volume rate. Recently it has also been shown to be sensitive to fast stimulus-evoked calcium responses mediated by spectral absorption changes in GCaMP proteins.^[8,21] However, accurate spectral unmixing of GCaMP-associated changes from the strong background absorption by hemoglobin remains challenging with fOA, especially when employing calcium-sensitive molecules with peak absorption in the visible spectral range. We have considered the possible confound where spectral unmixing of HbO and HbR with a linear model leads to quantitative errors in the sO₂ and HbT readings due to the wavelength-dependent light attenuation in the tissue. It has previously been shown that accurate estimation of oxygenation can be obtained by selecting a set of optical wavelengths that 1) minimize spectral coloring effects, 2) avoid ill-conditioning in the inverse problem, and 3) provide a sufficiently high signal-to-noise ratio (SNR).^[22] The wavelength combination in this work was carefully chosen to ensure optimal spectral unmixing of the fOA images. Specifically, we avoided the visible spectral range of 500–600 nm, where significant spectral coloring associated to the strongly wavelength-dependent blood absorption is produced. Furthermore, one of the excitation wavelengths was selected at the isosbestic point of hemoglobin (800 nm), another wavelength at the fingerprint spectral point of HbR (755 nm) and the remaining three wavelengths of 700, 730, and 850 nm exhibiting large absorption difference between HbR and HbO. This ensures a small condition number in the molar absorption coefficient matrix. NIR wavelengths are known to more efficiently penetrate biological tissues, thus providing a good balance between penetration depth and SNR for deep brain imaging.

While human applicability of fOA and FLOT is challenged by the limited penetration of light into living tissues as well as severe acoustic aberrations introduced by the thick skull bone,^[23] the approach can potentially be extended beyond the domain of murine brain imaging. A wearable optoacoustic system has recently been reported for imaging cerebral hemodynamics in awake-moving rats attaining lateral resolution down to 243 µm.^[24] Furthermore, the cerebral cortex of a nonhuman primate was accurately mapped using focused ultrasonic transducers through 2–4 mm thick skull bone.^[25] Fluorescence-based calcium imaging of neural activity has recently become commonplace both in rats^[26,27] and nonhuman primates.^[28,29] Finally, the feasibility of optoacoustic imaging through a whole adult human skull was demonstrated by employing a photon recycler to increase light transmittance, indicating that fOA can potentially be applied for in vivo human brain imaging.^[30]

Accurate estimation of the HRF shape is essential to preventing both false positives and false negatives that may result from a mismatch between the predicted and measured signals. Considering that no specific HRF model currently exists for the different hemodynamic readings provided by fOA, we employed a modified SPM canonical HRF with a TTP of 2.3 s and FWHM of 1.9 s to represent the rodent hemodynamic response.^[31,32] This simplified model was shown to effectively infer active voxels in the mouse brain. Note, however, that the shape of the HRF has been shown to be affected by the particular selection of tasks and brain regions, thus limiting its general applicability. Also, the hemodynamic response to a brief stimulus in the rodent brain is species-dependent and can significantly change during development.^[33] Ideally, the HRF parameters should be directly interpretable in terms of changes in neuronal activity and should be estimated so that statistical power is maximized. However, substantial differences among models in terms of power, bias, and parameter mismatches have been reported in fMRI studies,^[34] which are similarly expected to affect the fOA measurements. Hence, more sophisticated HRF models optimally fitting the different hemodynamic components are expected to bring more flexibility in estimating the functional parameters, which will be addressed in future work.

In summary, the unique capabilities of FLOT for simultaneous imaging of fast calcium signals associated to neural activity and the induced hemodynamic responses across the entire mouse cortex noninvasively can reveal new insights into the basic mechanisms and extent of stimulus-evoked responses. fOA is a versatile tool to resolve multiple hemodynamic parameters in real-time, thus providing additional information that may enhance or complement existing methods of studying the neurovascular and neurometabolic coupling mechanisms. Epifluorescence has high sensitivity and spatiotemporal resolution when it comes to tracking calcium activity in the cortical regions, which served to corroborate and complement the fOA observations. The fOA readings were in turn used to decouple the calcium responses from the slowly varying hemodynamic variations in the fluorescence recordings. The proposed approach does not require cranial windowing or any other invasive surgical procedures. It can also be used in a fully noninvasive manner in mice with nonpigmented scalps. All in all, FLOT is a powerful new brain imaging tool that can be broadly applied for investigations into neurovascular coupling, cerebrovascular and neurodegenerative conditions as well as monitoring of therapies.

Experimental Section Hybrid Epifluorescence and Optoacoustic Tomography (FLOT) System

Hybrid epifluorescence-optoacoustic imaging was achieved with a system that integrates a spherical matrix array transducer (Imasonic SaS, Voray, France) for volumetric fOA imaging with a custom-made fiberscope (Zibra Corporation, Westport, USA) for coregistered optical (epifluorescence) imaging (Figure 1a). The transducer array consists of 512 piezocomposite elements arranged on a hemispherical surface with a 150° angular coverage (1.48π solid angle). Individual elements have 7 MHz central frequency and >80% detection bandwidth. The array features a central 8 mm diameter aperture and three additional 4 mm diameter apertures located at 45° elevation angle and equally spaced in the azimuthal direction. The fiberscope consists of a 1.4 mm diameter optic image guide consisting of 100 000 fibers collecting the fluorescent responses and an illumination bundle composed of 19 fibers having 600 µm core diameter and 0.4 numerical aperture (NA) for optical excitation (Zibra Corporation, Westport, USA). Epifluorescence imaging was performed by inserting the tip of the fiberscope into the 8 mm aperture of the matrix array. At the output end of the image guide, two emission filters (FL514.5-10, Thorlabs, USA) were cascaded to isolate the fluorescence emission from GCaMP, which was subsequently captured with an EMCCD camera (iXon Life, Andor, UK) at 20 fps. The hybrid FLOT system provided a coregistered 12 mm diameter FOV with 40 µm lateral resolution for fluorescence and nearly isotropic 3D spatial resolution down to 113 µm for fOA imaging. Note that for characterizing the GCaMP6f impulse response with high temporal resolution in the millisecond range, the EMCCD camera having superior SNR was replaced with a high-speed CMOS camera (PCO.dimax S1, PCO, Germany) operating at 200 fps. Fluorescence excitation was provided with a continuous wave 488 nm laser (Sapphire LPX 488–500, Coherent, USA). On the other hand, a custom-made fiber bundle (CeramOptec GmbH, Germany) was used to guide a short-pulsed (<10 ns) beam generated with an optical parametric oscillator (OPO) laser (Spit-Light, Innolas Laser GmbH, Germany) through the three lateral apertures of the array. The laser wavelength was rapidly swept between five wavelengths (700, 730, 755, 800, and 850 nm) on a per-pulse basis at 100 Hz pulse repetition frequency (PRF). The pulse energy was ≈11 mJ at the output of the illumination fiber bundle. The generated signals were acquired with a custom-made data acquisition system (DAQ, Falkenstein Mikrosysteme GmbH, Germany) at 40 Megasamples per second (Msps), and recorded raw data transmitted to a PC via Ethernet. Synchronization of the excitation light pulses, the fluorescence and fOA data acquisition, and the electrical paw stimulation was achieved with an external trigger device (Pulse Pal V2, Sanworks, USA).

Animal Models

Female GCaMP6f mice (C57BL/6J-Tg(Thy1-GCaMP6f) GP5.17Dkim/J, the Jackson Laboratory, USA, 6–11 week-old, n = 6) were employed in this study. Animals were housed in individually ventilated, temperature-controlled cages under a 12 h dark/light cycle. Pelleted food (3437PXL15, CARGILL) and water were provided ad-libitum. All experiments were performed in accordance with the Swiss Federal Act on Animal Protection and were approved by the Cantonal Veterinary Office Zurich (license # ZH161/18).

In Vivo Imaging

Mice were anesthetized for the in vivo imaging experiments. Anesthesia was induced with intraperitoneal injection of ketamine (100 mg kg⁻¹ body weight, Pfizer) and xylazine (10 mg kg⁻¹ body weight, Bayer) cocktail. The bolus injections were separated into two injections with a 5 min gap to prevent cardiovascular complications. Maintenance injection was administered every 45 min with a mixture of ketamine (100 mg kg⁻¹ body weight) and xylazine (2.5 mg kg⁻¹ body weight). Both the scalp and skull of the mice were kept intact (n = 4) on the premise that no pigment in the scalp was presented, while the scalp was removed and skull was kept intact for mice with scalp pigmentation (n = 2). Scalp removal was performed after injection of buprenorphine (0.1 mg kg⁻¹ body weight, Temegesic, Indivior, Switzerland) together with hemostatic sponges (Gelfoam, Pfizer Pharmaceutical) to minimize bleeding. Imaging was performed by placing each mouse onto a 3D translation stage in a prone position that facilitated optimal positioning of the region of interest in the brain in the center of the FOV of both modalities. Ultrasound gel was applied on the mouse head to ensure optimal acoustic signal coupling. The mouse head was immobilized using a custom stereotactic frame (Narishige International Limited, London). During the experiment, an oxygen/air mixture (0.1/0.4 L min⁻¹) was provided through a breathing mask. Peripheral blood oxygen saturation, heart rate and body temperature were continuously monitored (PhysioSuite, Kent Scientific) during data acquisition. The body temperature was kept around 37 ℃ with a feedback-controlled heating pad.

The somatosensory electrical stimulation paradigm (i.e., pulse durations and intensity) was adopted from an earlier BOLD fMRI study.^[19] Unipolar rectangular electric pulses of 5 ms duration and 1.0 mA intensity were applied to the left hindpaw at 4 Hz stimulus frequency, 8 s onset time, and 82 s burst intervals (Figure 1b). The electric signals were generated using a stimulus isolator device (Model A365R, World Precision Instruments, USA) fed by an external trigger (Pulse Pal V2, Sanworks, USA). For each sequence, the total duration of the stimulation pulse train was 440 s with the individual stimuli synchronized with the excitation light and data acquisition system. The sequence was repeated 2–5 times per animal. After the experiments, the animals with scalp removed were euthanized, while still under anesthesia, whereas other mice with intact scalp emerged from the anesthesia and fully recovered.

Custom Data Analysis Pipeline for Both Fluorescence and fOA

A data analysis pipeline inspired by SPM12^[35,36] was developed for both fluorescence and fOA data analysis (Figure 1e). To this end, a custom MATLAB script was developed (MATLAB 2019b, MathWorks, USA), including data preprocessing, voxel-wise statistical parametric mapping with the general linear model (GLM) and an ROI-based activation curve plotting.

fOA volumetric data and fluorescence images were preprocessed before being fed into the GLM. Specifically, fOA raw signal matrices were first separated according to the corresponding excitation wavelengths and the signals bandpass filtered between 0.1 and 8 MHz. The images were reconstructed separately for each wavelength using a filtered back-projection algorithm^[37] (100×100×100 µm³ voxel size, 8×8×4 mm³ FOV). Generally, the fOA signal at an arbitrary point r within the imaged tissue volume is mainly resulting from the light absorption by HbO and HbR, which can be approximated as^[38][Image Omitted. See PDF]where Γ is the Grueneisen parameter, λ is the wavelength of the illumination laser beam, ϕ is the light fluence within the tissue, and ε and C are the molar extinction coefficients and concentrations of each absorbing substance (molecule), respectively. fOA signals were first normalized to the laser pulse energy at each wavelength. To reduce the spectral coloring effects,^[22] i.e., the spatially and wavelength dependent light fluence distribution within the tissue, an exponential light attenuation model was applied to compensate for the depth-dependent signal decay. The effective attenuation coefficient (μ_eff) was calculated as^[39][Image Omitted. See PDF]where μ_a is the absorption coefficient and was set to 0.1, assuming homogenous light absorption in the mouse brain at all the 5 wavelengths. $$\mu _{\rm{s}}^{\prime} = a{(\frac{\lambda }{{500( {{\rm{nm}}} )}})}^{ - b}$$ is the reduced scattering coefficient, whereas a = 24.2 cm⁻¹ and b = 1.611 were assumed based on the averaged scattering properties of the brain tissue.^[39]

After light fluence correction, the concentrations of HbO and HbR for a spatially independent Γ were estimated via linear spectral unmixing as [Image Omitted. See PDF]

Subsequently, HbT and sO₂ were calculated from the concentrations of HbO and HbR as [Image Omitted. See PDF][Image Omitted. See PDF]

For fluorescence image preprocessing, the image stack was rearranged to form a 4D volume with the third dimension equal to 1 to make it compatible to the overall processing pipeline. Motion estimation and correction were performed by using SPM12 Realign (estimate) function before the processed data were fed to the GLM for statistical analysis. The custom stereotactic frame helped mitigating the motion artifacts. The relative center of mass displacement in each trial typically remained below 0.2 voxel, thus ensuring negligible motion-induced errors (false activations). Spatial smoothing (Gaussian kernel FWHM = 0.3 mm) was further applied to reduce random fluctuations and reconstruction-related image artifacts. A high-pass filter of 1/135 Hz was then used to detrend the fOA data. For fluorescence data, a lowpass filter with 0.5 Hz cutoff frequency was applied to retrieve the slow hemodynamic response encoded in the fluorescence signal and a bandpass filter (3–8 Hz) was applied to extract the fast calcium transients.

In the GLM, the regressor was obtained by convolving the stimulation paradigm (Figure 1b) with either the measured GCaMP CRF (for fluorescence data) or a modified SPM canonical HRF (for fOA data) having a TTP of 2.3 s and FWHM of 1.9 s to represent the rodent hemodynamic response.^[31,32] The regressor together with a constant vector formed the design matrix. The GLM is expressed in the matrix format via^[40][Image Omitted. See PDF]where Y is the column vector of observations representing a signal voxel sampled at successive time points, ε is the column vector of error terms, β is the column vector of parameters, and X is the design matrix. The estimated parametric map is calculated as [Image Omitted. See PDF]

The statistical significance of the event-evoked responses in the observations was further evaluated with a contrast vector c = [1, 0]^T. The t statistic map (i.e., t-map) and probability levels (i.e., p-values) were calculated as the statistical inferences. Note that for rendering the activation map for fOA, the false discovery rate (FDR) controlling^[41] with p <0.05 was introduced to the t-map. For the fluorescence recordings, the parametric map (Betas in Equation (7)) was used directly to infer the activation.

The activation map and the corresponding structural image acquired at 850 nm excitation wavelength were overlaid to the Allen Mouse Brain Atlas 2015. The atlas volume has 320 × 456 pixels in each coronal plane and consists of 528 slices with 25 µm isotropic voxel size. Due to the lack of common features and utterly different intensity distribution, unbiased automated registration of fOA with MRI remains challenging.^[42] In this work, the commonly used manual landmark based coregistration method^[43] was employed for the alignment in Amira (version 5.4, Thermo Fisher Scientific, USA). First, the structural volumetric data of the mouse brain was reconstructed with the averaged signal matrix acquired at 850 nm. Second, this structural volume and the mouse brain atlas were loaded into Amira for prealignment. A coarse transformation with respect to translation and rotation was tuned according to the landmarks in both datasets. Finally, the position of the structural optoacoustic volumetric data was finetuned based on cross-sections along the three planes (perspectives) so that landmarks such as the brain curvature and vessel structures matched the MRI brain atlas. The calculated transform information in the transform editor in Amira facilitated registering all the other reconstructed parameters, such as HbO, HbR, HbT, and sO₂ from the same mouse to the brain atlas. The overlaid activation map was subsequently rendered and displayed along different planes.

The time courses of brain activation were retrieved from the raw 4D data, namely, 3D volume time series for fOA and 2D image series for fluorescence data. The calculated activation map served to infer responsive ROIs in both contralateral and ipsilateral hemispheres. A time window including 10 s prestimulation, 8 s onset, and 50 s poststimulation was selected. Baseline signals of each stimulation cycle were calculated by averaging the 10 s prestimulation time window. This was used for calculating fractional signal changes, which were averaged across different stimulation cycles. All the trials acquired at slightly different anesthesia depths and from different mice were averaged to reveal the stimulus-evoked brain activation in a more robust way. Statistics on activation intensity, TTP, FWHM as well as decay time (T_1/2) were further performed across all the trials. The statistical significance between the averaged values from the contralateral side and their ipsilateral counterparts was evaluated with two-tailed paired t-test.

Statistical Analysis

Statistical analysis was performed with MATLAB (2019b, MathWorks, USA). Results were presented as mean ± SEM for at least three independent measurements. Voxels in the activation map were considered statistically significant with respect to a null hypothesis of no activation with one-sample t-test after false discovery rate (FDR) correction with p < 0.05. Values outside the [1–91%] range of activation intensity, TTP, FWHM as well as decay time (T_1/2) were excluded before applying statistics. The significance of difference on functional parameters between the contralateral hemisphere and the ipsilateral counterpart was tested using two-tailed paired t-test with p < 0.05 across a total of 80 stimulation trials.

Acknowledgements

The authors acknowledge grant support from the US National Institutes of Health (No. UF1 NS107680), the European Research Council (No. ERC-2015-CoG-682379), and the Swiss National Science Foundation (No. 310030_192757).

Conflict of Interest

The authors declare no conflict of interest.

Author Contributions

Z.C. and D.R. conceived the concept and devised the study. Z.C. and Q.Z. performed the animal experiments. Z.C. and I.G. performed image reconstructions and data analysis. X.L.D.B., R.N., and M.R. helped with the experiments and provided helpful discussions. S.S. and D.R. supervised the work. All authors contributed to writing and revising the manuscript.

Data Availability Statement

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