Effects of phase regression on high-resolution functional MRI of the primary visual cortex

High-resolution functional MRI studies have become a powerful tool to non-invasively probe the sub-millimeter functional organization of the human cortex. Advances in MR hardware, imaging techniques and sophisticated post-processing methods have allowed high resolution fMRI to be used in both the clinical and academic neurosciences. However, consensus within the community regarding the use of gradient echo (GE) or spin echo (SE) based acquisition remains largely divided. On one hand, GE provides a high temporal signal-to-noise ratio (tSNR) technique sensitive to both the macro- and micro-vascular signal while SE based methods are more specific to microvasculature but suffer from lower tSNR and specific absorption rate limitations, especially at high field and short repetition times. Fortunately, the phase of the GE-EPI signal is sensitive to vessel size and this provides a potential avenue to reduce the macrovascular weighting of the signal (phase regression, Menon 2002). In order to determine the efficacy of this technique at high-resolution, phase regression was applied to GE-EPI timeseries and compared to SE-EPI to determine if GE-EPI's specificity to the microvascular compartment improved. To do this, functional data was collected from seven subjects on a neuro-optimized 7 T system at 800 μm isotropic resolution with both GE-EPI and SE-EPI while observing an 8Hz contrast reversing checkerboard. Phase data from the GE-EPI was used to create a microvasculature-weighted time series (GE-EPI-PR). Anatomical imaging (MP2RAGE) was also collected to allow for surface segmentation so that the functional results could be projected onto a surface. A multi-echo gradient echo sequence was collected and used to identify venous vasculature. The GE-EPI-PR surface activation maps showed a high qualitative similarity with SE-EPI and also produced laminar activity profiles similar to SE-EPI. When the GE-EPI and GE-EPI-PR distributions were compared to SE-EPI it was shown that GE-EPI-PR had similar distribution characteristics to SE-EPI (p<0.05) across the top 60% of cortex. Furthermore, it was shown that GE-EPI-PR has a higher contrast-to-noise ratio (0.5 ± 0.2, mean ± std. dev. across layers) than SE-EPI (0.27 ± 0.07) demonstrating the technique has higher sensitivity than SE-EPI. Taken together this evidence suggests phase regression is a useful method in low SNR studies such as high-resolution fMRI.


