Predicting the macrovascular contribution to resting-state fMRI functional connectivity at 3 Tesla: A model-informed approach

The use of resting-state blood-oxygenation level-dependent (BOLD) functional magnetic resonance imaging (rs-fMRI) has become increasingly valuable in assessing brain health and mapping brain connectivity (Srivastava et al., 2022; van den Heuvel & Hulshoff Pol, 2010). There has been extensive use of this technique in the study of numerous neurological conditions (Ovadia-Caro et al., 2014; Vemuri et al., 2012), as well as in the study of brain development and ageing processes (Andrews-Hanna et al., 2007; Geerligs et al., 2015). The most common method for calculating functional connectivity (FC) using resting-state magnetic resonance imaging (rs-fcMRI) is Pearson’s correlation analysis (Bandettini et al., 1993; Biswal et al., 1995), in which the BOLD signals from seed voxels are correlated with those from the rest of the brain. Commonly, correlation is interpreted as reflecting synchronous neuronal activity (Biswal et al., 1995; Pan et al., 2015). The contribution of the macrovasculature to functional connectivity (FC) in the grey matter (GM) has previously been demonstrated in human data (Huck et al., 2023; X. Z. Zhong et al., 2024). Specifically, our investigation indicated that, at 3 Tesla, such effects occurred not only within the confines of macrovasculature but also more than 10 mm away from individual microvascular blood vessels, consistent with findings by others (Bernier et al., 2021; Huck et al., 2023) including at 7 T. However, there is still limited systematic investigation into macrovascular contributions to rs-fMRI.

Unlike in rs-fMRI, the question of “vein versus brain” has been a recognized challenge for many years in task-based fMRI. Boxerman et al. (1995) demonstrated, through both rodent experiments and simulation studies, that macrovasculature affects gradient-recalled-echo (GRE) BOLD contrast more than microvasculature. Despite the fact that microvascular sensitivity increases with an increase in the main magnetic field (Menon & Goodyear, 1999), veins still account for the majority of the contribution to GRE BOLD signals at 7T (Menon, 2012; UIudag et al., 2009) and show strong association with both BOLD signal amplitude (1.5T) and timing (3T) (Provencher et al., 2018; Vigneau-Roy et al., 2014). Previous studies in rats showed that venous structure could be detected directly using fMRI (Hyde & Li, 2014), and that the BOLD response from veins is nearly double that from microvessels in brain tissue at 11.7 T (Yu et al., 2012). Moreover, several studies have suggested that penetrating arteries and ascending veins also do not strictly reflect local neural activity (Uludağ & Blinder, 2018), despite their much smaller size relative to the macrovasculature investigated to date.

Numerous approaches have been developed to remove the venous BOLD effect. The simplest method of correcting the bias appears to be to mask out the detectable macrovasculature using vascular masks, such as segmented from angiograms. However, as work by ourselves and others have shown (Bernier et al., 2021; Huck et al., 2023; X. Z. Zhong et al., 2024), macrovascular effects can extend well into the perivascular tissue, making masking-based correction inadequate. Spin-echo (SE) BOLD has also been proposed as an alternative to the conventional GRE BOLD to suppress the strong macrovascular susceptibility effects. Nevertheless, this option sacrifices signal-to-noise ratio (Menon, 2012), and the use of typical echo-planar imaging (EPI) readouts even with SE refocusing is still susceptible to large-vein T2* effects (Goense & Logothetis, 2006; Ragot & Chen, 2019). The BOLD image phase has also been proposed as a regressor to remove macrovascular effects in task-based fMRI (Menon, 2002; Stanley et al., 2021), assuming that the signal magnitude and phase scale with one another, and that the intravascular (IV) contribution is dwarfed by the extravascular (EV) contribution due to low voxel resolution at conventional field strengths. Phase regression has even been shown to be effective at reducing large-vein contributions even at high field and high resolution (Stanley et al., 2021). However, extraction of accurate phase values can be a challenge. As vascular diameter increases or as the voxel size decreases, it becomes less appropriate to ignore the IV phase offsets and vascular orientation effects, as discussed in our previous work (X. Zhong & Chen, 2022, 2023). Most recently, Huck et al. (2023) advocated the use of higher order polynomials for modelling the fall-off of extravascular fields over space and subsequently removing venous biases, both proximal and distal to the vasculature. However, the authors concluded that their model is not sufficient to correct venous bias in rs-fMRI data. In light of these recent findings, it may be necessary to consider a more comprehensive modelling approach for the vascular BOLD contributions as the initial step in correcting for them.

The modelling of macrovascular contributions to GM FC with a closed-form analytical model is a challenging undertaking. Earlier work by Ogawa, Menon, et al. (1993) demonstrated that the BOLD signal is dependent on the orientation and diameter of the vessels, as well as blood oxygenation, but some of these parameters are difficult to measure reliably in the human brain. The dependence of the BOLD signal on vascular geometry (mainly of large pial veins) was also evident in the biophysical simulations based on realistic vascular anatomy (Gagnon et al., 2015; X. Zhong & Chen, 2022) and in in-vivo fMRI measurements (Fracasso et al., 2018; Viessmann et al., 2019). Further simulations also revealed that the position of a large vessel position within an fMRI voxel also plays a significant role in modulating macrovascular BOLD (Sedlacik et al., 2007; X. Zhong & Chen, 2023), increasing the complexity of modelling macrovascular effects. Despite not having been investigated previously, other factors (such as patient positioning and spatial resolution) may also have a significant impact on modelling macrovascular effects.

Instead of using a closed-form analytical model, it is possible to use a numerical model for calculating macrovascular BOLD signals that may simplify these issues. A kernel-based approach, proposed in previous studies (Cheng et al., 2009; Salomir et al., 2003), involves convolving kernels (representing point-spread functions) with susceptibility maps to compute local field offsets and, consequently, to generate a BOLD signal based on the blood susceptibility. The approach allows us to directly simulate the BOLD signal arising from realistic macrovascular anatomical networks (macro-VANs: macrovascular-specific human VAN), although with the limitation that blood oxygenation remains unknown. It is important to note, however, that this numerical approach will result in approximation errors and require a significant amount of computational resources (as will be discussed in later sections). In this regard, it remains to be seen whether the numerical approach with realistic vasculature has an advantage over the analytical approach based on a cylinder model.

The current study is intended to provide a detailed analysis of the macrovascular contribution to FC using biophysical modelling. Moreover, two analytical models will be tested in addition to a numerical model in terms of establishing the feasibility of a theoretically driven framework for predicting the macrovascular BOLD signal in the resting state. Furthermore, the investigation would be extended to the perivascular tissue in order to gain a better understanding of the macrovascular contribution beyond the confines of the macrovasculature, which is of practical impact for connectivity mapping. As of the time of writing, this is the first study to investigate the prediction of the macrovascular contribution to resting-state FC using biophysical modelling.

