Data acquisition

Data were acquired on a whole-body Siemens Magnetom 7 T research scanner with a gradient system capable of maximum gradient amplitude of 70 mT/m and maximum slew rate of 200 T/m/s (Siemens Healthineers, Erlangen, Germany) and a 32-channel receive phased array head coil (Nova Medical, USA). The participant placement and preparatory procedure followed the protocol previously described in [34, 43]. In short, the eye centres were taken as iso-centre reference (instead of the eyebrows, as is typically done) and supplementary cushions were provided to the participants under the neck, to ensure that the large feeding arteries to the brain were parallel to the B₀. In addition, two 18x18x0.5 cm³ high-permittivity dielectric pads containing a 2.8:1 solution of calcium titanate (CaTiO₃) and heavy water (D₂O) by weight [52] were placed on either side of the head at the level of the participant’s temporal lobes to increase B₁+ (therefore, labelling) efficiency at 7 T [50]. This was done because our earlier work [34, 43] found that the relative orientation of large feeding arteries together with the use of dielectric pads gave the best perfusion results.

Data processing

The anatomical data were pre-processed in SPM12 r7487 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12) [60, 61] and FSL v6.0 (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) [62, 63] using the workflow illustrated in S4 Fig. This anatomical pre-processing workflow was developed particularly to work well for MP2RAGE data (https://github.com/srikash/presurfer). First, the second inversion image of the MP2RAGE was subjected to the automated segmentation in SPM12 [64]. The bias-corrected second inversion image was used to create a whole-brain mask using FSL BET [65]. The thresholded non-brain tissue classes from the SPM12 segmentation were summed together to create a mask of the non-brain tissue and large sinuses. The non-brain mask was manually curated in cases, in which the automatic masks were sub-optimal. The T₁-weighted MP2RAGE image (UNI) was bias-corrected using SPM12 and was stripped off the non-brain tissue and large sinuses using the mask obtained from the second inversion image. This pre-processed T₁-weighted MP2RAGE was supplied as input to the high-resolution recon-all pipeline of Freesurfer v6.0 (https://surfer.nmr.mgh.harvard.edu) [66]. Additionally, the MP2RAGE T₁ map was supplied as an additional input (proxy T₂-weighted image) for pial surface optimisation. The segmentation and surface construction were done in the native resolution and the segmentation quality in the occipital lobe was manually curated (S5 Fig). A probabilistic retinotopic atlas was applied to the Freesurfer reconstructed data using neuropythy (https://github.com/noahbenson/neuropythy) [67] to obtain participant-specific V1 and V2 regions-of-interest (ROIs) (S3b Fig). Following the automatic segmentation and reconstruction, the WM surface was extended into WM by 30% of the cortical thickness to account for any discrepancy of the GM-WM boundary when using T₁-weighted MP2RAGE images [68]. The first inversion image of the MP2RAGE was used to check the extended WM boundaries due to its sharp WM-GM contrast. We also extended the pial boundary by the same amount into the CSF to sample the signal away from the pial boundary. Then, we generated a total of twenty-one intermediate equi-volume surfaces within the GM using Surface tools (https://github.com/kwagstyl/surface_tools) [69] (S3c Fig). The functional datasets were pre-processed using Advanced Normalization Tools (ANTs) v2.3.1 (https://github.com/ANTsX/ANTs) [70, 71]. First, the functional runs were subjected to affine realignment. Next, the temporal mean of the functional run and the temporal mean of the opposite phase-encoded run were used to calculate an undistorted template image and the distortion-correction warps were saved. Lastly, a transformation matrix was calculated for each functional run to the T₁-weighted data using the visual alignment tools in ITK-SNAP v3.6 [72] and a final rigid alignment using ANTs. All transforms were concatenated and applied to the unprocessed functional datasets in a single resampling step using a 4th degree B-spline interpolation. This minimises resolution losses due to multiple interpolation steps while providing the required quality of registration accuracy in laminar fMRI studies [73, 74]. Statistical analyses of the functional data were carried out using the ‘Full Perfusion Signal Modelling’ pipeline [75] in FSL FEAT [76, 77]. Here, we modelled three regressors i.e., the stimulus design convolved with the canonical haemodynamic response function (HRF) representing the BOLD signal, the alternating label-control acquisition of the ASL sequence representing the baseline perfusion-weighting and the combination of these two regressors representing the perfusion activation [78, 79]. Due to the disparity in the spatial spreads of the BOLD and perfusion activation (Fig 2 left panel), a mask of the overlap between the BOLD and perfusion activation cluster thresholded masks from FEAT was created. This ensured that we sampled the BOLD and perfusion signals from the same voxels. Laminar analyses were carried out in Freesurfer by sampling the functional time-series signal from the ROIs using nearest-neighbour interpolation. No surface or intra-cortical smoothing was applied. The laminar time-courses sampled from V1 and V2 across all participants were imported into MATLAB R2016b (MathWorks, USA) for the time-series analyses. The BOLD and perfusion-weighted time-courses were obtained for each lamina by applying surround-averaging and surround-subtraction, respectively [80–82] and the event-related average time-courses were calculated. The event-related average BOLD time-course was subsequently rescaled to percent BOLD signal change relative to the pre-stimulus baseline (≈10 s). The analysis of the perfusion time-series followed several steps: First, the perfusion-weighted time-series is a measure of the modulation depth (i.e., the magnitude of the zig-zag) of the raw ASL time-course in MRI signal units (S6 Fig). It is important to note that these data are not scaled in physiological units and is representative of the perfusion SNR of the data. We then derived the following measures from perfusion-weighted time-course: absolute and relative perfusion change, and baseline perfusion. Absolute perfusion change was calculated by taking the change in the perfusion activation (i.e., by subtracting the pre-stimulus baseline) per lamina and then normalising the signal with the mean of the EPI (to account for transmit-receive biases and baseline T₂* effects). The absolute perfusion change, thus obtained, is in arbitrary units but proportional to the quantitative perfusion change. The absolute perfusion change can then be rescaled into physiological units, as typically done in perfusion quantification studies [38, 83]. Relative perfusion change is the percentage change in the perfusion signal due to activation per depth relative to its respective baseline. Note that the relative perfusion change does not need to be divided by the mean EPI image for scaling (as it appears both in the numerator and the denominator and thus cancels out). The baseline perfusion (Fig 1a) was calculated using simple subtraction of the label-control time-points during the baseline period (≈ 0-30 s at the beginning of the run) and pre-stimulus intervals (≈ 0-10 s before stimulus onset) of the stimulus blocks. The temporal signal-to-noise of the perfusion-weighted timeseries (perfusion tSNR) image (Fig 1b) was calculated as the ratio of the temporal mean to the temporal standard deviation of the perfusion-weighted timeseries. Laminar steady-state profiles of the BOLD signal, absolute and relative perfusion change signals were calculated by averaging the respective signals within the ≈ 14-28 s interval following stimulus onset. The baseline perfusion laminar profile (S7 Fig) was obtained by averaging within the entire ROI.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. The average baseline perfusion (a) and perfusion tSNR (b) maps from a superior (left) and inferior (right) slice of an example participant is shown overlaid on the corresponding T₁-weighted anatomical image (c).

https://doi.org/10.1371/journal.pone.0250504.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

(a, left-right) BOLD signal activation map of an example participant, averaged over all runs, overlaid on the inflated left hemisphere reconstruction from Freesurfer. Contours of the V1 (blue) and V2 (green) labels obtained from Neuropythy are also overlaid on the inflated surface. Cropped orthogonal views of the participant’s occipital lobe with the BOLD signal activation map overlaid in voxel space. Boundaries of the different laminar surfaces are also overlaid, colour-coded from cyan (pial)-to-magenta (white). The purple dotted square inset indicates the zoomed-out views presented below. Event-related average ASL time-course of highly BOLD activated voxels across runs for this participant is shown to the right. (b, left-right) The same are presented for perfusion activation in green. The error-bars indicate SEM across trials.

https://doi.org/10.1371/journal.pone.0250504.g002

Simulating the laminar BOLD signal from the measured perfusion profile

The experimentally measured laminar BOLD response profiles in V1 and V2 regions were compared to theoretical predictions of the dynamical laminar BOLD signal model of Havlicek and Uludag [13] using their publicly available MATLAB code (https://github.com/martinhavlicek/laminar_BOLD_model). In order to do this, some assumptions regarding baseline physiological parameters were made. That is, the laminar BOLD signal profile is fully determined given physiological values and measured CBF changes. These assumptions are within typical values of human vascular physiology (see [4] for a review, and references therein). Here, the measured absolute perfusion laminar profiles (both baseline and activation) were used to input the physiological parameters of the model. The total amount of venous baseline CBV (CBV₀) was set to 2 mL, of which 50% relates to microvasculature (ω_v=0.5) and 50% to the ascending veins (ω_d = 0.5) [84]. The baseline CBV (CBV₀) distribution was set to be constant across laminae in the microvasculature but increased linearly towards the surface in the ascending veins (slope, s_d = 0.4). Since the baseline perfusion (CBF₀) obtained from the ASL data was not in physiological units of mL/min/100 g, we rescaled the CBF₀ by a constant scaling factor x=5520 calculated using the following assumptions. First, that the CBV₀ in microvasculature (1 mL /100 g) is divided uniformly between all laminae and the mean transit time (t_v0) through microvasculature (averaged across all laminae) is 1 s (for complete list of other model parameters, see Table 2 in [13]), giving an average CBF₀ of 1 mL/s/100 g or 60 mL/min/100 g. This scaling factor allows a rescaling of the CBF from arbitrary MRI units to physiologically meaningful units for the purpose of reporting. However, the laminar BOLD model itself is driven by depth-dependent changes in relative CBF (in %), therefore, the exactness of the scaling in physiological units is not critical to the simulations. Next, by assuming a linear coupling between CBF and O₂ metabolism, depth-dependent changes in microvasculature CMRO₂ were obtained (n = 3) [85]. We also assume relative CBV changes in microvasculature and ascending veins are related to the relative CBF via Grubb exponents α_v = 0.35 and α_d = 0.1, respectively [86, 87]. Finally, we assume that 65% of the O₂ is extracted on arrival at the venules (E₀=0.35) with no further extraction along these vessels (i.e., O₂ saturation of haemoglobin (1-E₀) is assumed to be homogeneously distributed across laminae). All other parameters were defined as in the default scenario described in [13]. Please note that we did not fit the model to data but used experimentally obtained perfusion-weighted signal data and plausible biophysical parameters to generate a prediction of the laminar BOLD signal profile.

Baseline perfusion

Fig 1a and 1b show two representative slices (one superior, one inferior) of the average baseline perfusion map and the perfusion temporal signal-to-noise (tSNR) of a participant overlaid on the T₁-weighted anatomical image (Fig 1c). These maps show that the average perfusion signal is highly localised to the GM ribbon and demonstrates the quality of the co-registration between the acquisition slab with the anatomy as indicated by the absence of signal shifted into the ventricles and the clearly defined sulci (wherever resolvable). The perfusion-weighted data shown in Fig 1a is in arbitrary MRI signal units.

Functional activation

Robust statistical activation was obtained for all participants for both the BOLD (Fig 2a) and perfusion signals (Fig 2b). The BOLD activation envelopes a much larger swath of cortex than perfusion activation does (Fig 2, left panel). This is expected given the differences in the detection sensitivity (i.e. functional contrast-to-noise (fCNR)) between the BOLD and perfusion signals, and the presence of BOLD signal in pial veins. In addition, the BOLD activation obtained follows the characteristic localisation pattern observed with standard GE-EPI studies. That is, the largest BOLD activated voxels are always localised at the CSF-GM boundary (Fig 2a, purple zoomed-out boxes). In contrast, the perfusion activation was observed to be more spatially localised to the GM ribbon with the highest activated voxels localised mid- to deep-GM (Fig 2b, green zoomed-out boxes). The activation maps are shown in the three orthogonal views to highlight the consistency of the GM localisation of activation in 3D.

Finally, the ASL time-courses exhibit a zig-zag modulation that is characteristic of ASL sequences (due to the acquisition of alternating label and control volumes) demonstrating the high quality of the data. The modulation depth of this zig-zag represents the amount of labelled spins delivered to the tissue and is, therefore, proportional to tissue perfusion. The ASL time-course obtained from the highest BOLD signal activated voxels in grey matter shows the typically observed increase in the BOLD signal magnitude during activation with weaker zig-zag modulation (Fig 2a, right panel, purple). On the other hand, the ASL time-course obtained from the highest perfusion-activated voxels shows the strong zig-zag modulation throughout but with lower BOLD signal modulation (Fig 2b, right panel, green). Please note, that the zig-zag time-courses represent the event-related average response over all depths of grey matter from the ASL timeseries. All three key differences between the BOLD and perfusion activation signals were consistently observed in all the participants.

Laminar analysis

The group-average laminar time-courses of BOLD signal change, absolute perfusion change, and relative perfusion change are shown in Fig 3 for V1 and in Fig 4 for V2. The temporal behaviour of the three sampled signals across all laminae is presented as a heat-map in the top row with time along the X-axis, the cortical depth along the Y-axis, and the magnitude of the signal in colour code. We observed inter-regional differences with laminar responses of all three signals, with V2 having a lower amplitude than V1. The laminar profiles of the BOLD signal change exhibit positive slopes (Slope V1: 4.88 ± 0.129, Slope V2: 4.81 ± 0.195) with a strong linear trend (R² V1: 0.986, R² V2: 0.967). The laminar profiles of the relative perfusion change, on the other hand, exhibit negative slopes (Slope V1: -4.91 ± 0.27, Slope V2: -4.64 ± 0.111) with a strong linear trend (R² V1: 0.939, R² V2: 0.988). Interestingly, the absolute perfusion changes exhibit a moderately positive slope (Slope V1: 3.35 ± 0.68, Slope V2: 3.80 ± 0.485) albeit without a strong linear trend (R² V1: 0.537, R² V2: 0.745). In the perfusion signals, slight oscillatory behaviour is observed during the post-stimulus period.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Laminar BOLD and perfusion signal changes in human V1 ROI obtained using sub-millimetre 3D-EPI ASL at 7 T.

Top row: Heat-map representations of the group-average BOLD signal change, absolute perfusion change, and relative perfusion change with cortical depth along Y-axis and time along the X-axis. Middle row: Five out of the twenty-three total laminar time-courses for the respective sampled signals. Bottom row: Laminar profiles of the positive response during stimulation for the respective sampled signals with cortical depth along X-axis. Error-bars indicate SEM across participants. The black bar in the middle row indicates the stimulus duration.

https://doi.org/10.1371/journal.pone.0250504.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Same as Fig 3 for V2 ROI.

https://doi.org/10.1371/journal.pone.0250504.g004

Simulations of the laminar BOLD signal

Fig 5 shows the simulated laminar BOLD signal profile (solid blue lines) and the experimentally measured laminar BOLD signal profiles (dotted purple lines, see Figs 3 and 4). The measured and simulated profiles were highly congruent (Pearson’s correlation: r = 0.9984 for V1, r = 0.9977 for V2), demonstrating that, despite the discrepancy of the relative and absolute perfusion, and the BOLD signal profiles, they are in fact compatible with each other.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Laminar profiles of the BOLD signal for V1 and V2 ROIs.

The measured responses are the same as the BOLD signal profiles in Figs 3 and 4. The simulated profiles are obtained using the laminar BOLD model proposed in [13] (S7 Fig).

https://doi.org/10.1371/journal.pone.0250504.g005

Functional BOLD and perfusion activation

We obtained robust participant-specific, single-session activation maps for simultaneously acquired isotropic sub-millimetre spatial resolution BOLD and perfusion signals at 7 T. We observed a larger spread of activation for the BOLD signal (Fig 2a) compared to the perfusion signal (Fig 2b). This is expected because the detection sensitivity of the perfusion signal (fCNR) is much lower than that of the BOLD signal [34, 38, 39]. Additionally, this can also be explained by the higher spatial specificity of the perfusion signal compared to the BOLD signal, which is susceptible to non-local signal spread due to downstream venous bias away from the actual site of activation [15, 88]. This is also observed in high-resolution fMRI with the highest BOLD activated voxels located at the CSF-GM boundary (Fig 2a). On the other hand, the perfusion activation map exhibits a well-defined localisation to the cortical ribbon (Fig 2b), mostly located in cortical GM [18]. Importantly, given that perfusion signal has much lower fCNR than the BOLD signal in standard resolution studies (2-4 mm in each direction), it was not necessarily expected that ASL will have enough sensitivity at sub-millimetre resolution for detecting perfusion activation. One reason that with increasing resolution there is enough perfusion fCNR is that, not only image SNR, but also partial voluming with CSF and WM is decreased, i.e., thermal, and physiological noise coming from outside GM are reduced. This is different for the BOLD signal as pial vessels located in CSF (see Fig 2a) do contribute to the overall BOLD signal in low resolution studies and therefore, increases in spatial resolution decrease both image SNR and overall signal contribution. That is, going from low- to high-spatial resolution penalizes CNR of the BOLD signal more than of the perfusion signal. Recently, a novel fMRI approach called VAPER [89] has also been put forward as a contrast useful for perfusion-weighted high-resolution fMRI by mixing VASO and perfusion contrasts. Although the combination of two contrasts boosts VAPER’s sensitivity, it markedly complicates its ability to quantify perfusion and also, its physiological specificity. Thus, established ASL techniques remain the most feasible way to acquire in vivo perfusion-weighted images that can be straightforwardly validated using quantitative fMRI models, and can be expected to provide reproducible results across a wide range of sequence parameters and field strengths [34]. While the present study employs a PASL acquisition scheme, recent studies have explored the feasibility of using a pCASL acquisition scheme with a 3D-GRASE readout for perfusion-weighted laminar fMRI at 7T in the motor cortex [90, 91]. While the PASL laminar profile seems consistent with that of the present study, the pCASL profile shows a relatively higher baseline perfusion in the middle layers. The measured profiles are the outcome of a complex interplay between one or more of several factors such as readout approach, background suppression scheme, partial voluming, post-stimulus inhibitions in certain layers, residual pial arterial contributions despite long transit times [92]. The sources of this apparent discrepancy observed with the labelling schemes and their interplay with other ASL parameters (e.g., TE, spatial resolution, selection of voxels, labelling schemes and timings, background suppression etc.) needs further investigation. At this moment, there is very little in vivo data in terms of a full exploration of the impact of the parameter space that ASL offers for laminar fMRI. Thus, having demonstrated the feasibility of perfusion-weighted laminar fMRI at 7 T, future studies can consider systematic comparisons of the different labelling approaches, readout strategies and how they affect the laminar profiles of perfusion both during baseline and activation at UHF. A recent review used an exemplary dataset from the current study to illustrate the potential of perfusion-weighted imaging in the context of exciting avenues for non-BOLD laminar fMRI applications at UHF [93]. As highlighted in the review, the perfusion contrast has been highly desired for laminar fMRI as the perfusion signal is relatively unaffected by the venous compartments, both by the pial and ascending veins, and the large arterial compartments. In comparison, the BOLD signal is heavily weighted towards the venous compartments and the VASO signal can have contributions from both arterial and venous in addition to CBV changes in microvasculature [33, 94, 95]. The reason for the high perfusion localisation specificity is that the tagged arterial water is mostly exchanged with the tissue at the level of the capillaries. In addition, the transit delay for the labelled blood to arrive at the region-of-interest (in this case, occipital lobe) can be ≈1-1.3 s [96]. Together with the blood transit time within tissue on the order of ≈1-2.5 s [97], only little longitudinal magnetisation of the tag remains (due to T₁ decay), i.e. that almost no magnetisation of the label is present in venous blood (see S8 Fig), except for artefacts caused by labelling of venous blood superior to the imaging slab in some ASL schemes (see [83]). The transit time for the acquisition in the present study was optimised for the visual cortex [92], and is reflected in the inter-regional differences in the baseline perfusion signal and its temporal stability of the tissue (Fig 1). The absence of the venous bias and the signal being dominated by the capillary compartment implies that the perfusion contrast more closely follows both the spatial profile and the amplitude of cortical metabolism and neuronal activation. Another important aspect of ASL acquisitions is the possibility to obtain a quantitative estimate of the baseline signal across depths. The difference between the highly BOLD-activated or highly perfusion-activated voxels is readily visible in the ASL time-courses (Fig 2). The time-courses for perfusion activation show reduced amplitude of the signal envelope and larger difference between pairs of data points (i.e., the zig-zag modulation) indicating that these voxels contain signals from mostly the microvasculature and that observed responses are indeed capturing the changes in perfusion. In contrast, there are small zig-zag changes relative to the overall signal envelope in the time-courses for the highest BOLD activation, reflecting a smaller contribution from microvasculature. This means that the spatial non-overlap that we observe between the perfusion and BOLD signals is driven largely by differences in the underlying physiology and not the differences in SNR.

Laminar BOLD and perfusion responses

We replicated previous findings [21, 73, 98] that the event-related average BOLD signal amplitude (Figs 3 and 4, first column) increases towards the CSF-GM boundary (e.g., [99, 100]). The BOLD signal increase to the superficial layers is well understood and can be attributed to two signal biases: (a) increase in baseline CBV of the intra-cortical ascending veins, and (b) the non-local blooming effect from the pial veins ([101], and for overview see, [4]). The presence of these biases in the BOLD signal makes the interpretation of the measured laminar signal profile, specifically in the superficial layers, challenging [102]. One approach to deal with the issue of spatial bias in GE-BOLD signal is model-driven spatial ‘deconvolution’ [22, 24], which, however, has not yet been validated with another (simultaneously) acquired fMRI modality. The profile of the relative perfusion change (Figs 3 and 4, right column) exhibits the opposite behaviour (compared to the BOLD profile) [103] with the magnitude of the signal increasing towards the GM-WM boundary with a strong linear trend. Furthermore, QUIPSS II pulses were employed in the present study allowing clear-cut definition of the tagged bolus. This means that the observed patterns of laminar signal behaviour are unlikely to be due to undelivered tagged blood in the diving arterioles. Although the impact of the QUIPSS II pulse depends on the chosen parameters and the arrival times to the regions-of-interest, an increase in blood flow upon activation can result in a more complete delivery of the tagged spins to the tissue, including the deeper layers at the time of volume acquisition. This could yield a larger fractional perfusion change in the deeper layers relative to the baseline condition. While it is usually argued that for feed-forward stimuli the peak in activation must be in the middle layers, electrophysiological evidence, histology, and a previous BOLD signal study after spatial ‘deconvolution’ [22] support the view that V1 also receives increased input signal into layer VI in addition to layer IV. Please note that despite the high spatial resolution used in this study, we could not detect a fine-grained distinction between laminae. The perfusion spatial profile obtained, thus, represents a smoothed version of the underlying neuronal activity. For example, data shown in Fig 4b and 4d in [22] and in [104] (see Fig 9 in [22]) are compatible with the spatial profiles found in the current study. We find that the relative increase in the perfusion signal in the middle to deeper layers is also consistent with animal literature (see also [18, 105, 106]). The absolute perfusion signal change profile (Figs 3 and 4, middle column) exhibits a weak positive slope and non-linear behaviour across depths. However, both relative and absolute signal changes are derived from the same perfusion-weighted signal obtained after surround-subtraction and the difference stems from the spatial profile of the baseline perfusion (S5 Fig). Please note, that the increase of the absolute perfusion signal from WMB to CSFB (by ≈ 30-50%) is much smaller than that of the BOLD signal (by ≈ 100-120%). Additionally, in contrast to the BOLD signal, the absolute perfusion change drops beyond the CSF border. It is important to reiterate that the perfusion and BOLD signals were extracted from the exact same patch of cortex using a mask of the spatial overlap between the BOLD and perfusion activations. This was done due to the larger significantly activated volume using the BOLD compared to the perfusion contrast. By adopting this approach, we show that despite having the same underlying cortical architecture, the relative and absolute perfusion signal changes differ in their depth-dependent behaviour and both differ from the BOLD signal either in the sign of their slope or the relative increase of the profile towards the surface. Here, the BOLD signal profile obtained from ASL is consistent with profiles from previous studies acquiring GE-BOLD signal alone. In order to test if the apparent discrepancy between the relative and absolute perfusion profiles and BOLD profile can be reconciled, we simulated the BOLD signal profile from the measured perfusion profiles using the recent dynamical laminar BOLD signal model (for details, see S7 Fig). We show that the positive slope and the relative increase of the measured BOLD profile can be obtained from the laminar profile of the relative (having negative slope) and absolute (having much smaller increase towards the surface) perfusion signal by modelling the ascending vein bias, i.e., simulating the laminar BOLD response in a forward manner. Therefore, we conclude that despite their seemingly contrasting behaviours, the BOLD and perfusion signal profiles are, in fact, physiologically consistent with each other. Additionally, the BOLD time-courses exhibit a strong post-stimulus undershoot (PSU) consistent with previous studies [73, 107]. Interestingly, our perfusion measures also exhibit PSUs but with smaller amplitudes relative to the positive response. In contrast to the smooth recovery of the PSU to baseline in the BOLD signal, the perfusion PSU exhibits slight oscillatory behaviour. These post-stimulus oscillatory transients are consistent with previous reports of perfusion measurements in humans (e.g., [108]) and with optical imaging in rodents [109]. The oscillatory transients observed in the previous perfusion study in humans [108] could not resolve any depth-dependent modulations owing to its much lower spatial resolution (i.e. 2.65 x 2.65 x 5 mm³). The post-stimulus oscillations in our study near the WM boundary are smoother and evolve with a different oscillatory phase than near the CSF boundary, where the oscillations are more pronounced (Figs 3 and 4, middle panels). While this observation in itself is interesting, pin-pointing the exact vascular physiology that elicits this behaviour is beyond the scope of this study. Taken together, we believe that the current study presents a breakthrough in non-BOLD fMRI research with the development of sub-millimetre resolution perfusion fMRI using ASL and its feasibility for layer-specific investigations, which has hitherto been an uncharted territory in humans.

Data processing

We developed a novel workflow to pre-process anatomical images (S4 Fig) by using the second inversion image of the MP2RAGE and SPM12’s segmentation to automatically mask out the sagittal and transverse sinuses that are crucial for highly accurate pial surface delineation using Freesurfer’s recon-all. In some participants, the workflow required (albeit very little) manual corrections of the segmentation masks (S5 Fig). We used an open-source python package, neuropythy [67], to apply a probabilistic atlas of retinotopy in participant’s native space to generate automatic labels of V1 and V2 (S3b Fig). We qualitatively compared this atlas-based approach on a separate dataset of pRF mapping that was acquired using the same scanner, head coil and similar coverage, and found a high degree of overlap consistent with the findings of [67]. Please note, the focus of the present work is distinguishing the BOLD and perfusion signals, and does not rely on the perfect delineation of V1 and V2 borders. Cortical layering was done using the equi-volume approach [110] using Surface tools [111] as equi-volumetric layering is currently not natively supported in Freesurfer. Even though the layering is done on the whole cortical ribbon, we manually ensured that the delineations were accurate within the V1 and V2 masks in each participant. Nevertheless, for spatial resolutions such as the present study, the exact choices of the layering model [112] and the number of sampled depths [101] does not affect our main conclusions. As ASL measures the BOLD signal simultaneously with perfusion, the BOLD signal profile serves as an internal control and can be compared to other BOLD signal only studies, and the perfusion profile is compared to the BOLD signal profile that is concurrently acquired. The BOLD signal spatial profile for feed-forward stimuli (such as checkerboards used in this study) is well known and the BOLD signal derived from the ASL data reproduces this well-known amplitude increase towards the surface of the cortex (Figs 3 and 4), confirming the accuracy of the data processing. It is important to note that the TE in the present study was 15 ms owing to the EPI readout, and allowed us to minimise BOLD weighting of the acquisition. Furthermore, we obtained our perfusion-weighted timeseries using the surround-subtraction approach which has been demonstrated previously [80, 82] to significantly reduce BOLD contamination in the perfusion-weighted timeseries. We additionally calibrated the perfusion-weighted data using the mean EPI image in order to correct for any residual baseline T₂* effects. While minimal processing approaches for fMRI have been proposed [113], except for a few studies (such as [37, 73, 99, 114–116]), they are unfortunately not commonplace in laminar fMRI. We once again highlight the importance of minimal processing, and demonstrate its feasibility using the ANTs framework to estimate, combine and apply transformations of motion, distortion-correction and co-registration to the anatomical image in a single resampling step, thereby reducing the amount of smoothing resulting from the processing of the data [74]. Please note, in some ASL acquisitions there may be strong differences in image contrast between the label and control images, and the choice of realignment cost-function may impact the quality of correction. However, this was not the case for the present study (S9 Fig). Importantly, as can be seen in Fig 2, high values of perfusion were tightly confined to the GM ribbon illustrating the high accuracy of the segmentation and co-registration in the present study.

Limitations

The goal of the present study was to demonstrate the feasibility of using perfusion-weighted contrast with ASL for laminar fMRI and, to that end, we employed a block design with a strong feed-forward visual stimulus that is known to elicit widespread activation. Due to the lower SNR of the perfusion-signal, we averaged approximately 44 min worth of functional data. While there is the undoubted benefit in spatial specificity, ASL may not be well-suited for all laminar fMRI studies, particularly those with small effect sizes. In addition, GE-BOLD laminar fMRI data are routinely acquired with 0.6-0.8 mm isotropic resolutions, higher than the current ASL study. While the use of partial-Fourier acquisition [117] can reduce the effective spatial resolution along a dimension, the amount of blurring was reduced by using 8 iterations of the POCS reconstruction algorithm [59] instead of the default zero-filling [57]. At 7 T, large-coverage high-resolution single-TR acquisitions are unfeasible without substantial acceleration. Thus, a combination of GRAPPA 4 and phase partial-Fourier was used to shorten the echo time and gain perfusion SNR. The partial-Fourier in the slice direction may cause some SNR loss compared to an image with fully sampled slice dimension. The reason for the partial-Fourier sampling in the slice direction is to ensure optimal TI2 for the k-space centre and an adequate TR. Having the concurrently acquired BOLD signal information as an internal control, we do not expect this to have affected our main findings in the present study. We expect that any of the limitations from the partial-Fourier factors employed are offset by the gain in perfusion SNR (and CNR) due to the shorter TE. The lowest achievable TE in the present study was 15 ms owing to the EPI readout, which is still not ideally suited for perfusion imaging. It would be desirable to achieve even shorter TEs (e.g. ≈3 ms or less) for better perfusion-weighting, however it is currently not possible using conventional Cartesian EPI readouts. To this end, there has been recent progress in non-Cartesian (e.g. spiral readouts) ASL fMRI at ultra-high field [118, 119]. Dual-echo spiral acquisitions can be particularly useful for simultaneous perfusion and BOLD imaging achieving the first echo at ≈2 ms (perfusion-weighted) and the second echo at ≈25 ms (BOLD) at 7 T. However, these non-Cartesian acquisitions are prone to inaccuracies in the spiral trajectories due to gradient imperfections that require real-time monitoring and correction using specialised field-monitoring hardware [118]. However, research and development are still underway to address these technical challenges in non-Cartesian imaging and currently sub-millimetre fMRI acquisitions have not been demonstrated. Nevertheless, to the best of our knowledge, this study remains the highest spatial resolution functional ASL study in humans till date. Going forward, sub-millimetre resolution ASL can be invaluable to studies that are examining BOLD signal physiology, for validating existing models or for brain areas contaminated by close large pial veins. Future development work in both acquisition and reconstruction can push the boundaries of spatial and temporal resolutions of ASL for laminar fMRI at ultra-high fields. We show that high-resolution ASL at ultra-high field is possible using the standard commercial head-coil with single-channel transmit (NOVA Medical, USA). However, the B₁+ inhomogeneities remain a major hurdle [34]. While we were able to mitigate this to some extent using dielectric pads [50, 52], future studies will be able to take advantage of advances in parallel transmission (pTx) technology [120] or the use of dedicated labelling-only RF-coils [121–124] to potentially further optimise high-resolution ASL fMRI at ultra-high fields. Having demonstrated the feasibility of perfusion-weighted laminar fMRI using ASL at a sub-millimetre spatial resolution, future studies will be able to systematically evaluate different properties ASL and its impact on the perfusion signal evolution at ultra-high fields.