Introduction
The human cortex is organized into functionally distinct layers parallel to the pial surface and, in select areas, columns perpendicular to the surface. Cortical layers and columns are key functional units in understanding how the brain is organized. Similarly positioned layers perform similar tasks across different parts of the brain ( Douglas and Martin, 2004 ). Specifically, the neuronal inputs and outputs are contained in different cortical layers. Measurement of interactions between cortical features such as these could allow for a deeper understanding of intra-cortical and inter-cortical communication. Historically, investigat- ( Huber et al., 2017 ;Koopmans et al., 2010 ;Lawrence et al., 2019 ;Olman et al., 2012 ;Ress et al., 2007 ) and columns ( Cheng et al., 2001 ;Feinberg et al., 2018 ;Menon et al., 1997 ;Moerel et al., 2018b ) across the human brain.
High-resolution fMRI using the Blood Oxygenation Level Dependent (BOLD) technique struggles with various acquisition challenges that require consideration prior to data collection. These include lower SNR as resolution increases, macrovascular bias, and specific absorption rate (SAR) constraints at ultra-high fields. Several different acquisition approaches have been used to perform high-resolution fMRI although the two most commonly used are gradient echo EPI (GE-EPI) and spin echo EPI (SE-EPI). Conventional GE-EPI produces the largest signal changes but is not specific to microvasculature ( Budde et al., 2014 ;Koopmans and Yacoub, 2019 ;Menon, 2012 ;Rua et al., 2017 ;Yacoub et al., 2007 ). Unlike GE-EPI, SE-EPI is more specific due to the use of a refocusing pulse in order to suppress the macrovascular signal. Unfortunately, this sequence suffers from SNR and SAR penalties making it a less sensitive technique overall. Comparisons between the sensitivity and specificity of these techniques shows that GE-EPI has 1.29 times higher percent signal change in grey matter than SE-EPI and 5.33 times higher percent signal change in vessels at 7T ( Yacoub et al., 2005 ). These advantages have made GE-EPI the overwhelming choice for high resolution fMRI studies.
Several alternatives to GE-EPI and SE-EPI have been investigated such as vascular space occupancy (VASO) ( Huber et al., 2018 ;Lu et al., 2013 ), balanced steady state free procession (bSSFP) ( Báez-Yánez et al., 2017 ) and gradient and spin-echo imaging (GRASE) ( De Martino et al., 2013 ;Moerel et al., 2018b ). VASO focuses on imaging changes in cerebral blood volume which results in more specificity to the microvasculature ( Lu et al., 2013 ). bSSFP shows 2 -like weighting and SNR efficiency but is limited to a small slab to avoid excessively long acquisition times. GRASE reduces macrovascular signal by placing refocusing pulses throughout a GE-EPI sequence which lowers * 2 weighting ( De Martino et al., 2013 ) but also limits coverage to a specific region of interest to avoid reintroducing these * 2 effects. However, GE-EPI remains the most commonly used fMRI technique today due to its robustness and well understood signal properties.
Previous work by many groups has attempted to reduce the vascular bias from high-resolution GE-EPI while maintaining sensitivity. These methods include optical imaging to identify larger vessels ( Chen et al., 2013 ), susceptibility weighted imaging to identify veins ( Kay et al., 2019 ;Koopmans et al., 2010 ), removing venous bias through deconvolution with a vascular PSF ( Markuerkiaga et al., 2016 ), looking at the initial dip of the BOLD response ( Siero et al., 2013 ), contrast subtraction ( Cheng et al., 2001 ;Moerel et al., 2018a ) and removing the higher cortical layers where such veins are present from further analysis ( Polimeni et al., 2010 ). These techniques all rely on knowledge and/or assumptions of the vein's locations in the GE-EPI images which can require additional acquisitions or signal modelling, and this can complicate their use in high resolution fMRI. Another very recent approach that utilizes the temporal lag between microvascular BOLD signals and macrovascular BOLD signals shows promise in an initial report ( Kay et al., 2020 ). This paper proposes the use of the phase of the high-resolution GE-EPI images to estimate BOLD signal caused by large vessels and subtract it from the magnitude data. This data-driven method reduces macrovascular bias without using additional venous identification. fMRI phase is an intrinsic part of a conventional GE-EPI acquisition but is usually not reconstructed and saved as part of the fMRI pipeline. Phase regression has previously been used at low resolutions to reduce large vessel contributions in the magnitude images ( Caballero-Gaudes and Reynolds, 2017 ;Menon, 2002 ). This technique relies on the fact that although magnitude signal will contain BOLD changes from both large and small vasculature, phase data will primarily contain BOLD changes from large vessels ( Menon, 2002 ).
Some discussion of what constitutes a large vessel with respect to this technique is necessary. Cortical veins can be divided into three groups: pial veins ( > 280 m ), run along the cortical surface; intracortical penetrating veins (80-170 m ), run perpendicular to the cortical surface; and smaller intracortical tangential veins, which run at different depths parallel to the layers of the cortex ( Duvernoy et al., 1981 ). For the current experiments, it is unlikely that useful phase information can be obtained from vessels smaller than 150 m in diameter ( Klassen and Menon, 2005 ). Additionally, all phase related BOLD changes will increase in amplitude as the vessel size increases ( Ogawa et al., 1993 ) so larger vessels will dominate the phase time course. Thus, for the purposes of this paper we define the macrovasculature as vessels large enough to produce a detectable BOLD phase change, which will primarily be pial vessels and a few of the largest intracortical veins.
At low resolutions, BOLD related phase changes are primarily due to the intravascular BOLD signal ( Menon, 2002 ;Ogawa et al., 1993 ). The extravascular phase signal for the macrovasculature will be negligible at low resolutions due to the symmetric extravascular frequency profile. Therefore, it can be assumed any voxel with a high correlation between magnitude and phase contains signal from the intravascular component of macrovasculature. This assumption has yet to be tested in voxels near the size of pial vessels on the cortical surface which this paper seeks to investigate. Extravascular frequency shifts could produce a phase change in a sufficiently small voxel when the symmetry assumption is violated ( Vu and Gallant, 2015 ). This would result in phase changes and suppression of extravascular and intravascular BOLD signal, improving the reduction in macrovascular bias for high resolution data.
This paper investigates phase regression of high-resolution GE-EPI functional time series data as a method to reduce macrovascular bias. Laminar structures are evaluated in GE-EPI and SE-EPI functional acquisitions and compared with GE-EPI-PR (GE-EPI with phase regression) data. This paper examines the surface activation maps of the GE-EPI, SE-EPI and GE-EPI-PR as well as their activation distributions and contrastto-noise ratio (CNR). Furthermore, the laminar profiles of GE-EPI, SE-EPI and GE-EPI-PR are compared to determine the effect of phase regression on the laminar profile proximal to and distal from a vessel.

