1 Introduction

The brain’s vascular network plays a crucial role in delivering oxygen and glucose while clearing metabolic by-products to meet the high energy demands of neural information processing. Understanding the organization of the brain’s vasculature is vital for diagnosing and addressing clinical deficits related to stroke, vascular dementia, and neurological disorders with vascular components (Iadecola, 2013; Sweeney et al., 2018; Toledo et al., 2013). Furthermore, it is essential for advancing the applications of functional MRI (fMRI), as vascular density has implications for statistical power, and the arrangement of large vessels may impose limitations and biases on the spatial accuracy of functional localization. Despite its significance, our knowledge of the vascular network architecture in the primate cerebral cortex remains limited (Duvernoy et al., 1981; Schmid et al., 2019; Weber et al., 2008).

Anatomically, blood flows from the pial vessel network via feeding arteries and arterioles to capillary beds in each cortical layer, ultimately leading to draining veins and back to the pial vessel network. The capillary density varies with the rate of oxidative metabolism across cortical layers and exhibits sharp transitions between some cortical areas (Duvernoy et al., 1981; Ji et al., 2021; Zheng et al., 1991). Recent advances in immunolabeling and tissue clearing techniques have enhanced our understanding of brain vascularity in post-mortem mouse brains (Ji et al., 2021; Kirst et al., 2020). These studies have demonstrated heterogeneous vasculature varying in capillary length density 3-fold across brain regions and moderate variation across cortical layers. Still, methodological and analytical challenges have limited quantitative anatomical research in primates to a small number of cortical regions (Harrison et al., 2002; Lauwers et al., 2008; Weber et al., 2008) and quantitative anatomical research requires investigation across a broader range of cortical regions in primates.

Mapping the brain-wide vasculature using MRI faces several challenges due to the intricate nature of the vascular network. One crucial criterion for successful vascular mapping is arterio-venous density, which is necessary to delineate individual large-caliber vessels from microcapillary networks. The combined surface density of intra-cortical feeding arteries and draining veins is about 7 vessels/mm² (Weber et al., 2008). According to the sampling theorem, this implies that the minimal (spatial) sampling frequency is ≈14 voxels/mm² (≈0.26 mm isotropic) imposing stringent image acquisition requirements to critically sample cortical vasculature. Ferumoxytol contrast agent-weighted MRI offers a safe and indirect means to measure relative vascular volume and enhance the visibility of large vessels (Boxerman et al., 1995; Kim et al., 2013; Muehe et al., 2016; Yablonskiy and Haacke, 1994). Compared to gadolinium-based agents, ferumoxytol’s longer half-life allows for higher-resolution vascular volume measurements (Buch et al., 2022), albeit these methodologies are hampered by confounding factors such as vessel orientation relative to magnetic field (B₀) direction (Ogawa et al., 1993).

The macaque monkey is an excellent experimental non-human primate model to objectively investigate the MRI resolution and contrast requirements and their limitations for mapping arterio-capillary-venous networks. Quantitative vascular density data is available for a limited number of cortical areas (Weber et al., 2008), providing essential insights for determining vessel-density informed minimum image resolution requirements. Importantly, experiments in macaque monkeys can also help elucidate the neurometabolic factors that shape vascular network architecture. For instance, variations in the cellular composition (Collins et al., 2010), synaptic density (Elston, 2002), receptor distribution (Froudist-Walsh et al., 2023), neural connectivity (Felleman and Van Essen, 1991; N. T. Markov et al., 2014), myelination (Lewis and Essen, 2000), and oxidative metabolism (Sincich et al., 2003) are well-documented, but the connections between these factors and vascular architecture have only been investigated in a few cortical areas (Borowsky and Collins, 1989; Tsai et al., 2009; Weber et al., 2008).

In this study, we used ferumoxytol contrast agent-weighted 3D multi-echo gradient-echo MRI to investigate vascular heterogeneity in macaque monkey brains. We scaled image resolution to meet specific requirements, aiming to delineate cortical layers and individual vessels. Using this advanced MRI and laminar surface mapping, we then elucidate neuroanatomical factors underlying the heterogenous vasculature. Our analysis reveals insights into translaminar and regional heterogeneities, and signatures of neuroanatomical organization within the macaque cerebral cortex. By addressing these dual objectives— advancing vascular MRI technology and uncovering the neuroanatomical factors shaping cortical vascularity—we contribute to both methodological and conceptual advancements in the field. Our findings offer not only a framework for objectively validating cerebrovascular MRI but also a deeper understanding of how neural and vascular systems intricately interact in the primate cerebral cortices.

2 Methods

2.1 Data acquisition

Experiments were performed using a 3T MRI scanner (MAGNETOM Prisma, Siemens, Erlangen, Germany) equipped with 80 mT/m gradients (XR 80/200 gradient system with slew rate 200 T/m/s), a 2-channel B₁ transmit array (TimTX TrueForm) and a custom-made 24-channel coil for the macaque brain (Autio et al., 2020). The animal experiments were conducted in accordance with the institutional guidelines for animal experiments, and animals were maintained and handled in accordance with the policies for the conduct of animal experiments in research institution (MEXT, Japan, Tokyo) and the Guide for the Care and Use of Laboratory Animals of the Institute of Laboratory Animal Resources (ILAR; Washington, DC, USA). All animal procedures were approved by the Animal Care and Use Committee of the Kobe Institute of RIKEN (MA2008-03-11).

2.1.1 Anesthesia protocol

Macaque monkeys (Macaca mulatta, weight range 7.4 – 8.4 kg, age range 4-6 years, N = 4) were initially sedated with intramuscular injection of atropine sulfate (20 μg/kg), dexmedetomidine (4.5 µg/kg) and ketamine (6 mg/kg). A catheter was inserted into the caudal artery for blood-gas sampling, and endotracheal intubation was performed for steady controlled ventilation using an anesthetic ventilator (Cato, Drager, Germany). End-tidal carbon dioxide was monitored and used to adjust ventilation rate (0.2 to 0.3 Hz) and end-tidal volume. After the animal was fixed in an animal holder, anesthesia was maintained using 1.0% isoflurane via a calibrated vaporizer with a mixture of air 0.75 L/min and O2 0.1 L/min. Animals were warmed with a blanket and water circulation bed and their rectal temperature (1030, SA Instruments Inc., NY, USA), peripheral oxygen saturation and heart rate (7500FO, NONIN Medical Inc., MN, USA) were monitored throughout experiments.

2.1.2 Structural acquisition protocol

