Sub-millimetre resolution laminar fMRI using Arterial Spin Labelling in humans at 7 T

Laminar fMRI at ultra-high magnetic field strength is typically carried out using the Blood Oxygenation Level-Dependent (BOLD) contrast. Despite its unrivalled sensitivity to detecting activation, the BOLD contrast is limited in its spatial specificity due to signals stemming from intra-cortical ascending and pial veins. Alternatively, regional changes in perfusion (i.e., cerebral blood flow through tissue) are colocalised to neuronal activation, which can be non-invasively measured using Arterial Spin Labelling (ASL) MRI. In addition, ASL provides a quantitative marker of neuronal activation in terms of perfusion signal, which is simultaneously acquired along with the BOLD signal. However, ASL for laminar imaging is challenging due to the lower SNR of the perfusion signal and higher RF power deposition i.e., specific absorption rate (SAR) of ASL sequences. In the present study, we present for the first time in humans, isotropic sub-millimetre spatial resolution functional perfusion images using Flow-sensitive Alternating Inversion Recovery (FAIR) ASL with a 3D-EPI readout at 7 T. We show that robust statistical activation maps can be obtained with perfusion-weighting in a single session. We observed the characteristic BOLD amplitude increase towards the superficial laminae, and, in apparent discrepancy, the relative perfusion profile shows a decrease of the amplitude and the absolute perfusion profile a much smaller increase towards the cortical surface. Considering the draining vein effect on the BOLD signal using model-based spatial “convolution”, we show that the empirically measured perfusion and BOLD profiles are, in fact, consistent with each other. This study demonstrates that laminar perfusion fMRI in humans is feasible at 7 T and that caution must be exercised when interpreting BOLD signal laminar profiles as direct representation of the cortical distribution of neuronal activity.


Introduction
Neuronal activity in the brain is associated with an increased metabolic demand accompanied by changes in haemodynamics such as blood oxygenation, flow and volume (for reviews see [1][2][3][4]). Functional magnetic resonance imaging (fMRI) is a technique that can non-invasively measure these changes and allows inferring the spatial pattern of neuronal activity while performing a task or at rest. Improvements in MRI technology over the past decades, such as higher magnetic field strengths, novel sequences, optimised pulse designs, and parallel imaging, have pushed the spatial and temporal limits to an extent wherein MRI at ultra-high magnetic field (UHF, �7 T) can routinely achieve sub-millimetre spatial resolution voxels in humans, for both structural and functional imaging (see Special Issues [5,6] and reviews therein). While fMRI investigations have yielded robust, reproducible functional parcellation [7] of different brain areas consistent with previous ex vivo cyto-and myelo-architectural studies [8,9], the advantages of UHF fMRI have enabled neuroscientists to investigate the mesoscopic circuitry within regions across cortical depths and, to a lesser extent, columns in humans (see Special Issue [10] and reviews therein). A vast majority of standard-resolution and laminar fMRI studies have been performed using the Blood Oxygenation Level-Dependent (BOLD) contrast [11,12]. While the BOLD contrast excels in its sensitivity to detect signal changes due to its high signal-to-noise (SNR), it is inherently limited in its spatial specificity relative to site of neuronal activation because of strong signal bias introduced via the intra-cortical ascending veins [13] and by the non-local signal spread (drainage effect) through pial veins [14,15]. Studies investigating the specificity of the laminar BOLD response in humans and animals [16][17][18][19][20] have consistently observed the largest signal change in the BOLD signal at the superficial layers and pial surface despite the fact that the peak of the neuronal activity is expected in the input layers (layer IV in human V1) for feed-forward stimuli [21,22]. Some earlier studies have investigated the leakage of the signal between laminae during steady-state [22][23][24]. Recently, a fully dynamical model of the laminar BOLD signal has been developed [13] that enables model-driven 'deconvolution' (i.e. removal of the intra-cortical ascending venous signal) of the measured BOLD signal profiles to unravel the underlying neuronally-driven signal profiles. However, theoretical assumptions of these model-driven approaches have not yet been subjected to experimental validation. The versatility of MRI provides the means to also measure other (non-BOLD) haemodynamic response parameters such as cerebral blood volume (CBV) using vascular space occupancy (VASO) [25][26][27] or cerebral blood flow (CBF) through tissue (perfusion) using arterial spin labelling (ASL) [28][29][30]. Most studies using these non-BOLD approaches have been carried out in animal models [1,31,32] and have only been applied to high-resolution human studies with the advent of UHF fMRI [33][34][35]. From the perspective of laminar fMRI, animal studies have shown that perfusionweighting is a highly desirable contrast, even more so than total CBV, due to its spatial proximity to neuronal activation [18,36]. While CBV-weighted imaging using VASO has seen a resurgence for laminar fMRI applications [37], perfusion-weighted fMRI using ASL has been mostly limited to relatively low spatial resolution (� 2-4 mm) studies [38] (but see [39]). Achieving higher spatial resolutions, let alone sub-millimetre resolution, with perfusionweighting and adequate brain coverage is challenging. This is due to the relatively lower SNR of the perfusion-weighted signal owing to the lower microvascular density relative to tissue volume (� 1-2%) and T 1 recovery of the labelled arterial water signal, and the higher RF power deposition of ASL sequences in general. The SNR limitation can be addressed to some extent by moving to UHF. The gain in SNR due to increased field strength [40] and the prolonged longitudinal relaxation times (T 1 ) [41,42] allows longer post-labelling delays, thereby, improving the perfusion SNR [34]. Recent developments using ASL at 7 T [34, 35, 43-45] have enabled pushing the spatial resolution for perfusion-mapping to the sub-millimetre regime [46,47] by overcoming several technical challenges; i.e. optimisation of sequence and pulse design [33,35,45,48], using dielectric pads [49,50] in order to improve the labelling efficiency [34], and the utilisation of a 3D-EPI readout [51]. Taking together the advantages of UHF, the spatial specificity of the perfusion signal and the fact that ASL acquires both BOLD and perfusion-weighted images simultaneously, makes ASL a very attractive tool for laminar fMRI. In the present study, we build on our previous work [46,47] to acquire, for the first time, submillimetre resolution simultaneous BOLD and perfusion-weighted fMRI of the human visual cortex at 7 T. We demonstrate that robust, participant-specific, single-session, high-resolution perfusion activation maps can be obtained for laminar fMRI in humans at 7 T. We probe the cortical depth-dependence of BOLD and perfusion-weighted signals in response to visual stimulation in humans and reconcile our experimental findings using the recently proposed dynamic model of the laminar BOLD signal.