Imaging protocol
Data from seven subjects was acquired (5 male, 2 female, 25.8 ± 4.0 years). Each individual was positioned supine on the MRI bed with a mirror placed over the eyes for viewing a rear-projection screen 28 cm away producing a left-right visual angle of 27.5 degrees. Foam cushions were placed around the head for comfort and immobilization as well as medical tape across the forehead for haptic feedback to reduce motion. Informed consent of all participants was collected in accordance with and approved by the Human Subjects Research Ethics Board at the University of Western Ontario.
Imaging was performed using a 680 mm neuro-optimized 7 T MRI (Siemens Magnetom Step 2.3, Erlangen, Germany) equipped with an AC84 Mark II head gradient coil. An 8-channel Tx, 32-channel Rx radiofrequency coil optimized for occipital-parietal imaging with no visual obstruction over the face was chosen for data collection ( Gilbert et al., 2017 ). The actual flip-angle imaging (AFI) technique ( Yarnykh, 2007 ), augmented with an RF and gradient spoiling scheme ( Nehrke, 2009 ), was used to map the transmit field. In addition, 8 images with Fourier B 1 + encoding were acquired to map relative transmit profiles. RF shimming was subsequently performed, which consisted of setting the phase and magnitude of each transmit channel using a leastsquares optimization that balanced transmit efficiency and uniformity ( Curtis et al., 2012 ). The B 1 + shim solution was optimized over the region of interest relevant to the BOLD measurements. In order to ensure phase was not inappropriately filtered zero-filling partial Fourier was used. The scanning protocol consisted of GE-EPI, SE-EPI, a multi-echo gradient echo sequence for venous localization and an MP2RAGE sequence with high gray-white contrast to extract functional surfaces. Parameters for all imaging sequences can be found in Table 1 . Phase reconstruction of the GE-EPI data was completed using the fitted SVD sensitivities method ( Stanley et al., 2018 ) to prevent destructive interference. Coil sensitivity estimates are obtained by utilizing the multi-image prescan collected for B 1 + mapping and performing a singular value decomposition. These estimates are fit to a functional basis in order to allow for their interpolation to other fields of view during the imaging session. The fitted SVD derived coil sensitivities are multiplied with the uncombined coil data to align it prior to a complex sum. This multiplication and complex sum were completed as part of the Siemens reconstruction chain of the CMRR multiband sequence through insertion of a custom functor. Maxwell correction is turned off to prevent any spatial translation differences between the magnitude and phase images after combination. This method of combination is memory efficient due to the low resolution prescan and through the use of custom functors during the normal Siemens reconstruction pipeline, typically a few hundred megabytes. For non-Siemens sites, phase combination can use the fitted SVD sensitivities method or additional methods such as COMPOSER ( Bollmann et al., 2018 ;Robinson et al., 2017 ), the virtual reference coil method ( Parker et al., 2014 ) or the voxelwise SVD method ( Walsh et al., 2000 ). We recommend an online combination for computational efficiency when dealing with large datasets ( Hansen and Sørensen, 2013 ). A sum-of squares combination from all Rx channels was used to reconstruct the magnitude data.
The magnitude of the complex point spread function of the EPI acquisitions was calculated for the phase encode direction, by simulating a real uniform EPI echo train in k -space in order to estimate the effective image resolution ( Table 1 ). The effects of acceleration were applied, lines were skipped from GRAPPA, * 2 and/or 2 signal decay added and the echo was zero filled consistent with the partial Fourier technique used during acquisition. The phase encode profile was inverse Fourier transformed and the full width half maximum of the magnitude was measured and reported as the point spread function (PSF) of each EPI sequence. The protocols for the SE-EPI and GE-EPI acquisitions were matched as closely as possible resulting in the SE-EPI having a slightly wider PSF than the GE-EPI. This estimate was performed to compare the two EPI acquisitions and may not be entirely representative of the true resolution ( Chaimow and Shmuel, 2017 ).

Functional Stimulus
The visual stimulus was an 8 Hz contrast reversing checkerboard created using Pyschtoolbox (3.0.11 ( Brainard, 1997;Kleiner et al., 2007 )) in Matlab (2015a( MATLAB, 2015). This was delivered in a restactivation paradigm of 15 s off, 15 s on lasting for 8 repetitions and ending on a rest block. In order to help maintain attention, a button press task was used where participants were asked to respond when a central fixation cross changed orientation by 45 degrees. Three runs were acquired for each participant and for each sequence type: GE-EPI and SE-EPI.

Data preprocessing 2.2.1. Data analysis software
All imaging data was converted to the brain imaging data structure (BIDS) format using in-house conversion tools wrapped around heudiconv (Heuristic Dicom Conversion). Analysis was completed using Nipype pipelines (1.1.8 ( Gorgolewski et al., 2017 )) including several custom interfaces for phase analysis. The in-house software used can be found at: https://github.com/ostanley/phaseprep . An overview of reconstruction and preprocessing is available in Fig. 1 .

Functional data preprocessing
Both GE-EPI and SE-EPI underwent the same magnitude preprocessing. Each functional imaging run (GE-EPI or SE-EPI) was motion corrected and aligned to the first volume of the first run using AFNI (18.1.24 ( Cox, 1996 )). Brain extraction was completed on the same first functional volume using FSL's BET tool and the mask was applied to all functional runs (FSL version 5.0.10). The preprocessed magnitude data was then used as the magnitude input for phase regression. After phase regression but prior to general linear modelling all data was scaled to a mean of 10000 and high-pass filtered with a window of 100 s (identical to conventional FSL FEAT preprocessing ( Jenkinson et al., 2012 )).
Preprocessing of the phase data was performed using the in-house Nipype workflow, preproc_phase_wf.py . The workflow consists of conversion of the magnitude and phase data to real and imaginary. Motion correction was performed in real and imaginary space since it is spatially smooth and interpolatable. The transformations from the magnitude images were applied to the real and imaginary data and then converted back to magnitude and phase images. The phase data was further processed by performing first volume subtraction, temporal unwrapping, and linear detrending. Voxelwise detrending was performed to remove systematic linear frequency drift and B 0 field variations over time. The result was a motion corrected phase timeseries, free of temporal and spatial wraps which also accounts for linear system and B 0 field variation.

Phase regression
Previous work on phase regression has shown BOLD related phase changes will correlate with an associated BOLD related change in magnitude and this can be used to estimate signal originating from macrovasculature ( Menon, 2002 ). This estimated signal can then be subtracted from the magnitude signal in order to reduce signal from large vessels. This method relies on two assumptions: (1) the temporal correlation of magnitude and phase is different in the microvasculature than in macrovasculature, which prevents complete suppression of the tissue signal in a voxel with a large vessel ( Jorge et al., 2018 ), (2) that all large vessels produce a phase change, which may not be true for vessels at certain orientations ( Ogawa et al., 1993 ).
Phase regression was performed using voxelwise orthogonal distance regression (ODR) in the in-house Nipype gadget, PhaseFitODR.py . ODR uses residuals perpendicular to the line of best fit and was selected due to the noise present in both magnitude and phase data. The regression was completed to solve the following equation:

= +
where M is the magnitude signal, is the phase signal and A and B are the fit coefficients. ODR requires inputs to estimate error ellipses prior to fitting. In order to estimate these errors for magnitude and phase, each time course is high pass filtered at 0.15 Hz (above the task frequency). The temporal standard deviation of these filtered signals was then used as the inputs to the ODR for uncertainties and the unfiltered signals are used as input to the fits. ODR is then used to estimate the component of the signal with high magnitude and phase correlation which is assumed to be macrovascular signal. Subtraction of this estimated macrovascular signal results in a signal weighted towards microvasculature ( or GE-EPI-PR): The effect of this on the timeseries of both a tissue and venous voxel is shown in Fig. 2 . Both the estimated macrovasculature and GE-EPI-PR timeseries underwent the same preprocessing steps as the GE-EPI and SE-EPI time courses (scaling to 10000 and high-pass filtering with a filter window of 100 s).

Functional Data Fitting
In order to perform physiological noise correction, regressors were created using anatomical Compcor ( Behzadi et al., 2007 ) from each subject's eroded white matter mask transformed to native EPI space. Compcor masks were generated with two erosions in fslmaths using a 3 × 3 × 3 kernel. Six compcor regressors and the six motion regressors were included as regressors of no interest in the GLM to account for noise caused by physiology and/or motion. All four time series (GE-EPI, SE-EPI, GE-EPI-PR, and the estimated macrovascular timeseries) were analyzed using the FSL film_gls tool. The output was converted to % BOLD signal change through normalization to the mean intensity of the timeseries. CNR was calculated by dividing the amplitude of the fit signal by the standard deviation of the residuals ( Welvaert and Rosseel, 2013 ).

Structural image analysis
The MP2RAGE image was run through the Freesurfer high resolution recon-all pipeline with two modifications (6.0.0 Fischl et al., 1999 )). First, the Talairach registration was turned off as the structural image was limited to the posterior part of the brain (notalariach , due to coil construction). Second, the corpus collosum and pons were manually seeded to ensure proper initialization. The cortical segmentations were manually inspected for agreement with the borders in the region of interest and brain mask corrections were performed, if necessary. The white matter surfaces were equidistantly expanded to allow for depth analysis using Freesurfer's mris_expand tool ( Kay et al., 2019 ). All results were calculated at 10% cortical depth intervals from 0 (pial surface) to 100% (white matter surface). This does not represent the expected anatomical distribution of the cortical layers but allows for investigation across surfaces and depths. All results were presented across the flattened surfaces by sampling voxel results onto the vertices that make up the surface at each depth.
In order to restrict analysis to a reasonable area an occipital patch was cut from the rest of the cortex and flattened using Freesurfer's mris_flatten tool. A patch over the calcarine sulcus was selected by manual delineation from the white matter curvature and the field of view of all acquisitions projected onto the occipital flat patch ( Fig. 3 ). This selected patch is expected to be within primary visual area V1, but more importantly this cortical patch would contain a variety of activation levels and contain vessels required for this investigation. This area was 1000 ± 300 mm 2 per subject leading to on average 1600 voxels analyzed for each hemisphere's surface patch.

