A mechanistic computational framework to investigate the hemodynamic fingerprint of the blood oxygenation level‐dependent signal

Blood oxygenation level‐dependent (BOLD) functional magnetic resonance imaging (fMRI) is one of the most used imaging techniques to map brain activity or to obtain clinical information about human cortical vasculature, in both healthy and disease conditions. Nevertheless, BOLD fMRI is an indirect measurement of brain functioning triggered by neurovascular coupling. The origin of the BOLD signal is quite complex, and the signal formation thus depends, among other factors, on the topology of the cortical vasculature and the associated hemodynamic changes. To understand the hemodynamic evolution of the BOLD signal response in humans, it is beneficial to have a computational framework available that virtually resembles the human cortical vasculature, and simulates hemodynamic changes and corresponding MRI signal changes via interactions of intrinsic biophysical and magnetic properties of the tissues. To this end, we have developed a mechanistic computational framework that simulates the hemodynamic fingerprint of the BOLD signal based on a statistically defined, three‐dimensional, vascular model that approaches the human cortical vascular architecture. The microvasculature is approximated through a Voronoi tessellation method and the macrovasculature is adapted from two‐photon microscopy mice data. Using this computational framework, we simulated hemodynamic changes—cerebral blood flow, cerebral blood volume, and blood oxygen saturation—induced by virtual arterial dilation. Then we computed local magnetic field disturbances generated by the vascular topology and the corresponding blood oxygen saturation changes. This mechanistic computational framework also considers the intrinsic biophysical and magnetic properties of nearby tissue, such as water diffusion and relaxation properties, resulting in a dynamic BOLD signal response. The proposed mechanistic computational framework provides an integrated biophysical model that can offer better insights regarding the spatial and temporal properties of the BOLD signal changes.


| INTRODUCTION
Blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) is a powerful noninvasive tool that employs susceptibility-sensitive imaging techniques to map in vivo brain functioning. 1However, the BOLD signal is only an indirect measurement of neuronal activity.0][21] Furthermore, computational approaches have provided better understanding of MRI signal features related to the biophysical interactions between magnetized water molecules that freely diffuse within the tissue and the magnetic susceptibility-induced effects produced by oxygenation changes in the vasculature at the mesoscopic level. 22[25][26][27][28][29][30] As a first attempt, numerical simulations of the BOLD signal using nonrealistic vascular representations, such as spheres or cylinders, have described static BOLD signal characteristics.[33][34] On the other hand, computational simulations using realistic 3D vascular models obtained from the parietal cortex of mice, acquired with two-photon microscopy, have demonstrated that the vascular architecture and the vessel/cortex orientation with respect to the main magnetic field produce a significant effect on the amplitude of the measured BOLD signal. 35,36The latter has also been observed in experimental data. 37,38dent models, however, may not be reflecting the human brain because of differences in the vascular architecture between the species.Angioarchitectural differences appear in the artery/vein ratio that feed/drain a specific volumetric region: reports based on histology show an artery/vein ratio of about 3:1 in humans as opposed to 1:3 in mice. 30,39,40These differences might lead to quantitative errors on the computed BOLD signal for the human brain, and thus to a misinterpretation of the data acquired in humans.
To achieve a comprehensive understanding of the spatial and temporal hemodynamic changes of the BOLD signal in the human brain, it is necessary to adopt a computational framework that (1) enables one to simulate the human cortical vascular architecture; (2) mimics hemodynamic changes within this simulated vascular network; and (3) considers the intrinsic biophysical and magnetic properties of the tissues and the influence of the MRI pulse sequence parameters.
Here, we present a computational framework that permits simulating the hemodynamic fingerprint of the BOLD signal for gradient-echo (GE) and spin-echo (SE) based on a statisticallydefined 3D model approximating the human cortical vasculature and associated hemodynamics.

| MATERIALS AND METHODS
The mechanistic computational framework consists of three modules: (1) generation of a statistically defined 3D vascular model (SVM) approximating the human cortical vasculature; (2) simulation of hemodynamic changes, for example, cerebral blood flow (ΔCBF), cerebral blood volume (ΔCBV), and blood oxygen saturation (ΔSO 2 ), which occur in the simulated vascular network in response to virtual arterial dilation; and (3) simulation of the biophysical effects induced by SO 2 changes and intrinsic properties of the tissue (diffusion and magnetic susceptibility changes), and the MR pulse sequence.