T1w images were acquired using a 3D Magnetization Prepared Rapid Acquisition Gradient Echo (MPRAGE) sequence (0.5 mm isotropic, FOV 128×128×112 mm, matrix 256×256, slices per slab 224, coronal orientation, readout direction of inferior (I) to superior (S), phase oversampling 15%, averages 3, TR 2200 ms, TE 2.2 ms, TI 900 ms, flip-angle 8.3°, bandwidth 270 Hz/pixel, no fat suppression, GRAPPA 2, turbo factor 224 and pre-scan normalization). T2w images were acquired using a Sampling Perfection with Application optimized Contrast using different angle Evolutions (SPACE) sequence (0.5 mm isotropic, FOV 128×128×112 mm, matrix 256×256, slice per slab 224, coronal orientation, readout direction I to S, phase oversampling 15%, TR 3200 ms, TE 562 ms, bandwidth 723 Hz/pixel, no fat suppression, GRAPPA 2, turbo factor 314, echo train length 1201 ms and pre-scan normalization) (Autio et al., 2021, 2020)

In a separate imaging session, additional high-resolution structural images were acquired (Autio et al., 2022). T1w images were acquired using a 3D Magnetization Prepared Rapid Acquisition Gradient Echo (MPRAGE) sequence (0.32mm isotropic, FOV 123×123×123 mm, matrix 384×384, slices per slab 256, sagittal orientation, readout direction FH, averages 12-15, TR 2200 ms, TE 3 ms, TI 900 ms, flip-angle 8°, bandwidth 200 Hz/pixel, no fat suppression, GRAPPA 2, reference lines PE 32, turbo factor 224, averages 12-15, and pre-scan normalization). T2w images were acquired using a Sampling Perfection with Application optimized Contrast using different angle Evolutions and Fluid-Attenuated Inversion Recovery (SPACE-FLAIR) sequence (0.32mm isotropic, FOV 123×123×123 mm, matrix 384×384, slice per slab 256, sagittal orientation, readout direction FH, TR 5000ms, TE 397 ms, TI 1800 ms, bandwidth 420 Hz/pixel, no fat suppression, GRAPPA 2, reference lines PE 32, turbo factor 188, echo train duration 933 ms, averages 6-7 and pre-scan normalization). The total acquisition time for structural scans was ≈3h.

2.1.3 Quantitative transverse relaxation rate acquisition protocol

Data was acquired before and after (12 mg/kg) the intravascular ferumoxytol (Feraheme, ferumoxytol AMAG Pharmaceuticals Inc, Waltham, MA, USA) injection using gradient- and RF-spoiled 3D multi-echo gradient-echo acquisition (0.32 mm isotropic, FOV 103×103×82 mm, matrix 320×320, slices per slab 256, sagittal orientation, bipolar read-out mode, elliptical scanning, no partial Fourier, ten equidistant TEs, first TE (TE1)=3.4 ms, time between echoes (ΔTE) 2.4 ms, TR 33 ms, FA 13° (corresponding to Ernst angle of gray matter; median T₁=1370 ms), bandwidth 500 Hz/pixel (fat-water shift one voxel), scan duration 20 min, GRAPPA 2, reference lines 32, and pre-scan normalization). The total acquisition time before and after ferumoxytol injection were 40 and 100 min, respectively.

2.1.4 Vessel-density informed data acquisition protocol

To investigate the periodicity of the penetrating large vessel network, we performed auxiliary ferumoxytol-weighted experiments using image resolution adjusted to satisfy critical (spatial) sampling frequency (14 voxels/mm² ≈0.26 mm isotropic) of intra-cortical vessels (7 vessels/mm²; Weber et al., 2008). The original gradient- and RF-spoiled 3D multi-echo gradient-echo product sequence, however, did not allow sufficient matrix size to satisfy the critical sampling frequency of penetrating vessels. To achieve the target resolution, the sequence was customized by easing the matrix size limitations (by Y.U.). Using the customized sequence, we performed experiments at 0.25 (N=1) and 0.23 mm (N=2) isotropic spatial resolution. Scan #1: (FOV 104×104×80 mm, matrix 416×416, slices per slab 320, sagittal orientation, bipolar read-out mode, elliptical scanning, no partial Fourier, three TEs 6, 10 and 14 ms, TR 22 ms, FA 11°, bandwidth 260 Hz/pixel (fat-water shift 1.6 voxels), scan duration 21 min, GRAPPA 2, reference lines 32 and pre-scan normalization). Scans #2-3: (FOV 103×103×81 mm, matrix 448×448, slices per slab 352, sagittal orientation, bipolar read-out mode, elliptical scanning, no partial Fourier, three TEs 5, 9 and 13 ms, TR 23 ms, FA 11°, bandwidth 340 Hz/pixel (fat-water shift 1.2 voxels), scan duration 25 min, GRAPPA 2, reference lines 32 and pre-scan normalization). The total acquisition time was 150 min.

2.2 Data analysis