Venous maps
As a simple, robust method for identifying venous vasculature, the product of * 2 and the initial magnetization ( 0 ) from multi-echo susceptibility weighted imaging was used. This was done because * 2 is a physical value and is therefore expected to be consistent across subjects except for the presence of a low frequency background field ( Fernández-Seara and Wehrli, 2000 ). In addition, 0 can be calculated from the same fit and does not require free parameters. In order to calculate these parameters, the multi-echo GRE data was run through the qsm_sstv pipeline ( https://github.com/AlanKuurstra/qsm _ sstv/releases , 1.0.0). Briefly, this BIDS app extracted the brain from the multi-echo data and performed complex fitting to calculate * 2 , frequency, and 0 . These maps were then registered to the T 1 image using the same methods as the EPIs (described below) to transfer them to surface space for functional analysis.

Registration to structural data
Registration of the functional maps to structural space was completed using ANTs (2.2.0 ( Avants et al., 2011 )). After initialization using the center of mass, a rigid transform was completed followed by two affine transformations, one general and one targeted at the region of interest. Mutual information was used as the target metric and all interpolation was completed using order 3 splines. In most laminar studies it is common to bring the anatomic surfaces into functional space. This was not done in this study due to the different fields of view of the GE-EPI and SE-EPI. As an alternative all transformations prior to the GLM were kept to a minimum (one spatial transform per volume) and then performed one single transform of each result to T 1 space. These results were then transferred onto the cortical surface ROI using mri_vol2surf and allowed for surface comparison between both pulse sequences. The T 1 transform and sampling to surface space will result in some effective blurring of the data, however these effects are minimized by doing phase regression and GLMs prior to transformation. Fig. 3. Visual demonstration of the surface pipeline. 1: Quality assurance figures for registration and segmentation. FSL Fast segmentation of a) GE-EPI and b) SE-EPI overlaid on the T 1 weighted MP2RAGE image, c) Freesurfer surfaces overlaid on the T 1 weighted MP2RAGE image. 2: Surface processing pipeline. First the MP2RAGE is used to generate surfaces (a). From the generated surfaces (b) a calcarine patch is extracted (c) and flattened (d). The ROI (blue area) for analysis is then manually delineated using tksurfer and the curvature map as well as slice coverage from the functional scans (e).