| Generation of the SVM
The SVM was generated in two different steps, reflecting the microvascular (Figure 1A) and macrovascular (Figure 1B) networks.2][43][44] The macrovasculature was mathematically adapted from mice to reflect the human artery/vein ratio.This was carried out to simulate the impact of a realistic macroarchitecture on the magnetic field disturbances.The two networks were spatially superimposed to yield the SVM (Figure 1C), but they were not physically connected for hemodynamic simulations (details of this assumption are provided in Section 2.2.3).The generation of each vessel subpopulation is described below.

| Generation of the microvascular network by means of Voronoi tessellation
We defined an isotropic simulation space that resembles the MR imaging voxel.For convenience, we constrained the size of the simulation space to a 1 mm 3 isotropic space, as normally used in high-resolution fMRI experiments.For simplicity, we assumed that within this voxel a large portion of the cortical vasculature, extending from the pial vessels to the gray matter (GM)/white matter (WM) boundary, was covered.
The microvasculature was simulated by means of a 3D Voronoi tessellation.Theoretically, it has been demonstrated that Voronoi tessellations give a good representation of the regions supplied by the capillary bed. 44A 3D Voronoi tessellation was generated by fragmentation of the volume-or voxel-into subvolumes that encompass a given set of seed points.Specifically, a Voronoi tessellation divides the simulation volume into spatially optimized polyhedra (Figure 1A).The Voronoi tessellation was computed using an adapted version of the polytope-bounded Voronoi diagram algorithm, 45 which uses linear inequalities formed with perpendicular bisectors between any two connected points in the Delaunay triangulation.Here, the number of polyhedra inside the simulation voxel was determined by the resemblance between statistical density distributions of the generated microvascular network-using the number of connections between the vessel's joints, the vessel's length and volume fraction-and the corresponding statistical density distributions of the human microvasculature obtained from histological data [41][42][43][44] (see also Figure 1A 39 The macrovascular compartment has an artery/vein ratio (3:1) obtained from anatomical studies in human cortex. 41(C) Representative SVM approximating the human cortical vasculature.
The edges of the polyhedra simulated the microvascular network.Each edge was simulated as a cylinder and the length of each edge simulated the vessel's length; the vessel length was determined as the Euclidean distance between joint connections (i.e., the joint connections generated by the polyhedra [small green dots in Figure 1A]).Further, a vessel radius was imposed on each polyhedron edge.The value of the vessel radius was randomly selected from a normal distribution with mean = 6.2 μm and std = 1.2 μm, which describes the range of microvessel radii observed in the human brain. 42,43The overall microvascular network occupied $4% of the total volume fraction of the simulation volume.Simulation parameters of a representative microvascular SVM compartment are summarized in Table 1.
The complete structure of the simulated microvascular network was defined by a set of vectors and an adjacency matrix.An m Â 3 vector describes the spatial coordinates of the microvessel joints inside the voxel space.The connectivity of the microvascular network was defined by an m Â m matrix M, where each element M ij is equal to one if the node i is connected to the node j, or is otherwise zero.The spatial information of each microvessel was stored in an n Â 6 information vector (spatial position and orientation direction of the cylinder).The vessel radius was stored as a column-vector and the vessel length value was stored as an additional column-vector.

| Spatial superposition of the macrovascular architecture
In the present work, pial vasculature as well as penetrating arteries and ascending veins were extracted from the 3D vascular architecture obtained by Blinder et al. 39 acquired with two-photon microscopy from the parietal cortex of mice.The mouse macrovascular compartment was segmented based on a vessel radius threshold.Vessels larger than 10 μm in radius were selected and labeled as the macrovascular compartment (Figure 1B).Then it was classified by randomly labeling three-quarters of the vessels as arteries and one-quarter of the vessels as veins, to approach human brain characteristics according to the artery/vein ratio in humans of 3:1, 41 as opposed to 1:3 in mice. 30,39Each vessel was then simulated as a cylinder.Vessel radius and length are parameters given by the mouse vascular model (≥ 10-60 and ≥ 70-500 μm, respectively; see Table 1 39 ).The modified macrovascular compartment consisted of $2% of the total volume fraction in the simulation volume.
The topology of the resulting macrovascular network was defined by vectors, similarly to the microvasculature.The spatial information corresponding to each artery or vein segment was stored in a p Â 6 information vector (spatial position and orientation/direction of the cylinder), and vessel radius and length were stored as column-vectors.
Finally, the labeled macrovascular compartment was spatially superimposed on the generated Voronoi microvasculature (as depicted in Figure 1C) to yield a representative SVM approximating the human cortical vasculature.

| Hemodynamic simulation using the SVM
Hemodynamics were simulated based on the principle of conservation of energy.We computed ΔCBF, ΔCBV, and ΔSO 2 considering the rheological properties of oxygenated blood and the anatomical properties of the vessels, such as geometry, resistance, and compliance.

| Microvascular hemodynamics at basal state
The microvascular network was converted into an electric circuit representation-resistor-capacitor circuit (RC circuit)-to simulate the transient hemodynamics in the network resulting from the similarities between an electric circuit and fluid dynamics. 7,46A B L E 1 Simulation parameters used to yield a representative statistically defined 3D vascular model approximating the human cortical vasculature.

Microvessels Macrovessels
Vessel radius Gaussian distribution mean: 6.2 μm std: 1.2 μm We considered that the Hagen-Poiseuille law describes the microvessel resistance (R) which depends on the radius of the microvessel (r), the length of the microvessel (l), and the viscosity of blood (η).The compliance describes the vessel's capacity to expand and the delayed return of the vessel's diameter to its basal state (i.e., the damping effect 46 ).In this work, the compliance of each vessel was modeled as a capacitance in the electric circuit and was randomly selected from a normal distribution (mean = 8 μF, std = 2 μF).The combined vessel compliance and resistance describe the dilation or constriction behavior of the vessel (i.e., the balloon effect). 3e simulation took account of the Fahraeus-Lindqvist effect that relates the viscosity of a fluid to the radius of the microvessel and the hematocrit (HcT) level. 18,47We implemented a vessel-dependent HcT level subject to the anatomical properties of the microvascular network. 48T ranged from 0.25 to 0.5, dependent on the vessel's size. 49ygenated blood was supplied through the microvascular network by three virtual arteries (input blood source) and was drained by one virtual vein (output blood source).The spatial position of the virtual input/output blood sources in the microvascular network was chosen arbitrarily, that is, the virtual blood sources were connected in an arbitrary node, which was randomly chosen from the microvascular network.It is worth noting that the virtual inlet/outlet blood sources are not directly connected to the spatially superimposed macrovascular compartment (see details in on 2.2.3).The consequences of this gross assumption are commented upon in the Discussion section.The blood pressure of the three virtual arteries was set to 75 mmHg, and for the virtual vein it was set to 25 mmHg. 50e basal-state cerebral blood flow (CBF) and blood pressure (BP) for each microvessel across the network were calculated based on the conservation of energy using: which describes the relation between the blood pressure in the inlet/outlet source points (ΔBP) and the microvascular resistance (for specific details, see 7,49,51,52 see also Table 2).
Note that in the current SVM implementation, the hemodynamic changes only occur in the microvascular network.

| Blood volume, blood flow, and oxygenation changes as a result of virtual arterial dilation
The three virtual arteries and one virtual vein were spatially connected at random to the microvascular SVM network.The vessel properties of these input/output blood sources were defined similarly as Boas et al. 7 input/output blood sources.
The virtual arterial dilation, which effectively changes the virtual arterial resistance, was modeled as a Gaussian function with properties of timeto peak (TTP) at 4 s and a full width at half maximum (FWHM) of 1.2 s. 53 We simulated a 10% change dilation compared with basal state. 7e computational framework enables tracking of the differential hemodynamic state for each vessel (Figure 2).Simulation results are shown for a global hemodynamic evolution within the voxel reflected by the dilation of all three virtual feeding arteries.
The ΔSO 2 for each microvessel was computed relating the balance between oxygen flow into the i-th microvessel via blood flow and oxygen diffusion into the tissue: T A B L E 2 Simulation parameters used to compute the hemodynamic basal state and activation transition.where SO 2,in,iÀvessel and SO 2,out,iÀvessel are each the oxygenation content (in μmol) carried in and out of the i-th vessel, respectively, and t is the simulation time step (200 μs).The oxygenation parameter in the virtual arterial sources was set constant to 97%.OEF is the oxygen extraction fraction from the i-th vascular segment into the tissue.For a detailed description of the oxygen dynamics, see the Oxygen dynamics section in Boas et al., 7 that is, the vascular anatomical network (VAN) model.
We assumed that tissue maintained a spatially homogeneous cerebral metabolic rate of oxygen (CMRO 2 ) inside the simulation volume that depends on its partial pressure of oxygen.We assumed a fixed partial pressure of oxygen in tissue of 15 mmHg. 50

| Blood volume and oxygenation changes in the macrovascular compartment
As described earlier, the virtual inlet/outlet blood sources only supply blood to the microvascular network.However, the values of the virtual inlet/outlet blood sources control the macrovascular compartment's CBV and SO 2 behavior, influencing the MRI signal.The spatial superposition of the macrovascular compartment only has an impact on the local magnetic field calculation, and not directly on the hemodynamic simulation.
Thus, the ΔSO 2 computed at the microvascular network's virtual inlets (three arteries) was used to control the oxygen saturation level of the entire arterial macrovascular compartment (red vasculature in Figure 1B).Moreover, the virtual arterial dilation controlled the vessel diameter of the entire arterial macrovascular compartment and the microvascular network, which corresponds to the arterial ΔCBV.We assumed that the arterial macrovascular compartment presents a constant oxygen saturation level (SO 2 = 97%), given that arteries do not show a significant change in blood oxygenation. 54milarly, the ΔSO 2 computed at the virtual outlet (one vein) was used to control the oxygen saturation level of the entire venous macrovascular compartment (blue vasculature in Figure 1B).We assumed that veins do not present volume changes (venous ΔCBV = 0), but only the changes in oxygen saturation computed at the outlet blood source (venous ΔSO 2 ).

| Simulation of the BOLD signal formation
The BOLD signal was then computed by the interaction of moving spins within the local magnetic field distortions. 32For each hemodynamic time step (200 μs), we computed local magnetic field changes created by the SO 2 levels of the macrovascular and microvascular compartments.Magnetic field disturbance caused by a segment of the macrovascular or microvascular compartment was modeled as the dipolar magnetic response of a finite cylinder, presuming negligible effects on the cylinder's extremities. 35The local magnetic field offset for each vessel segment was computed using: where γ is the gyromagnetic ratio, Þis the susceptibility difference produced by the SO 2 t ð Þ (Equation 3) and the hematocrit level (HcT), B 0 is the main magnetic field (7 T), R is the vessel radius, d is the Euclidean distance from the center line of the cylinder to a particular spatial position in the simulation volume, θ is the angle between the cylinder and the spatial position, and φ is the angle between the orientation of the cylinder and the main magnetic field.
The dephasing experienced by a bulk of diffusing spins was simulated using a Monte Carlo approach (1 x 10 10 spins).An ensemble of spins senses the extravascular local magnetic field with a diffusion coefficient of 1 μm 2 /ms.The calculation of the dephasing was obtained in the volume using: where ϑ t ð Þ is the phase acquired during the echo time (TE) and ΔB 0 x, t ð Þ is the local magnetic field offset at position x at time t.The dephasing of spins was then integrated over the TE.We assumed isotropic motion of spins. 32Further, we considered that exchange of spins between vascular compartments was prohibited (i.e., an impermeable vascular network).
To contain all spins inside the simulation space, voxel boundary conditions were established as an infinite space: if a spin moves out of the voxel, it is then mirrored on the opposite side of the simulation space, preserving the magnetization history.
Simulations shown here were computed for GE (TE = 27 ms) and SE (TE = 45 ms) at 7 T with a time step of 200 μs for a total simulation time of 20 s with an angular orientation of the normal vector of the SVM parallel to the main magnetic field (φ = 0 degrees).We used tissue relaxation times according to the nonlinear relationship given by Khajehim and Nasiraei Moghaddam 55 and Uluda g et al. 20 for the GM at 7 T (Table 3).
Additionally, to demonstrate the impact of changes in SO 2 and CBV on the BOLD signal, and the contribution of the microvasculature, we simulated different conditions where either ΔSO 2 or ΔCBV was constant for a representative SVM.
The computational pipeline was implemented in Matlab/Julia.The biophysical parameters used to compute the representative dynamic BOLD response are summarized in Table 3.
The computational framework allows the positioning of the main magnetic field at any spatial orientation.We defined the angular orientation (φ = 0 degrees) of the SVM to be parallel to the main magnetic field (z-axis) and the normal vector that rises from the cortical surface.The cortical surface does not present any curvature along the simulation volume, that is, pial arteries and veins are orthogonal to the main magnetic field. 37,38 demonstrate the impact of the angular orientation of the human vasculature with respect to the main magnetic field on the amplitude of the BOLD signal, additionally, we computed the local magnetic field distortions at two exemplary time points, that is, a representative "active" (virtual dilation) and basal state.
These local magnetic field calculations assumed, in the basal state, an SO 2 level in the arteries of 95%, microvessels of 75%, and veins of 68%; and in the "active" state an SO 2 level of 98%, 85%, and 75%, respectively.The MRI signal was calculated using the inverse Fourier transform of the computed local frequency components for both oxygenation states. 56This type of calculation assumes a GE pulse sequence and no diffusion effects (i.e., a static diffusion regime).
To emphasize the effect of a particular vascular compartment, we computed the local magnetic field distortion for: 1. (I + II) the total microvasculature: pial arteries, penetrating arteries, pial veins, and ascending veins.
2. (II) the cortical penetrating arteries and ascending veins.the input assumptions of our model.The macrovascular compartment (Figure 1B) is spatially superimposed on the microvascular network, resulting in a representative SVM consisting of $2500 microvessels (Figure 1C).

| Hemodynamic simulation of virtual arterial dilation and simulation of the BOLD signal using the SVM
Figure 2 shows the hemodynamic response of the microvascular network to virtual arterial dilation.Upon arterial dilation, the spatial hemodynamics in the microvascular network are affected by the pattern, depending on the location of the virtual inlets/outlet blood sources and the layout of the microvascular network.
We show eight different time points (A-H), where the change in blood flow on the network is depicted on the left side of each panel, and the blood pressure on the right side of each panel.The increase in blood flow follows the path of minimum resistance from the virtual feeding arteries to the virtual draining vein, determined by vessel resistance, compliance, and topology.The blood pressure remained locally constant across the microvascular network at vessel joints; it also depended on the vessel resistance, compliance, and topology, as well as the positioning of the inlets/outlet sources.
In Figure 3, we show the hemodynamic response, that is, ΔBP (A), ΔCBF (B), and ΔCBV (C), of the representative SVM network in response to virtual arterial dilation.The virtual arterial dilation modifies the vessel resistance of the vessel network (D).Using this hemodynamic behavior, a dynamic BOLD response is computed using Monte Carlo approaches.GE BOLD and SE BOLD responses are displayed in (E) and (F), respectively.
In (A) (ΔBP), as the virtual arterial dilation takes place between the second and the fourth second, the microvascular network linearly increases the volume of the vessels according to this dilation (e.g., there is no lag behavior upon activation in these simulations in the vessels), thus maintaining the blood pressure equally distributed within the microvascular network.Then the principle of conservation of mass is demonstrated by the shaded region, which represents the effect on the blood pressure changes, as induced by the differences in microvascular compliance.The exact configuration of this shape is influenced by the spatial arrangement of the inlet and outlet sources and the vascular properties of the network.The blood pressure reverts to its baseline state once the entire capillary network has also returned to its baseline state, that is, the balloon effects return to baseline state (D: the effect after the fourth second of simulation).
Given the simulation parameters already described, we obtain an $5% BOLD signal difference in contrast to the basal state in a GE pulse sequence; and an $2% signal change under a SE pulse sequence.This difference is related to the refocus dephasing of a particular vessel's contribution to the MRI signal in a SE pulse sequence.
In Figure 4 the large pial veins.Microvascular structures are poorly discernable in the maps, although their oxygen saturation does contribute to the local magnetic field changes in the order of $10-30 Hz.
In Figure 5, we show the dynamic BOLD signal response for GE and SE pulse sequences.The average over the simulated signals (n = 6) is depicted by the emphasized dots in the BOLD signal response, while the standard deviation is illustrated by the shaded region.To demonstrate the effect of changes in SO 2 and CBV, and the contribution of the microvasculature, we simulate different conditions: (1) a complete SVM (A-F), that is, microvascular and macrovascular compartments; and (2) a SVM composed only of the microvascular network (G-L).
We observe that the change in oxygen saturation drives the positive amplitude of the BOLD signal change for both the complete SVM (C and F) and the microvasculature (I and L).The amplitude of the BOLD signal change is dependent on the vascular contribution, being larger for a complete SVM compared with only microvessels.It is important to highlight that GE and SE exhibit distinct standard deviations, with SE demonstrating lower variability in the computed signals.
The negative %BOLD signal corresponds to CBV changes (B, E, H, and K) in the simulated signals.Moreover, the contribution of the large arteries in (B) has a large impact on the signal.It is important to mention that, in these simulations, the veins were not dilating, based upon the notion that veins do not significantly expand during functional hyperemia. 30re, the TTP depends on the vascular resistance and compliance and appears approximately some milliseconds later than the virtual arterial dilation.The FWHM of the BOLD signal in the positive part depends not only on the FWHM of the arterial dilation, but also on the compliance of the simulated vasculature.The FWHM of the simulated BOLD response in the microvascular-only model is narrower than that for the complete model.
To demonstrate the effect of the macrovascular artery/vein ratio on the local magnetic field distortion, we show the macrovascular magnetic field distortion induced for an artery/vein ratio of 3:1 (humans) and 1:3 (mice) in Figure 6.Here, a 95% and 65% SO 2 level was used in the arterial and venous compartment, respectively.We obtained an $21% difference in the distribution of magnetic field changes between the two artery/ vein ratio scenarios, with the 1:3 ratio (mice) yielding the largest magnetic field changes.
In Figure 7 we show the angular orientation dependence of the BOLD signal change obtained with the SVM with respect to the main magnetic field using four different macrovascular topologies. 39The previously reported cos 2 (φ) behavior [35][36][37][38] is demonstrated: when the pial arteries and veins are included (panels I + II and I + II + III), the reported behavior is in phase with the simulations; when the pial arteries and veins are not included (panels II, III, and II + III), the reported behavior shows a -π/2 shift in phase compared with the simulations.Moreover, the %BOLD of the microvasculature shows a negligible angular dependence (III).The contribution to the BOLD signal of the different vascular compartments can be an order of magnitude different.For instance, the computed %BOLD signal change when the pial vessels are included in the simulation is approximately four times larger compared with when only penetrating vessels contribute to BOLD signal formation.Importantly, the microvasculature shows a very small BOLD signal change compared with the signal change induced by the macrovessels (panels II and II + III), almost by a factor of $20.

| General discussion
To understand the hemodynamic evolution of the BOLD signal response in humans, it is beneficial to have a computational framework available that virtually resembles the human cortical vasculature, and simulates hemodynamic changes and corresponding MRI signal change via interactions of intrinsic biophysical and magnetic properties of the tissues.We have presented a computational framework that permits simulations of the hemodynamic fingerprint of the BOLD signal for GE and SE based on a statistically defined 3D model that approximates the human cortical vasculature, and we simulate the associated hemodynamics in response to virtual arterial dilation, assuming the intrinsic biophysical properties of the tissues.Our work extends previous computational frameworks by enabling simulations of hemodynamics and the BOLD signal evolution in a realistic 3D model approximating the human cortical vasculature.These points are discussed below.

| On the SVM
The proposed SVM approximates the geometrical and topological characteristics of the human cortical vasculature.In addition, this framework enables mimicking of the temporal evolution of hemodynamic changes in a typical voxel used in high-resolution fMRI acquisitions. 57Moreover, it simulates the intrinsic biophysical and magnetic properties of the tissues, such as diffusion of water molecules and local magnetic field distortions created by oxygen saturation changes, and the MRI pulse sequence.
The use of a 3D model gives valuable insights regarding the spatial contribution of the vascular topology on the MRI signal formation, which cannot be reflected with two-dimensional (2D) approaches. 10,11,26Further, the use of a statistical vascular model, as implemented here, circumvents assumptions often adopted in previous models, such as a uniform angular distribution of vessels and/or a mono-sized vessel radius population. 19,20uda g and Blinder 30 and Schmid et al. 58 made an extensive analysis of vascular architecture differences between species.An important difference between species lies in the arterial/venous ratio, which was the main difference between mice and humans (1:3 in mice vs. 3:1 in humans).However, the general macrovascular morphology and topology appear to be similar across species. 59This topological analysis supports our approach of spatially superpositioning the experimental macrovascular data of Blinder et al. 39 on top of a random microvascular bed after relabeling larger vessels to fulfill the human artery/vein ratio.A difference of $21% on the frequency components was obtained by simulation of the local magnetic field distortions between the human and mouse models for the same constant blood properties (arteries = 95% oxygenation level, veins = 65% oxygenation level, hematocrit = 40%).
Our approach provides a solution to mitigate the need for an experimental 3D reconstruction of the human vasculature, which is currently difficult to obtain, because the vasculature is statistically defined.[43]

| On the hemodynamics
The proposed mechanistic computational framework enables simulation of hemodynamic changes-CBV, CBF, and SO 2 -over time within the 3D vascular network.The SVM allows manipulation of the number of inlet/outlet blood sources and their spatial position inside the network.Moreover, the proposed computational framework permits simulations of different hemodynamic conditions, such as the effect of the compliance of the vessels and the HcT level dependence on vessel parameters that directly influence the CBF, CBV, and SO 2 changes upon dilatation/ activation.
Several models have been proposed to simulate the hemodynamic behavior of CBF and CBV changes in a 2D vascular network. 3,7,11,12Our SVM comprises all three compartments-arteries, capillaries, and veins-assuming morphological and spatial distributions described by statistical density characteristics.Thus, it offers the possibility to study hemodynamic behavior in a 3D network and the direct impact on the BOLD signal evolution.In general, this computational framework represents the initial stage towards a complete statistical model.Using a statistical approach provides an advantage when attempting to replicate variations in angioarchitecture and intrasubject/intersubject variability.Additionally, the presented framework offers a second benefit, as it allows for the independent study of microvessels and their hemodynamic/MR behavior, including the simulation of effects such as angulation with respect to B0 and the microvascular contribution to GE BOLD.However, a drawback of this approach is that the microvessels are not connected to larger vessels, which creates discontinuity in hemodynamics: inflow coming through large penetrating arteries, as well as draining effects in veins.
Hemodynamic simulations on the microvascular network are shown for three arteries (input blood source) and one vein (output blood source).
The simulations were performed using four different macrovascular topologies, where black dots represent the mean behavior, and the shaded areas indicate the standard deviation of the computed BOLD contrast.
The computation of the BOLD signal presented here assumes only the contribution of the extravascular space, given by the interaction of the diffusing spins and the susceptibility effects induced by the vascular compartments and the differential oxygenation level.The intravascular BOLD contribution is negligible because of the relatively short T2 relaxation time constant of oxygenated blood at ultrahigh magnetic fields. 20Further improvements of the simulation framework could entail inclusion of the intravascular contribution, which might (slightly) affect the hemodynamic signature of the BOLD signal. 63nally, an important issue to tackle is the validation of our model and results, which of course is an important next step to probe the reliability of our computational framework and the accuracy of the simulated BOLD signals compared with those obtained from the human brain.However, our objective is to showcase the feasibility of a more robust and flexible computational framework that can contribute to understanding the various interplays involved in the formation of the MRI signal and the evolution of the hemodynamic signature in the BOLD contrast.Validation can be achieved qualitatively through comparisons between simulated and experimental signals.For instance, within our research group, Siero et al. 64 reported the hemodynamic response function (HRF) characteristics of GE BOLD and SE BOLD at 7 T based on short duration evoked stimuli in the visual cortex, and the signals obtained in that study share BOLD amplitudes and dynamic shapes with the ones simulated.The angular amplitude-dependence of the BOLD signal with respect to B0 reported by Fracasso et al. 37 and Viessmann et al. 38 provides an additional qualitative validation of our computational framework.These references serve as support for the notion that our framework is heading in the right direction.However, we are actively working on conducting a quantitative validation that will help us narrow down the selection of broad parameters and provide substantial information about BOLD signal evolution.

| CONCLUSIONS
We present a mechanistic computational framework that enables simulations of the hemodynamic fingerprint of the BOLD signal to be performed, based on a SVM approximating the human cortical vasculature and associated hemodynamics.Further, this framework comprises an integrated computational method assuming intrinsic biophysical effects induced by specific oxygen saturation changes and tissue magnetic properties that results in a dynamic BOLD signal response.All these biophysical interactions provide a view of the relation between the vascular topology of the human brain and the associated hemodynamics that results in a dynamic BOLD signal response.This computational framework may help with the interpretation-quantitative and qualitative-of the BOLD signal, and provide better insights regarding the hemodynamic changes of the signal evolution.
.1-A.3).F I G U R E 1 Representative statisticallydefined three-dimensional vascular model (SVM).(A) The microvascular compartment generated through Voronoi tessellation.The coloring is chosen at random and does not represent any particular feature.Yellow dots represent the center of mass for each polyhedron.Blue dots represent the vessel joints in the vascular model.A statistical description of the microvascular compartment is illustrated in A.1-A.3: Density distributions of the vessel's (A.1) Radius, (A.2) Length, and (A.3) Orientation with respect to the cortical surface.(B) Spatial distribution of the macrovascular compartment (arteries in red and veins in blue) adapted from the mouse model acquired by Blinder et al.

