Characterizing BOLD activation patterns in the human hippocampus with laminar fMRI

The hippocampus is known to play a vital role in core cognitive functions, including learning and memory. It is subdivided into several subfields, including the dentate gyrus (DG), cornu ammonis (CA) 1-4, and the subiculum. Subfields are arranged sequentially along its transverse axis, forming a characteristic crescent shape. Within each subfield, layers are stacked radially, meaning they are arranged perpendicular to the surface of their respective subfield like the layers of an onion, following its curve. This radial organization of layers is maintained as each subfield curves along both the transverse and longitudinal axes of the hippocampus, creating a continuous three-dimensional structure.

While the human hippocampus has been extensively studied at the macroscale using functional magnetic resonance imaging (fMRI), the underlying microcircuits at the mesoscale (e.g., laminar level) have barely been investigated (Maass et al., 2014). In contrast to the six-layered neocortex, the hippocampus contains three (partially subdivided) layers. Furthermore, animal research has shown that in stark difference to the classical canonical microcircuits in the neocortex, differentiating feedforward and feedback pathways (Douglas et al., 1989), the hippocampal projections are primarily feedforward. Nevertheless, inputs to hippocampal pyramidal cells can be differentiated based on their position along the dendritic tree (a schematic depiction is given in Supplementary Fig. S1): While synaptic connections along the “trisynaptic pathway” (Cheng, 2013; Jing et al., 2021; Schultz & Rolls, 1999) (entorhinal cortex (EC) → DG → CA3 → CA1) terminate at perisomatic CA3 and CA1 layers (i.e., stratum (St.) pyramidale and St. radiatum), direct EC projections terminate at distal apical dendrites (in St. lacunosum-moleculare) of CA1 and CA3; CA3 receives additional input to proximal apical dendrites via its abundant recurrent collaterals (Sammons et al., 2024). Importantly, each layer maintains its identity and specific connectivity patterns within its respective subfield throughout the entire hippocampus. This means that a layer’s function and inputs are always constrained by the subfield it belongs to.

Recent advances in ultra-high field (UHF) fMRI made it possible to obtain structural (Lüsebrink et al., 2021) and physiological (Haast et al., 2024; Maass et al., 2014) information from the human hippocampus with sub-millimeter resolution. This high spatial resolution allows fMRI responses to be probed as a function of depth (also known as laminar fMRI). This methodology holds the key to investigating the above-mentioned projections non-invasively in humans. However, methodological limitations and physiological constraints pose challenges for the design and interpretation of cortical depth-dependent fMRI experiments of the hippocampus. The most-widely used contrast for laminar fMRI is the gradient echo (GRE) blood oxygenation level dependent (BOLD) contrast due to its high sensitivity (Boxerman et al., 1995; Koopmans & Yacoub, 2019; Weisskoff, 1996) and efficient sequence implementation with echo planar imaging (EPI) readouts. However, GRE-BOLD responses at the laminar level are weighted toward large venous vessels, the well-known draining vein bias. While the main alternative to GRE-BOLD, vascular space occupancy (VASO), can measure laminar responses free of venous contamination in the neocortex (Finn et al., 2019; Huber et al., 2014, 2017), it suffers from an intrinsically low temporal efficiency and reduced sensitivity to CBV changes in inferior brain structures which renders VASO currently unviable in the hippocampus (Huber et al., 2024). In the neocortex, the venous drainage pattern and its effect on laminar fMRI has been thoroughly described (Havlicek & Uludağ, 2020; Polimeni et al., 2010; Turner, 2002; Uludağ et al., 2009). However, because both the anatomy and the vascularization of the hippocampus are considerably different from those of the neocortex, it is currently not known how the venous bias would affect hippocampal laminar BOLD responses. Hence, it is imperative to characterize the distribution of large blood vessels and their influence on hippocampal laminar fMRI. Furthermore, even though UHF fMRI provides a higher baseline signal-to-noise ratio (SNR) and an improved BOLD contrast (Gati et al., 1997; Ogawa et al., 1993), it is also more prone to artifacts associated with main magnetic field (B₀) and transmission field (B₁) inhomogeneities (Farahani et al., 1990; Van De Moortele et al., 2005). The anatomical location of the hippocampus renders it particularly susceptible to these inhomogeneities and hence a challenging region for imaging. Therefore, it is necessary to characterize high-resolution functional imaging of the human hippocampus in terms of depth-dependent mapping of the T₂* relaxation time, signal stability, and noise characteristics to allow for a better interpretation of high-resolution GRE-BOLD activation patterns.

As a prerequisite for transferring neuroscientific fMRI experiments from low to high spatial resolutions and for validating other fMRI contrasts and/or sequences, a benchmark experiment is necessary. In the context of autobiographical memory (AM), neuropsychological and fMRI studies have shown that several subfields of the hippocampus are involved in AM. The subfields include CA1 (Bartsch et al., 2011), CA3 (Bonnici et al., 2013; Chadwick et al., 2014; Miller et al., 2020), and the subiculum (Barry et al., 2021; Leelaarporn et al., 2024), making an AM experiment an ideal candidate for benchmarking.

In this work, we leverage sub-millimeter GRE-BOLD fMRI performed at 7 T to answer two fundamental questions on laminar fMRI of the human hippocampus:

• 1)

  What is the laminar distribution of venous draining patterns in different hippocampal subregions and how does it affect the interpretation of hippocampal laminar responses?

• 2)

  Can we perform a benchmark laminar fMRI experiment which robustly elicits single-subject hippocampal activation using the GRE-BOLD contrast? Can we obtain interpretable results despite the venous bias?

To answer the first question, we invited nine participants to perform a breath-holding experiment while acquiring high-resolution multi-echo T₂*-weighted images. Additionally, we acquired high-resolution structural T₁- and susceptibility-weighted images as well as time-of-flight (TOF) images. The hypercapnic challenge led to a ubiquitous increase in blood flow and thus to a widespread BOLD response. Sampling at different echo times allowed us to differentiate signal changes inside venous vessels from signal changes in the parenchyma. To address the second aim, the same cohort was scanned while performing an autobiographical memory task on a separate day. This task was shown to elicit reliable single subject activity at high-spatial resolution using GRE-BOLD (Leelaarporn et al., 2024). By modeling the hippocampus as a folded surface, we investigated unique depth-dependent features of the hippocampal subfields during memory processing contrasted against a control condition, and during the construction (i.e., initial search phase for a specific memory) versus the elaboration (i.e., reliving the perceptual details of the memory with autonoetic consciousness) of autobiographical memories (Daviddi et al., 2023). We linked the results of the two experiments to interpret the resulting profiles during AM with regard to neural activity and to obtain MRI parameters relevant for future neuroimaging studies of the hippocampus. Altogether, the presented methods and results will facilitate interpretation of hippocampal laminar fMRI responses and validation of other fMRI contrasts. Furthermore, our findings provide mechanistic insights into hippocampal function at the mesoscale.