Fig. 4. Surface visualization demonstration. a)
The vertices for a cortical patch plotted over a grey background (used to calculate all laminar profiles and distributions). b) A triangular mesh to make the data more contiguous and easier to visualize when examining surface activation maps qualitatively.

Surface visualization
Surface activation maps become distorted during flattening resulting in uneven vertex placement across flat space. The vertices for each layer were converted to a three-dimensional mesh and the laminar surface activation maps were plotted as a triangular mesh in order to reduce this effect. By doing so, it becomes easier to view as it does not involve varying amounts of dead space. The effect this has on visualization is displayed in Fig. 4 . All laminar profiles and signal distributions were calculated across vertices and did not use an interpolation.

Vessel segmentation
Manual segmentation was performed in order to delineate visible vessels from tissue. The product of the * 2 and 0 surface maps was selected for manual segmentation as it showed reduced noise compared to the * 2 surface map. Each laminar surface map was manually segmented for every subject and every cortical depth. Hyperintensities were outlined as polygons on top of the mesh using matplotlib (3.1.0 ( Hunter, 2007 )). All vertices in these hyperintense region polygons were then labelled as vessels. To control for bias, no indication of cortical depth or subject was given when each map was presented, each map was presented with an identical colour bar, and the maps were presented in randomized order. After manual segmentation was completed, two forms of continuity clustering were used. First, vertices which were labelled as a vessel across two adjacent depths were included in the final vessel map to include penetrating vessels. Second, marked vertices greater than 0.3 mm away from another marked vertex were excluded, this provides for the possibility of a vessel running along the surface of a single layer. These thresholds were both applied to improve vessel continuity and reduce the presence of single noisy vertices causing mis-labelling. An overview of this process can be seen in Fig. 5 .

Laminar profile generation
Laminar profiles were plotted by averaging across all vertices in the calcarine patch of interest for each of the nine depths from 10% to 90%. Vertices were also classified as proximal to or distal from a vessel based on their minimum Euclidean distance to a vessel vertex thresholded at 2.4 mm. At this distance, an activation based frequency shift of 24 Hz is expected, compared to 220 Hz at the surface of a 0.8 mm vessel (calculated from ( Vu and Gallant, 2015 )). This was considered sufficiently out of the influence of large veins for this study.

Results
The temporal SNR in the field of view of interest is uniform ( Fig. 6 ). A poor B 1 + shim in one subject's hemisphere was observed and verified on the actual flip angle map. This hemisphere was excluded from the group  metrics reported below and from all further analysis. Temporal SNR across the cortical ribbon of all subjects was 10.2 ± 1.2 (mean ± standard deviation) for the GE-EPI and 8.36 ± 0.83 for the SE-EPI data. The tSNR of the GE-EPI and SE-EPI was significantly different in a Welch's t-test ( p = 0.015) and this is an important note for later surface activation map comparison. Finally, the phase standard deviation of the timeseries was 0.21 ± 0.12 radians across all subjects.
To investigate the changes in the laminar surface activation maps due to phase regression the BOLD % signal change was projected onto surfaces at various cortical depths (an example subject is shown in Fig. 7 ). Two additional examples are shown in supplemental information. The equidistantly projected data shows surface veins that are clearly visible in the higher layers of cortex (towards the pial surface). This is to be expected even in the SE-EPI case as the purely T 2 weighting only applies for the central measurement of k-space. The lower tSNR in the SE-EPI does affect the laminar surface activation maps as they appear noisier, but it is still clearly less sensitive to large vessels. The hyperintense venous regions in the GE-EPI data exhibit the largest signal suppression after phase regression compared to surrounding areas. The spatial distribution of GE-EPI-PR appears to more closely match the SE-EPI case. In order to validate the areas of high activation in the GE-EPI are truly large vessels the data was examined in conjunction with the structurally derived, vessel sensitive * 2 and M 0 data ( Fig. 8 ). Fig. 8 shows the * 2 and M 0 product maps projected onto the cortical surfaces indicating the vessel locations from independent anatomy without the functional data. Also shown are the two metrics that illustrate the performance of the phase regression. These are the correlation between the fitted phase and magnitude (R 2 ), and the activation resulting from the fitted phase time series (estimated macrovascular activation). Sup-plemental information also provides two additional examples. The * 2 and M 0 map shows vessel like structures where the largest reduction in GE-EPI-PR % signal change occurred. The estimated macrovascular activation also shows areas of hyperintensity at these locations indicating that phase regression is suppressing venous signal. This can be further quantitatively investigated through examining the distributions of the different functional imaging methods.
SE-EPI exhibits specificity to the microvasculature making it an appealing method for BOLD imaging of cortical substructures like columns ( Yacoub et al., 2007 ). By comparing the GE-EPI and GE-EPI-PR distributions to SE-EPI, a direct comparison to a microvascular control can be evaluated ( Fig. 9 ). A group of Kolmogorov-Smirnov tests with Bonferroni multiple comparisons correction were used in order to investigate similarities between distributions of the imaging methods. These tests show the distributions are all significantly different ( p < 0.05) except the distributions of GE-EPI-PR and SE-EPI from depths of 10 to 60% demonstrating that the distribution of GE-EPI-PR is more characteristic of SE-EPI than GE-EPI in the higher layers of cortex. This supports the hypothesis that GE-EPI-PR is suppressing pial vessel signal and producing a SE-EPI-like activation map.
One concern in using phase regression in fMRI processing is the reduced contrast-to-noise ratio. GE-EPI-PR shows signal suppression relative to GE-EPI ( Fig. 10 ). CNR was calculated by dividing the amplitude of the activation by the standard deviation of the residuals ( Welvaert and Rosseel, 2013 ). This was done to investigate whether phase regression is introducing any noise through the fit subtraction process which could potentially reduce the method's efficacy. The average CNR across layers of the GE-EPI data is 0.9 ± 0.3 (mean ± std dev. across layers), for GE-EPI-PR the CNR is 0.5 ± 0.2 and finally SE-EPI has a CNR of  0.27 ± 0.07. This means the CNR of the SE-EPI data is only 30% of the GE-EPI data compared to the CNR of GE-EPI-PR which is 60% of the GE-EPI data. This shows that the phase regression method reduces GE-EPI CNR as expected, however it has higher CNR than SE-EPI. These findings suggest that although some noise may be introduced GE-EPI-PR is still an advantageous method to use over SE-EPI as it has more statistical power.
Large venous vessels exhibit both an intra-and extra-vascular BOLD response. Removal of the extravascular bloom is an important component in reducing the signal bias from these large draining veins. The laminar profiles distal from all vessel vertices were examined in order to determine if this extravascular bloom was being successfully reduced. Two bins of vertices were created, one proximal to and one distal from a vessel vertex. Fig. 11 shows laminar profiles over all vertices as well as for vertices proximal ( < 2.4 mm) and distal to a vein ( > 2.4 mm). The GE-EPI-PR data shows activation similar in profile to the GE-EPI data but with a lower percent signal change when distal from vasculature. The difference between the GE-EPI-PR and GE-EPI is most prominent in the higher depth vertices proximal to veins where the GE-EPI-PR laminar profiles have a lower slope than GE-EPI. This would indicate that the phase regression is reducing contribution from pial veins to a higher degree than tissue as we expect. Fig. 9. Test statistic of the two-sided Kolmogorov-Smirnov test between distributions as a function of depth. The dashed line represents the significance threshold ( p < 0.05) after Bonferroni comparisons correction across all depths. Fig. 10. Laminar CNR profiles across subjects. Error bars represent the standard error of the mean across subjects. All vertices were used for this calculation.