F I G U R E 2
Representative hemodynamic transient on the microvasculature in response to virtual arterial dilation of 10%.The left-side cube for each panel displays the blood flow changes computed for each microvascular segment (dark = low cerebral blood flow [CBF], bright = high CBF).The right-side cube of each panel shows the blood pressure state across the microvascular network resulting from the conservation of energy and the position of the inlets/outlet sources (red = low pressure, yellow = high pressure).At the maximum dilation state (10%, panel D), the vascular network depicts higher blood flow changes, in contrast to the other states.

3 | RESULTS 3 . 1 |
Figure 1A displays a representative microvasculature generated through Voronoi tessellation.The coloring is chosen at random to enhance visualization and does not represent any particular feature.The yellow dots represent the center of mass for each polyhedron.The blue dots-vertices of the polyhedral-represent the microvessel joints.Using vector analysis and graph theory (ion 2.1.1),it is straightforward to characterize the geometric properties of the microvascular compartment, as shown in the density distributions in Figure 1A.1-A.3, which are in agreement with

F I G U R E 3
Hemodynamic response of the microvascular network in response to virtual arterial dilation.(A) Blood pressure changes at vessel nodes, (B) Blood flow changes, (C) Blood volume changes, and (D) Vessel resistance changes.The black line in (A-D) represents the mean value and the shaded area the standard deviation.(E and F) Blood oxygenation level-dependent (BOLD) signal computed assuming hemodynamic changes, magnetic tissue properties, and MR pulse sequence (E: gradient echo; F: spin echo).F I G U R E 4 Local magnetic field distortions produced by a representative statistically defined three-dimensional vascular model (SVM) (arteries in pink, veins in light blue, microvessels in yellow) for an exemplary hemodynamic time point: (F) in Figure 2. Magnetic field maps are shown for three different locations inside the simulated voxel: top, middle, and bottom planes.The main magnetic field was simulated to be orthogonal to the cortical surface and thus the depicted planes.F I G U R E 5 Dynamic blood oxygenation level-dependent (BOLD) signal response computed for a gradient-echo (GE) (echo time [TE] = 27 ms) and a spin-echo (SE) (TE = 45 ms) pulse sequence at 7 T for different hemodynamic conditions.(A-F) BOLD signal change produced by a complete statistically defined three-dimensional vascular model (SVM): macrovasculature and microvasculature.(G-L) BOLD signal change computed using only microvessels.The bold dots in the BOLD signal response represent the average over the simulated signals (n = 6), and the shaded area shows the standard deviation.The hemodynamic parameters can be found in

F I G U R E 6
Representative macrovascular architecture between (A) Human (simulated data), and (B) Mouse (experimental data), and (C and D) Frequency line-shape comparison.The arterial compartment is shown in red and the venous compartment in blue.(C) Frequency line-shape computed at 7 T for each model, and (D) Normalized frequency line-shape difference between human and mouse macrovascular compartments.

F I G U R E 7
Impact of the angular orientation of the statistically defined three-dimensional vascular model (SVM) on the blood oxygenation level-dependent (BOLD) signal/contrast change.We compute the %BOLD signal/contrast change for several scenarios: (I + II) Pial arteries and pial veins plus penetrating arteries and ascending veins (i.e., complete macrovasculature); (II) Only penetrating arteries and ascending veins; (III) Microvasculature; (I + II + III) The complete SVM; and (II + III) Penetrating arteries and ascending veins plus microvasculature.The x-axis of all plots shows the angular orientation (φ) between the main magnetic field and the normal vector of the cortical surface.Note that the scale of the BOLD signal change (y-axis) is notably different for the different vascular combinations.The light blue line represents a scaled reference to facilitate the visualization of the cos2(φ) obtained in the experimental data.