Abstract
A major challenge for human neuroimaging using functional MRI is the differentiation of neuronal excitation and inhibition which may induce positive and negative BOLD responses. Here we present an innovative multi-contrast laminar functional MRI technique that offers comprehensive and quantitative imaging of neurovascular (CBF, CBV, BOLD) and metabolic (CMRO2) responses across cortical layers at 7 Tesla. This technique was first validated through a finger-tapping experiment, revealing ‘double-peak’ laminar activation patterns within the primary motor cortex. By employing a ring-shaped visual stimulus that elicited positive and negative BOLD responses, we further observed distinct neurovascular and metabolic responses across cortical layers and eccentricities in the primary visual cortex. This not only indicates feedback inhibition of neuronal activities in both superficial and deep cortical layers underlying the negative BOLD signals in the fovea, but also illustrates the neuronal activities in visual areas adjacent to the activated eccentricities.
1. Introduction
High resolution functional magnetic resonance imaging (fMRI) at ultrahigh field has significantly enhanced our capability to study neurovascular and metabolic activities in-vivo at the mesoscopic scale. A major challenge for fMRI, however, is the differentiation of neuronal excitation and inhibition which may induce positive and negative BOLD responses, depending on the complex interplay between cerebral blood flow (CBF), cerebral blood volume (CBV), and cerebral metabolic rate of oxygen (CMRO2) (Huber et al., 2019; Shmuel et al., 2002). The relationship between changes in neuronal activities and changes in fMRI signals has long been a focus of debate, especially regarding negative BOLD responses (NBR). Shmuel et al. hypothesized two possibilities for the formation of NBR: 1) from the reduction of neuronal activity; and 2) from hemodynamic changes independent of neuronal activity (Shmuel et al., 2002). NBR has been shown to be associated with reduction in CBF and neuronal activity (Boorman et al., 2010; Devor et al., 2007; Northoff et al., 2007; Shmuel et al., 2006). However, other studies have found that NBR may occur with heightened neuronal activity (Goense et al., 2012; Mullinger et al., 2014; Schridde et al., 2008). These contradicting findings highlight the complexity of the NBR and BOLD signals. As illustrated in Figure 1, different scenarios for changes in CBF, CBV and CMRO2 could lead to positive or negative BOLD responses respectively. This is a major challenge to precisely interpret the underlying neural mechanisms using existing laminar fMRI techniques that measure a single hemodynamic parameter (Fracasso et al., 2018; Han et al., 2021; Huber et al., 2017; Shao et al., 2021).
In vivo imaging of the full set of hemodynamic parameters including CBF, CBV and BOLD as well as metabolic parameters such as CMRO2 provides comprehensive characterization of neurovascular responses, and precise inference of the neuronal activities. In particular, changes in CMRO2 are closely related to changes in neuronal activity measured by electrophysiology across cortical layers (Herman et al., 2013). Yang et al first proposed a technique at 3 Tesla to simultaneously measure CBF, CBV and BOLD contrasts by combining pulsed arterial spin labeling (ASL) with vascular space occupancy (VASO) (Yang et al., 2004), which was further developed by Cheng et al (Cheng et al., 2017) at 3T and adapted for 7T by Krieger et al (Krieger et al., 2015). However, pulsed ASL is susceptible to vascular contaminations due to arterial transit effects (Shao et al., 2019). In addition, the spatial resolution and/or image coverage of the above techniques does not allow laminar fMRI. Optical imaging in animals provides evidence that CBF can be precisely regulated on an extremely fine scale at the level of capillaries (Chaigneau et al., 2003; Hall et al., 2014; Hamilton et al., 2010; Peppiatt et al., 2006). Therefore, it is highly desirable to develop high resolution multi-contrast fMRI at ultrahigh field for simultaneous hemodynamic and metabolic imaging at laminar level.
Here we present an innovative multi-contrast laminar fMRI technique at 7 Tesla, which enables simultaneous measurement of CBF, CBV, T2-weighted BOLD, and CMRO2. The key innovation of the proposed method lies in the use of pseudo-continuous ASL (pCASL) with intracranial labeling and a zoomed 3D gradient-and-spin echo (GRASE) readout, coupled with a non-selective inversion pulse at an optimized post-labeling delay (PLD), to generate concurrent ASL perfusion and VASO contrast for CBV. T2-weighted BOLD signals, which has improved specificity than T2*-weighted BOLD, were acquired after the ASL/VASO acquisition with GRASE readout (Fig. 2, details in Section 4.1). Compared to existing laminar fMRI methods, our technique offers several unique advantages: 1) CBF, CBV, T2-weighted BOLD contrasts are simultaneously measured to enable detailed and precise characterization of hemodynamic changes across cortical layers at 7T; 2) Concurrent CBF, CBV, T2-weighted BOLD contrasts enable accurate estimation of the laminar profile of CMRO2 using calibrated fMRI with a respiratory challenge (breath-hold) and an innovative algorithm to estimate laminar profiles of calibration parameters; 3) The use of pCASL with optimized PLD offers higher SNR than pulsed ASL, while minimizing pial arterial signals to ensure that the measured perfusion signal arises from capillaries and brain tissue (Shao et al., 2021).
In this study, we performed 4 experiments to demonstrate the capability of the proposed multi-contrast fMRI for comprehensive neurovascular and metabolic imaging to understand the underlying neuronal mechanisms associated with positive and negative BOLD responses: 1) Validating the specificity of our proposed technique by measuring activations in both deep and superficial layers of the primary motor cortex (M1) induced by finger-tapping (FT) tasks; 2) Estimating calibrated fMRI parameters through a breath-hold (BH) hypercapnia method, enabling the calculation of CMRO2 (Davis et al., 1998) from simultaneously measured CBF, CBV, and BOLD signals; 3) Conducting eccentricity mapping using the population receptive field (pRF) approach (Greene et al., 2014), a technique that quantifies visual field representations in the brain to understand visual processing and cortical organization; 4) Performing a ring-shaped visual stimulation that elicited positive and negative BOLD responses at the projected visual area and the fovea respectively, while investigating layer-dependent CBF, CBV, T2-weighted BOLD, and CMRO2 activities and their variations across eccentricities, aiming to understand the underlying neuronal activities associated with positive and negative BOLD responses.
2. Results
2.1. Double-peak activation in M1
FT task elicits neuronal activities across both superficial and deep cortical layers: the superficial layers are engaged by sensory inputs from the somatosensory (S1) and premotor cortices, while the deep layers are involved in motor outputs (Supplemental Fig. S1). The FT task, commonly utilized in laminar fMRI studies for its efficacy in probing cortical layer-specific activities, served as a key component in our study to validate the sensitivity of our proposed technique. Fig. 3 shows FT-induced signal changes and laminar profiles for CBF, BOLD, VASO and CBV (CBV without BOLD correction was also shown for comparison). Employing this technique, we observed increases in CBF and BOLD signals, alongside decreases in VASO signals along the omega-shaped M1 hand-knob, with double-peak patterns observable upon visual inspection (Fig. 3 Middle). The right panels of Fig. 3 display average laminar profiles from four participants, showing up to 126±22% increase in CBF and a 28±2% increase in CBV in the deep layers, and up to 111±31% increase in CBF and a 46±7% increase in CBV in the superficial layers. The observed double-peak activation patterns were statistically significant for both CBF (R2diff = 0.79, P < 0.001 for %CBF and R2diff = 0.54, P = 0.005 for ΔCBF) and CBV (R2diff = 0.75, P < 0.001 for CBV without BOLD correction, R2diff = 0.35, P = 0.026 for BOLD corrected CBV), while the BOLD response predominantly peaked in superficial layers, attributed to significant venous contamination (R2diff = 0.11, P = 0.14). These findings align with our hypothesis and verify the intricate layer-specific neurovascular dynamics during motor execution tasks.
2.2. Calibrated fMRI with breath-hold (BH)-induced hypercapnia
Using the calibrated fMRI approach, we employed BH hypercapnia to estimate M and β parameters and calculate CMRO2, which provides a surrogate index of the underlying neuronal activity (detailed in Method section 4.2.2). By assuming that CMRO2 remains unchanged, the simplified Davis model enables the calculation of M and β from measured time series of CBF, CBV and BOLD signals across cortical layers.
Fig. 4 shows neurovascular responses induced by BH hypercapnia. The observed global increases in both BOLD and CBF signals, accompanied by a decrease in VASO signals (indicative of CBV increase), reveal the dynamic nature of blood supply and its regulation in response to elevated CO2 levels (Fig. 4 top). The laminar profiles further show an up to 45±16% increase in CBF peaking in the deep layers and up to a 16±3% increase in CBV peaking in the superficial layers (Fig. 4 middle). This layer-specific vascular response may reflect the complex architecture of cortical vasculature probed by VASO and ASL. CBV changes reflect hypercapnia-induced vasodilation in both macro- and micro-vasculature, where CBV sampled from macro-vasculature (arterioles) are hypothesized to increase from deep to superficial layer but not for the micro-vasculature compartment (capillaries) (Schellekens et al., 2023). ASL signal primarily originates from capillary vessels, and peak of CBF increase in the deep layer is consistent with the higher microvascular density in middle to deep cortical layers (Shao et al., 2021). The observed patterns of BOLD, CBF, and CBV signals indicate a layer-specific vascular response to hypercapnia, implying that the profiles of M and β are also likely to vary across cortical layers. The bottom panel of Fig. 4 shows the laminar profiles of M and β. Our findings provide a detailed laminar profile for both M and β, which enables the calculation of CMRO2 across cortical layers in subsequent experiments.
2.3. Distinct neurovascular and metabolic responses across cortical layers and eccentricities in the primary visual cortex
Our visual stimulation experiment utilized a ring-shaped stimulus within the 4-6 degree of eccentricity range (Supplemental Fig. S2) while the participants were asked to perform a color detection task on the central fixation point. We measured layer dependent CBF, CBV, BOLD and CMRO2 responses for each eccentricity from 0 to 8 degrees. The definition of ROI is from the pRF experiment (Supplemental Fig. S3). Fig. 5 shows our hypothesized neural circuits for the stimulated area, fovea and their adjacent regions of V1 and associated layer-dependent response patterns.
Fig. 6 shows the activation map of a representative participant and the schematic diagram of the ROI definition. Fig. 7A shows the average laminar profile of the changes of CBF, CBV, BOLD and CMRO2 responses of nine participants across cortical layers and eccentricity ROIs associated with stimuli, fovea, peripheral and transition regions. For different eccentricity areas, we performed multi-contrast fMRI analysis separately to derive the laminar profiles of CBF, CBV, BOLD and CMRO2 for different eccentricity ranges respectively.
2.3.1. Feed-forward visual stimuli (4 - 6 degree of eccentricity)
From the activation maps (Fig. 6), we observed a significant increase in CBF and BOLD signals, and a decrease in VASO within the 4 - 6 degree of eccentricity range. This matches well with the areas that our stimulation activates. We observed that the multi-contrast fMRI responses of the 4 - 6 degree of eccentricity range (Fig. 7A) were significantly greater than the respective responses in other eccentricities (P<0.001 for BOLD, CBF, CBV and CMRO2). Although the peak activation from feed-forward stimuli is commonly expected to occur in the middle towards deep layers (Fig. 5A&B) (Lamme et al., 1998), we observed a slight shift between the profiles of CBF and CBV responses (Fig. 7 top). The peak of the laminar profile in the CBF signal is located in the deeper middle layer, while the peak of the CBV signal is in the more superficial part of the middle layer. Despite the shifted peaks observed in CBV and CBF profiles, CMRO2 increase was relatively consistent across cortical layers, with a large peak in the middle to deep layers and two subtle bumps in middle and superficial layers.
2.3.2. Negative BOLD responses in fovea (0 – 1 degree eccentricity)
Previous studies reported negative BOLD responses in the V1 fovea region when visual stimuli were presented in the peripheral visual field (de la Rosa et al., 2021). In this study, participants were asked to perform a color judgment task on the fixation point. Previous research has indicated that when an irrelevant distractor is presented with a target stimulus, the capture of attention by the distractor impairs processing of the target as reflected by slower reaction times to the target (Schreij et al., 2008; Theeuwes, 1992).
Moreover, a study by Theeuwes et al (1998) has shown that the attention can also be captured by a new-appearing irrelevant distractor (Theeuwes et al., 1998). There is also evidence from fMRI and electrophysiological experiments that interfering stimuli (destractors) can automatically capture attention and thus affect the processing of task stimuli (Hickey et al., 2006; Spinks et al., 2004). However, the exact neural pathways underlying this interference remain unknown.
Interfering stimuli will engage participants’ bottom-up attention processes and thus interfere with participants’ performance on task stimuli. A saliency map is the bottom-up contribution to the deployment of exogenous attention. The mainstream view is that a salience map should integrate various visual information and be formed in the high-level cortex (Itti & Koch, 2001). In contrast, Li proposed that V1 creates a salience map (Li, 2002). Previous studies have indeed found that V1 (Zhang et al., 2012) and superior colliculus (SC) (White et al., 2017) can spontaneously generate a salience map. In view of these two perspectives, we propose two hypotheses (Fig. 5C). One possibility is that attention decreases through feedback connections from upstream cortex, cortico-cortical feedback (Felleman & Van Essen, 1991) (Hypothesis 1 in Fig. 5) targeting the superficial and deeper layers. Another possibility is a SC->TRN->LGN->V1 feedforward pathway (Kolmac & Mitrofanis, 1998; McAlonan et al., 2008) (Hypothesis 2 in Fig. 5) which involves decreased neuronal activity in the middle layer.
For the fovea region (0-1 degree eccentricity range), we found a decrease in CBF and BOLD signals, and an increase in VASO signal in the activation map (Fig.6), consistent with findings of the previous study (de la Rosa et al., 2021). We performed two-peak activation tests on the laminar profile in the fovea region (black lines in Fig. 7A). We found significant two-peak activation patterns in CBF (R2diff = 1.09, P < 0.001) and CMRO2 (R2diff = 0.78, P = 0.003) signals, while BOLD decreases were relatively flat across cortical layers (R2diff = -0.38, P = 0.25). In other words, the decrease in CBF and CMRO2 is greater in superficial and deep layers than in middle layer.
2.3.3. Adjacent region (7 – 8 and 2 – 3 degree eccentricity)
In addition to retinal topographic projection region and fovea, the proposed technique allowed investigation of neuronal activities in peripheral (7-8 degree) and transitional (2-3 degree) regions. Since the fMRI activation patterns of these two areas are relatively similar, we pooled the two together as the adjacent regions to study the activation pattern. There is a distance between the adjacent regions and the stimulus activation region that cannot be ignored, we therefore expect that they will be minimally affected by lateral inhibition when activated by visual stimulation. Previous study used injections of the CVS-11 strain of rabies virus to examine disynaptic long-range horizontal connections within macaque monkey V1 and demonstrated that long-range horizontal connections mainly occur on layer 6 and are responsible for integration across a large region of the visual field (Liu et al., 2014). A recent study found neurons corresponding to this integration in layer 6 of V1 mainly function as GABA inhibitory neurons (Srivastava et al., 2022). Therefore, we predicted that these two areas would exhibit suppressed activation in deep layer when activated by visual stimulation (Fig. 5D).
For adjacent regions, it can be clearly seen from the activation map in Fig. 6 that these areas show a decrease in BOLD and CBF signals and an increase in VASO signal (indicative of CBV decrease), which is opposite to the stimulus activation pattern. Furthermore, the laminar profiles of adjacent regions shown in Fig. 7A confirmed our prediction of lateral inhibition. The CBF and CMRO2 signals show a trend of lower signal intensity in the deep layers than in the superficial and middle layers in the eccentricity range of 2-3 degrees and 7-8 degrees.
2.4. The CBF signal has the closest trend to the change of the CMRO2 signal
Fig. 7 shows that CBF signal changes have the closest trend to the CMRO2 signal changes compared to those for CBV and BOLD. We performed correlation analyses between CMRO2 and CBF/CBV/BOLD signals across eccentricities and cortical layers (Fig. 8). Compared to CBV (r = 0.156, P = 0.008) and BOLD (r = 0.147, P = 0.013), CBF (r = 0.782, P < 0.001) showed the highest correlation with CMRO2 signal.
3. Discussion
In this study, we introduced an innovative laminar multi-contrast fMRI technique to study the neuronal mechanisms underlying the positive and negative BOLD responses at 7T. We demonstrated the high specificity in detecting distinct CBF, CBV, and BOLD signal profiles across cortical layers utilizing the FT experiment. The FT task induced double-peak activation patterns in CBF and CBV signals within the M1 hand-knob area, reflecting the underlying neurocircuits involved in receiving and transmitting signals (Huber et al., 2017). Such high specificity is crucial for our study to accurately estimate laminar CMRO2 profiles and examine the complex neurovascular and metabolic responses associated with both positive and negative BOLD signals.
There is a limited number of studies utilizing concurrent high resolution CBF, CBV and BOLD measurements to investigate the neuronal activities (Cheng et al., 2017; Krieger et al., 2015; Yang et al., 2004). Moreover, previous calibrated fMRI studies (Blockley et al., 2013; Griffeth & Buxton, 2011; Krieger et al., 2014) generally assumed that changes in CBV and CBF follow Grubb’s power law (CBV ∝ CBFα, α = 0.38 (Grubb Jr et al., 1974)). The Grubb’s power law allows CMRO2 estimation from BOLD signals using either CBF or CBV data. However, given the complexity of the neurovasculature, the relationship between CBV and CBF becomes complex, particularly when examining neuronal activity in different brain regions (Chen & Pike, 2009). By applying the Grubb’s power law to our concurrent CBF and CBV measurements, we found an average α = 0.36±0.10 across all eccentricities and cortical layers, aligning well with the literature (Grubb Jr et al., 1974). However, notable variations of the Grubb’s constant were observed across cortical layers as well as during specific tasks (Fig. S4 A). Within the stimulation region, α peaks towards superficial layer (∼0.37) while lower α values were found in deep layers (∼0.21), which reflects the differences of the underlying microvasculature (Chen & Pike, 2009). Difference profiles of α values were also observed between the stimulation region and the fovea, which further emphasizes that assuming a constant α value may be inaccurate for assessing the underlying neuronal activities for mesoscopic fMRI studies. In our study, when comparing CMRO2 values measured from concurrent CBF, CBV, and BOLD measurements, we observed significant underestimations of CMRO2 within the stimulation region by 68.0±40.6% when using BOLD and CBV, or by 27.2±18.5% when using BOLD with only CBF data (Fig. S4 B).
Several approaches were developed to estimate the CMRO2 (Lu et al., 2004; Rodgers et al., 2016; Xu et al., 2009), while the calibrated fMRI technique based on the Davis model has been commonly utilized (Davis et al., 1998). This model demonstrates the correlation between CMRO2 changes and variations in CBF, CBV and BOLD signals via two parameters, M and β (Davis et al., 1998). M denotes the maximal change in the BOLD signal achievable based on baseline deoxyhemoglobin levels, while β accounts for the impact of localized magnetic susceptibility on R2* differences, emphasizing the vascular distinctions between capillaries and larger vessels. Commonly, M and β can be estimated under the assumption that the CMRO2 remains constant during hypercapnia. However, the slow dilation of vessels induced by breath-holding (BH) may lead to a mixed physiological state between BH and resting-state measurements. We adjusted acquisition timing based on BH-signal delay measured from a pre-scan, however, resting state CBF, CBV and BOLD signals can still be over-estimated due to residual hypercapnia. To mitigate this delay and potential fluctuations of BH hypercapnia levels, we considered all scans as BH condition (i.e. resting-state signal is zero to low-level BH) and simultaneously estimated M, β, as well as baseline CBF, CBV and BOLD values from all measurements. Supplemental Fig. S5 shows an example of dynamic CBF, CBV, and BOLD signals from experiment 2. Despite signal fluctuations across the measurements, our method can effectively predict BOLD signals using measured CBF and CBV, along with estimated M and β values and RSME between measured and predicted BOLD signals was 0.0050±0.0019 across all participants and cortical layers. BH-induced CMRO2 changes were averaged to be 0.3% (ranging from -0.96% to 2.39%) across eight cortical layers (Eq. [5]), which indicates minimal CMRO2 variation due to BH and further validates the M and β measurements.
In this study, we were able to measure detailed laminar profiles for M and β (Fig. 4). Since GRASE acquisition was utilized to generate T2-weighted BOLD contrast, we observed a plateau from middle to superficial layers with relatively lower amplitudes (∼5.5%) compared to EPI BOLD signals (∼14%) (Guidi et al., 2016; Schellekens et al., 2023). To the best of our knowledge, no study has reported detailed laminar profiles of β at 7T, while β is usually assumed to be 1 for UHF studies when the venous signal becomes negligible in extravascular BOLD signals (Griffeth & Buxton, 2011; Kida et al., 2000; Krieger et al., 2014). We found β gradually increases from deep (0.97) to superficial layer (1.22) with an average value of 1.10±0.09 across layers. Measuring laminar profiles of M and β provides insights into the underlying neurovascular coupling mechanisms and is crucial for accurate CMRO2 estimations in laminar fMRI studies.
For the visual stimulation region in V1, activation map and laminar profile across participants show the strongest response and a pattern where middle layer has a stronger response than superficial and deep layers. It is known that this neural activity mainly arises from the feedforward pathway of the LGN. Sublamina 4Cα receives mostly magnocellular input from the LGN, while layer 4Cβ receives input from parvocellular pathways. After processing, the signal is transmitted to layers 4A/B and 2/3. Layer 2/3 then sends afferents to layer 5/6 and the higher-level visual cortex (Fig. 5A., Accordingly, we predicted that this neural activation should show a pattern of greater response in the middle towards deep layer (4) with moderate increase in superficial (4A/B and 2/3) and deep layers (2/3) (Fig. 5B). While our results matched this hypothesis very well, we observed a slight shift between the profiles of CBF and CBV responses (Fig. 7 top). Specifically, for CBV measured by VASO, the signal changes reflect contributions from all vasodilation vessels associated with neuronal activity, including pial vessels, penetrating arterioles and capillaries (Huber et al., 2019). We observed that CBV increase peaks towards the mid-to-superficial layers, suggesting that pial arteries and penetrating arterioles in the mid-to-superficial cortex may exhibit greater vasodilatory responses compared to those in deeper layers (Duvernoy et al., 1981). The signal of labeled blood in ASL is considered to primarily originate from capillaries. We employed pCASL with higher SNR than pulsed ASL (Huber et al., 2019; Kashyap et al., 2021). We have measured the arterial transit time (ATT), which was found to rang from 0.9 to 1 sec from superficial to deep layers in V1, and applied sufficiently long PLD to minimize pial arterial artifacts (Alsop & Detre, 1996; Alsop et al., 2015), thereby offering a “clean” depiction of CBF changes within the capillaries that are closely related to neuronal activation. The laminar profile of CBF activation indicates greater increase within the mid-to-deep layers compared to baseline, which is consistent with previous laminar perfusion fMRI studies (Kashyap et al., 2021; Shao et al., 2021) and findings from animal research on microvascular density (Jin & Kim, 2008). However, CMRO2 increase was relatively consistent across cortical layers, with one large peak in middle to deep layers and two subtle bumps in middle and superficial layers. The laminar profile of CMRO2 suggests the feed-forward visual stimuli involve the primary visual signal from the LGN predominantly targeting layer IV, with slight extensions into layers V and VI (Clay Reid & Alonso, 1995) and the potential projections from layer IV to superficial layers (Callaway, 2004) (Fig. 5A, B).
For the fovea region, we found that the decrease in CBF and CMRO2 is significantly greater in superficial and deep layers than in middle layer, while BOLD signal shows no difference between layers. The reason for the lack of such results for the BOLD signal may be due to its lower layer specificity (Huber et al., 2019). With the estimation of the laminar profile of CMRO2, which is considered a surrogate of neuronal activity, we conclude that the negative BOLD response in the fovea of V1 is due to neuronal inhibition through a feedback pathway from higher visual areas. Specifically, the presented interference stimulus integrates a salience map in high-level brain areas, and then inhibits the original attention response of the primary visual cortex task area through the feedback pathway.
For the adjacent regions, it can be clearly seen from the activation map (Fig. 6) and laminar profile (Fig. 7) that deep layer signals are more inhibited than superficial layer. This is consistent with previous research showing that stimulating neurons in the activated region inhibits activation of neurons in both regions via deep long-range neurons (Liu et al., 2014; Srivastava et al., 2022). In contrast, the BOLD and CBV signals do not show this obvious trend. The high resolution of the proposed technique allowed investigation of hemodynamic and metabolic activities within the transitional and peripheral areas of positive and negative BOLD activations. To the best of our knowledge, such investigation hasn’t appeared in literature.
In this study, we found that the CBF (r = 0.782) signal has the highest correlation with the CMRO2 signal compared to CBV (r = 0.156) and BOLD (r = 0.147) signals. Despite the significant correlations between CBV/BOLD and CMRO2, the directions of changes in CBV/BOLD to positive and negative CMRO2 remains exclusive as indicated by large spread along y-axis in Fig. 7. Synchronous increase or decrease in CBF and CMRO2 were found in 81.6% of all ROIs across eccentricities and cortical layers, while only 53.1% and 47.9% of ROIs had the same direction of CBV and BOLD changes along with CMRO2. These results suggest CBF can be the best proxy to CMRO2 and the underlying neuronal metabolic activity (Gsell et al., 2000).
There are several limitations in this study. The TR was relatively long due to the use of 1280 msec pCASL labeling with a 1160 msec PLD, limiting the capability to capture dynamic responses of CBF, CBV, and BOLD signals (Kashyap et al., 2021; Kim & Kim, 2010). Despite this, pCASL offers advantages including enhanced SNR and less contamination from pial arterial signals compared to pulsed ASL (Alsop et al., 2015), providing higher specificity in laminar CBF profiles (Shao et al., 2021). Future studies might benefit from exploring a dithered acquisition scheme to improve temporal resolution (Kim et al., 2020). Additionally, we employed a single-shot GRASE acquisition method to improve acquisition efficiency, however, relatively thick slices (2.2 mm) were used to provide sufficient spatial coverage. For visual experiments, imaging planes were placed nearly parallel to the calcarine sulcus to preserve high spatial resolution (1 mm) across cortical layers, which allows to detect ∼2-3 independent voxels across V1 cortical layers, respectively (Fischl & Dale, 2000). The specificity of the multi-contrast fMRI can be improved with higher spatial resolution utilizing accelerated acquisition and constrained reconstructions (Chang et al., 2017; Park et al., 2021; Spann et al., 2020; Zhao et al., 2023). Another limitation is the lack of direct validation of our observations, such as the behaviors within the fovea region through methods such as positron emission tomography or electrophysiological recordings (Hickey et al., 2006; McGregor et al., 2020; Spinks et al., 2004). The distinctive patterns observed in CBF, CBV, BOLD, and CMRO2 across cortical layers and eccentricities need to be validated in the future studies to associate with underlying neuronal activities.
In conclusion, we presented an innovative multi-contrast laminar functional MRI technique that offers comprehensive and quantitative imaging of neurovascular and metabolic responses across cortical layers at 7 Tesla. We demonstrated that the proposed sequence provides quantitative CBF, CBV, BOLD and CMRO2 measurements with high resolution and specificity to delineate distinct neurovascular and metabolic responses across cortical layers and eccentricities in the primary visual cortex. This novel technique allows the differentiation of neuronal excitation and inhibition underlying positive and negative BOLD responses, illustrating the associated neural circuit activities across cortical layers.
4. Methods
4.1. MRI pulse sequence
The diagram of MRI pulse sequence is shown in Fig. 2A. Perfusion contrast was generated by pCASL, and a non-selective hyperbolic secant (HS)-inversion pulse was added after pCASL pulse train to suppress background signal and generate VASO contrast (Shao et al., 2021). Due to limited B1+ at bottom of brain at 7T, we assume fresh in-flow blood only experience one inversion and control images acquired at PLD = 1160 msec, which nulls in-flow blood signal (T1,blood = 2.1sec (Rane & Gore, 2013)), exhibit VASO contrast. The PLD = 1160 msec is slightly longer than ATT in M1 and V1, which provides desirable SNR for functional ASL study (Shao et al., 2021). Concurrent T2-BOLD signals were acquired 500 msec after the ASL/VASO acquisition (Fig. 2A). Intracranial vasculature was revealed by maximal-intensity-projection (MIP) of the MP2RAGE (0.7mm3) images, and the pCASL labeling plane was placed above the circle of Willis (CoW) and simultaneously perpendicular to the anterior cerebral artery (ACA), middle cerebral artery (MCA) and posterior cerebral artery (PCA) (Fig. 2B). This labeling location took advantage of the higher B1 and more homogeneous B0 field around the center of the brain. Imaging parameters were: FOV = 80×40 mm2, 14 slices, 1-mm2 in-plane resolution, 2.2 mm slice thickness, 3 segments along phase direction, TE = 17.2 msec, labeling duration = 1280 msec, TR = 4000 msec. Volume TR was 24 sec for one CBF and VASO image, and 12 sec for one BOLD image (Fig. 2C).
4.2. MRI experiments
We conducted four experiments on a 7T Siemens Terra MRI scanner with 1Tx/32Rx head coil to validate the high sensitivity and specificity of the proposed technique in detecting neuronal activity and reveal the distinct neurovascular and metabolic responses across cortical layers and eccentricities under positive and negative BOLD signals induced by a ring-shape visual stimulus. Experiment details are shown as following:
4.2.1. Experiment 1 – Finger tapping
Four runs with interleaved blocks of rest and sequential dominant-hand FT were acquired. Each run consisted of five blocks of interleaved rest and FT tasks (48 sec per block, 8 min each acquisition block). Total scan time was 32 min. As shown in Supplemental Fig. S1, we expect FT to induce neuronal activity in both superficial layers (sensory input from primary sensory cortex (S1) and premotor cortex) and deep layers (spinal output), and increased CBF and decreased VASO signals in both deep and superficial layers.
4.2.2. Experiment 2 - Laminar profiles of M and β (Davis model) estimated by breath-hold (BH) induced hypercapnia
M and β are two unknown parameters which are likely to change across participants and cortical layers. Multi-contrast fMRI data were acquired with 12 cycles of interleaved 24 sec rest and 24 sec BH with visually presented instruction, with two runs collected in a total of 19 min 20 sec. Under the assumption that CMRO2 does not change during hypercapnia, we can directly estimate M and β using concurrently measured CBF, CBV, and BOLD signals between resting state and BH using the following simplified Davis model (Eq. [1]):
4.2.3. Experiment 3 - Eccentricity mapping with a pRF stimulus
Eccentricity mapping was conducted through pRF analysis using a flickering checkerboard pattern that sweeps (30 steps, TR=1.2s) across the screen, sequentially stimulating different visual field regions with a 15-degree rotation after each sweep (Supplemental Fig. S2 top). We collected two pRF runs using 1.6 mm³ isotropic EPI (SMS-2, in-plane GRAPPA-2, TR=1.2s), and acquired 435 volumes in 8 min 42 sec for each run. Eccentricity masks ranging from 0 to 8 degrees within V1 were then reconstructed by estimating each voxel’s receptive field from the BOLD signal changes.
4.2.4. Experiment 4 – Partial visual field activation stimulus
Finally, we conducted fMRI experiments using visual stimulus consisted of a mean-gray background with a ring-shaped sector subdivided into 12 subsectors, covering an eccentricity range of 4° to 6° (Supplemental Fig. S2 top). Each subsector contained a high-contrast (100%) radial grating (1 cycle per degree, cpd) that reversed contrast at a rate of 4Hz. Each run comprised six blocks of 48 sec of visual stimulus followed by 48 sec of rest, and four runs were acquired in a total of 39 minutes.
4.3. Human participants
All participants in this study were healthy without known neurological disorders and refrained from caffeine 5 hours before the scan. Experiment 1 included four participants (3 male, age = 25.3±3.2years). A total of nine participants were recruited for experiment 2 to 4 (3 male, age = 26.0±6.3 years). All participants provided written informed consents according to a protocol approved by the Institutional Review Board (IRB) of the University of Southern California. Head motion was minimized by placing cushions on top and two sides of head and taping participant’s chin to coil.
4.4. Data processing
4.4.1. ASL data prep-processing
A separate proton density weighted M0 image with the same zoomed GRASE readout was acquired for ASL calibration. Dynamic ASL images were realigned and coregistered to M0 using SPM12 (Functional Imaging Laboratory, University College London, UK). Control and label images were pairwise subtracted to obtain raw perfusion images. Signal decay in k-space was prominent in partition direction due to long echo train length and short arterial blood T2 at 7T (68 msec (Krishnamurthy et al., 2014)). We implemented a singular value decomposition (SVD) based deblurring algorithm to reduce blurring along slice direction without amplifying noise (Shao et al., 2021). CBF was calculated from the relative perfusion signal (normalized by M0) using the one compartment model (Alsop et al., 2015; Buxton et al., 1998) with the following model parameters: Arterial blood T1 = 2.1 sec (Rane & Gore, 2013), tissue-to-blood partition coefficient λ = 0.9 mL/g (Alsop et al., 2015), labeling efficiency of 82.1% and 95% background suppression inversion efficiency (Shao et al., 2021). Task induced absolute CBF (ΔCBF) and relative CBF (%CBF) changes were obtained by subtraction and division between task CBF and resting state CBF.
4.4.2. VASO and BOLD data prep-processing
Task induced BOLD signal changes were calculated by: VASO signals were divided by BOLD to minimize BOLD effect and vein contamination (Huber et al., 2014). Relative CBV (%) changes were calculated assuming baseline CBVrest = 0.055 ml/ml (Lin et al., 2008; Lu et al., 2005):
4.4.3. Calculation of CMRO2
Concurrently measuring CBF, CBV and BOLD changes allow us to calculate CMRO2 according to the Davis model (Davis et al., 1998) using M and β estimated from experiment 2 for each participant:
4.4.4. Dynamic signal processing and model fitting for BH hypercapnia (Experiment 2)
Dynamic VASO and BOLD signals were detrended (repeat for each run) to remove low-frequency physiological signal drift (Yan et al., 2009). Detrending was not performed for ASL due to pair-wise subtraction of label and control signals. Time points were identified as outliers if they fell beyond the 15th and 85th percentiles of the BOLD and VASO signals and 25th and 75th percentiles of ASL signals across resting-state and task runs.
To estimate the delay between the onset of BH and vessel dilation, we conducted a 3-min pre-scan using low-resolution BOLD (2.5 mm-iso, TR = 1000 ms). This delay was incorporate into the timing of BH cues during the experiment for each participant (e.g. participants were cued to start BH a few seconds ahead of acquiring BH data). Due to the inherently slow process of BH induced vessel dilation, a mixed physiological state of BH and resting conditions could potentially present throughout our measurements. This overlap might result in elevated resting state signals due to residual hypercapnia effects. To address this issue, a least-square curve fitting algorithm (MATLAB lsqcurvefit function) was employed to fit M, β as well as baseline CBF, VASO and BOLD simultaneously from all time points including all rest and BH scans (5 unknowns and 54 time points excluding outliers).
4.5. Segmentation of cortical layers
ASL images were upscaled to a finer grid of 0.25×0.25×0.5 mm3 resolution to avoid singularities at the edges in angular voxel space using the AFNI ‘3dresample’ program with linear interpolation (Cox, 1996). The borderlines of CSF/GM and GM/WM in M1 hand-knob area were manually drawn on ASL control images. Fifteen cortical layers of M1 were segmented by LAYNII (Huber et al., 2021) using the equi-volume layering approach (Waehnert et al., 2014). Segmentation of cortical layers in V1 was done on reconstructed surface. A boundary-based algorithm was used to co-register structural MP2RAGE UNI images to ASL control images (Saad et al., 2009). MP2RAGE volume was skull removed and segmented into WM, GM and CSF using FreeSurfer (7.1.1). AFNI/SUMA and custom python codes were used to generate the equi-volume surfaces (https://github.com/herrlich10/mripy) between WM and pial surfaces. Eight cortical layers (1 WM, 6 GM and 1 CSF) in V1 were then projected back to volume space for analysis using the AFNI ‘3dSurf2Vol’ program with mapping function ‘max’ (https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dSurf2Vol.html).
4.6. Statistical analysis
Paired t-tests were performed to test the significant increase of BOLD, CBF, CBV and CMRO2 in the visual stimulus region (4 – 6 degree of eccentricity) compared adjacent eccentricities. Two-peak patterns were observed in FT activation in M1 (CBF and CBV) and fovea inhibition in V1 (CBF and CMRO2). To test the significance of the two-peak pattern, we calculated R2diff score for each laminar profile (Shao et al., 2021). R2diff was defined as the difference between the adjust R2 estimated from curve fitting assuming two or single Gaussian distributions. Larger R2diff indicates the profile can be better described as two peaks instead of a single peak. Statistical significance of the R2diff was estimated by 5000 Boot strapping runs with random noise in each cortical layer estimated by inter-participant variation. P value was calculated as the probability that R2diff score of profile can be explained by noise only. Linear regressions were performed to assess the correlations between CBF, CBV, BOLD and CMRO2 signals of the visual stimuli experiment. P < 0.05 was considered as significant for both t-test and two-peak detection.
4.7. pRF calculation
First, we sampled the visual field into a 21*21 grid (both the x-axis and y-axis are from - 10 degrees to 10 degrees, with a step of 0.5 degrees). As shown in Supplemental Fig. S3, we obtained the time series of the visual stimulus sweep of each node in the grid. Then we selected a voxel to extract its time series. The time series of this voxel was deconvolved and combined with the sequences of all nodes in the grid to calculate the response of the neurons in this voxel to the visual stimulation at each node. In this way we could get the receptive field map of the neurons in this voxel, then find the center position of the receptive field and calculate its position in the visual field (eccentricity and angle). Finally, we obtained the centrifugation corresponding to this voxel, and calculated the eccentricity map by repeating this step for all voxels.
Data Availability
The data presented in this manuscript can be requested by contacting the corresponding author. The pulse sequence will be available after we establish Material Transfer Agreement (MTA) between user's institute and University of Southern California. MATLAB code for two-peak pattern significance test can be downloaded from our lab website (http://loft-lab.org/index-5.html).
6. Data and code availability
The data presented in this manuscript can be requested by contacting the corresponding author. The pulse sequence will be available after we establish Material Transfer Agreement (MTA) between user’s institute and University of Southern California. MATLAB code for two-peak pattern significance test can be downloaded from our lab website (http://loft-lab.org/index-5.html).
7. Supplemental information
5 Acknowledgement
This work was supported by National Institute of Health (NIH) grant UF1-NS100614, S10-OD025312, R01-NS114382, R01-EB032169, RF1AG084072, R01-EB028297 and R01-NS121040.