Discussion
In this study, we investigated the use of phase regression on highresolution GE-EPI data to assess the feasibility of the technique for use in intracortical BOLD fMRI at the laminar and/or columnar level. GE-EPI is an attractive sequence for use in high-resolution fMRI as it has an inherently higher contrast-to-noise ratio per unit time (sensitivity) compared to other popular intracortical fMRI approaches (i.e. SE-EPI, VASO, GRASE, or bSSFP), as well as lower SAR requirements and higher spatial coverage making it easier to achieve high temporal resolutions and shorter imaging times ( Beckett et al., 2020 ;Koopmans and Yacoub, 2019 ). Unfortunately, GE-EPI suffers from macrovascular contamination leading to low specificity to the capillary bed microvasculature ( Menon, 2012 ). Phase regression of the GE-EPI images was investigated to determine if the specificity could be improved without sacrificing microvascular sensitivity, improving GE-EPI utility in high resolution studies. GE-EPI and GE-EPI-PR were compared to SE-EPI, a sequence that has been well studied and provides functional signal with specificity to the microvasculature ( Yacoub et al., 2007 ). We demonstrated the utility of phase regression for intracortical fMRI by showing GE-EPI-PR (1) has specificity across cortical surfaces comparable to SE- Fig. 11. Laminar activation profiles across subjects. Error bars represent the standard error of the mean across subjects. a) Profile across all vertices, b) Profile across vertices proximal to a vein (thresholded at a Euclidean distance of 2.4 mm to a vessel vertex) and c) Profile across vertices distal to all veins. EPI, (2) has higher sensitivity than SE-EPI across all cortical layers, and (3) reduces the extravascular and intravascular functional contributions from pial veins compared to GE-EPI. With these advantages in mind, GE-EPI-PR is a useful addition to a laminar imaging toolkit as it improves specificity of GE-EPI with only a minor reduction in sensitivity.

GE-EPI: Specificity and Sensitivity
GE-EPI is the workhorse sequence for fMRI studies and has advantages over T 2 based methods such as SE-EPI as it requires less RF power (lower SAR) and has higher SNR efficiency. As a result of the high SNR efficiency, GE-EPI has higher contrast-to-noise per unit time than SE-EPI, from 2 to 2.9 experimentally ( Moerel et al., 2018a ;Rua et al., 2017 ;Schumacher et al., 2011 ) and this has beneficial effects when voxels are evaluated for activation models such as tuning or encoding as the fits are more robust ( Moerel et al., 2018a ;Olman et al., 2010 ). However, GE-EPI profiles as a function of cortical depth have a positive slope towards the cortical surface, indicative of large BOLD changes due to large pial vessels ( Budde et al., 2014 ). The tradeoff of sensitivity for specificity between GE and SE is further complicated by the high SAR requirements of SE-EPI which lengthen acquisition time and limit coverage ( Kemper et al., 2015 ). Alternative sequences such as GRASE and VASO have been developed that have improved specificity compared to GE-EPI but also have reduced CNR, spatial coverage limitations and SAR restrictions ( Huber et al., 2015 ;Kemper et al., 2015 ). Thus GE-EPI remains the most popular fMRI sequence to date and is widely used in the high-resolution fMRI field.