In this work, we compare experimental rs-fMRI data with simulated rs-fMRI data on a voxel-wise basis. We used data selected from the Midnight Scan Club (MSC) data set (Gordon et al., 2017). The study protocol was approved by the Human Studies Committee and Institute Review Board at Washington University School of Medicine in accordance with the Declaration of Helsinki. These data can be obtained from the OpenNeuro database, with accession number ds000224. Due to the high computational demand of our modelling approach, we balanced computational cost with the generalizability of the results, selecting data from four healthy right-handed young participants, two males and two females, ages 28–34 years.

Each subject underwent 12 imaging sessions on a Siemens TRIO 3T MRI scanner (Siemens Healthcare GmbH, Erlangen, Germany) on separate days. Here, only the protocols related to this study are listed.

In total, 12 of 2D time-of-flight (TOF) angiograms were acquired, including 4 ascending (transverse, 0.6 × 0.6 × 1.0 mm³, 44 slices, TR = 25 ms, TE = 3.34 ms, flip angle (⁠$α$⁠) = 20^o), 4 left-right encoded (sagittal, 0.8 × 0.8 × 2.0 mm³ thickness, 120 slices, TR = 27 ms, TE = 7.05 ms, $α$ = 60^o), and 4 anterior-posterior encoded (coronal, 0.7 × 0.7 × 2.5 mm³ thickness, 128 slices, TR = 28 ms, TE = 7.18 ms, $α$ = 60 degrees) data sets. Resting-state fMRI data were acquired with a gradient-recalled-echo BOLD-weighted EPI sequence (TR = 2.2 s, TE = 27 ms, $α$ = 90^o, 4 mm isotropic resolution, 36 slices, scan time = 30 min). T1-weighted anatomical MPRAGE data were included as well (sagittal, 224 slices, 0.8 mm isotropic resolution, TE = 3.74 ms, TR = 2,400 ms, TI = 1,000 ms, $α$ = 8 degrees). The participants were instructed to fixate on a white crosshair on a black background during the rs-fMRI scans. An EyeLink 1000 eye-tracking system (SR-Research, Ottawa, Canada, http://www.sr-research.com) was used to monitor participants to determine whether they had any prolonged eye closures, which may indicate sleepiness. There was only one participant who showed prolonged eye closures during the scans.

A summary of the rs-fMRI processing procedures can be found in Figure 1a. fMRI preprocessing pipeline was implemented with tools from FSL (Jenkinson et al., 2012), AFNI (Cox, 1996), and FreeSurfer (Fischl, 2012). The following steps were included in the preprocessing steps: (a) 3D motion correction (FSL MCFLIRT), (b) slice-timing correction (FSL slicetimer), (c) brain extraction (FSL bet2 and FreeSurfer mri_watershed), (d) rigid body coregistration of functional data to the individual T1 image (FSL FLIRT), (e) regression of the mean signals from white-matter (WM) and cerebrospinal fluid (CSF) regions (fsl_glm), (f) bandpass filtering to obtain frequency band 0.01–0.1 Hz (AFNI 3dBandpass), and (g) spatial smoothing with a 6 mm full-width half-maximum (FWHM) Gaussian kernel (FSL fslmaths). For each participant, after preprocessing, we calculated the voxel-wise venous–venous, arterial–arterial, and arterial–venous FC metrics based on the maximum cross-correlation coefficients for the voxels with positive correlation coefficients and the minimum cross-correlation coefficients for the voxels with negative correlation coefficients. In spite of the fact that most fMRI connectivity analysis is based on Pearson’s correlations, we used cross-correlations in order to compare the theoretically predicted correlation coefficients directly with the experimentally measured ones without taking into account time lags.

The strategies for macrovasculature segmentation and processing are summarized in Figure 1b. TOF images were registered to T1 space (FSL MCFLIRT) and segmented using the Braincharter Toolbox (https://github.com/braincharter/vasculature) (Bernier et al., 2018) (Fig. 2a). Visual inspection was performed to ensure the absence of major artefacts. The centreline image encoded with diameter information (“centredia” output from the Braincharter Toolbox) obtained from the raw images (at 0.8 mm isotropic voxel resolution) was pooled from all TOF images across all encoding directions, and where vessel overlaps were detected, the highest diameter estimate is assumed for the overlapping vessel. Accordingly, the vessels used in the subsequent analysis are primarily major cerebral vessels with a minimum diameter of close to 0.8 mm (e.g., the superior sagittal sinus and the Circle of Willis). This combinatorial approach maximized the completeness and signal-to-noise ratio of the resulting macro-VANs. The combined centrelines were then upsampled to a resolution of 0.2 mm using AFNI (Cox, 1996) with nearest-neighbour interpolation (3dresample). Based on each vascular centreline, a line was constructed that connects two vascular voxels that both neighbour the central voxel in a 9 × 9 × 9 voxel matrix. Vascular orientation was calculated as the angle between this line and the z-axis (zenith) and x-axis (azimuth). The zenith angle is further corrected by applying position metrics in MRI header files to compensate for the difference between the participants’ slice direction and the B₀ direction (scanner z-axis).

Vascular blood-volume fraction (fBV) and orientation maps were manually registered to BOLD space using coordinates from AFNI’s volume selection function (3dAutobox). fBV was estimated by counting the number of high-resolution voxels occupied by vessels in each fMRI voxel (4 mm isotropic voxel resolution), and macrovascular maps were derived from binarized fBV maps. In order to ensure the accuracy of registration, necessary quality control measures have been added. Due to the fact that all TOF images contained both arteries and veins, each arterial and venous map was manually separated based on an anatomical atlas (Tortora & Derrickson, 2018) (Fig. 2b, c). Lastly, the measured voxel-wise values of orientation and fBV are used to create different simulated voxels containing blood vessels, described in the following.

All simulation parameters and the values used in all models (Fig. 3) are listed in Table 1. In the remainder of this work, the tabulated parameter values are assumed unless otherwise stated. We do not expect the results to be affected by the choice of specific Y for arteries, veins, and tissues. As the spatial extents of B₀ inhomogeneities generated by large vessels are much greater than the diffusivity of water molecules at body temperature (Jensen et al., 2005), we assumed negligible diffusion effects. The simulated BOLD scan parameter values (e.g., TR and TE) were adjusted to match the parameter values used in the rs-fMRI scans. At these settings, the arterial inflow effect is deemed negligible (Gao et al., 1996). The Y-dependent T₂ of tissue and blood was calculated according to previous relaxometry studies (UIudag et al., 2009).

In this model, magnetic susceptibility maps are based on experimentally acquired macro-VANs (Fig. 4), based on this acquisition, all resolutions used in the simulation pipeline are listed in Table 2. A mask of the vasculature was generated by upsampling the macro-VANs to 0.2 mm isotropic resolution (The upsampling aims to allow at least four voxels to define the minimum vascular diameter (~0.8 mm, based on the spatial resolution of the TOF images), in order to minimize the ringing artefact due to undersampling), and necessary zero-padding by the size of a full field of view on each side of the matrix to avoid cycle convolution wraparound. The experimental data and simulated data were first aligned with the coordinates extracted from the processed macrovascular segmentation, as described earlier. To minimize computational complexity, the susceptibility map was derived using a Fourier-based model-independent convolution approach (Fig. 4) (Eqs. 1 and 2) (Cheng et al., 2009; Salomir et al., 2003), as follows.

FT denotes the Fourier transform, χ the local susceptibility, k_z the distance in the k space along the z-axis, and k the distance in k-space (⁠$k2=kx2+ky2+kz2$⁠). R₂’ is then calculated through the magnitude of the complex-valued mean magnetization of the dephasing spins resulting from the B₀ offset.

We further assumed the tissue type outside macrovasculature was grey matter (GM) to calculate the appropriate susceptibility difference between blood and tissue. The mean BOLD signal was calculated as defined by Equations 3 and 4 as

where $γ$ is the gyromagnetic ratio, and the operator | $μ$(·)| represents the magnitude of the mean of a complex number.

On the other end of the complexity spectrum is the single-voxel two-dimensional (2D) infinite-cylinder model, in which field offset is estimated in an analytical approach (Fig. 5a) (Eqs. 5 and 6) (Ogawa, Menon, et al., 1993). Here, instead of generating a volumetric mask of the macro-VAN, we calculated at each location of the macro-VAN the diameter and orientation of the best-fitting cylinder and used this information (combined with the blood-tissue susceptibility difference) to analytically calculate the magnetic field offset. This calculation is based on a closed-form expression for the extravascular and intravascular field offsets for an infinite cylinder that are a function of the vessel geometry and blood susceptibility, that is,

where ω₀ is the main magnetic field in terms of angular frequency; ∆χ, magnetic susceptibility of fully deoxygenated blood; Hct, Hematocrit; Y, the percent blood oxygenation (Y_a for artery and Y_v for vein); (Y_tissue−Y), the degree of deoxygenation difference between blood and tissue; θ, the angle between the B₀ and the vascular axis; a, the cylinder radius; r, the distance between the point of interest and the centre of vascular cross-section; φ, the angle between vector <r> and B₀ projection in the plane. This model assumes that the shape of the blood vessel is an infinite cylinder with zenith θ, azimuth φ, and radius a derived from fBV. These input parameters can all be obtained from the arterial and venous segmentations based on the TOF data. This 2D Cylinder Mode is suitable when the magnetic susceptibility does not vary along the length of the cylinder, and facilitates our understanding of the dependence of IV and EV rs-fMRI signals on a single large vessel. As in the case in the work of Bandettini and Wong (1995), we constructed a 2D voxel for each vascular orientation-fBV combination derived from the experimental data, divided into 4,000 × 4,000 1-μm isotropic sub-voxels, with each voxel defined by a set of vascular (arterial or venous) orientation and fBV values derived from the TOF data in the previous step. The 2D model is more conducive to analytical calculations of susceptibility, and the dimensionality reduction can be traded off for higher spatial resolution.

This was a 3D version of the infinite-cylinder model, with the exception that the vessel is not assumed to penetrate the voxel at a right angle, such that fBV will vary with θ, and that the susceptibility pattern can vary along the length of the vessel. We constructed a 3D voxel containing 100 × 100 × 100 40-μm isotropic sub-voxels, in which the infinite cylinder is defined by the collection of measured parameters as was done in the 2D Cylinder Model. The field-offset estimation procedure is the same numerical Fourier-transform-based method as used in the macro-VAN Model. This model serves as an intermediate between the macro-VAN Model and the 2D Cylinder Model, and is evaluated as a less computationally intensive approximation of the macro-VAN Model. The 2D and 3D Cylinder Models are both single voxel and thus do not take into account interactions between neighbouring vessels that may influence EV field offsets.

Vascular dynamics were simulated as sinusoidal variations at 0.1 Hz (Chu et al., 2018) for arterial fBV and venous Y (Uludağ & Blinder, 2018; Vazquez et al., 2013). We used a sinusoid with ±10% variations (peak-to-peak) for Y and fBV fluctuations in the 2D and 3D Cylinder Models. For macro-VAN Model, the macro-VAN was upsampled to 0.2 mm isotropic resolution, and the smallest detectable vessels had approximate diameters of 0.8 mm. Due to the coarser sample, the circular cross-section of each 0.8 mm vessel, assuming it is at the centre of its voxel, would be upsampled in 2D into a 13-voxel neighbourhood of 0.2-mm isotropic voxels. Thus, it was necessary to expand the fBV variations to ±30% to allow the quantized vascular fBV to fluctuate by one voxel in the diameter dimension. The simulated BOLD signal can be calculated as Equation 4.

From these simulated BOLD signals, we computed both resting-state fluctuation amplitudes (RSFA) and FC (i.e., the correlation between voxel time courses). For computing the RSFA, both the simulated and experimental BOLD signals were normalized by their respective means, with outlier time points (defined as lying beyond 1.5 quartiles of the mean) removed. Simulated FC metrics were calculated for voxel-wise arterial–arterial, venous–venous, and arterial–venous correlations with the simulated signals from all three models. As the simulated data are noiseless but are to be compared with noisy experimental results, we added noise for additional realism. Thus, as suggested by Golestani and Goodyear (2011), we used the temporal standard deviations of simulated time courses as a surrogate for SNR; that is, with noise being similar across voxels, a higher baseline signal fluctuation would translate into a higher SNR. We then multiplied the cross-correlation coefficients in the FC metrics with the σ for each simulated voxel to obtain SNR-independent FC measurements.

The simulations were conducted using our servers equipped with 14 cores of Intel Xeon X5687 CPU (at 3.6 GHz) (Intel Corporation, Santa Clara, CA, United States) and 180 GB of memory running Red Hat Enterprise Linux Server 7.7 (Red Hat Inc., Raleigh, NC, United States). A customized simulation script was written in MatLab 2019b (MathWorks Inc., Natick, MA, United States).

This is done for the macro-VAN case only. Predicted and measured rs-fMRI signals were compared in the perivascular tissue as well. Figure 1c illustrates the methods used to assess perivascular effects of the macrovasculature on the rs-fMRI signal, whereby the perivascular ROIs were generated as follows. First, the arterial and venous segmentations were dilated by one voxel in 3D. Second, the original vascular ROIs were subtracted from the dilated vascular ROIs. This produced the extravascular ROI that is centred at one voxel distance from the centre of the vasculature. This process was repeated for perivascular distances of two and three voxels (4 mm isotropic voxels). As discussed in the previous section, cross-correlation-based connectivity values were calculated separately for BOLD signals in arterial and peri-arterial ROIs, and between venous and peri-venous ROIs. This was done for both simulated and in-vivo experimental signals.

Each simulated voxel results in the simulated macrovascular-based BOLD signal specific to a vascular location in the experimental data, from which a corresponding experimental BOLD signal can be extracted. In this way, to quantify model prediction accuracy, the simulated RSFA and vascular-driven correlation coefficients (FC values) predicted by all three models were regressed against the experimental cross-correlation coefficients (FC value). All vasculature-related correlation coefficients (simulated and experimental) were then binned. Given that there are more voxels for venous ROIs than arterial ROIs, the bin sizes were adapted approximately proportionately, that is, 100 for arterial–arterial correlations, 1,000 for arterial–venous correlations, and 10,000 for venous–venous correlations (as there are more venous voxels than arterial voxels), to ensure a similar number of fitted data sets for all three categories. This was done separately for each data set. The R² values were calculated to assess the quality of the prediction in terms of the simulated FC values. For perivascular BOLD, a similar regression procedure was conducted for the perivascular FC metrics using the macro-VAN Model. The correlations were grouped by distance from the edge of the vasculature (i.e., one, two, or three voxels away). Similarly, the procedure described above was applied to assessing the predictability of BOLD signal fluctuation amplitudes.

The analytical equations (Eqs. 5 and 6) allowed us to gain an intuitive understanding of the interactions between the BOLD signal and various vascular parameters. When simulated based entirely on the default parameters listed in Table 1, arterial R₂’ increased with fBV but only up to fBV = 0.4, and decreased as fBV increased further (Fig. 6a). This is consistent with decreasing phase coherence as fBV increases in a tissue-dominated voxel followed by an increasing phase coherence once the voxel becomes blood dominated. However, as the vascular orientation varied, so did the behaviour of R₂’ (Fig. 6b). The peak R₂’ was generally not found at the θ extrema (corresponding to when the vessel is parallel or perpendicular to B₀). On the venous side, R₂’ declined with increasing Y_v, but this behaviour also varied with fBV. Venous R₂’ was maximized at an intermediate value of fBV. As shown in Figure 6c, it was higher for fBV = 0.4 than for lower and higher fBVs (0.1 and 0.7, respectively). Lastly, for the default fBV, like arterial R₂’, the venous R₂’ dependence on Y_v varied with θ, with R₂’ higher at an intermediate value of 3π/4 than at θ of π/2 or π (Fig. 6d). BOLD signal intensity was disproportionate to R₂’ (Fig. 6e–h): the R₂’ was maximized at intermediate values of fBV, the baseline BOLD signal intensity was minimized at intermediate values of fBV; when the R₂’ was maximized at intermediate theta values, the baseline BOLD signal intensity minimized at intermediate theta values.

In rs-fMRI, fluctuations in R₂’, rather than the mean R₂’, is more relevant. The quantity ∆R₂’, computed as the difference between the maximum and minimum R₂’ (associated with the maxima and minima in Y and fBV), can be viewed as the derivative of R₂’ relative to the oscillations in Y and fBV (Fig. 6a–d). Arterial ∆R₂’ initially rose with increasing fBV, then decreased as fBV surpassed ~0.2, eventually increasing again for fBV > 0.4 (Fig. 6i). Arterial ∆R₂’ decreased with θ as θ increases to 3π/4 and then increased for θ up to π (minimum at θ = 3π/4) (Fig. 6j). Venous ∆R₂’ in Figure 6k showed the opposite trend as arterial with fBV > 0.4 (Fig. 6i), which decreased with increasing fBVs. Lastly, like arterial |∆R₂’|, venous |∆R₂’| also decreased with θ until θ = 3π/4, then increased for greater values of θ. Moreover, ∆R₂’ became more stable while θ approached π. ∆BOLD was in turn computed as the difference in the BOLD signal corresponding to the two extremes of the input signals (fBV for artery and Y for vein), as derived from the associated R₂’ values. As seen in Figure 6m–p, the magnitude of ∆BOLD can be viewed as proportional to |∆R₂’|, whereas the polarity of ∆BOLD will influence the correlation between ∆BOLD time courses. For instance, the venous ∆BOLD in Figure 6o showed the opposite trend as arterial ∆BOLD (Fig. 6m), decreasing with increasing fBV until fBV ~0.2, increasing for higher fBVs and remaining positive throughout all fBVs. Thus, an equivalent increase in fBV for both an artery and a vein would lead to their BOLD time courses remaining negatively correlated. Also, arterial and venous ∆BOLD showed the same magnitude variation trend until arterial ∆BOLD changed polarity at fBV ~0.4, which suggests a negative correlation between arterial and venous BOLD time courses for fBV > 0.4. Lastly, arterial ∆BOLD increased with θ until θ = 3π/4, then increased for greater values of θ. This is contrary to the trend observed for venous ∆BOLD (peaking at θ = 3π/4) (Fig. 6p).

As shown in Figure 7a and d for a representative data set, arterial RSFA is less well predicted than venous RSFA by the macro-VAN Model. For the 2D Cylinder Model, both arterial and venous RSFA are poorly predicted (Fig. 7c, d), possibly attributable to a large cluster of measured BOLD RSFA values in voxels that exhibit a low simulated RSFA. The 3D Cylinder Model is an improvement over the 2D Cylinder Model especially for predicting venous RSFA, but the performance remains lower than for the macro-VAN Model (Fig. 7e, f). These trends are generalizable to all subjects (see Supplementary Materials). See Figure S2 for results for all subjects.

The simulated correlation coefficients based on the macro-VAN Model show the highest predictiveness for venous–venous correlation coefficients (Fig. 8g), followed by arterial–venous correlation coefficients (Fig. 8d), and with the lowest predictability for arterial–arterial correlation coefficients (Fig. 8a). For the 2D Cylinder Model, as expected, while significant positive associations between simulated and in-vivo correlation coefficients dominated, the goodness of fit was low (Fig. 8b, e, h), with the arterial–arterial correlation R² being the lowest (<0.01). We also noted that the data points appear to sub-divide into two groups, one spanning correlation coefficients of >0.04, and the other one spanning the area between 0 and 0.04. For the 3D Cylinder Model, the predictability of venous–venous correlations remains remarkably high (Fig. 8i), followed by arterial–venous correlations (Fig. 8f). The accuracy of the prediction of arterial–arterial correlations was relatively low (Fig. 8c), like in the cylinder models. See Figures S3–S5 in Supplementary Materials for results for all subjects.

Only the macro-VAN provided the possibility to examine the perivascular effects in a meaningful way, as it conveniently incorporates perivascular contributions from neighbour vessels. We found the majority of vascular contributions to the surrounding tissue to be limited to one voxel away (4 mm isotropic) from the original vascular masks (Fig. S6 in Supplementary Materials). As shown in Figure 9, both arterial and venous perivascular RSFA are poorly predicted by the macro-VAN Model, irrespective of the distance from the vascular voxels.

Unlike for venous RSFA, venous-driven correlation remains quite well predicted by the macro-VAN Model at one voxel away (voxels are 4 mm isotropic) from the vascular voxel, as shown by the R² value (Fig. 10f). The regression remains significant at two and three voxels away, however, R² decreased substantially (Fig. 10g, h). The predictability of arterial-driven correlations remains poor for all perivascular distances (Fig. 10a–d). See Figures S7 and S8 in Supplementary Materials for results for all subjects.

Comparisons of FC-prediction R² values from the three models based on data from all four subjects are summarized in Figure 11. In general, the macro-VAN Model provided the highest FC prediction accuracy (i.e., R² values), followed by the 3D Cylinder Model, and lastly the 2D Cylinder Model (Fig. 11a–c). Additionally, the simulated venous–venous FC values (Fig. 11c) are associated with the highest R² values (all <0.8), followed by arterial–venous FC values (Fig. 11b), and then arterial–arterial FC values (Fig. 11a). Moreover, perivascular prediction differences are also summarized in Figure 11 across all subjects and distances. Generally, the R² values of arterial–arterial FC were low, and increased with distance away from the original vascular mask (Fig. 11e). As with the FC simulation, the RSFA simulation shows a similar pattern, but with generally lower R². Both arterial and venous signals were better predicted by the macro-VAN Model than by the other models (Fig. 11f, g). Moreover, there is a trend of decreasing goodness of fit with increasing distance from the macrovascular voxel (Fig. 11h, i).

Despite the strong influence of macrovascular contributions to BOLD fMRI (Huck et al., 2023; Menon, 2012; Tong et al., 2015; X. Z. Zhong et al., 2024), previous research demonstrated that removing this bias is not an easy task (Huck et al., 2023; Menon, 2002; Ragot & Chen, 2019; Stanley et al., 2021). Using biophysical simulations informed by angiography data, the present study examines the feasibility of predicting the macrovascular BOLD contributions proximal and distal to the vasculature. We explored the use of a whole-brain macro-VAN Model as well as single infinite-cylinder approximations (2D and 3D). We found that:

1. With the availability of angiographic data, it is feasible to model macrovascular BOLD FC using both the macro-VAN-based model and 3D infinite-cylinder models, though the former performed better.

2. Biophysical modelling can accurately predict the BOLD FC (in terms of correlations) for veins, but not for arteries.

3. Biophysical modelling provided less accurate predictions for RSFA than for FC.

4. Modelling of perivascular BOLD connectivity was feasible at distances close to veins but not arteries, with performance deteriorating with increasing distance.

The macrovascular contribution to the rs-fMRI BOLD signal has many physiological origins. The concept of systematic low-frequency oscillations (sLFOs) was first proposed in an earlier review article and defined as a low-frequency vasogenic BOLD signal travelling through the brain (Tong et al., 2015). This phenomenon does not appear to originate within the brain (Frederick & Tong, 2010; Li et al., 2018; Tong & Frederick, 2010) and is likely associated with heart rate variability (Thayer et al., 2012), gastric oscillations (Mohamed Yacin et al., 2011; Rebollo et al., 2018), vasomotion (Hundley et al., 1988; Mayhew et al., 1996; Rivadulla et al., 2011), respiratory volume variability (Birn et al., 2006; Chang et al., 2009), and/or variations in carbon dioxide levels (Sassaroli et al., 2012; Wise et al., 2004), among others. In previous studies using near-infrared spectroscopy, it was demonstrated that such oscillations travel through the peripheral to the cerebral vasculature (Frederick et al., 2013; Tong & Frederick, 2012, 2014; Tong, Bergethon, et al., 2011; Tong, Hocke, et al., 2011; Tong, Lindsey, et al., 2011; Tong et al., 2014, 2015). As venous blood is predominantly responsible for the BOLD signal (Ogawa, Lee, et al., 1993; Segebarth et al., 1994), the presence of sLFOs may suggest that certain findings from rs-fMRI experiments may not reflect neural activity but rather venous bias (Aso et al., 2020).

To understand the influence of the macrovasculature on rs-fMRI, it is essential to have a basic biophysical understanding of the macrovascular BOLD effect. The macrovascular contribution is likely to be strongest within the confines of both macrovascular and perivascular tissues, as discussed in our previous study (X. Z. Zhong et al., 2024), but a full understanding of the large variability in the perivascular BOLD effect across brain regions and acquisitions (Huck et al., 2023; X. Z. Zhong et al., 2024) cannot be obtained from in-vivo experiments alone. We hope that our simulated results can contribute to a better understanding of this phenomenon.

The three simulation models compared in this study range in complexity and realism from low to high, namely from the 2D infinite-cylinder Model, to the 3D infinite-cylinder Model, to the empirical macro-VAN Model. Unlike microvascular BOLD simulations (Berman et al., 2018), macrovascular simulations cannot ignore fBV changes caused by vessel obliquity, as macrovascular fBV is more substantial and arises from much larger “cylinders.” The 3D Cylinder Model is expected to be more accurate than the 2D Cylinder Model due to its ability to capture the effect of orientations (azimuth and zenith) on fBV. As an example, a vessel with an azimuth angle of 45° or zenith angle of 45° would have an fBV of approximately 1.41 times that of a vessel with both angles at 0. Thus, the 2D Cylinder Model is expected to be outperformed by the 3D Cylinder Model, which also models vessels as cylinders but more accurately captures fBV and vessel orientation. While the 2D Cylinder Model was expected to provide low accuracy, it was unclear a priori how these simplifying assumptions would affect our BOLD predictions, and given its practical advantages it was evaluated in our study. Our conclusion is that, for modelling large vessels, the 2D model is inadequate and the disadvantages outweigh the advantages.

Likewise, we expected the macro-VAN Model to produce the most realistic BOLD simulations since it incorporates the shape and position effects of real macrovasculature that cannot be captured by the 2D and 3D Cylinder Models. The macro-VAN Model can capture details missing from the 3D Cylinder Model, such as the geometry of blood vessels (with bends, branching, etc.), partial volume effects, and positional effects (when the vessel is not centred in a voxel). Despite the advantages of the macro-VAN Model, it still has some limitations. Because it is based on the Fourier-based discrete Green’s function, the accuracy of this approach depends on the grid resolution. More specifically, to ensure accuracy, the diameter of the smallest blood vessels should span at least four voxels (Tang et al., 1993). Coarse spatial sampling of the vascular susceptibility gradient and boundaries leads to truncation or “ringing” effects, which can be alleviated by spatially upsampling, thus smoothing of the kernel function. This is in contrast to the 2D Cylinder Model, which can be constructed at an arbitrarily high spatial resolution (by using smaller sub-voxels) to provide a more accurate representation of fBV variations, which is especially relevant for small fluctuations in the resting state. Furthermore, given the complexity of the macro-VAN Model, a more intuitive understanding of the BOLD signal dependence on vascular geometry is challenging. Such an understanding can be more easily achieved with the 2D Cylinder Model, which allows for the effect of different geometric parameters to be examined in isolation (Ogawa, Menon, et al., 1993). Furthermore, the computational complexity of the macro-VAN Model is also the highest, as it takes a considerable amount of computational memory and processing time to simulate the entire macrovasculature, as shown in the next section. The 3D Cylinder Model provides a trade-off between geometrical accuracy and computational requirements, and was investigated as a more practical alternative to the macro-VAN Model, to potentially make such modelling methods more practical.

We also note that high accuracy comes at a cost. The macro-VAN Model can consume as much as 160 GB of memory and 75 h of computation time for arterial simulations alone. The 3D Cylinder Model, though much simpler (2 GB memory and 108 h computation time), remains much more computationally demanding than the 2D Cylinder Model (1 GB memory and 7.5 h computation time) (Table 3). According to the discussion above, the 2D Cylinder Model could perform very similarly to the macro-VAN Model in certain scenarios, and in fact, most S_v,model3:S_v,model1 values are lower than 3 (zoomed window of Fig. S1b in Supplementary Materials). It is the same for the prediction of correlation coefficients. Thus, in order to select the appropriate simulation model, it is important to consider the error tolerance and the availability of computational resources. Our study is not focused on numerical algorithms, and the computational resources listed here are only meant to provide a comparison and may not be based on the most efficient algorithms or reflect all hardware architectures.

According to the classic biophysical model proposed by Ogawa, Menon, et al. (1993), the BOLD signal is dependent on the vascular orientation, size, and oxygenation. θ and fBV are two parameters widely investigated in relation to BOLD signal dependencies (Ogawa, Menon, et al., 1993; Viessmann et al., 2019; X. Zhong & Chen, 2022). The ΔBOLD metric, calculated based on the difference in the BOLD signal corresponding to the two simplified cases of the input signals (variations in fBV for arteries and in Y for veins), was used to demonstrate the dependency of each parameter on the simulated BOLD fluctuations, whose magnitude should be proportional to the RSFA and whose polarity will directly affect the correlation with other BOLD signals from other locations.

In contrast with previous experimental observations at 7 T based primarily on mesoscopic pial veins and the associated extravascular fields (Viessmann et al., 2019), the magnitude of macrovascular ΔBOLD in this work does not consistently decrease with increasing θ, but is rather minimized at θ ~ 3π/4 (Fig. 6m, n), which, according to the original biophysical model (Ogawa, Menon, et al., 1993), leads to zero intravascular magnetic field offset. The macrovasculature, in contrast to microvasculature, has a higher fraction of intravascular volume, such that its IV contribution cannot be neglected. It is intriguing to note that, although venous ΔBOLD remains positive for all fBV, arterial ΔBOLD changes polarity from negative to positive with an fBV of ~0.4 (Fig. 6o, p). It is possible that a decrease in the extravascular compartment associated with an fBV increase could lead to a dominance of intravascular contribution to the measured BOLD fluctuations. (In the case of arteries, an increase in the intravascular compartment would be expected to increase phase homogeneity, since all intravascular spins in a given voxel would be in phase.) As a result of the change in polarity of the simulated signal, arterial–arterial correlations could potentially be negative with specific sets of fBV, just as the arterial–venous correlation could be positive. These complicated dependencies suggest that multiple factors must be considered simultaneously when modelling macrovascular contribution to both BOLD fluctuation magnitude and correlation, and that simply considering diameter would not suffice for voxel-based modelling (Huck et al., 2023).

However, these are not the only parameters that could affect the simulated BOLD signals. Our previous study demonstrated that even unexpected parameters such as the position of a vessel within a voxel can have an impact on macrovascular BOLD (Sedlacik et al., 2007; X. Zhong & Chen, 2023) and this can be extended to the impact of partial-volume effects. It is, however, difficult to measure the exact position of vessels within the fMRI voxels in our in-vivo experimental data, so the simulation was not able to include this effect. It is also possible that our results may be affected by other parameters (such as TE, spatial resolution, and the main magnetic field), which are not addressed in this initial study but could be considered in the future.

A general conclusion can be drawn that the macro-VAN Model provides better predictions than the 2D Cylinder and 3D Cylinder Models for RSFA (Figs. 7 and 11). Interestingly, the macro-VAN Model was able to predict venous-driven RSFA reasonably well. As discussed in the theory section, the macro-VAN Model has the most realistic geometric representation of macrovasculature, and while the 3D Cylinder Model is less realistic, it still embodies the interactions between vascular orientation and fBV, whereas the 2D Cylinder Model does not. This is also reflected in the comparison of the simulated baseline BOLD signal intensities (see Fig. S1 in Supplementary Materials), whereby the unrealistic geometry of the 2D Cylinder Model leads to predicted signal fluctuation amplitudes that deviate (up to 60%) from that of the macro-VAN Model. Nevertheless, none of the three models is capable of providing satisfactory prediction for arterial RSFA. In the 2D Cylinder and the 3D Cylinder Models, we observed a cluster of low simulated values that do not fit well into a linear relationship, which suggests unreliable prediction for some of the voxels.

Moreover, the RSFA prediction accuracy is highly participant dependent as shown in Figure 11, which may in turn be dependent on the parameters of the acquisition and on the distinct anatomical structures of each participant (e.g., participant positioning and geometry of the macrovasculature). Moreover, different participants may also be affected by different levels of physiological noise (e.g., heart rate variability, respiratory variability), which further increases the level of inter-subject variability. As well, the effect of motion artefacts cannot be ignored. Different participants may have different levels of motion (either gross head motion or local physiological pulsation) which result in different levels of noise for both rs-fMRI and TOF imaging.

In general, unlike for venous RSFA, none of our models provided satisfactory predictions of arterial RSFA in any of our subjects (Fig. 11). Since the diameter of the arteries is constantly changing, it is almost impossible to establish an accurate baseline diameter for use in our simulations. Additionally, the duration of the TOF angiogram acquisition is significantly longer than the arterial variation (the MRA acquisition in the MSC data set lasted approximately 5 min) (Gordon et al., 2017; Reimers et al., 2018), resulting in a further blurring of the arterial vasculature boundaries. A further challenge is emulating arterial movement in simulation, particularly for the macro-VAN Model, which requires extremely high resolution and is, therefore, not practicable. While some state-of-the-art acquisition techniques allow dynamic imaging of arteries, they still only cover one blood vessel rather than the entire vasculature (Varadarajan et al., 2023). Further research to improve the temporal resolution of angiograms will be necessary.

The discrepancy between the predicted and the experimental FC could be, in part, due to the simulated signal having a much higher signal-to-noise ratio (SNR) than the in-vivo signal. To account for this, we weight the simulated connectivity with the simulated RSFA in order to emulate the in-vivo experimental signals better (Golestani & Goodyear, 2011), as described in the methods section. This approach prevented the FC maps from being dominated by large vessels.

Like in the case of the RSFA, the macro-VAN Model resulted in higher FC prediction accuracy than the 2D and 3D Cylinder Model (including arterial–arterial, arterial–venous, and venous–venous FC). The 3D Cylinder Model in turn performed better than the 2D Cylinder Model in all cases (Fig. 7). Venous–venous correlation was remarkably well predicted by the macro-VAN Model, followed by arterial–venous correlation. Both the macro-VAN Model and the 3D Cylinder Model generated satisfactory predictions for venous–venous correlation, which suggests that the 3D Cylinder Model may be used as a substitute for the macro-VAN Model. However, for the 2D Cylinder Model, we noted in cases with poor prediction that the simulated values appear to be sub-divided into two groups, the first spanning correlation coefficients greater than 0.04 and the second spanning correlation coefficients between 0 and 0.04. The existence of two clusters of correlation coefficients may result from signal intensity outliers as discussed above (see Model Theoretical Comparison). In light of the results shown in Figure S1 in Supplementary Materials, we suggest that the 2D Cylinder Model is likely to produce spurious BOLD signal fluctuation amplitudes with specific vascular orientation and fBV values that may distort the predicted time courses, which translate into extreme RSFA values and thus widely fluctuating FC values.

Similar to the case of RSFA, we were not able to predict arterial–arterial correlation with any of the models. Interestingly, predicted arterial–venous FC values were anti-correlated with the experimental values (Fig. 8d–f). Arterial BOLD contributions are generally weaker than venous contributions, whereas the angiograms are more sensitive to the fast-flowing arterial blood than the slower moving blood of the neighbouring veins (such as the cavernous sinus, which is in proximity to the Circle of Willis). Thus, when a single voxel contains both an artery and a vein, one dominating the segmentation and the other dominating the susceptibility effect, the model prediction may be particularly challenging. Certainly, given the poor accuracy of arterial predictions seen in this study, it is clearly more feasible to pivot towards estimating large-vein contributions using our approach.

Our R² values are lower than those of the binned predictions reported by Huck et al. (2023) (R² > 0.98), but much higher than those of their voxel-wise predictions (R² < 0.05). The results of our study, although still based on binning, are based only on predicted correlation coefficients from first principles, which should allow better representation of the voxel-wise macrovascular effects. Since Huck et al. (2023) also examined other connectivity metrics (low-frequency fluctuation amplitude, regional connectivity inhomogeneity, Hurst exponent, and eigenvector centrality), a direct comparison of these two studies is difficult.

Interestingly, our FC predictions appear to be more accurate than our prediction of RSFA, especially when comparing venous RSFA and venous–venous correlation. There may be several reasons for this. First, measured RSFA likely includes factors that are not incorporated in our model such as multiplicative scaling from receive coil sensitivities, which may scale amplitudes of the fluctuations but not their timing. Correlation coefficients are normalized quantities and thus are less influenced by these regional scalings. Second, following similar logic, the RSFA and correlation coefficient are also differentially sensitive to partial-volume effects. Third, factors such as vascular position in a voxel also affect signal intensity fluctuations but not necessarily the correlation (X. Zhong & Chen, 2023). The basic concept stems from the fact that a given voxel may not be able to sample the entire dipole, but the details are beyond the scope of this paper. The same argument can be extended to the relative immunity of correlation coefficients to errors in estimating fBV and orientation from the segmented macrovasculature because of how they are normalized. Such errors will also have a greater impact on the fluctuation magnitude than the temporal features of the fluctuations.

In both arteries and veins, our model was unable to predict the in-vivo experimental RSFA at locations farther than one voxel away (~4 mm) from the macrovasculature. This may correspond to different number of voxels depending on the voxel size. A potential explanation for the unsatisfactory predictions could be found in the RSFA maps (Fig. S6 in Supplementary Materials), which suggest that the macrovascular susceptibility effect may only extend, or be predictable, over a limited distance from the macrovasculature. Consequently, in ROIs with limited macrovascular susceptibility effects, the in-vivo signal may be dominated by other sources of fluctuations (e.g., neural activity and physiological noise), which cannot be reflected by the predicted RSFA.

From the perspective of correlations, we found that in-vivo venous–venous correlations, while extremely well predicted by the macro-VAN Model, also became substantially less well predicted at one voxel away from the macrovasculature. This finding is similar to that of a previous study that found that the contribution of the veins to connectivity decreased with distance from the venous vasculature (Huck et al., 2023). Arterial–arterial correlations were not well predicted by the simulations even in the voxels containing the vessels.

One of the important questions raised by our previous studies is whether perivascular FC is related to EV magnetic susceptibility or to contributions of undetected vessels. As discussed in the previous section, the magnetic susceptibility difference between blood and tissue would generate an extravenous dipolar magnetic field offset around the macrovasculature that could extend well beyond the voxel containing the vessel (Ogawa, Menon, et al., 1993). We observed the susceptibility effect become attenuated rapidly beyond the macrovascular voxel (Fig. S6 in Supplementary Materials), unlikely to persist more than 4 mm away from the macrovasculature (even when thresholded a 1% of the maximum RSFA). Our regression results, for both RSFA and FC, support this claim. Such would be the case whether the macrovascular R₂’ perturbations are due to neuronal activity or other physiological processes. Of course, some physiological noise effects are global (Birn et al., 2006; Chang et al., 2009; Tong et al., 2015); however, they would not cause the model performance to deteriorate with increasing distance from the macrovasculature. Since the latter is what we see, we propose that EV susceptibility effects are not the sole contributor to perivascular BOLD. For one thing, there are numerous medium-to-small blood vessels in the extra-macrovascular space that are not detected by the TOF acquisition due to their smaller size and lower flow velocity. Their unmodelled BOLD contributions, which are likely also modulating the BOLD signal behaviour more distal to the microvasculature, may be the main reason that our modelling accuracy deteriorates with increasing distance from the macrovasculature. Our future work will target the improvement of macrovascular detection sensitivity, perhaps through the use of contrast agents (Bernier et al., 2020).

This study aims to demonstrate that modelling the macrovascular BOLD FC (in terms of correlation coefficients) using a biophysical model is feasible, and can provide a post-acquisition means to reduce the macrovascular bias during data analysis and interpretation in rs-fMRI. This study demonstrates, in concept, that TOF data acquired within rs-fMRI sessions can be used to predict the macrovascular bias, as we suggested in our previous study (X. Z. Zhong et al., 2024). The use of a full macrovascular VAN (the macro-VAN Model in this work) through numerical field-offset calculations demonstrated the highest performance in predicting venous-driven macrovascular FC, and thus should be the starting point for future efforts to eliminate macrovascular contribution on FC. To address the high computational cost of such an approach, pre-computed look-up tables could be developed that encapsulate all possible simulated signals, thereby increasing the accessibility of our approach. However, the prediction of arterial macrovascular BOLD is poor, potentially calling for better methods of capturing arterial dynamics. Lastly, perivascular BOLD contributions by large veins could be modelled reasonably accurately at ~4 mm beyond the vascular voxel, but not more. This could speak of the need to use more sensitive TOF imaging and higher spatial resolution to model perivascular effects from smaller vessels.

Based on the current results, adjusting FC values based on simulated FC would be the most promising way to reduce the macrovascular effect on resting-state fMRI. That is, macrovascular contributions to FC could be accurately predicted by our macro-VAN Model. It is, however, important to realize that such an FC correction is currently only realizable on the cross-correlation coefficients (instead of Pearson’s correlation, e.g.) between voxels in the vicinity of large veins. This is because our models do not account for lags between times series and their effects on correlation coefficients. Besides correlations, there may also be other, lag-insensitive approaches for calculating FC, such as regression and independent component analysis (ICA), but adapting our investigation to these other metrics will be the subject of our future work.

This study has several limitations. This report will not discuss the limitations of in-vivo experiments, which are discussed in our previous work (X. Z. Zhong et al., 2024). We will instead focus on the limitations most relevant to the simulations.

First, the comprehensiveness of the simulations is limited by the sensitivity of the TOF. Thus, we are unable to model the entire macrovasculature, which makes the input macrovasculature imperfect for fully accurate simulations and will reduce prediction performance. Moreover, the accuracy of vascular fBV estimates can be biased by differences in flow rates across vessels. Additionally, the simulation is based on simple sinusoidal signals as inputs to mimic resting-state fluctuations, which is likely not realistic (for either fBV or Y). Also, certain parameters, particularly for the 2D and 3D Cylinder Models, were shown to be important but were difficult to extract from TOF images, including the position of the vessel within the fMRI voxel (X. Zhong & Chen, 2023). Moreover, there have been various definitions of “macrovasculature” in the literature. Recognizing that not all venous BOLD effects are artefactual, we will adhere to the definition of “macrovasculature” outlined in Menon (2012) In our specific case, “macrovascular” refers to vessels detectable from our TOF MRA data, which are far larger in diameter than the intracortical veins. Furthermore, the TOF preprocessing procedures, such as brain extraction, may have led to some loss of macrovascular voxels, such as at the periphery of the prefrontal cortex. Future analysis will target ways to maximize vascular segmentation completeness using improved tools.

Secondly, due to the nature of TOF contrast, it is impossible to rely on the directionality of flow to distinguish between arteries and veins. Moreover, due to the proximity of some arteries and veins (such as the Circle of Willis and the Cavernous Sinus), partial-volume effects cannot be completely avoided. This limitation is mentioned in the Limitations section. While this adds to the challenge in identifying macrovasculature, we wish to point out that (1) the majority of venous data points are from identifiable veins such as the sagittal sinus and transverse sinus; (2) the venous BOLD effect is much stronger than the arterial BOLD effect, likely contributing to the high predictability from “venous voxels” compared with arterial ones in this work despite potential spatial overlaps between arteries and veins in some regions.

Thirdly, the availability of computational resources also limits the feasibility of the simulation. In spite of having more than 180 GB of memory, our server can only support 0.2-mm resolution for the macro-VAN Model simulation, which is not sufficient for some narrow segments of the macrovasculature. The long computational times for the simulations are another factor, which only allowed us to include four participants within a reasonable time frame. There is a possibility that this problem may be alleviated by the addition of computational resources in the near future.

Lastly, the applicability of these results likely depends on scan parameters (including spatial resolution, main magnetic field, and temporal resolution) of the MSC data set; therefore, comparisons with previous studies, especially those with higher main magnetic fields, should be made with caution. In practice, predictions should be tailored to different scanning protocols.

To conclude, we found that the macrovascular FC could be predicted using biophysical modelling based on vascular anatomical information, and that the venous contributions could be more accurately predicted than the arterial contributions. In addition, we found that the models were less reliable in predicting macrovascular RSFA than FC. Moreover, we found that modelling perivascular BOLD effects using simulations is feasible at a distance close to macrovasculature, also more accurately for veins than for arteries. This study paves a path for model-based correction of the macrovascular bias in resting-state and other types of fMRI, the validation of which will be part of our future work.

All data described in this manuscript were obtained from the Midnight Scan Club (MSC) data set. The data are available publicly at OpenNeuro. The code will be made available upon request.

Xiaole (Zachary) Zhong: Conceptualization, methodology, software, investigation, writing—original draft preparation. Jonathan R. Polimeni: writing—reviewing and editing. J. Jean Chen: Conceptualization, methodology, software, investigation, writing–original draft preparation, supervision, writing—reviewing and editing.

The study was approved by Baycrest REB (#11-47).

None.

The authors would like to acknowledge financial support from Canadian Institutes of Health Research and the Canada Research Chairs Program (JJC), funding support from NIH NIBIB grant P41-EB030006, NIMH grant R01-MH111419, and NINDS grant U19-NS123717 (JRP), and funding support from Ydessa Hendeles Graduate Scholarship (XZZ).

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