All data were acquired on a 7 T system (MAGNETOM Terra, Siemens Healthcare, Erlangen, Germany) equipped with a 1 channel transmit, 32 channel receive head coil (NOVA Medical Inc.). Nine healthy participants (3 female, 26.4 ± 6.0 years) were scanned in both experiments after giving informed consent according to the guidelines of the local ethics committee. All processing and analysis code can be found in the public Github repository (https://github.com/viktor-pfaffenrot/hippocampus_laminarfMRI_code).

Whole-brain T₁-weighted anatomical reference data were acquired with an MP2RAGE sequence (Marques et al., 2010) with fat navigators (Gallichan et al., 2016) at 0.75 mm isotropic resolution. To obtain veno- and angiograms, vendor-provided susceptibility-weighted (SWI) and time-of-flight (TOF) sequences were used to acquire images at 0.3 x 0.3 x 0.6 mm resolution (interpolated to 0.15 mm in-plane). The field-of-view (FOV) was positioned such that the long axis of the FOV was parallel to the long axis of the hippocampus (see Supplementary Fig. S2A). Sequence details are shown in Table 1.

For the breath-hold experiment, data were acquired at 0.8 mm isotropic resolution using a 2D multi-echo FLASH sequence. The FOV was positioned similarity to that of the SWI and TOF acquisitions. Echo times were TE = 2, 6, 10, 14, 17, and 21 ms. An additional phase stabilization echo (Hu & Kim, 1994) was acquired at the start of each echo train. This echo serves the purpose of capturing phase fluctuations in the raw MR signals which are corrected in the subsequent offline reconstruction (see section 2.1.2.2). The effect the phase stabilization has on the data quality is shown in Supplementary Figure S2B-E. Shot-to-shot TR was 976 ms and 40 slices were acquired, resulting in a volume TR of TR_vol = 40992 ms. For each participant, two runs of 19 volumes were acquired with each run lasting 14:36 min.

The maximum BOLD response depends on the breath-hold duration and, for the volume TR used in this study, peaks about 10 s after the participant started to breathe again (Haller & Bartsch, 2009). To account for this delay, the paradigm was designed as follows: One volume of normocapnia (labeled as ‘rest’) was followed by a transition volume where the participant was cued to hold their breath in 5 s. The participants were instructed to slowly inhale into their abdomen once the breath-holding command appeared and to hold their breath until the command disappeared. Ten seconds before the breath-hold period ended, the acquisition of the active volume started. The temporal shift between paradigm and image acquisition was designed such that the peak of the BOLD response (10 s after hypercapnia) coincided with the acquisition of the k-space center of the ‘active’ volume. The participants were instructed to slowly exhale after the breath-hold command disappeared. The block of ‘rest,’ ‘transition’, ‘active’ was repeated six times and a run ended with a rest volume.

The image corresponding to the first inversion time (INV1) and the unified image (UNI), that is, the image calculated from both INV1 and the 2^nd redout (INV2), of the MP2RAGE sequence were pre-processed using presurfer (Kashyap, 2021). Both processed images were denoised using the BM4D filter (Lüsebrink et al., 2021; Maggioni et al., 2013). The image quality improvement is shown in Supplementary Figure S3, and filter parameters deviating from default can be found in our github repository. To obtain surface boundaries and hippocampal subfield segmentations, the snakemake (Köster & Rahmann, 2012) implementation of hippunfold (DeKraker et al., 2018, 2022) (https://hippunfold.readthedocs.io) was used. In particular, the denoised INV1 image was used as a proxy for a T₂-weighted input due its T₂-like image contrast (see Fig. 1 and Supplementary Fig. S3). Figure 1 shows the 3D input of a representative subject to hippunfold (1A, B) and zoomed sections of the hippocampus together with the generated segmentations and surfaces (1C, D). The orientation shown in D is identical to that of the schematic drawing in Supplementary Figure S1. In order to present subject-averaged surface representations, we created a ‘canonical’, that is, average surface for all subjects using the connectome workbench (Marcus et al., 2011) command surface-average with the hippunfold surfaces of all subjects as input. Figure 1E shows the 3D canonical surface model of the outer surface (blue surface in C and D) of the left hippocampus in native, that is, folded space, while 1F depicts the unfolded 2D representation of that surface.

Both SWI and TOF images were denoised using BM4D with a fixed noise estimate σ of 30. SWI was receiver bias-field corrected and filtered for veins using a Frangi vesselness filter (Frangi et al., 1998) (α = β = 0.5, Fig. 2B). TOF images were filtered for arteries using a custom filter based on a Haar 2-D wavelet transform. Both vessel-filtered images were subsequently registered to the T₁-weighted image using ITK-SNAP Version 3.6.0 (Yushkevich et al., 2006) and SPM12 Version 7487 (Welcome Trust Centre for Neuroimaging, London). Finally, the resampled images were projected onto the inner and outer surface of the hippocampus as given by hippunfold using the connectome workbench’s (Marcus et al., 2011) command volume-to-surface-mapping. For subject-averaged results, these surfaces were averaged over subjects and plotted on the canonical surface.

Breath-hold scans were reconstructed offline using BART (v. 0.8.00) (Uecker et al., 2015). In particular, ESPIRiT (Uecker et al., 2014) was used to obtain coil sensitivity maps from the autocalibration signal (ACS) before running an iterative SENSE (Pruessmann et al., 1999, 2001) reconstruction with a regularization factor of 0.005 (pics -r 0.005). To account for potentially low SNR, all volumes were denoised using the BM4D filter with the same settings as used for the MP2RAGE.

Denoised FLASH images were pre-processed similar to previous work (Pfaffenrot & Koopmans, 2022; Pfaffenrot et al., 2021) using MATLAB code based on SPM12. In short, after slice-time correction, data were realigned within each run using a brain mask as weights. Motion between runs was accounted for by realigning a biasfield-corrected volume average of the third TE (10 ms) of each run. Registration to the anatomy was performed using ITK-SNAP and a rigid transformation, guided by a mask which excluded areas of strong signal dropout in the sphenoid sinus and maxillary sinuses. All transformation matrices were concatenated, and a single resampling operation using a 4th-order spline was applied.

The functional data were sampled as a function of cortical depth using custom MATLAB code (https://github.com/viktor-pfaffenrot/hippocampus_laminarfMRI_code/blob/main/layerfication/VPF_create_hippocampus_layers.m). In particular, we used the surfaces generated by hippunfold, labeled inner, midthickness, and outer. These, in turn, are formed first, by solving the Laplace equation for two perpendicular directions to obtain a single surface (DeKraker et al., 2022). This is then expanded along the depth of hippocampal GM while obeying the equivolume principle (Waehnert et al., 2014) to form the midthickness and the outer surface. For each hippocampal subfield, signals were sampled between the surface vertices of the inner and outer hippocampal boundary, that is, the inner and outer hippunfold surface (see Supplementary Figs. S1 and S4A, B). In particular, the distance between the vertices defining the inner and the vertices defining the outer surface was equidistantly split into 20 coordinates, i.e., bins. Here, the outer surface refers to the transition region between stratum oriens and parts of St. pyramidale closest to St. oriens. The inner surface corresponds to the transition between St. radiatum and parts of St. pyramidale closest to St. radiatum. The voxel signal was resampled at these coordinates using a 3rd-order polynomial. To account for the high curvature of the hippocampus, the midthickness depth, provided by hippunfold and obeying the equivolume principle, was used as a landmark in the sampling process. That is, the sampling algorithm is forced to include the midthickness depth into the set of bins. The effect on the laminar profiles is shown in Supplementary Figure S4C-F. In case of CA1 and CA2, the sampling was extended by 10 bins beyond the inner surface, that is, further toward the hippocampal sulcus, to also sample parts of the St. lacunosum and St. moleculare, i.e., the innermost layers corresponding to the location of distal apical dendrites. In case of CA4 and the DG, only two bins were defined as there is no clear definition of ‘inner’ and ‘outer’ surface. Hence, given the anatomical structure as shown in Figure 2A (mod. From Duvernoy et al., 2013, Fig. 3.7) and the segmentation of hippunfold, the entire DG was labeled ‘inner’ and the entire CA4 was labeled ‘outer’. This procedure was carried out for each functional volume and allowed us to obtain, separately for each subfield, a time course at each cortical depth that follows both the unique anatomical and laminar characteristics of each subfield.

The projected veno- and angiograms were rescaled to be within the interval [-100,0] and [0,100], respectively, and vessel density was estimated by dividing the number of vertices passing a threshold over the total number of vertices within each subfield. The threshold was -3 a. u. and 1 a .u. for veins and arteries, respectively, and was chosen to reduce the number of falsely classified vessels.

The laminar data of one subject were disregarded due to too excessive motion (1.5 mm maximum displacement due to translation or rotation in any run). For the remaining eight subjects, data of each run were high-pass filtered with a cutoff period of 1/(246 s) and the first and the last volumes of each run were removed to account for edge effects at the boundaries of each run’s time series as a result of high-pass filtering (Stoewer et al., 2012). Finally, the profiles were averaged over runs and over the rest and the active condition, respectively.

For the rest condition, cortical depth-dependent T₂* was estimated by fitting a mono-exponential signal decay function to each bin and subfield using MATLAB’s Levenberg-Marquardt nonlinear least-squares algorithm.

Breath-hold induced signal change was calculated as breath-hold − normal breathing. To relate signal change to venous vessel density, we averaged over echoes and fitted a linear line to each signal change profile for each subject and subfield to estimate the slope between the outer and inner surface. Although phase stabilization accounts for a high degree of unwanted artifacts, some nuisance signal fluctuations, scaling with echo time, were still present. Using a similar weighted echo approach as described in Poser et al. (2006), we tried to reduce their effect by weighting each echo average and penalizing longer TE as given by

The venous vessel density was compared between the inner and outer surface for each subfield using MATLAB and its build-in Wilcoxon signed-rank test corrected for multiple comparisons with Hochberg’s step-up Bonferroni procedure (Hochberg, 1988). Similarly, the Wilcoxon signed-rank test was used to test whether the estimated slopes of the signal change profiles are significantly different from zero.

The relationship between the slope of the signal change as a measure for the bias in laminar fMRI, and the vessel density difference was modeled using a linear mixed-effects model in R (version 4.3.3) utilizing the lme4 package. To investigate whether the difference in vessel density was a significant main effect, we compared a model with the density change as a predictor with one including only the random effects. These were modeled as random intercepts for subject and subfield. Both models were compared using an approximated F-test based on the Kenward-Roger approach for small sample sizes (Kenward & Roger, 1997) as implemented in the pbkrtest package. Marginal and Conditional R² were obtained with the performance package (Lüdecke et al., 2021).

The same participants who performed the breath-hold task were invited for the autobiographical memories (AM) task, hence anatomical data were not re-acquired. Functional AM data were obtained using the GRE-BOLD contrast and a 3D-EPI readout (Stirnberg & Stöcker, 2021) at 0.9 mm isotropic resolution. The key parameters were: TE/TR/TR_vol = 22.2/58.3/2210 ms, FOV = 162 x 178 x 34 mm. The remaining sequence details are given in Table 1. For each participant, three runs of 489 volumes were acquired with each run lasting 18:09 min. In order to investigate data reliability, two subjects were scanned twice on different days. For distortion correction, the same sequence with reversed-phase encoding was used and 13 volumes were acquired.

The functional task was an autobiographical memory paradigm adapted from Leelaarporn et al., 2024 (see also McCormick et al., 2015). In brief, the task consisted of AM trials as a test condition that were randomly interleaved with mental arithmetic (math) trials as control. In the AM trials, participants were asked to retrieve autobiographical episodic memories following the presentation of general cues (e.g., ‘party’, ‘pet’, etc.). They had to press a button to indicate when they had found a memory related to the cue, and to vividly imagine the event with autonoetic consciousness with as much detail as possible for the remainder of the trial. The before and after button press phases of each trial, thus, correspond to the construction and elaboration phases of autobiographical memory retrieval, respectively (Daviddi et al., 2023). Participants were additionally asked to retrieve recent autobiographical memories (no older than 2 years) to ensure hippocampal activation as remote autobiographical memories have been suggested to evoke less hippocampal activity compared to recent ones (Gilmore et al., 2021).

The AM trials were randomly interleaved with math trials in which subjects had to perform a simple arithmetic operation (e.g., 30 + 15). Participants were instructed to press a button once they found the result, and to recursively add 3 to the result (e.g., 45 + 3 + 3 + 3…) until the end of the trial in order to keep them engaged for the whole trial duration.

Each AM and mental arithmetic trial lasted 17.6 s. The inter-trial interval was a fixation cross with a duration of either 2.2 or 4.4 s (in order to be a multiple TR_vol), which was randomly determined. Each run consisted of 44 trials. Unique cue words and arithmetic operations were presented during each run.

A brief training was undergone by each participant prior to the experiment as well as a post-experiment interview to ensure that participants understood the instructions well, and deemed the task’s difficulty manageable. On average, participants took 3.99 s ± 1.78 s to find an autobiographical memory, and 3.19 s ± 1.72 s to resolve the arithmetic operations.

The first three volumes of all runs were disregarded to allow magnetization to reach a steady state. The remaining, online reconstructed magnitude and phase images were used in an offline partial Fourier reconstruction utilizing the POCS algorithm (Haacke et al., 1991; Völker, n.d.). POCS-reconstructed magnitude data were subsequently processing with a custom pipeline based on Advanced Normalization Tools (ANTs) (Avants, Tustison, Wu, et al., 2011) (https://github.com/ANTsX/ANTs). In short, both functional and reversed-phase encoded images of each run were subject to rigid-body correction, followed by a rigid alignment of the temporal mean of the functional and the reversed-phase encoded data. For distortion correction, an undistorted template image was calculated using the script antsMultivariateTemplateConstruction2 (Avants, Tustison, Song, et al., 2011) employing ANTs’ symmetric normalization (SyN) algorithm. The distortion-corrected template image of each run was subsequently rigidly registered to the anatomy using ITK-SNAP and the same masking principle as described in section 2.1.2.3. All transformations and distortion-correction warps were concatenated, and only one single resampling operation using Lanczos-windowed sinc interpolation was applied to the unprocessed functional data to minimize the resolution loss due to multiple interpolation steps.

After pre-processing, functional data were sampled as described in section 2.1.2.4 to obtain depth-dependent time courses. In addition to the procedure described above, we sought to reduce the venous bias in GRE-BOLD laminar fMRI by masking out venous vessels. To this end, the projected SWI surface, filtered for venous vessels (see section 2.1.2.1 and Fig. 2B), was sampled and used as a venous mask. In particular, surface vertices exhibiting a value below -3 a.u. were labeled as belonging to a vessel and disregarded.

Functional data were analyzed in a two-step procedure using robust weighted least squares (Diedrichsen & Shadmehr, 2005) within the GLM framework in SPM12. In a first pass, data were analyzed in voxel space. A single design matrix was used with three regressors of interest: The AM condition, split into construction phase and elaboration phase, and the math condition. The six motion parameters were used as nuisance regressors, and the data were highpass filtered with a 1/(128 s) cutoff before fitting the model.

It is well known that physiological noise, that is, respiration and heartbeat, is a major confounding factor in cortical areas near large vessels and ventricles (Krüger & Glover, 2001; Murphy et al., 2013) especially at ultrahigh fields (Triantafyllou et al., 2005). To account for the pronounced physiological noise at the level of the hippocampus, we applied acompCor (Behzadi et al., 2007). AcompCor (anatomical component based noise correction) estimates significant principal components from noise regions of interest (ROIs) which are used as nuisance regressors in the GLM. In particular, we used two regions of interest (ROIs). For one region, we created a mask of high residuals as given by the residuals map of the first GLM fit. A residual was classified as high if its value was larger than three times the SD of all residuals. This mask should reflect signal fluctuations associated with respiratory and cardiac cycles. A second ROI was manually drawn on the white matter of the brain stem close to the hippocampus. Care was taken to not include GM voxels of the hippocampus in any noise ROI. For the first ROI, the number of nuisance regressors was chosen such that 50 % of the variance was explained. For the second ROI, the number of components was fixed to five. In total, this resulted in approximately 11 to 14 additional acompCor regressors for each run. To increase the predictability of the model, the acompCor regressors were orthogonalized with respect to the motion regressors. After obtaining the acompCor regressors, the model was re-fitted at the voxel level. An evaluation of this procedure can be found in the Supplementary Material (Supplementary Fig. S5). Finally, two contrasts were generated, one as the difference between construction and elaboration phase (labeled pre > post) and one as the difference between the average of the two AM phases and the control condition (labeled memory > math).

Prior to applying the same GLM model on depth-dependent signals, the layer time courses were baseline z-transformed. That is, the z-transformation was performed using the mean and the SD calculated from the volumes corresponding to the math condition.

Temporal SNR was calculated by dividing the temporal mean over the temporal SD where only the volumes corresponding to the math condition were included.

To investigate cortical depth-dependent physiological noise variations in the hippocampus, we performed a Weisskoff test in a depth-specific manner. The Weisskoff test (Weisskoff, 1996) was originally devised as a method to assess scanner stability. Originally, the test was performed by continuously increasing the number of voxels in a spatial average before calculating the temporal SD. If the temporal noise is thermal, that is, of random nature, the noise should decrease with $N$⁠, where N is the number of voxels used in the average. But, at some level of N, any additional voxel in the average does not lead to a noticeable reduction in the noise level. At this point, physiological noise, that is, structured noise, dominates the overall noise profile. On the laminar level, the Weisskoff test was used to estimate the minimum extend of an ROI that was needed to be in the physiological noise dominated regime (Koopmans et al., 2011; Markuerkiaga et al., 2021).

In this work, we averaged over hippocampal hemispheres and performed the Weisskoff test for each subject, each subfield, and each depth separately. In particular, we selected profiles at random and varied the number of profiles included in the average from 1, that is, no averaging, to the maximum number of vertices within each subfield before calculating the temporal noise level (as defined by the temporal SD over time for volumes in the math condition). Finally, we averaged over subjects. This procedure helped us to investigate the laminar variation in physiological noise and to answer the question of whether our current laminar fMRI experiment is performed in the physiological noise dominated regime.

To confirm that the used autobiographical memory task elicits activation in the same areas as previously reported, we conducted a group-level GLM analysis on the memory > math contrast smoothed with a 2.7 mm Gaussian kernel and registered to MNI space. Permutation testing as implemented in the SnPM toolbox (Nichols & Holmes, 2002) was used to control for the FWE rate. Specifically, 2⁹ permutations were used, and the variance of the model was smoothed with the same kernel as recommended by the developer of the toolbox for small sample sizes. The resulting pseudo t-statistic was thresholded at .05 significance level.

All group-level laminar statistical analyses were performed in R (version 4.3.3) for both functional contrasts (memory > math and pre > post). To assess whether there is a significant main effect of layer, that is, whether the profiles resemble each other, we performed a linear mixed model comparison using the lme4 package for each subfield. We compared three models: A baseline model of the form z ~ (1|subject), where z is the baseline-transformed z-value and subject is the random effect, a linear model (z ~ depth+(1|subject)), where depth is the cortical depth as the main effect, and a quadratic model (z ~ depth+depth²+(1|subject)). In particular, we compared the baseline with the linear model and the linear with the quadratic using an approximated F-test based on the Kenward-Roger approach for small sample sizes (Kenward & Roger, 1997). This was preferred over other methods, such as the Likelihood Test, as it has been shown to produce acceptable Type 1 Error Rate, particularly for small sample sizes (Luke, 2017). Multiple comparisons were FDR-corrected using the Benjamini & Hochberg method (Benjamini & Hochberg, 1995).

Functional contrasts (memory > math vs. pre > post) were compared by performing an LME analysis for each subfield (except CA4/DG) with depth, contrast and their interaction as main effects (z ~ depth+contrast+depth:contrast+(1|subject)) followed by a Type III analysis of variance with Kenward-Roger approximation for degrees of freedom and an FDR correction for multiple comparisons.

To test whether venous vessel masking had an effect on the comparison between cortical depths, we performed a similar analysis as described in the previous paragraph and modeled the response as z ~ depth+masking+depth:masking+(1|subject) followed by testing for significant main and interaction effects for each subfield.

Results related to the breath-hold experiment (T₂*, vessel density and breath-hold induced signal change) as well as results related to the autobiographical memory experiment (tSNR) were mapped on the inner and outer surface of the folded and unfolded hippocampus using modified code from the hippunfold toolbox (https://github.com/jordandekraker/hippunfold_toolbox). In addition to the presented visualizations, surface maps and laminar profiles can be interactively explored in our online app (https://viktor-pfaffenrot.shinyapps.io/hippocampus_data_viewer/).

To better understand the venous drainage pattern in laminar fMRI of the human hippocampus (Fig. 2A), we mapped the hippocampal vascular anatomy using high-resolution SWI and TOF angiography. By applying filter techniques (Frangi et al., 1998) (Fig. 2B) followed by mapping on the folded (Fig. 2C, D) and unfolded (Fig. 2E) surface representation of the hippocampus (see Fig. 1), we obtained a distinct vascular pattern. Figure 2C, D shows the subject-averaged vascular distribution of large vessels at the inner (Fig. 2C) and outer (Fig. 2D) surface (veins and arteries in blue and red, respectively). The outer surface putatively corresponds to St. oriens and parts of St. pyramidale closest to St. oriens. The inner surface putatively corresponds to the transition between St. radiatum and parts of St. pyramidale closest to St. radiatum. Figure 2E shows the same data plotted on the unfolded hippocampal surfaces (red arrows indicate exemplary points where E shows the same data as C and D). Comparing between surfaces, Figure 2C and E show an increased venous vessel density at the inner surface of subiculum and CA1. Both regions show lower densities of large veins at the outer surface (D). In contrast, regions CA2-CA3 (Fig. 2D) exhibit a higher venous density at the outer surface compared to the inner surface. These results appear consistent with the schematic drawings depicted in the Duvernoy atlas (Fig. 2A).

We performed a breath-holding experiment to investigate how these subfield-specific vascular differences translate into biases in the laminar fMRI responses. Because the signal change during breath-holding is completely driven by vascular density distribution (rather than a cognitive task), the signal change distribution should provide information on the expected venous bias. Figure 2F shows unfolded surface maps of the breath-hold induced, subject- and echo-averaged signal change with vascular weighting, respectively. In the echo average, the vascular weighting penalizes signals acquired at longer echo times which are more contaminated by B₀-related artifacts not showing any vascular origin (see section 2.1.3 for details). The maps of signal change largely resemble the pattern seen in the venous vessel density distribution, that is, a high signal change at the inner layer of subiculum and CA1 compared to the outer layer and the reverse pattern, although lower in magnitude, for CA2-4.

To quantitatively link venous vessel density ρ to vascular bias, we first statistically compared the differences in vessel density between inner and outer surfaces for each subfield. Figure 2G shows ρ obtained by aggregating the vertices as shown in Figure 2E for each subfield. A lighter color represents the inner, a darker color the outer surface, respectively. The venous vessel density is significantly higher at the inner surface for subiculum and CA1 (Wilcoxon signed-rank test W = 36, p_up = .031). In case of CA3, the density is significantly higher at the outer surface (W = 1, p_up = .047). CA4/DG shows a trend after multiple comparisons correction (W =3, p_up = .078). Figure 2H shows laminar profiles of the subject- and echo-averaged signal change during the breath-hold experiment with vascular-weighting. Single-subject and subject-averaged multi-echo laminar profiles can be found in Supplementary Figure S6. Shaded areas correspond to the SEM across subjects. The profiles show a strong positive slope from the outer toward the inner surface for subiculum (blue profile, Wilcoxon signed-rank test W = 0, p_up = .023) and CA1 (orange profile, W = 0, p_up = .023). The slopes of CA2 and CA3 (green and red profiles) do not show a trend in either direction (p_up = .945 for both subfields).The DG exhibits a stronger signal change compared to CA4. Figure 2I shows the slope of the signal change (corresponding to Fig. 2H) as a function of vessel density difference Δρ computed from Figure 2G. We conducted a comparison of LME analyses using an approximated F-test (see section 2.2.4). The test revealed a significant effect of vessel density difference (Kenward-Roger approximated F-test, F(1,33.5) = 5.28, p = .028, marginal R² = 0.18, conditional R² = 0.48). The estimated model, plotted as a black line in Figure 2I (shaded area is the model’s SE), indicates a positive relationship between the slope of the breath-hold induced signal change, that is, the venous bias, and the vessel density.

In summary, these results show subfield- and layer-dependent differences in venous vasculature and in breath-hold induced signal changes as a measure of physiological bias. They are strongly related and shed light on the venous drainage pattern and how it affects laminar fMRI of the human hippocampus. In particular, subiculum and CA1 show very prominent layer differences in venous drainage, both exhibiting a strong bias from outer to inner layers. CA3 exhibits a bias in the opposite direction, although smaller in magnitude. Layer differences in CA2 and between DG and CA4 appeared less pronounced and did not reach significance after correction for multiple comparisons.

Next, we sought to establish a benchmark laminar fMRI experiment of the human hippocampus utilizing the most-widely used contrast and sequence for laminar fMRI: GRE-BOLD 3D EPI. To this end, inspired by the results of Leelaarporn et al. (2024), we performed an autobiographical memory experiment (see section 2.2.1.1) with the same participants as in the breath-holding experiment.

To verify that the AM experiment elicited similar activation patterns as previously described (Leelaarporn et al., 2024), we performed a non-parametric group-level analysis at the voxel level (Supplementary Fig. S7) showing activations at the expected locations (primarily anterior body of the hippocampus).

Figure 3 shows the cortical depth-dependent profiles for the baseline z-transformed and venous vessel masked memory > math contrast for all subjects and subfields along with the subject average (shaded area correspond to SEM). CA4 and DG again were not subdivided into layers. For all remaining subfields, we tested for significant linear and quadratic main effect of depth by performing an approximated F-test with Kenward-Roger approximation for degrees of freedom and restricted maximum likelihood (for details, see section 2.2.4 and Luke, 2017). Figure 3F shows the fitting coefficients obtained (errorbars indicate SE of the coefficient). The asterisk indicates whether a higher-order model describes the data statistically better than a lower-order model. The LME-based analysis revealed a significant linear effect of depth in all subfields (F (1,260)> 76, p_FDR < .001). A quadratic model fitted the data significantly better than a linear model for the subfields CA1 (F(1,259) = 76, p_FDR = 1.5e-15) and CA3 (F(1,178) = 9.4, p_FDR = .006). In case of CA3, the quadratic model suggests a higher response In middle layers which can be tentatively interpreted as autobiographical memory being predominantly driven by recurrent collateral inputs to CA3, terminating at proximal apical dendrites, rather than EC inputs to apical distal dendrites or DG inputs to perisomatic layers. CA1 exhibits a prominent peak at inner layers which suggests a dominant input from the trisynaptic pathway, that is, from CA3, rather than from EC inputs to distal apical dendrites. Obviously, this interpretation should be considered with caution since we could not directly identify anatomical layers and since these anatomical layers show pronounced physiological differences (see section 4). It is worth noting though that the peak at inner layers of CA1 cannot be explained by the venous bias as our hyperemia experiment would predict a monotonic increase toward St. lacunosum-moleculare (dashed line in Fig. 3B) which is absent in the functional profiles.

As autobiographical memory recall exhibits a temporal dynamic (Daviddi et al., 2023), we investigated whether possible laminar differences between memory construction versus memory elaboration exist. Figure 4 shows the laminar profiles of this contrast (pre > post). Again, LME analysis showed a highly significant linear main effect of depth (F(1,260) > 104, p_FDR < .001). In case of CA3, LME analysis showed a significant quadratic effect of depth (F(1,178) = 37, p_FDR = 3e-8), suggesting that recurrent connections within CA3 play a more important role during memory construction than elaboration.

We compared the two functional contrasts (memory > math vs. pre > post) by performing an LME analysis for each subfield (except CA4/DG) with depth, contrast, and their interaction as main effects (see section 2.2.4 for details). All subfields show a significant difference between contrasts (F(1,528) > 18, p_FDR < .001 for all subfields) and a significant interaction between depth and contrast (F(1,528) > 4, p_FDR < .05 for all subfields). The strongest effects of contrast and interaction are seen for CA1 (F(1,528) = 101, p_FDR = 3e-21 for contrast, F(1,528) = 20.5, p_FDR = 3e-5 for the interaction) where the laminar profiles in Figure 4B do not exhibit a peak at inner layers as shown in Figure 3B but rather a maximum in mid to outer layers. Importantly, this indicates that laminar profiles depend on the exact task demands (i.e., contrasts) rather than reflecting unspecific differences of apparent overall recruitment that may be confounded by vascular (or other) factors.

In summary, these results show that we can obtain consistent laminar fMRI responses from all subfields of the human hippocampus using GRE-BOLD. Specifically, the subiculum shows, on average, the strongest BOLD activity in the memory > math condition, in line with Leelaarporn et al. (2024). Furthermore, CA1 shows laminar responses that differ strongly between contrasts and are not explainable by a venous bias alone. These results can serve as a benchmark for future laminar fMRI studies of the human hippocampus at UHF.

We were interested in whether using a mask of venous vessels as obtained from the SWI image (see section 3.1) would significantly reduce the venous bias in high-resolution GRE-BOLD fMRI. Figure 5 displays the difference between laminar profiles without and with venous masking for both contrasts. To assess whether venous masking has an effect on the laminar profiles, we performed an LME analysis similar to that used in the previous section (see section 2.2.4). For the memory > math contrast (Fig. 5A), LME analysis revealed a significant main effect of masking only for the subiculum (F(1,366) = 9.76, p_FDR = .005). The effect is the strongest at the inner surface as predicted by the hyperemia experiment (Fig. 2H). In case of the pre > post contrast (Fig. 5B), masking had a weak effect on CA2 (F(1,528) = 5.47, p_FDR = .049). CA3 showed a significant interaction between depth and masking (F(1,366) = 9.05, p_FDR = .007) but no main effect of masking (F(1,366) = 4.36, p_FDR = .06).

In summary, these results indicate that vein masking has a noticeable effect on the subiculum which is shown to be most affected by the venous bias. Furthermore, venous masking has a stronger effect during the memory > math contrast compared to the pre > post contrast.

In addition to the laminar fMRI results, the experiments performed in this study allowed us to obtain novel insights into depth-dependent MRI metrics of the human hippocampus relevant for the design and evaluation of other methods than GRE-BOLD.

By fitting a mono-exponential decay function to the multi-echo data during normocapnia (sections 3.1, 2.1.3), we were able to obtain baseline T₂* relaxation times for each subfield as a function of cortical depth. Figure 6A-D shows the fitted T₂* on the unfolded inner (Fig. 6A) and outer (Fig. 6B) surface as well as on the folded outer surface (Fig. 6C). Aggregating vertices results in the depth-dependent plots shown in Figure 6D. Interestingly, the subiculum exhibits a comparably short T₂* relative to the remaining subfields. This can be explained by a high myelin content in the subiculum as noted with histology (Ding & Van Hoesen, 2015).

The outer surface as shown in Figure 6B,C exhibits a noteworthy area of long T₂* at the demarcation between CA1 and CA2. One possible cause for this could be a partial voluming effect with CSF in the proximity of the hippocampus (cf. Supplementary Fig. S4).

As a measure of overall signal stability, temporal SNR is often reported. In this study, we obtained tSNR on the voxel-level (tSNR of two subjects scanned twice on different days shown in Supplementary Fig. S8) and as a function of cortical depth by calculating the temporal mean and standard deviation over volumes corresponding to the control condition of the AM experiment. In Figure 6E-G, the tSNR is shown as a projection on the unfolded inner and outer hippocampal surface, as well as the folded outer surface. Here, no masking for venous vessels was performed. The same structure of long T₂* is now seen as high tSNR at the demarcation between CA1 and CA2, resulting in an overall high tSNR at the outer surface of CA1 as shown in Figure 6H. This result supports the notion that partial voluming with CSF is the most likely cause.

Masking venous vessels (Fig. 6I-L) led to an overall reduction in tSNR as the masking process also excluded GM vertices which are partial voluming with veins. However, when averaging over vertices to obtain layer profiles (Fig. 6L), the reduction in tSNR is negligible (i.e., comparing Fig. 6H with 6L).

An important question to answer is whether the noise in the laminar profiles is dominated by physiological noise, that is, whether averaging does not further improve the tSNR (Krüger & Glover, 2001). To investigate the laminar noise characteristics, we performed a Weisskoff test for each subfield (Weisskoff, 1996) (see section 2.2.3.2). Figure 7A-E shows the subject-averaged temporal noise as a function of cortical depth and amount of averaged vertices (⁠$N$⁠). For all subfields, the noise does not vary considerably with averaging when approaching the maximum number of available vertices. In particular, the physiological noise regime is reached when approximately 2300, 5700, 1100, 1300, and 3400 vertices are averaged for the subiculum and CA1-4, respectively. The number of vertices is well below the total number of vertices available for averaging (for both hemispheres together: 3324, 7139, 1401, 1979, and 4257 for the subiculum and CA1-4), highlighting that in our current experimental setup, we can safely assume to be in the physiological noise-dominated regime.

The noise profiles in the physiological noise regime (Fig. 7F) exhibit a subfield-specific laminar structure. Koopmans et al. (2011) argued that the physiological noise profiles follow the vascular density and hence give estimates of blood volume (Guidi et al., 2020). We also sampled an image of an India Ink stain (Supplementary Fig. S9C, adapted from Duvernoy et al., 2013), to compare profiles of the vascular density with the obtained physiological noise profiles. Results can be found in the Supplementary Material (Supplementary Fig. S9).

In this study, we show that is it possible to image the human hippocampus at a precision high enough to robustly sample BOLD responses and other physiological and MR-related parameters as a function of cortical depth in each individual hippocampal subfield. Based on our findings, we derive the following recommendations, discussed in turn below:

1. Venous bias in high-resolution fMRI and signal interpretation:

   • GRE-BOLD signal changes in inner layers of CA1 and subiculum should be interpreted with caution as they may be amplified due to blood draining from outer to inner layers, that is, the signal change at inner layers of CA1 and subiculum may stem from effects unrelated to the experimental task.

   • For CA1: Compare experimental effects against vascular bias baseline—strong effects different from the here reported bias patterns may indicate true neural activity.

   • For subiculum: Exercise particular caution in interpreting inner layer activity due to strong vascular bias.

   • Consider experimental designs which modulate the strength of different hippocampal inputs rather than contrasting ‘task’ vs ‘rest’ condition to allow for additional control.

2. Reduction Attempts:

   • Consider including a short, high-resolution SWI scan to map large venous vessels for each subject.

   • Use the map of venous vessels to reduce venous bias, in particular for the subiculum.

3. Analysis Strategies:

   • Account for the complex and variable anatomy of hippocampal subfields in the analysis pipelines by using manual segmentation of subfields or highly reliable automatic segmentation tools (e.g., HippUnfold).

   • Correct for physiological noise, especially when imaging at ultra-high fields.

The results of the breath-hold experiment indicate that among all subfields, the subiculum, closely followed by CA1, has the strongest venous bias, blurring activity toward inner layers (Fig. 2H). By mapping venous vasculature (Fig. 2A-G), we could show that the bias is driven by a heterogenous distribution of large venous vessels across layers of the hippocampus. This result was to be expected: The most-widely used GRE-BOLD contrast exhibits a high sensitivity toward brain activity-induced signal changes. However, its sensitivity stems from magnetic field effects around small and large venous vessels. While a sensitivity toward the small vessels (i.e., capillaries) is desirable due to their proximity to the neuronal source of activity, large venous vessels drain blood out of the parenchyma and blur local activity over a large spatial scale. In other words, the GRE-BOLD contrast scales with vessel size (Ogawa et al., 1993) and follows the underlying vascular architecture (Havlicek & Uludağ, 2020; Uludağ et al., 2009). Considering hippocampus-specific features, a recent study by Haast et al. (2024) additionally showed that the subiculum is proximal to arterial macrovessels. Therefore, we advise interpreting strong laminar fMRI responses in inner layers of the subiculum, but also CA1, with caution.

The anatomical location of the hippocampus renders it prone to artifacts associated with B₀-inhomogeneity. Although the acquisition of phase stabilization echoes strongly improved image quality and stability (Supplementary Fig. S2), performing a breath-hold experiment by itself inevitably amplifies B₀-inhomogeneity, especially at longer TE. We tried to reduce these artifacts by down-weighting longer echoes prior to echo averaging. In this regard, attempts to reduce venous bias by performing a breath-hold calibration scan at typical echo times when using an EPI readout might be less successful in the hippocampus compared to the neocortex (Cohen et al., 2004; Eun et al., 2022).

By acquiring 3D-EPI images at 0.9 mm isotropic resolution, we could obtain GRE-BOLD responses as a function of cortical depth during an autobiographical memory paradigm (Figs. 3 and 4). Comparing these results with those of the hyperemia experiment (Fig. 2), CA1 shows a clearly different behavior than what would be predicted by a venous bias alone. If the venous vasculature distribution was the sole driver of these effects, a peak in St. lacunosum-moleculare of CA1 would be expected in both contrasts (dashed orange line in Figs. 3B and Fig. 4B). However, this is not what we see. In the memory > math contrast (Fig. 3B), the activity increases from outer to inner layers and peaks at the inner layer (around the putative transition between St. pyramidale and St. radiatum) before leveling off when approaching the most inner layers, St. lacunosum-moleculare. This peak suggests a predominant relevance of inputs via the trisynaptic pathway, that is, from CA3, rather than of direct EC inputs during autobiographical memory retrieval. In the pre > post contrast, the profiles of CA1 (Fig. 4B) show a significantly different behavior compared to the memory > math contrast. Here, no difference between construction and elaboration phase exists in the innermost layers (putatively St. lacunosum-moleculare). Instead, the contrast increases from inner to outer layers, suggesting that trisynaptic path inputs from CA3 are more relevant during memory construction than elaboration.

Laminar differences of hippocampal BOLD responses may depend on multiple different factors, which should be disentangled in the future, and the current results can serve as a baseline the community can refer to. Some factors to consider would be differences between the laminar distribution of “inputs”, that is, synaptic inputs (which may arrive at all layers) versus “outputs”, that is, action potentials (which are generated perisomatically, i.e., St. pyramidale) or differences between different types of synaptic inputs (largely entorhinal inputs vs. perforant path inputs). These could be tested via functional connectivity analyses, for example, via PPI (Sharoh et al., 2019). In addition, differences in laminar distribution of inhibitory interneurons might also play a role. The latter could be tested pharmacologically by combining laminar fMRI with an intervention such as application of benzodiazepines (which are GABA-A receptor agonists) or by studying epilepsy patients with altered function in specific populations of interneurons.

Accompanying the laminar fMRI results (a supplemantary laminar fMRI analysis along the longitudinal axis of the hippocampus can be found in Supplementary Figure S10), our study also provides insights into MRI and physiological parameters at the mesoscale. Our T₂* profiles are in line with previous reported T₂* values obtained at a lower spatial resolution at 7 T (Cho et al., 2010), and the ability to obtain depth-dependent profiles adds a new dimension to investigate MRI parameter differences between healthy controls and patients (Cho et al., 2010). From an fMRI acquisition perspective, our results can help to guide researchers in selecting echo times optimal for their specific study. In particular, Figure 6D shows that a longer TE than used in this study would be theoretically more optimal in terms of maximizing BOLD sensitivity. However, the breath-hold results indicate what artifact level can be expected at longer TE. Because we could obtain robust BOLD responses at a TE of 22 ms, we advise scanning at a shorter than theoretically optimal TE to improve temporal efficiency and signal stability.

The obtained tSNR maps (Fig. 6, Supplementary Fig. S8) and profiles indicate that we could obtain robust laminar fMRI responses using 3D-EPI. Excluding vertices shown to be likely co-located with veins reduced tSNR but the effect masking has on the tSNR of aggregated laminar profiles is negligible and it reduces signal changes at surfaces known to be prone to venous bias (Fig. 5). However, it must be noted that the measures to reduce venous bias used in this study can only serve as a 1st-order correction as the BOLD effect of large veins affects all depths of the parenchyma (Bause et al., 2020; Viessmann et al., 2019) and the convolutional effect of draining veins (Markuerkiaga et al., 2016) cannot be accounted for with masking. Hence, other fMRI contrasts like SE-BOLD (Beckett et al., 2020) or non-BOLD alternatives like VASO (Faes et al., 2023; Huber et al., 2014, 2024) are needed to be developed further to obtain laminar fMRI responses with an as low as possible effect of large vessels and our results can serve as a benchmark to compare these contrasts to. Alternatively, other post-processing methods like phase regression (Stanley et al., 2021) could also be employed to reduce venous bias although we were not able to reliably reduce effects of large veins using phase regression, most likely due to differences in SNR between the hippocampus and the neocortex where this method has previously been used (data not shown). In addition to the methods mentioned here, the design of the neuroscientific experiment can also help to reduce effects of blood drainage (Olman et al., 2012). In this study, we contrasted two phases that putatively both involve the hippocampus (construction versus elaboration phase). Alternatively, modulating the strength of different hippocampal inputs, similar to work performed in sensory areas (Lawrence et al., 2019), might be a possibility to control the vascular bias in GRE-BOLD.

The noise curves in Figure 7 all reach a plateau, showing that when aggregating over vertices, physiological noise is the dominant noise source. This is in line with previous studies (Koopmans et al., 2011; Markuerkiaga et al., 2021), indicating that the effect of thermal noise can be neglected. However, given the hippocampus’ susceptibility to signal fluctuations driven by physiology, some measures of noise reduction need to be applied to increase data reproducibility, especially in inferior brain regions. We demonstrated this by assessing model residuals before and after using acompCor (Supplementary Fig. S5) showing that modeling physiological noise sources in the GLM is advisable as it increases the model’s predictability. Future work could include a more systematic evaluation of acompCor and other methods, for example, RETROICOR (Glover et al., 2000). Our own attempts (Pfaffenrot et al., 2024) to compare acompCor and RETROICOR show a superior performance of acompCor, possibly due to a more robust estimation of physiological noise sources when using acompCor.

In this study, we used the default canonical HRF from SPM12. It has been shown that subcortical areas exhibit different HRF nonlinearities compared to the neocortex (Lewis et al., 2018) and that for the neocortex, the HRF is not constant along the cortical depth (Dresbach et al., 2023; Puckett et al., 2016) but exhibits a higher BOLD response at superficial layers and a longer tail compared to deep layers. In future work, a hippocampus- and depth-specific HRF could be mapped analogous to the study of Dresbach et al. (2023) to further characterize the BOLD response of the human hippocampus along the temporal dimension.

An estimate of the subject-averaged cortical thickness for each subfield as shown in Supplementary Table ST1 indicates that while an isotropic spatial resolution of 0.9 mm is sufficient to separate signals from the inner layers (i.e., St. lacunosum-moleculare) and outer layers (i.e., St. pyramidale), resolving more depths without significant partial voluming requires a higher spatial resolution. This, in turn, needs to be traded for temporal resolution, spatial coverage, or SNR.

In conclusion, sampling BOLD responses as a function of cortical depth makes investigation of the underlying functional microcircuits of the human hippocampus possible, and has potential for non-invasively investigating the directionality of the projections to and from other brain regions. Here, we present relevant insights on the vascular bias of laminar fMRI of the human hippocampus which will help the community to interpret hippocampal responses when utilizing the most widely used fMRI contrast, GRE-BOLD. Finally, the presented AM results and methods will help researchers to test new hypotheses on a laminar level and to validate other fMRI contrasts to improve the specificity of laminar fMRI of the human hippocampus.

All processing and analysis code can be found in the public Github repository (https://github.com/viktor-pfaffenrot/hippocampus_laminarfMRI_code). Final surface data can be downloaded from our online app (https://viktor-pfaffenrot.shinyapps.io/hippocampus_data_viewer/). Raw data can be obtained from one of the corresponding authors on reasonable request.

V.P.: Conceptualization, Methodology, Software, Formal analysis, Investigation, Writing—Original Draft, Writing—Review & Editing, and Visualization. A.B.: Conceptualization, Methodology, and Writing—Review & Editing. C.A.G.: Software, Writing—Review & Editing. S.K.: Software, Writing—Review & Editing. N.A.: Conceptualization, Writing—Review & Editing, Supervision, and Funding Acquisition. D.G.N.: Conceptualization, Writing—Review & Editing, Supervision, and Funding Acquisition.

None.

The authors thank Roy Haast for his helpful comments on hippocampal surface generation and for providing code for the online app. The authors would also like to thank the Khan Lab for their support in using hippunfold. The authors thank Rüdiger Stirnberg and the team at DZNE Bonn for providing the 3D-EPI sequence. The research leading to these results has received funding from MERCUR grant number Ko-2021-0010. The MAGNETOM Terra used in the study was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), grant number 432647511. N.A. received funding via ERC-CoG “GridRepresentations” (grant number 864164) and via the DFG SFB 1280 “Extinction Learning” (project number 316803389, subprojects F02 and A02).

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00532.