GE-EPI-PR: Specificity and Sensitivity
Large venous vessel BOLD signal reduction through phase regression produces activation maps ( Fig. 7 ) and laminar profiles ( Fig. 11 ) comparable to SE-EPI. Activation map similarity was quantified through two-sided Kolmogorov-Smirnov tests which show GE-EPI-PR and SE-EPI activation map distributions are not significantly different on the upper laminar surfaces of cortex (10-60% cortical depth). This demonstrates that phase regression produces a SE-EPI-like signal which will have higher specificity to microvasculature without incurring the conventional penalties for that specificity such as higher SAR and longer imaging times. Using GE-EPI-PR also offsets one of the major problems with SE-EPI, namely reduced sensitivity.
The low BOLD sensitivity of SE-EPI ( Koopmans and Yacoub, 2019 ;Yacoub et al., 2005 ) produces data with low contrast-to-noise efficiency and often requires repeated acquisitions to increase the statistical power. Consistent with prior studies, we found SE-EPI has 30% of the CNR of GE-EPI data averaged across layers. Our approach demonstrates GE-EPI-PR doubles the CNR compared to SE-EPI across all layers which will make imaging using this technique more statistically powerful ( Fig. 10 ). This technique also shows GE-EPI-PR has 60% of the CNR of GE-EPI, which is the same as the CNR of VASO ( Huber et al., 2015 ). Utilizing GE-EPI-PR will therefore create a more microvasculature-weighted signal with increased sensitivity and some practical acquisition advantages over alternative fMRI sequences.

GE-EPI-PR: Venous signal suppression
Our hypothesis that phase regression decreases macrovascular signal in GE-EPI-PR activation maps predicts a lower activation in the areas that correspond to veins. Indeed, in this study, the GE-EPI-PR activation maps ( Fig. 7 ) show spatially varying suppression compared to GE-EPI with the largest suppression in the 'vessel' regions as identified from the multi-echo GRE scan ( Fig. 8 ). This observation is further supported in Fig. 11 showing that the GE-EPI-PR laminar profiles in vertices proximal to vessels show increased signal suppression compared to vertices distal to vessels. These results are a promising indication that at high resolution, phase regression has the ability to also suppress extravascular signal ( Vu and Gallant, 2015 ), which was not observed in previous studies at low spatial resolution ( Barry and Gore, 2014 ;Menon, 2002 ). Extravascular signal is the dominant BOLD producer at 7T ( Duong et al., 2003 ) and the specificity improvement from extravascular signal removal is particularly significant for pial veins as they have been shown to impact the signal distribution across the entire cortical ribbon of the visual cortex ( Fracasso et al., 2018 ). In addition to the extravascular suppression, this study observed intravascular suppression as expected by phase regression ( Menon, 2002 ). This has a similar effect to applying a diffusion gradient to a GE sequence to suppress intravascular BOLD effects ( Boxerman et al., 1995 ;Duong et al., 2003 ). Both extra-and intra-vascular signal suppression was greatest at the pial surface supporting the theory that phase regression exhibits the highest suppression effects near large vessels. All of the tangential and penetrating vasculature in cortex combines to confound BOLD signal distant from the capillary bed but for GE-EPI the effects from pial veins are dominant ( Bause et al., 2020 ) as proximity to a vessel affects the amplitude of the BOLD response to a greater degree than cortical depth. Phase regression results in functional maps with higher microvascular specificity to the capillary bed. This is important as previous studies have shown that performing venous removal on GE BOLD data results in laminar profiles more closely matched to the expected laminar profiles from histology ( Huber et al., 2015 ;Koopmans et al., 2010 ).

Venous signal removal from GE-EPI in literature
Removing GE-EPI signal contributions from large venous vessels to increase the specificity to the microvasculature remains one of the open problems in high resolution GE-EPI fMRI research today. Several studies have demonstrated reducing large vessel signal contributions from GE-EPI BOLD data using masking, profile correction, experiment setups or selective analysis. One such approach, masking, can be performed using additional acquisitions such as multi-echo GREs to identify and mask venous vessels ( Chen et al., 2013 ;Moerel et al., 2018a ) but suffers from poor localization of the venous voxels after registration of the multi-echo scan to the distorted EPI space ( Polimeni et al., 2018 ). Phase regression is performed in native EPI space so will not suffer from these potential registration errors. Alternatively, it is possible to mask vessels by determining cutoff thresholds of EPI intensity or percent BOLD change in order to separate venous voxels from non-venous voxels in native EPI space ( Kay et al., 2019 ;Koopmans et al., 2010 ) but hard cutoffs may not be able to separate venous and non-venous signal completely and may require manual segmentation (as was done in this study) or additional filtering ( Koopmans et al., 2010 ). Additionally, hard cutoffs fail to account for the gradual distance dependent reduction in extravascular effects. Fortunately, phase regression requires no cutoffs and GE-EPI-PR also is useful at removing extravascular effects from pial veins proximal to vessels. Laminar profile correction can be completed spatially through PSF estimation and deconvolution to remove bias from penetrating vessels ( Markuerkiaga et al., 2016 ) but it does not consider pial vein effects ( Koopmans and Yacoub, 2019 ) unlike the phase regression technique. It is also possible to correct the profiles temporally by estimating an early and late response across an area through temporal decomposition through manifold fitting ( Kay et al., 2020 ). This technique, like phase regression is agnostic to venous size or orientation but does require finite impulse response modelling as an initial step which can be challenging in resting state or naturalistic paradigms. Using phase and magnitude data separately could add additional power to the manifold fitting approach. Phase regression uses the correlation between magnitude and phase and therefore works across functional paradigms. Some experiments can reduce vessel bias using their experimental design, such as ocular dominance columns, where contrast subtraction removes most of the large venous effects ( Cheng et al., 2001 ;Moerel et al., 2018a ) but this assumes linearity in the BOLD response and still shows some venous contamination ( Yacoub et al., 2007 ). It also eliminates the desirability of using single-condition maps. Other forms of selective analysis can be performed, such as focusing analysis on the initial dip of the BOLD signal as it is more spatially specific to the active microvascular blood pool, unfortunately it is smaller and requires additional modelling in order to determine the HRF voxel by voxel ( Siero et al., 2013 ). Alternatively, one can deliberately remove upper layers from further analysis ( Polimeni et al., 2010 ) which limits the utility of intracortical fMRI to neuroscience problems fully described in the lower cortical depths. These experimental restrictions are not required by phase regression. Phase regression is an additional viable tool for this venous reduction literature as phase data is already available for many gradient echo sequences and only requires a robust phase sensitive method of coil combination.