Data analysis utilized a version of the HCP pipelines customized specifically for use with non-human primates (https://github.com/Washington-University/NHPPipelines) (Autio et al., 2020; Glasser et al., 2013; Hayashi et al., 2021).

2.2.1 Structural image processing

PreFreeSurfer pipeline registered structural T1w and T2w images into an anterior-posterior commissural (AC-PC) alignment using a rigid body transformation, non-brain structures were removed, T2w and T1w images were aligned using boundary based registration (Greve and Fischl, 2009), and corrected for signal intensity inhomogeneity using B₁-bias field estimate (Glasser et al., 2013). Next, data was transformed into a standard macaque atlas by 12-parameter affine and nonlinear volume registration using FLIRT and FNIRT FSL tools (Jenkinson et al., 2002).

FreeSurferNHP pipeline was used to reconstruct the cortical surfaces using FreeSurfer v6.0.0-HCP. This process included conversion of data in native AC-PC space to a ‘fake’ space with 1-mm isotropic resolution in volume with a matrix of 256 in all directions, intensity correction, segmentation of the brain into cortex and subcortical structures, reconstruction of the white and pial surfaces and estimation of cortical folding maps and thickness. The intensity correction was performed using FMRIB’s Automated Segmentation Tool (FAST) (Zhang et al., 2001). The white matter segmentation was fine-tuned by filling a white matter skeleton to accurately estimate white surface around the blade-like thin white matter particularly in the anterior temporal and occipital lobe (Autio et al., 2020). After the white surface was estimated, the pial surface was initially estimated by using intensity normalized T1w image and then estimated using the T2w image to help exclude dura (Glasser et al., 2013).

The PostFreeSurfer pipeline transformed anatomical volumes and cortical surfaces into the Yerkes19 standard space, performed surface registration using folding information via MSMSulc (Robinson et al., 2018, 2014), generated mid-thickness, inflated and very inflated surfaces, as well as the myelin map from the T1w/T2w ratio on the mid-thickness surface. The volume to surface mapping of the T1w/T2w ratio was carried out using a ‘myelin-style’ mapping (Glasser and Van Essen, 2011), in which a cortical ribbon mask and a metric of cortical thickness were used, weighting voxels closer to the midthickness surface. Voxel weighting was done with a Gaussian kernel of 2 mm FWHM, corresponding to the mean cortical thickness of macaque. For quality control, the myelin maps were visualized and potential FreeSurfer errors in pial or WM surface placement were identified. The errors were manually corrected by editing wm.mgz and by repositioning and smoothing the surfaces using FreeSurfer 7.1, the curvature, thickness, and surface area on each vertex were recalculated, and then PostFreeSurfer pipeline was applied again and T1w/T2w was visually inspected for quality control.

Twelve cortical laminar surfaces were generated based on equivolume model (Autio et al., 2022; Van Essen and Maunsell, 1980) using the Workbench command ‘-surface-cortex-layer’ and the native pial and white surface meshes in subject’s AC-PC space. Throughout the text, the equivolumetric layers (ELs) are referred to as EL1a (adjacent to the pial surface), EL1b, EL2a, EL2b,… and EL6b (adjacent to the white matter surface). This nomenclature is intended to ease but also distinguish comparison between anatomically determined cortical layers which vary in thickness. Anatomical layers are referred to using roman numerals (e.g. Ia, Ib, Ic, IIa,… and VIb). To assess the vascularity on the white matter and pial surfaces, additional layers were generated underneath and just above the gray matter in the superficial white matter and pial surface, respectively. Surface models and data were resampled to a high-resolution 164k mesh (per hemisphere).

2.2.2 Quantitative multi-echo gradient-echo data processing

The original 3D multi-echo gradient-echo images were upsampled to 0.25 or 0.15 mm spatial resolution for the data with spatial resolution 0.32 or 0.23 and 0.26, respectively; and transformed using cubic-spline to the subject’s AC-PC space using a rigid body transformation. Pre- and post-ferumoxytol runs (two and six, respectively) were averaged and R₂*-fitting procedure was performed on multi-TE images with ordinary least squares method in the in vivo histology using MRI (hMRI) Toolbox (Tabelow et al., 2019). The baseline (pre-ferumoxytol) R₂* was subtracted from the post-ferumoxytol R₂* maps to calculate ferumoxytol induced change in ΔR₂* (Boxerman et al., 1995). Subcortical region-of-interests (ROIs; thalamus, striatum, cerebellum, hippocampus, inferior colliculus and corpus callosum) were manually drawn whilst avoiding large vessels using T1w image as a reference. The quantitative R₂* and ΔR₂* maps were mapped in the twelve native laminar mesh surfaces using the Workbench command ‘-volume-to-surface-mapping’ using a ribbon-constrained algorithm. MSMSulc surface registration was applied, the data was resampled to Mac25Rhesus reference sulcus template using ADAP_BARY_AREA with vertex area correction, and left and right hemispheres were combined into a CIFTI file.

Since large-caliber pial vessels run along cortical surface, large penetrating vessels are mainly oriented normal to the cortical surface and the capillary network may be orientated more random to the cortical surface (Ji et al., 2021; Reina-De La Torre et al., 1998), the orientation of the cerebral cortex relative to the direction of static magnetic field (B₀) may bias the assessment of R₂* and ΔR₂* (Bolan et al., 2006; Lee et al., 2011; Ogawa et al., 1993; Viessmann et al., 2019; Yablonskiy and Haacke, 1994). Because the brain R₂* measures are primarily determined by extravascular MR signal, we may assume that

where Θ is the angle between normal of the cortex relative to the direction of B₀. Each vertices Θ was determined in the subject’s original MRI space. The (Pearson’s) correlation coefficient between Eq. 1 and with R₂* and ΔR₂* was estimated in each EL. To remove orientation bias, in Eq. 1 was regressed out from each laminar R₂* and ΔR₂* surface map.

To examine repetitive patterns in the vascularity, the R₂* and ΔR₂* laminar profiles were parcellated using the M132 Lyon Macaque brain atlas (N. T. Markov et al., 2014). Because some of the M132 atlas cortical parcels exhibited a degree of laminar inhomogeneity due to artifacts (e.g., areas adjacent to major sinuses, large vessels penetrating to white matter and FreeSurfer errors in surface placement), median values were assigned to each parcel in each EL. The effect of blood accumulation in the superficial layers was estimated using linear model and regressed out from the parcellated ΔR₂* maps. Subjects (including a single test-retest data set), ELs and hemispheres were combined (5×12×2=120) and hierarchical clustering was applied to parcels using Ward’s method. A dendrogram was used to determine the number of clusters.

To explore sharp transitions in cortical vascularization, each B₀ orientation corrected ΔR₂* EL was smoothed using a factor of 1.2 mm. The smoothing factor was twice the average distance of draining veins (Weber et al., 2008). Then, gradient-ridges were calculated using maps using wb_command -cifti-gradient for each EL. The resulting gradients were then cross-referenced with potential FreeSurfer surface error displacements and when required manual corrections (wm.mgz and using reposition surface in the FreeView 3.0) were performed to the FreeSurfer segmentations.

2.2.3 Vessel detection

To improve contrast-to-noise ratio for vessel detection, multi-echo gradient-echo images were aligned in the native AC-PC space and averaged across runs. Vessels were identified using the Frangi “vesselness” filter which enhances the vessel/ridge-like structures in 3D image using hessian eigenvalues (Avadiappan et al., 2020; Frangi et al., 1998). Volume images in the native AC-PC space were also non-linearly transformed into a standard “SpecMac25Rhesus” atlas (Hayashi et al., 2021).

To facilitate visualization of the pial vessel network, low-frequency fluctuations were removed by subtracting extensively smoothed versions of the post-ferumoxytol TE-averaged EL surface maps. To detect vessels running parallel to the cortical surface, continuous signal drop-outs were clustered along ELs using wb_command -cifti-find-clusters with a criterion of a 0.5 mm² minimum cluster area and visually determined intensity threshold.

To enable surface detection of penetrating vessels, an ultra high-resolution 656k cortical surface mesh was generated using wb_command -surface-create-sphere resulting in an average 0.022 vertex surface area, approximately half the isotropic 0.23 mm voxel (face) surface area (0.053 mm²). Then, TE-averaged multi-echo gradient-echo images and Frangi-filtered vessel images were mapped to twelve native mesh laminar surfaces in the subject’s physical space. To detect penetrating vessels oriented perpendicular to the laminar surfaces, localized signal drop-outs were visualized using wb_command -cifti-gradient and their central locations were identified by detecting the local minima using wb_command -find-extrema.

The number of vessels in V1 was estimated using the M132 parcellation as a reference (Markov et al., 2014). In the native space, the surface area of each vertex was determined using wb_command -surface-vertex-areas. Then, the surface-area map was transferred into atlas space and the area of V1 was determined using M132 areal atlas. The V1 vessel density was determined by dividing the number of vessels (by a conservative estimate using Frangi-filtering and a liberal estimate using local minima in gradient) by surface-area. These MRI estimates were then compared with histological large vessel densities in the V1 (Zheng et al., 1991, Weber et al., 2008).

To determine the periodicity of the cortical arterio-venous networks, non-uniformly sampled Lomb-Scargle geodesic periodogram analysis (Matlab Signal Processing Toolbox, The MathWorks Inc., US) was performed on the spatially low-frequency filtered twelve native mesh ELs in the subject’s physical space. The analysis was limited to the closest 2000 geodesic vertices within 20 mm geodesic distance of manually selected vertices in V1. Geodesic distance between vertices was calculated with wb_command -surface-geodesic-distance in each EL. Periodograms were binarized with an equidistant interval (=0.05 1/mm) up to 10 1/mm, and then 95% confidence interval of the mean of magnitude was estimated from bootstrap and compared across cortical laminae.

2.2.4 Comparison with histological datasets

In V1, translaminar ΔR₂* was compared with cytochrome oxidase (CO) activity, capillary and large vessel volume fractions (Weber et al., 2008) These ground-truth measures were estimated from Weber and colleagues Figure 4. Each measure was peak normalized, so that values ranged between zero and one, to compare different contrasts across the cortical layers.

To investigate the correspondence between regional variation in cerebrovascular volume and heterogeneous neuron density, we used the 42 Vanderbilt tissue sections covering the entire macaque cerebral cortex (Collins et al. 2010) available from the Brain Analysis Library of Spatial Maps and Atlases (BALSA) database (Froudist-Walsh et al., 2023, 2021; Van Essen et al., 2017). These sections were processed using the isotropic fractionator method to estimate neuron densities (Collins et al., 2010). The sections were used to parcellate R₂* and ΔR₂* and these were then compared with neuron density using Pearson’s correlation coefficient. To compare neuron and total receptor densities with R₂* and ΔR₂*, we also applied the Julich Macaque Brain Atlas for parcellation (Froudist-Walsh et al., 2023). The parcellated neuron and total receptor densities were used in linear regression model to predict R₂* and ΔR₂* across ELs and the resulting T-values were then threshold at significance level (p<0.05, Bonferroni-corrected).

Relation between cerebrovascular volume and parvalbumin and calretinin positive interneurons, collated from multiple studies and ascribed to M132 macaque atlas by Burt and colleagues (Burt et al., 2018; Condé et al., 1994; Gabbott and Bacon, 1996; Kondo et al., 1999), were compared across ELs using Pearson’s correlation coefficient. The resulting p-values were Bonferroni-corrected for the number of layers and contrasts.

The number of dendritic spines (putative excitatory inputs) and dendrite tree length size were also obtained from BALSA (Elston, 2007; Froudist-Walsh et al., 2023, 2021). The region-of-interests, described in Elston 2007 and plotted on the Yerkes surface by Froudist-Walsh and colleagues, were used to obtain median value of R₂* and ΔR₂* and these were then compared with the number of dendritic spines and dendrite tree length size using Pearson’s correlation coefficient.

3 Results

3.1 Laminar R₂* and ferumoxytol induced ΔR₂* MRI in macaque cerebral cortex

Figure 1 displays representative gradient-echo images before (Fig. 1A) and after (Fig. 1B) intravascular ferumoxytol injection (N=1). The ferumoxytol effectively reversed the signal-intensity contrast between gray matter and white matter while enhancing the visibility of large vessels, as expected. For quantitative assessment, R₂* values were estimated from multi-echo gradient-echo images acquired both before and after the administration of ferumoxytol contrast agent (Table 1). Subsequently, the R₂* and ΔR₂*, an indirect proxy measure of CBV (Boxerman et al., 1995), volume maps for each subject were mapped onto the twelve native equivolumetric layers (ELs) (Fig. 1C). Each vertex was then corrected for normal of the cortex relative to B₀ direction (Supp. Fig. 1A-C). Surface maps for each subject were registered onto a Mac25Rhesus average surface using cortical curvature landmarks and then averaged across the subjects (Fig. 1D, E). Around cortical midthickness, the distribution of R₂* resembled the spatial pattern of ΔR₂* vascular volume. However, in translaminar domain they exhibited reversed patterns: R₂* increased whereas ΔR₂ decreased toward the WM surface (Fig. 1D, G).

Figure 1.

Table 1.

To explore heterogeneous brain vascularity, we investigated ΔR₂* in selected subcortical regions (Fig. 1F). We found the highest CBV in the inferior colliculus, an early auditory nucleus, and the lowest in the corpus callosum. Overall, the relative blood volume variations among the investigated subcortical regions were comparable to those reported in mice (Kirst et al., 2020).

Adjacent to the pial surface, large vessels exhibited notable signal-loss in the superficial gray matter. To visualize the pial vessel network, we removed very low-frequency components from TE-averaged post-ferumoxytol signal-intensity maps and identified continuous signal drop-outs along ELs using clustering. Within the most superficial layers (e.g. EL1a-2a), this analysis revealed an extensive arterio-venous pial vessel network spanning almost the entire cortical surface (Fig. 2A). We attempted to delineate pial arteries and veins using pre-contrast R₂*-values; however, due to the ‘blooming’ effect of ferumoxytol (Buch et al., 2022) distinguishing adjacent large-caliber vessels was difficult to differentiate with high confidence. Additionally, the continuity of the pial vessel network may also have been influenced by veins crossing the sulci (Duverney et al., 1981). Pial vessel network was consistently observed across subjects (Supp. Fig. 2A), although the precise locations of vessels did vary across subjects. In contrast, in the middle and deep cortical layers (EL2b-6b) continuous signal-dropout clusters were largely absent demonstrating that the influence of the large pial vessels was minimal in these layers (Supp. Fig. 2B).

Figure 2.

To visualize the intra-cortical vessel network, we next performed ferumoxytol-weighted experiments with isotropic image resolution of 0.23 mm adjusted below to the critical (spatial) sampling frequency of large penetrating vessels (19 voxels/mm² vs 7 vessels/mm²) (Zheng et al., 1991; Weber et al., 2008). The post-ferumoxytol signal-intensity maps (Supp. Fig. 3A) were used to identify vessels in volume space using the Frangi filter (Supp. Fig. 3B) and in the cortical surface by calculating sharp gradients and determining their local minima (Supp. Fig. 3C). Local minima, however, by mathematical definition can capture 1 vessel per 7 vertices (each vertex contains six neighbors). To address this limitation, we generated an ultra-high cortical surface mesh (656k) with an average vertex area of 0.022 ± 0.012 mm² (1st-99th percentile range: 0.006-0.065 mm²) (Fig. 2B). Within the cortical gray matter, we identified an average of 24,000 ± 2,000 penetrating vessels per hemisphere (Fig. 2B, C; for all ELs see Supp. Fig. 4). In V1, we found 1.9-2.2 vessels/mm² using Frangi filter and surface vessel detection, respectively (Fig. 2E). This vessel density corresponds to about 30% of the anatomical ground-truth (Zheng et al., 1991, Weber et al., 2008).

To corroborate the periodicity of the cerebrovascular network, we next applied a non-uniformly sampled lomb-scargle geodesic periodogram analysis on the signal-intensity averaged native ELs (Fig. 2D). The periodograms revealed dominant periodicity at spatial frequency of ≈0.6 1/mm in the most superficial ELs, likely reflecting the presence of large-caliber vessels in the pial network. In the middle ELs, bimodal distribution was observed with peaks at around 0.6 and 1.2 1/mm. In the deep ELs, peak power occurred at a shorter distance (≈1 1/mm), potentially indicative of the large arteries supplying the white matter. These findings underscore the substantial variation in vascular organization across the cortical layers.

To explore areal differences in translaminar features, we next parcellated the dense R₂* and ΔR₂* maps using M132 cortical atlas (Fig. 2A-D). To mitigate confounding effects from the large-caliber vessels, median parcel values were used for parcellation, ΔR₂* profiles were detrended across ELs and then averaged across subjects. In the EL4b, which approximately corresponds to the location of histologically defined thalamic input L4c, the primary visual area exhibited larger vascular volume in comparison to surrounding cortical areas (Fig. 3B, D). Within the visual system, inspection of laminar profiles revealed distinct features.

Figure 3.

Since the ΔR₂* is an indirect proxy measure of vascular volume (Boxerman et al., 1995), we next sought to validate the noninvasive laminar ΔR₂* maps with respect to quantitative histological assessment of vascular properties in the macaque visual cortex (Weber et al., 2008). In V1 we found that ΔR₂* more closely resembled microcapillary and oxidative metabolism rather than large vessel volume fraction (Fig. 3F), albeit we could not identify the vascularity peak in L6 potentially due resolution limitations. Moreover, the V2/V1 ΔR₂* ratio in EL4b (79% ± 5%) (Fig. 3E) was also in excellent agreement with the previous reports of capillary volume (Zheng et al., 1991; Weber et al., 2008). Finally, the average ΔR₂* ratio between V1 gray matter and the underlying white matter (2.2 ± 0.1) was also close to the histological assessment (1.8-2.1). Taken together, we found comparable relative variations in vascular volume with anatomical ‘ground-truth’ substantiating the validity of our noninvasive methodology.

3.2 Variations in cerebrovascular network architecture reveal inter-areal boundaries

Since cellular composition (Collins et al., 2010) and oxidative metabolism (Sincich et al., 2003) are known to exhibit sharp transitions between cortical areas, we next tested the hypothesis whether the variations in vascular network architecture may also reveal inter-areal boundaries (Zheng et al., 1991). To address this question, we calculated the gradient ridges of ΔR₂* in each EL (Fig. 1E). Due to the strong cortical contrast (ΔR₂* = 39 ± 2 ms), the resulting gradients were notably strong and revealed several sharp transitions (Fig. 5A, B). A particularly strong gradient was observed at the boundary between V1/V2 at EL4b (Fig. 3D and Fig. 4A, B), attributable to the relatively large capillary density difference between the areas (Fig. 3E) (Duverney 1981; Zheng et al., 1991; Weber et al., 2008). We also found a sharp ΔR₂* transition between the primary sensory cortex (area 3) and the primary motor cortex (area 4), in line with histological evaluation of capillary density in humans (Duverney 1981). Moreover, we discovered a sharp ΔR₂* transition between area 3 and the secondary sensory area (Brodmann area 2). The estimated area boundary locations were supported by comparison of cortical area boundaries as defined in the M132 atlas (Fig. 4A,B) (N. T. Markov et al., 2014). The auditory cortex exhibited also relatively high ΔR₂*, however, the gradient-ridges were less distinguishable in this region. Multiple layer-specific gradient-ridges were also observed, albeit these were weaker in magnitude and more challenging to delineate.

Figure 4.

Figure 5.

Because terminals of myelinated axons often overlap with high oxidative metabolism (Horton, 1984; Rockoff et al., 2014), we also examined the association between ΔR₂* and T1w/T2w-FLAIR, an indirect proxy measure of cortical myelin density (Fig. 4A, C) (Glasser and Van Essen 2011, Autio et al., 2022). We found that M132 atlas parcellated ΔR₂* was positively correlated with intra-cortical T1w/T2w-FLAIR myelin (R = 0.49 ± 0.16), and negatively correlated with cortical thickness (R = -0.37 ± 0.10) (Fig. 4D).

3.3 Translaminar vascular volume variations link with neuroanatomical organization

Across the cortical areas and layers, average ΔR₂* profiles exhibited moderate variability (Fig. 5A, Supp. Fig. 5). To search for repetitive patterns in translaminar vascularity, we applied agglomerative clustering to concatenated group data. This analysis revealed distinct groups of vascularity arranged between eulaminate and agranular regions (Fig. 5B, C, D). A unique vascular profile was identified in V1, characterized by very dense vascularity and prominent peak density in EL4b (Fig. 5A, E). Average cluster profiles demonstrate that translaminar ΔR₂* was relatively high in isocortical areas, and low in agranular areas (Fig. 5E). The cluster boundaries (Fig. 4D) typically occurred in vicinity of the strong ΔR₂* gradients (Fig. 5B). Given the clustering between agranular and granular eulaminate cortices, we further corroborated whether the heterogeneous vascularization is associated with local microcircuit specialization of the cortex. For this objective, we used the designation of the cytoarchitectonic classification mapped onto M132 atlas (Burt et al., 2018; Hilgetag et al., 2016). This analysis confirmed that CBV, indeed, varies along the cytoarchitectonic types (Kendall’s tau τ = 0.69, p < 10^-5) (Fig. 5F). We also found a close association between baseline R₂* and cytoarchitecture (τ = 0.73, p < 10^-6). In the isocortex, the majority of the areas exhibited a distinctively high ΔR₂* in EL4 (Fig. 4A, E). The primary input layer, approximated as EL4a/b, exhibited systematically higher vascularity than in the primary output layer (p < 10^{- 25}), approximated as EL5a/b. In contrast, the majority of the agranular and dysgranular areas (cluster3; Fig. 3) exhibited weak laminar differentiation and a modest vascular density peak in EL1-2 (Fig. 5A).

Given neurons and receptors collectively constitute approximately 75-80% of the brain’s total energy budget (Howarth et al., 2012; Hyder et al., 2013), we next asked whether the regional variation in cerebrovascular network architecture (Fig. 1E) is associated with heterogeneous cellular and receptor densities. To address this question, we applied linear regression model, utilizing quantitative neuron (Collins et al., 2010) and receptor density maps (Froudist-Walsh, et al., 2023), to predict variation in ΔR₂* (Fig. 6A, B, D). This analysis revealed a positive correlation between CBV and neuron density in the middle cortical layers, where neuron density is typically highest, while revealing a negative correlation with receptor density in the superficial layers where synaptic density is highest (Fig. 6F). Additionally, we observed that baseline R₂* exhibited positive correlation with neuron density and negative correlation with receptor density (Fig. 6E).

Figure 6.

Because the Julich cortical area atlas covers only a section of the cerebral cortex, and the neuron density estimates are interpolated maps, we extended our analysis using the original Collins sample borders encompassing the entire cerebral cortex (Supp. Fig. 6A-C). This analysis reaffirmed the positive correlation with ΔR₂* (peak R = 0.80, p < 10^-20) and baseline R₂* (peak R = 0.77, p < 10^-30), yielding linear coefficients of ΔR₂* = 110 × 10³ neurons/s and R₂* = 28 × 10³ neurons/s (Supp. Fig. 6D-G). This suggests that the sensitivity of quantitative MRI in detecting neuronal loss is relatively low, and the introduction of the Ferumoxytol contrast agent has the potential to enhance this sensitivity by a factor of four.

Having established that vascular volume is associated with fundamental building units of cortical microcircuitry (Fig. 5F, Fig. 6E, F), our subsequent inquiry aimed to explore connection with interneurons that govern the neuroenergetics of local neural networks (Buzsáki et al., 2007). By utilizing the interneuron densities mapped onto M132 atlas (Burt et al., 2018), we identified a positive correlation between ΔR₂* and parvalbumin interneuron density (peak at EL5b, R = 0.72, p < 10^-6; Bonferroni corrected). In contrast, ΔR₂* showed negative correlation with the density of calretinin-expressing slow-spiking interneurons which preferentially target distal dendrites (peak EL1b, R = -0.61, p < 0.001; Bonferroni corrected). In conjunction, we found that ΔR₂* was also negatively correlated with dendritic tree size (R = -0.46, p < 0.01) (Supp. Fig. 7A, D, E) , the number of spines in L3 pyramidal cells (R = -0.37, p = 0.06), and also with R₂* (R = -0.69, p < 0.001) (Supp. Fig. 7B, D, F) (Elston 2002, Froudist-Walsh et al., 2023). These findings establish the intricate relationship between vascular density and the regulatory mechanisms governing diverse neural circuitry within the cerebral cortex.

4 Discussion

We present a noninvasive methodology to evaluate layer variations in vascular network architecture in the primate cerebral cortex. The quantitative cortical layer thickness adjusted ferumoxytol-weighted MRI enables exploration of systematic variations in cortical energy supply architecture and vessel-frequency informed image acquisition enables benchmarking penetrating vessel density measures relative to the anatomical ‘ground-truth’. These advances enabled us to unravel the systematic relation between vascularity and neurometabolic factors such as neuron and synaptic densities. Altogether, our study provides methodological and conceptual advancements in the field of cerebrovascular imaging.

4.1 Methodological considerations - vessel density informed MRI

To gain insights into the organization of cerebrovascular networks, it is important to critically sample the large irrigating arteries and draining veins while preserving adequate SNR in gray matter. However, MRI of the dense vascularity in the large and highly convoluted primate brain presents formidable challenges. Here, we used a combination of ferumoxytol contrast agent and laminar-resolution 3D GRE MRI to map cerebrovascular architecture in macaque monkeys. These methods allowed us to delineate large vessels and indirectly estimate translaminar variations in cortical microvasculature.

This methodology, however, has known limitations. First, gradient-echo imaging is more sensitized toward large pial vessels running along the cortical surface and large penetrating vessels, which could differentially bias the estimation of ΔR₂* across cortical layers (Fig. 2A, 2B) (Boxerman et al., 1995; Zhao et al., 2006). Additionally, vessel orientation relative to the B₀ direction introduce strong layer-specific biases in quantitative ΔR₂* measurements (Supp. Fig. 1C) (Lauwers et al., 2008; Ogawa et al., 1993; Viessmann et al., 2019). To address these concerns, we conducted necessary corrections for B₀-orientation, obtained parcel median values and regressed linear-trend thereby mitigating the effect of large-caliber vessels across ELs (Fig. 2C, Supp. Fig. 1). These analytical solutions yielded ΔR₂* V1 translaminar profiles that more closely resembled capillary rather than large-vessel volume profiles thus substantiating the validity of our methodology (Fig. 3F; Weber et al., 2008).

Ferumoxytol-weighted MRI of macaque cerebral cortex also enables the benchmarking of the methodological strengths and limitations to noninvasively measure vessel and vessel network length densities relative to the ‘ground-truth’. In macaque V1, large vessel length density (threshold at 8 μm diameter) is 138 mm/mm³ and mean area irrigated or drained by the vessels are about 0.26 and 0.4 mm² for arteries and veins, respectively (Weber et al., 2008). Combined, the total vessel density (= artery/0.26 + vein/0.40 mm²) is 6.3 vessels/mm² (Zheng et al., 1991; Weber et al., 2008), but see (Adams et al., 2015; Keller et al., 2011). Based on the former literature estimates, we hypothesized that isotropic voxel of 0.23 mm (19 voxels/mm²) may enable critical sampling of large vessels in accordance with sampling theorem (the sampling frequency equals to or is greater than twice the spatial frequency of the underlying anatomical detail in the image). In V1, we found an average vessel density of 2.2 ± 0.7 vessels/mm² (Fig. 2E) which corresponds to ≈30% of the ‘ground-truth’ estimate (Weber et al., 2008). Using cortical thickness as a reference, we estimate that the vessel length density is ≈8 mm/mm³ which corresponds to a modest ≈10% of the ‘ground-truth’ (Weber et al., 2008). The latter underestimate may be attributed to under-sampling of the branching arteriole and venule networks. Indeed, anatomical studies accounting for branching patterns have reported much higher vessel densities up to 30 vessels/mm² (Keller et al., 2011; Adams et al., 2015). Further investigations are warranted, taking into account critical sampling frequencies associated with vessel branching patterns (Duverney 1981) and achieving higher SNR through ultra-high B₀ MRI (Bolan et al., 2006; Harel et al., 2010; Kim et al., 2013). Such advancements hold promise for improving vessel classifications for veins and arteries and constructing detailed cortical surface maps of the vascular networks which may have diagnostic and neurosurgical utilities (Fig. 2A, B) (Iadecola, 2013; Qi and Roper, 2021; Sweeney et al., 2018).

4.2 Sharp transitions in microvasculature are indicative of cortical area- and layer-specific energy requirements

These methodological advances enabled us to unveil variations in vascular density within the primate cerebral cortex. Primary sensory cortices, known for their high energy demands (Kennedy et al., 1978), exhibit distinctive vascularization patterns (Fig. 1C, D, Fig. 3B). Notably, V1, area 3, auditory cortex, and also MT, all demonstrate elevated levels of CO enzymatic activity compared to surrounding cortical regions (Hackett et al., 1998; Horton, 1984; Huntley and Jones, 1991; Krubitzer et al., 2004; Matelli et al., 1985; Morel et al., 1993; Sincich et al., 2003). CO staining often reveals sharp transitions historically employed to delineate cortical area boundaries and modular features of the cortex. Our results corroborate the over three-decade-old, yet previously untested, hypothesis that some inter-areal boundaries may be determined by their microcapillary density using contrast-agent-weighted MRI (Zheng et al., 1991). Dense vascularity in these areas, the sharp gradient-ridges observed between surrounding areas and co-alignment with existing areal atlases further support this hypothesis (Fig. 4A, B).

Beyond the primary sensory areas, our observations extend to various smaller layer-specific vascular transitions (Fig. 5A, B). Specifically, the lateral intraparietal (LIP) area exhibits high CBV and gradient-ridges relative to ventral intraparietal (VIP) area and associative Brodmann area 7 (BA7). In EL5b, a strong gradient-ridge was observed distinguishing areas F4 and 44 from 45B and 8L/m. We also note weaker vascular transitions between areas such as 3a vs 3b, 5 vs anterior intraparietal (AIP) area, supplementary motor cortex (SII) vs insula, A1 vs medial belt and 46v vs 46d. However, our ability to confidently determine these borders is constrained by the presence of large vessels, as well as potential surface placement errors, and validating these areal boundaries would benefit from utilizing multimodal approaches.

In their work, Zheng and colleagues also proposed that modular features of the cortex, characterized by greater vascular density in the V1 CO blobs (42%) than interblobs, could be delineated using contrast-agent-weighted MRI (Zheng et al., 1991). Such a large vascular density difference should be well-within our contrast-to-noise and spatial sampling limitations. However, our results do not support this hypothesis, as we do not observe a distinct vascular peak at the spatial frequency of CO blobs (≈2.2 1/mm; Fig. 2D). Our results align with anatomical studies challenging the existence of high capillary density in the V1 blobs (Kernel et al., 2011, Adams et al., 2015).

The variation in vascular volume also has implications for statistical power in fMRI. For example, the pallidum exhibits lowest vascular volume (Fig. 1C) while in V1 EL4, the vascular volume is approximately twice as high (Fig. 3D). According to the classical single-tissue compartment model, CNR is optimized when TE is matched with T₂*. Consequently, there is no single optimal TE for cerebral blood volume weighted fMRI. Multi-echo EPI acquisition (Kuroiwa et al., 2014; Poser and Norris, 2009) may provide a more balanced comparison for statistical power across different brain regions and cortical layers.

4.3 The vascular network architecture is intricately connected to the neuroanatomical organization within cerebral cortex

Given the fivefold variability in neuron density across the cortex (Cohen et al., 2010), one might expect that there is a corresponding variation in cerebral blood flow (CBF) and CBV (Fig. 6) (Tsai et al., 2007). In the cerebral cortex, neurons account for a significant portion (≈80-90%) of energy demand, with most of this energy allocated to signaling (≈80%) and maintaining membrane resting potentials (≈20%) (Attwell and Laughlin, 2001; Howarth et al., 2012). Since firing frequency is modulatory and the neural networks utilize distributed coding, the maintenance of resting-state membrane potential determines the minimal energy budget and the lower-limit for cerebral perfusion. Based on neuronal variability and energy dedicated to maintaining surface potential, this suggest an approximate (4 × 20% ≈) 80% variation in CBF and a resultant 25% variation in CBV across the cortex, in line with Grubbs’ law (CBV = 0.80 × CBF^0.38) (Grubb et al., 1974). However, this is a simplified estimation, and a more comprehensive assessment would need to consider an aggregate of biophysical factors such as neuron types, neuron membrane surface area, firing rates, dendritic and synaptic densities (Fig. 6F-G), neurotransmitter recycling, and other cell types (Kageyama 1982; Elston and Rose 1997; Perge et al., 2009; Harris et al., 2012). Indeed, the majority of the mitochondria reside in the dendrites and synaptic transmission is widely acknowledged to drive the majority of the energy consumption and blood flow (Wong-Riley, 1989; Attwell et al., 2001).

Contrary to our initial expectations, we observed a relatively smaller CBV in regions and layers with high receptor density (Fig. 6B, D, F). This unexpected relationship extends to other factors, such as number of spines (putative excitatory inputs) and dendrite tree size across the entire cerebral cortex (Supp. Fig. 7) (Froudist-Walsh et al., 2023, Elston 2007). These results align with the work of Weber and colleagues, who reported a similar negative correlation between vascular length density and synaptic density, as well as a positive correlation with neuron density in macaque V1 across cortical layers (Weber et al., 2008). This relation is also compatible with the opposing relation between CBV and GABAergic (GABA, γ-aminobutyric acid) interneuron subtypes (Fig. 7G): parvalbumin-expressing fast-firing interneurons target perisomatic parts of pyramidal neurons, whereas calretinin-expressing slow-firing interneurons target distal-dendrites. The interneurons are also well-positioned to play an important role in integrating activity of large numbers of excitatory principal cells and translating this into neurovascular regulation of local microcirculation via subcortical pathways (Cauli et al., 2004). Another perspective on our results considers that a relative smaller faction of synapses may be simultaneously active in larger neurons, and neurons with large dendrites may exhibit lower excitability, supporting pattern separation necessary for higher-level functions (Chavlis et al., 2017; Hawkins and Ahmad, 2016). To comprehensively understand the factors contributing to the vascular organization of the brain, experimental disentanglement through integrated analysis of laminar cell types and receptor densities is needed (Hayashi et al., 2021, Froudist-Walsh et al., 2023). Moreover, employing more advanced statistical modeling, including considerations for synapse-neuron interactions, may be important for refined evaluations.

Additionally, our investigation also uncovered an overlap between vascular volume and myelin (Fig. 3A, C). Myelin may indirectly contribute to increased energy consumption in the cortical gray matter by facilitating higher frequency firing in comparison to unmyelinated axons (Perge et al., 2012; Saab et al., 2016). A large fraction of cortical myelin enwraps the axons of parvalbumin-expressing fast-spiking interneurons (Stedehouder et al., 2017). Although interneurons represent a minority of the neuronal population (10-15%), the parvalbumin-expressing interneurons’ ability to sustain high-frequency gamma oscillations (30-100 Hz) may match the sparse firing of the majority principal cells (Buzsáki et al., 2007). Indeed, parvalbumin-expressing interneurons exhibit higher mitochondrial volume compared to other cells in the brain and a 3-fold higher CO activity than principal neurons (Gulyás et al., 2006; Kageyama and Wong-Riley, 1982; Kann et al., 2014; Nie and Wong-Riley, 1995). Thus, the high metabolic load of parvalbumin-expressing interneurons makes them potentially vulnerable to failures in the vascular network due to aging, Alzheimer’s disease as well as stroke (Kann, 2016).

Given the distinct pre- and postsynaptic metabolic requirements, heterogenous translaminar vascularization may also indicate distinct cortical layers each characterized by anatomically and physiologically distinct feedforward and feedback pathways (Bastos et al., 2015; Borowsky and Collins, 1989; Kageyama and Wong-Riley, 1982; Takahata, 2016). A prime example is V1 where the primary input layer 4 (L4) has more dense vascularization and 50% higher CO activity in comparison with primary output L5 (Fig. 3F) (Keller et al., 2011; Livingstone and Hubel, 1982). This elevated energy demand in L4 may arise, in part, from spontaneous supra-threshold gamma-frequency oscillations between the retina→lateral geniculate nucleus→L4 (Castelo-Branco et al., 1998) along with recurrent amplification of local and distant inputs (Douglas and Martin, 2007; Shu et al., 2003). When viewed in terms of information flow, CBV appear to decrease along the canonical circuit pathway (e.g., L4→L2/3→L5) in the primary visual cortex (Douglas and Martin, 2007) and as one ascends the hierarchy (e.g., V1→V2→V3&4→MT→7A) from primary sensory areas (Fig. 3F, Supp. Fig. 8) (Felleman and Van Essen, 1991; Nikola T. Markov et al., 2014; N. T. Markov et al., 2014). In the agranular and dysgranular regions, characterized by a lack of distinct L4, the translaminar CBV profiles did not exhibit a distinct peak at around EL4 nor signatures of canonical circuitry (Fig. 5). These results demonstrate a greater allocation of the energy budget to early stages of feedforward processing in primary cortical areas characterized by strong sensory inputs, as opposed to higher-level cortical areas characterized by high synaptic densities supporting cognitive and behavioral functions (Fig. 6B) (Yokoyama et al., 2021).

The anatomical uniformity of the primate neocortex is thought to reflect extensive replication of a few specialized microcircuits varying along the brain’s hierarchical organization (Douglas and Martin, 2007; Hilgetag et al., 2016; Markram et al., 2004). Analogous to cortical circuitry, our study reveals the large-scale replication of translaminar vascular network motifs in primates (Fig. 4B, C, D). Since the cerebrovascular system evolved to support the high energy demands of neural information processing, this raises questions about whether the large-scale replication of anatomical and vascular circuits is evolutionarily coupled (Carmeliet and Tessier-Lavigne, 2005). In mice, comprehensive analysis of the cerebrovascular system has revealed distinct translaminar types between sensory and motor-integrative areas (Kirst et al., 2020). In macaque, we found the strongest distinction between isocortex and allocortex and its adjacent regions (Fig. 6B). The species difference may reflect evolutionary expansion and emergence of new cortical layers. For instance, in mice the primary somatosensory cortex exhibits highest vascularization (Kirst et al., 2020) whereas our results show that in macaque the highest vascularization is in the V1 (Fig. 5A, E) (Duverney 1981). According to the theory that sensory systems, behavior, and habitat choice are all influenced by evolutionary processes (Endler, 1992; Ikeda et al., 2023), this may reflect an evolutionary adaptation to an environment in which the visual landscape is an ecologically more important sensory domain in primates.

Author contributions

Conceptualization: JAA; Methodology: JAA; Software: JAA, IK; Formal analysis: JAA, IK, TI; Investigation: JAA, TO, YM, MO; Writing - original draft: JAA; Writing - review & editing: JAA, IK, MFG, DVCE, TH; Visualization: JAA, IK, TI; Resources: JAA, YU, TH; Project administration: TO; Funding acquisition: JAA, TH.

Acknowledgements

The authors appreciate discussions and technical contributions from Akiko Uematsu, Timothy S. Coalson, Katsutoshi Murata, and Reiko Kobayashi. This research is partially supported by JSPS KAKENHI Grant Number (JP20K15945, JAA), by the program for Brain/MINDS and Brain/MINDS-beyond from Japan Agency for Medical Research and development, AMED (JP18dm037006, TH) and by NIH R01MH60974 (DCVE, MFG). The authors have no conflicts of interest to declare.