Materials and methods
Seven healthy volunteers (median age = 28 years) participated in the study following screening and having given written informed consent. The study was approved by the Ethics Review Committee for Psychology and Neuroscience (ERCPN) at Maastricht University and all procedures followed the principles expressed in the Declaration of Helsinki.

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 isocentre 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 0 . In addition, two 18x18x0.5 cm 3 high-permittivity dielectric pads containing a 2.8:1 solution of calcium titanate (CaTiO 3 ) and heavy water (D 2 O) by weight [52] were placed on either side of the head at the level of the participant's temporal lobes to increase B 1 + (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.
Stimulus paradigm. Full contrast black-and-white radial flickering checkerboard was presented using PsychoPy v1.90.0 [53] for 20 s (stimulus on) followed by 40 s of an iso-luminant grey background (stimulus off). Each functional run was � 11 min long consisting of a 30 s initial baseline period and ten stimulus on-off blocks. The participants were instructed to remain motionless and fixate on a central fixation dot throughout each of the four functional runs.
Functional MRI. Functional data were acquired at a nominal resolution of 0.9 mm isotropic using a pulsed ASL (PASL) sequence [

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 1 -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 1 -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 1 map was supplied as an additional input (proxy T 2weighted 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 1 -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 1 -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 intracortical 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][81][82] and the eventrelated 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 prestimulus baseline) per lamina and then normalising the signal with the mean of the EPI (to account for transmit-receive biases and baseline T 2 � 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-tonoise 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.

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 0 ) was set to 2 mL, of which 50% relates to microvasculature (ω v =0.5) and

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 1 -weighted anatomical image (c).
https://doi.org/10.1371/journal.pone.0250504.g001

PLOS ONE
50% to the ascending veins (ω d = 0.5) [84]. The baseline CBV (CBV 0 ) 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 0 ) obtained from the ASL data was not in physiological units of mL/min/100 g, we rescaled the CBF 0 by a constant scaling factor x=5520 calculated using the following assumptions. First, that the CBV 0 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 0 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 2 metabolism, depth-dependent changes in microvasculature CMRO 2 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 2 is extracted on arrival at the venules (E 0 =0.35) with no further extraction along these vessels (i.e., O 2 saturation of haemoglobin (1-E 0 ) 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. 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 1 -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 timecourse 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 zigzag 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 2 V1: 0.986, R 2 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 2 V1: 0.939, R 2 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 2 V1: 0.537, R 2 V2: 0.745). In the perfusion signals, slight oscillatory behaviour is observed during the post-stimulus period.

Discussion
The present study, is the first demonstration in humans of the improved spatial specificity of the perfusion signal compared to the BOLD signal using isotropic sub-millimetre spatial resolution ASL at 7 T. We found incongruent cortical depth profiles between the BOLD signal and perfusion changes, which, however, turned out to be physiologically consistent with each other after employing a dynamical BOLD signal model.

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

PLOS ONE
(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 perfusionweighted 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 1 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 3 ). 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 2 � effects. While minimal processing approaches for fMRI have been proposed [113], except for a few studies (such as [37,73,99,[114][115][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 costfunction 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, largecoverage 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 perfusionweighting, 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 realtime 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 ultrahigh field is possible using the standard commercial head-coil with single-channel transmit (NOVA Medical, USA). However, the B 1 + 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][122][123][124] to potentially further optimise highresolution ASL fMRI at ultra-high fields. Having demonstrated the feasibility of perfusionweighted 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.
Supporting information S1 Fig. Point spread function of the 3D-EPI acquisition along the principal phase encoding and readout directions. Partial-Fourier [117] was employed along the principal phase encoding direction. Please note, this simulation used zero-filling with partial-Fourier representing the default image reconstruction scenario. However, in the present study, a POCS reconstruction [59] with 8 iterations was carried out which minimises partial-Fourier blurring, and consequently improves the PSF Tissue classes C3-C5 were thresholded on a subject-by-subject basis to include as much of the sagittal and transverse sinuses as possible. In most of our subjects, this procedure did not require any manual correction of the combined mask. ITK-SNAP v3.6 [72] was used to make any manual corrections when required. The Freesurfer T 1 -weighted input is presented as a 3D render showing an intact GM surface at the occipital lobe with the green arrows indicating locations of the now automatically stripped sagittal and transverse sinuses.