Study limitations
It is important to note, for studies requiring voxel sizes of less than a millimeter, appropriate echo times may only be achievable by using acceleration such as GRAPPA ( Griswold et al., 2002 ) or SENSE ( Pruessmann et al., 1999 ) possibly in combination with partial Fourier acquisition ( Feinberg et al., 1986 ) which is ubiquitous in high-resolution EPI fMRI acquisitions in order to obtain short echo train lengths. The phase regression technique works regardless of partial Fourier as long as more than half of the k-space is collected. Although we expect the phase to be affected by the partial Fourier we do not expect, nor do we observe, its complete destruction. Our data did not exhibit artifacts such as Gibbs ringing in the phase data which we could expect from the use of zero-filled partial Fourier. This may be due to the relatively low SNR this data was collected with obscuring the expected ringing. This evaluation demonstrates that partial Fourier will produce data with a lower effective resolution but without contributing significant additional artifact in the images and/or resultant functional maps. Future work on phase regression will have to perform similar quality assurance in order to determine that the GE sequence and acceleration parameters used are appropriate for phase data.
To examine the effects of sequence parameters on resolution, the magnitude of the complex PSF was reported for both GE-EPI and SE-EPI.
The PSF provides additional acquisition information to the commonly reported nominal resolution and allows for an improved understanding of the effect acceleration has on our data. Several existing studies have attempted to compare PSFs this way in order to better explain the effects that different sequences and acceleration parameters have on their data ( Feinberg et al., 2018 ;Kemper et al., 2015 ). This method calculating the magnitude of the complex PSF will not represent the physiological PSF ( Koopmans and Yacoub, 2019 ;Markuerkiaga et al., 2016 ;Menon and Goodyear, 1999 ) and will not provide results targeted at resolving a specific pattern such as ODCs ( Chaimow and Shmuel, 2017 ). However, our reported PSFs still provide a direct comparator between sequence parameters. This magnitude of the complex PSF represents the level of influence neighboring voxels have on each other in the phase-encode direction, the worst blurring case in our sequences ( Chaimow and Shmuel, 2017 ). Despite the blurring due to acceleration limitations these voxel sizes were sufficient to study reductions in macrovascular signal and still showed the phase regression effect at high resolution for the first time.

Future work
Future work is needed to explore the properties of phase regression at high resolution, This study was conducted with limited spatial coverage, 22 mm in the slice direction limiting the ability to assess other aspects such as the relationship between phase regression and cortical orientation ( Stanley et al., 2019 ). This was a deliberate choice due to the nature of the visual stimulation used. More study is needed to assess whether the phase regression effects reduce the orientation dependence of GE-EPI which is driven by pial vessels ( Viessmann et al., 2019 ). One additional area of interest would be the inclusion of phase regression into laminar modelling ( Havlicek and Uludag, 2019 ;Markuerkiaga et al., 2016 ) as these methods focus on removal of penetrating vasculature and not correction for vessels on the pial surface. Finally, extension of phase regression to other GE sequences such as 3D-EPI could allow for wider adoption of this technique ( Hendriks et al., 2020 ). These proposed studies would help to expand the utility of phase regression beyond the investigation performed in this study.

Conclusions
This study has demonstrated that phase regression can be applied to reduce large vessel bias in high resolution functional acquisitions with complex data. Applying phase regression to GE-EPI data results in a similar activation map to SE-EPI while maintaining a higher contrastto-noise ratio. Phase regression may be a useful tool in the laminar fMRI toolkit. This valuable technique can be used without additional acquisitions or equipment and requires only a method to combine phase data. Phase regressed GE-EPI is a powerful technique to reduce venous bias considered to be an important confounding factor at high fields and thus allowing GE-EPI imaging to have increased utility in laminar fMRI studies.

Data availability statement
All tools developed for this paper are shared in the phaseprep repository ( https://github.com/ostanley/phaseprep ). Any code not provided in the phaseprep repository that supports the findings of this study are available from the corresponding authors upon request.
Raw data from four subjects in this study have their complete datasets uploaded here in BIDS format: https://openneuro.org/ datasets/ds003427 . The remaining subjects did not consent to open data sharing and as such are not available.