|Home | About | Journals | Submit | Contact Us | Français|
Neuronal activity-induced changes in vascular tone and oxygen consumption result in a dynamic evolution of blood flow, volume, and oxygenation. Functional neuroimaging techniques, such as functional magnetic resonance imaging, optical imaging, and PET, provide indirect measures of the neural-induced vascular dynamics driving the blood parameters. Models connecting changes in vascular tone and oxygen consumption to observed changes in the blood parameters are needed to guide more quantitative physiological interpretation of these functional neuroimaging modalities. Effective lumped-parameter vascular balloon and Windkessel models have been developed for this purpose, but the lumping of the complex vascular network into a series of arterioles, capillaries, and venules allows only qualitative interpretation. We have therefore developed a parallel vascular anatomical network (VAN) model based on microscopically measurable properties to improve quantitative interpretation of the vascular response. The model, derived from measured physical properties, predicts baseline blood pressure and oxygen saturation distributions and dynamic responses consistent with literature. Furthermore, the VAN model allows investigation of spatial features of the dynamic vascular and oxygen response to neuronal activity. We find that a passive surround negative vascular response (“negative BOLD”) is predicted, but that it underestimates recently observed surround negativity suggesting that additional active surround vasoconstriction is required to explain the experimental data.
The hemodynamic response to neuronal activation involves a complex interplay between processes, which leads to changes in cerebral perfusion and to changes in glycolytic metabolism and oxygen consumption. Vasoactive agents, released during neuroglial activation, interact with the vascular smooth muscle of arterioles, inducing changes in muscle tone and vascular diameter (Devor et al., 2007; Hamel, 2006; Iadecola, 2004; Toda and Okamura, 1998). Capillary diameter can also be controlled through activation of pericytes (Peppiatt et al., 2006). Arteriolar dilation, a common denominator in evoked hemodynamic response, alters the vascular resistance and results in an increased blood flow into the vascular network and increased pressure in the downstream vascular segments. The compliant capillaries and venules subsequently expand in response to this increased pressure with some debate as to whether dominant volume changes occur on the arteriole side (Hillman et al., 2007; Kim et al., 2007; Lee et al., 2001) or venous side (Buxton et al., 1998; Kong et al., 2004; Mandeville et al., 1999). Consequently, increased flow results in a greater delivery of oxygenated hemoglobin. This flow-induced dilution of deoxy-hemoglobin overbalances increases in oxygen consumption. Thus, oxygen extraction fraction is found to decrease (Buxton et al., 2004; Fox and Raichle, 1986; Hoge et al., 1999b).
Physiological interpretation of functional MRI (fMRI) or optical imaging measures of the hemodynamic response to brain activation requires a detailed understanding of the competing effects of increased perfusion and metabolism. Biologically plausible models of the cerebral vascular response have been described not only for determining functional changes in the cerebral metabolic rate of oxygen (CMRO2) (Boas et al., 2003; Hoge et al., 1999a; Kida et al., 2000; Zheng et al., 2002), but also for aiding in the interpretation of the characteristics of the fMRI blood oxygen level-dependent (BOLD) response (Mandeville et al., 1999; Obata et al., 2004) as reviewed in Buxton et al. (2004). The first such models to detail these interactions in the context of fMRI data were the vascular balloon (Buxton et al., 1998) and Windkessel (Mandeville et al., 1999) models, coupled with models of oxygen delivery including the diffusion limited delivery model (Buxton and Frank, 1997; Gjedde et al., 1999) and later extended to include the interstitial oxygen-buffered models (Herman et al., 2006; Zheng et al., 2002).
These initial models treated the cerebral vascular network as a serial pipe with a dilating arteriole and a passively responding venule. While this is a gross oversimplification, it does make predictions that agree reasonably well with experimental fMRI data (Deneux and Faugeras, 2006; Obata et al., 2004) and provides an interpretation of response characteristics such as the initial “dip” and post-stimulus under-shoot of the BOLD signal (Buxton et al., 2004) and the delayed recovery of blood volume (Mandeville et al., 1999). However, when compared to the immensely more complex biological vascular network, it is not surprising to find systematic biases in these oversimplifications, especially in the context of the rich information available from invasive functional imaging studies. Recent work by Zheng et al. (2005) has highlighted discrepancies between a simplified two-compartment model and experimental data acquired at high temporal resolution with invasive optical measurements of hemoglobin and blood flow changes. In order to better explain temporal differences observed between oxy- and deoxy-hemoglobin changes, Zheng et al. (2005) extended the earlier vascular and oxygen transport models by Mandeville et al. (1999) and Zheng et al. (2002) to a three-compartment model of an actively dilating arteriole and a passively responding capillary and venule. The compartmentalized model was reported to predict the oxygenation response as measured by optical imaging data (Zheng et al., 2005) and more recently was confirmed by Huppert et al. (2007), who showed that a three-compartment model statistically improved the fit to experimental data and provided more robust estimates of CMRO2 changes. These models provide predictions of the hemodynamic response to brain activation and, inversely, precise estimates of the oxygen consumption response from experimental data (Huppert et al., 2007; Zheng et al., 2005) and offer opportunities for better analysis of multimodal datasets (Riera et al., 2005).
While the goal is to develop the simplest model that effectively and accurately captures the correct behavior of the complex system, care must be taken in interpreting the lumped parameters of the simplistic model that describe the effective properties of the brain’s vasculature. The relation between these effective parameters and the true physical properties of the hemodynamic response is unknown, and thus the accuracy of any estimates obtained by these models is uncertain. In addition, the serial pipe models may not fully describe spatial characteristics of the hemodynamic response, such as the surround negativity response observed in whisker barrel cortex with optical imaging (Cox et al., 1993; Devor et al., 2007; Devor et al., 2005; Woolsey et al., 1996) and negative BOLD signals in the visual cortex (Harel et al., 2002; Shmuel et al., 2002). Although the recent development of multi-compartment models has improved the interpretation of the temporal dynamics of hemodynamic changes, the spatial characteristics of signals have not been explored.
In order to obtain a more detailed understanding of the spatio-temporal characteristics of the hemodynamic response to brain activation, and ultimately to obtain more accurate estimates of the underlying physiological factors driving the hemodynamic response, we have developed a parallel vascular anatomical network (VAN) model based on the physically measurable properties of the vessels, hemoglobin, and oxygen. Our VAN model is an extension of our previous three-compartment serial pipe Windkessel model (Huppert et al., 2007) to now include parallel arterioles, capillaries, and venules and enable model parameters to derive from physical properties of the vascular network including vessel segment diameters, lengths, and viscosity. Given this VAN model, we then predict baseline blood pressure and hemoglobin oxygen saturation consistent with literature values. Examining the dynamic response to local arteriole dilation, we observe a hemodynamic response consistent with experiments and the serial pipe models. Finally, we find that spatial surround negativity in the hemodynamic response (i.e., decreased blood flow, volume, and oxy-hemoglobin with increased deoxy-hemoglobin) arises naturally from the passive redistribution of pressure gradients in the vascular network, but does not fully account for surround negativity recently reported in the literature (Devor et al., 2007; Devor et al., 2005), suggesting a combination of passive and active processes.
This section describes our model of the microscopic vascular anatomy including the blood flow and volume dynamics and the oxygen delivery. We first cover the flow and volume aspects of the model, including parameter selection, and then address the oxygen delivery aspects of the model, including parameter selection.
The vascular anatomical network (VAN) model is built up from a branching series of individual arterioles, through the capillaries, and then on to a converging series of venules, each with its own characteristic properties. The simplest version of this is a series of a single arteriole, capillary, and venule as depicted in Fig. 1. This simple serial vascular pipe model has been developed and described by the Balloon (Buxton et al., 1998) and Windkessel (Mandeville et al., 1999) models and more recently extended to multiple vascular compartments (Huppert et al., 2007; Zheng et al., 2005).
As stated in the works of Huppert et al. (2007), Mandeville et al. (1999), and Zheng et al. (2005) for the serial pipe Windkessel models, the resistance of each vascular segment is related to the vessel diameter, d, length, l, and fluid viscosity, η, according to Poiseuille’s law (Washburn, 1921)
The arterial and venous pressures are fixed at the boundaries as PA and PV, respectively. Given the resistance of each vascular segment, the pressure and flow along the vascular network can be determined using the flow–pressure circuit relationship
where P1 is the pressure at the first segment, P2 is the pressure at the second segment, F12 is the flow from the first to the second segment, and R12 is the resistance between the two segments.
Neuronal activation gives rise to arteriolar dilation. Dilation results in decreased arteriolar resistance, as determined by Eq. (1), and increased flow, as determined by solving the circuit equations. The change in blood volume within a vessel segment per unit time is the differential of blood flow into and out of the vessel segment
where the subscript ‘i’ refers to the ith vessel segment.
Vessels generally have a non-linear compliance that decreases with increasing pressure. Following the example of Mandeville et al. (1999), this non-linear compliance produces a pressure–volume relationship of the form
for the arterioles, capillaries, and veins and the subscript ‘i’ refers to the ith vessel segment. PIC is the intra-cranial pressure, which we assume is constant. The non-linearity in the pressure–volume relationship is given by compliance parameter β ≥ 1 for diminished volume reserve at higher pressures. The constant Ao is given by the initial volume, Vo, and pressure, Po, of the ith segment
Extension of the serial pipe model to a network model requires adding branch nodes in which vessel segments diverge or converge. These branch nodes have no blood volume and are simply modeled by a conservation of flow equation, i.e., net flow out must equal net flow in, as illustrated in Fig. 1.
We are modeling the network of cortical vessels from the pial arterioles through the capillaries to the pial venules. We choose parameters to resemble the rat’s cerebral cortex. We use a simple balanced branching network from a single arteriole through 2N capillaries back to a single venule, with N =6 branches of arterioles and venules and a single branch of capillaries. The length of each arteriole and venule segment is fixed at 100 μm, and the capillaries have a length of 250 μm (Turner, 2002; Vovenko, 1999).
The capillary diameter is set to 8 μm (Kleinfeld et al., 1998; Lipowsky, 2005; Villringer et al., 1994). Arteriole diameter decreases by 20% at each branch, starting at 30.5 μm and decreasing to 10.0 μm just before the capillaries. While this is a simplistic model of the distribution of arteriole diameters, it was chosen to match the literature on the reduction of velocity from arteriole to capillary. Red blood cell velocity decreases from ~13 mm/s in 30 μm arterioles to ~2.5 mm/s in the capillaries, to a minimum of ~2 mm/s in the post-capillary venule, as measured in cat mesentery (Lipowsky, 2005; Zweifach, 1974). This requires that the net cross-sectional area of the capillaries be greater than the net cross-sectional area of the initial feeding arteriole. The velocity at each branch will decrease only if the diameter decreases by no more than 30%, otherwise the net cross sectional area decreases, leading to a velocity increase. A decrease from 13 mm/s to 2 mm/s over 6 branches requires an average diameter reduction of 20% per branch.
Venule diameters are set to increase by 20% per branch, with the post-capillary venules given a diameter of 12.0 μm, resulting in a final venule diameter of 36.6 μm, as venules are typically larger than arterioles (Vovenko, 1999).
The resistance of each vascular segment is dependent on its diameter, length, and apparent viscosity, as indicated by Poiseuille’s law (Eq. (1)). Viscosity increases monotonically with increasing vessel diameter and linearly with hematocrit, as summarized by Pries et al. (1992), and we assume a plasma viscosity of 2 cP (Haidekker et al., 2002). Microvessel hematocrit decreases with vessel diameter and ranges from ~14% in 30 μm arterioles to ~8% in capillaries (Lipowsky et al., 1980). All of these vascular segment-dependent values are summarized in Table 1.
The baseline parameters not associated with individual vascular segments are summarized in Table 2. The non-linear vascular compliance is given by β = 2 as representative of the range 1 < β <3 given in the literature (Dunn et al., 2005; Jones et al., 2001, 2002; Mandeville et al., 1999; Sheth et al., 2004b). Finally, PA=60 mmHg for the 30-μm arteriole blood pressure, PV =25 mmHg for the 35-μm venule blood pressure (Lipowsky, 2005), and the intracranial pressure, PIC =10 mmHg (Albeck et al., 1991; Alperin et al., 2000).
As illustrated in Fig. 2, the oxygen content in a vascular segment is a balance between oxygen flow into a segment via blood flow and out of the segment via blood flow and oxygen diffusion into the tissue. That is,
where SO2,in and SO2,out is the hemoglobin oxygen saturation flowing in and out of the segment, respectively, OE is the rate of oxygen extraction from the segment into the tissue in units of mol/s, and NHbO2 is the moles of oxygenated hemoglobin in the segment. The factor γB converts from microliters of blood to moles of hemoglobin, while the factor of 0.25 accounts for the fact that each molecule of hemoglobin can carry 4 molecules of oxygen. The molecular weight of hemoglobin is 64,500 g/mol. With a typical blood hemoglobin concentration of 150 g/l (Habler and Messmer, 1997), γB = 2.3 × 10−9 mol/μl. We assume the average oxygen saturation of each segment as
The quantity of oxygenated hemoglobin is related to the blood volume (units of μl) by
where NHbT is the amount of hemoglobin (oxygenated and deoxygenated) in the segment. Note that
The oxygen flux across the vessel wall can be expressed using Fick’s Law (Popel et al., 1989)
where PO2,tis,i and PO2,i are the partial pressure of oxygen in the tissue around and in the ith vascular segment, respectively; w is the thickness of the vascular wall; and KW is the oxygen permeability of the vascular wall. The oxygen extraction is just the flux times the surface area of the vessel wall
where the last term on the right converts the volume of consumed oxygen (at standard pressure and temperature of 1 atm and 273 K, respectively) to moles. Note that the relationship between PO2 and SO2 is given by the oxygen disassociation curve. For our modeling, we use the empirical fit given by Kelman (1966) and Severinghaus (1979). We combine Eqs. (11) and (12) to lump all of the coefficients into a single parameter that we call an oxygen extraction coefficient kw,i with units of mol/mmHg/s.
Finally, changes in the tissue partial pressure of oxygen are caused by differences in the delivery of oxygen to the tissue by OE and extraction by the cerebral metabolic rate of oxygen CMRO2, i.e.,
where α is the Bunsen solubility coefficient for oxygen, and VCMRO2 = πliρ2cmro2 is the volume in which the extracted oxygen is consumed.
Vovenko (1999) has made detailed measurements of the hemoglobin oxygenation in arterioles, capillaries, and venules of the rat cortex. We use Vovenko’s numbers to guide the initial condition for our model. Vovenko’s results show a significant drop in SO2 from the systemic arterial oxygen saturation of 94% to a pre-capillary oxygen saturation of 84%, indicating significant oxygen flux from the arterioles. A drop of ~25% is observed across the capillaries. Negligible change is observed throughout the venules.
Hemoglobin oxygenation in each vascular branch is a function of the blood flow, blood volume, and oxygen extraction. The flow and volume of each branch are set as described above. We thus must define the initial oxygen extraction from each segment to obtain the initial hemoglobin oxygenation. This is achieved by adjusting the effective oxygen extraction coefficient kw(Eq. (13)). To simplify, we use the same kw for all vessel segments independent of vessel diameter, length, and thickness. Given a vessel length of 250 μm, diameter of 8 μm, wall thickness of 2 μm, and permeability Kw = 5 × 10−8 μl/mm/s/mmHg (Popel et al., 1989), we find kw = 7 × 10−15 mol/mmHg/s. We use this value for all vessel segments.
The oxygen solubility coefficient and radial extent of oxygen consumption are defined in Table 2. As described in the works of Hudetz (1997) and Turner (2002), the volume fraction of capillaries is approximately 2% and thus each capillary supplies about 50 × its cross-sectional area with nutrients. Thus, for an 8-μm-diameter capillary, oxygen is delivered radially outward to a distance of approximately 30 μm.
Given the initial resistance of each vascular segment and the input and output pressure, the steady-state pressure distribution throughout the parallel vascular network is obtained by simultaneously solving the coupled linear equations arising from Eq. (2) with the boundary condition that flow in and out of a node point must be conserved. The steady-state flow is then calculated from Eq. (2). The steady-state solution for oxygen saturation is obtained from Eqs. (10) and (13) given the steady-state flow, input oxygen saturation, and partial pressure of oxygen in the tissue. These coupled differential equations are solved iteratively until convergence. Specifically, the time derivatives in Eq. (10) are set to 0. The oxygen saturation is set to the arteriole oxygen saturation at the input node and to 0 everywhere else. The partial pressure of oxygen is estimated at each node point—Eq. (10) is solved for SO2,out at each node point, SO2,in is estimated from flow conservation—and this process was iterated until convergence.
The solution with arteriole resistance and oxygen consumption varying in time is similar. At each time step, the arteriole diameter at specific nodes can vary, causing a variation in resistance. The pressure distribution is recalculated at each time point, considering the incremental arteriole volume increase and further considering compliant volume changes in capillaries and venules given the pressure at the previous time instance. Flow in and out of each node point is then calculated from the pressure distribution and used to calculate the time derivate in SO2 utilizing Eqs. (10–14) with the possibility of changes in CMRO2. We ran with a time step of 0.1 ms to ensure stability. Stability was tested by running with smaller time steps and obtaining the same solution. More efficient adaptive time step solvers could be implemented in the future.
Given the vessel lengths and diameters summarized in Table 1, we find blood volumes of 887, 804, and 1279 pl in the arterioles, capillaries, and venules, respectively. This is consistent with the volume fraction partitioning given in Duong and Kim (2000).
The resistance of each vascular segment derives from length, diameter, and viscosity. These values are summarized in Table 1 and result in cumulative resistances of 11,000, 39,000, and 5400 mmHg s/μl in the arterioles, capillaries, and venules, respectively, placing 20% and 70% of the resistance in the arterioles and capillaries, respectively.
Given the pressure drop from 60 mmHg to 25 mmHg from feeding arteriole to draining venule, we obtain an estimate of the pressure and red blood cell velocity in each segment (flow divided by cross-sectional area, see Table 1). These pressures and velocities in each vascular segment are consistent with those given in Lipowsky (2005), as compared in Fig. 3. The volume divided by flow gives the vascular transit time τflow. This transit time is given for each segment in Table 1, and sums to 0.090, 0.082, and 0.131 s in the arterioles, capillaries, and venules, respectively, for a total of 0.303 s. We note that the apparently short transit time arises from the fact that we are only considering the transit of blood over a distance of ~1.5 mm from a 30-μm arteriole to a 40-μm venule, and not the transit time of blood through a localized volume of tissue, which has the longer transit time of 1.5 to 3 s.
The hemoglobin oxygen saturation, SO2, in each branch is determined by the blood flow and the flux of oxygen across the vascular membrane as driven by the oxygen pressure gradient (Eq. (13)). The baseline SO2 distribution is given in Table 1 and is in qualitative agreement with Vovenko (1999) with SO2 dropping to ~75% just before the capillaries and to ~55% just after the capillaries, and slowly decreasing a few percent thereafter (see Fig. 3d).
The sensitivity of baseline velocity, pressure, and SO2 to vascular segment length and diameter is explored in Fig. 4. Increasing the lengths increases resistance and thus leads to a reduced velocity. This also results in an increased transit time, enabling a greater extraction of oxygen as evidenced by the decrease in SO2. Decreasing the diameter of the vascular segments increases resistance and decreases flow. Velocity is flow divided by the cross-sectional area. Both flow and cross-sectional area decrease with a smaller vascular diameter, with the result of an increase in velocity for our network. The increased velocity with constant vascular length indicates a smaller vascular transit time, suggesting less extraction of oxygen. There is also less oxygen in the blood because of the reduced diameter, with the result that there is an increased oxygen extraction fraction evidenced by the significant reduction in SO2. Increasing the vascular diameter significantly decreases velocity with a small increase in SO2. Taken together, we observe that our model is in agreement with literature values for velocity, pressure, and SO2, and is good for vascular lengths ranging ± 50% and vascular diameters ranging less than ± 20%.
To examine the dynamic response to arteriolar dilation, pre-capillary arteriolar resistance was locally decreased as
for two neighboring 4° arterioles and their downstream 5° and 6° branches with ΔRa =0.05. Recall that resistance is proportional to the diameter as d−4, and thus the 5% resistance decrease corresponds to a 1.3% increase in arteriole diameter. This results in a rapid local increase in flow (2.3%) and volume and a slower local increase in oxygen saturation as indicated by the red lines in Fig. 5. The surround flow is seen to decrease as flow prefers to follow the path of least resistance through the dilated arterioles. Surround volume initially decreases in the compliant capillaries and venules, but subsequently increases downstream due to backpressure as the venules downstream from the dilating arterioles increase in volume. Oxygen saturation is decreased in the surround as reduced flow results in increased oxygen extraction fraction. Note that oxygen consumption remained constant for these simulations.
To verify that our model is not limited to small changes only, we also simulated diameter changes of 5.5%, 16%, and 44%, which produced flow responses of 9.6%, 28%, and 65%, respectively. No convergence or stability issues were observed for these larger changes. The only notable difference in the relative response of the different vascular compartments was that the fractional changes in total hemoglobin are weighted higher in the arteriolar compartment as the driving arteriolar diameter changes are increased.
As observed in Fig. 5, there is a surround negativity in flow, volume, and oxygen saturation. The surround behavior of these hemo-dynamic parameters compared to the center region of interest is compared for the arterioles, capillaries, and venules in Fig. 6. Flow exhibits a surround decrease that is approximately 15% of the center increase. Blood volume does not exhibit a noteworthy increase or decrease. Upon closer inspection, we observe that surround volume decreases in the capillaries due to reduced pressure, but increases in the venules due to increased backpressure. Deoxy-hemoglobin, however, exhibits a noticeable increase in the surround in the capillaries and venules, while oxy-hemoglobin exhibits a modest decrease. The magnitude of the surround negativity relative to the center response is only weakly diminished as the absolute magnitude of the center response increases.
In Fig. 8, we plot the modulation of the hemoglobin concentration response by changes in the rate of oxygen consumption. The change in oxygen consumption was given the same dynamics as the arteriole dilation, i.e.,
Results were compared for 0%, 0.5%, 1%, and 1.5% increases in CMRO2. Note that, as the percent change in CMRO2 approaches the percent change in blood flow, the deoxy-hemoglobin response goes from negative towards positive as the increase in oxygen consumption exceeds the increase in oxygen delivery.
We created a transient arteriole dilation to explore the relative time-to-peaks and magnitudes of the different hemodynamic parameters. The form of the dilation was given by (recall that resistance is proportional to diameter−4)
where Rao is the initial arteriole resistance, ΔRa is the magnitude of the resistance change caused by the dilation, Toffset is the temporal offset to the onset of the dilation, σwid is the temporal width of the dilation, H(t) is a Heaviside step function equal to 0 when t < 0 and 1 when t > 0, and the denominator normalizes the numerator to a maximum of 1. The hemodynamic responses to this transient arteriole dilation are shown in Fig. 9 for ΔRa =0.05, Toffset =0.1 s, and σwid = 1.0 s. All other baseline parameters are given in Tables 1 and and22.
The fourth column shows the results of summing over the arteriole, capillary, and venule compartments represented in the first 3 columns. Units are in micromolars.
We created a transient increase in CMRO2 in addition to the arteriole dilation depicted above. The form of the increased oxygen consumption was given by
where CMRO2,o is the initial oxygen consumption, ΔCMRO2 is the magnitude of the resistance change caused by the dilation, Tc,offset is the temporal offset to the onset of the dilation, σc,wid is the temporal width of the dilation, H(t) is a Heaviside step function equal to 0 when t < 0 and 1 when t > 0, and the denominator normalizes the numerator to a maximum of 1. The hemodynamic responses to this transient increase in CMRO2 without and with accompanying arteriole dilation are shown in Fig. 9 compared with the case of arteriole dilation only. Parameters are ΔCMRO2 = 0.008, Tc,offset =0.05 s, σc,wid =0.2 s, and ΔRa =0.05, Toffset =0.1 s, σwid =1.0 s. All other baseline parameters are given in Tables 1 and and2.2. The magnitude of the CMRO2 change was chosen to be half of the flow response to agree with literature observations of a flow:consumption ratio of 2 (Hoge et al., 1999b, 2005; Kastrup et al., 2002; Leontiev and Buxton, 2007). Under the assumption that major metabolic costs are associated with oxidative phosphorylation and Na/K ATPase operation in the neuronal compartment, the time course of CMRO2 might follow that of the local field potential. Based on the work of Devor et al. (2005) on a single whisker, we chose a brief 200-ms width of the CMRO2 response to coincide with the brief electrophysiological response. Note that an early HbR increase arises from the early CMRO2 increase but that the peak HbR decrease is not altered because it occurs 1 s later, allowing the vascular transit time of 0.3 s to attenuate any HbR modulation.
In Fig. 10, the transient hemoglobin responses from the model are compared with experimental measurements from the rat whisker barrel cortex published previously in Devor et al. (2005). The model results were obtained with ΔCMRO2 = 0.10, Tc,offset =0.05 s, σc,wid =0.2, ΔRa =0.45, Toffset =0.1 s, σwid =1.7 s, and the remaining parameters are the same as used for Fig. 9. We observe good qualitative agreement in the center with the relative magnitudes of HbO, HbR, and HbT in agreement between the model and experiment. The relative magnitude of the early HbR increase and the later decrease is also in agreement between model and experiment. The experiment does differ in exhibiting a greater delay in time-to-peak for HbO and HbR relative to the model. This likely arises from the experiment being sensitive to more downstream veins having a longer transit time, which our model does not include. Both the model and the experiment show that the HbR time-to-peak is delayed with respect to HbO. In the surround, the time-to-peak of HbO and HbR is about the same in both the model and the experiment. However, in the surround, the relative magnitudes of HbO, HbR, and HbT are not in agreement between the model and the experiment. This suggests that the experimental observation of a surround decrease in HbO and increase in HbR does not arise solely from a passive redistribution of blood flow, but is actively controlled by vasoconstrictive mediators (Hamel, 2006).
The velocity in different branches of the vascular network depends on the topology of the network, the resistance of each vascular segment, and the pressure difference between the feeding input arteriole and the draining output venule. The resistance of each vascular segment depends on the diameter, length, and blood viscosity. The capillary diameter was set to 8 μm, consistent with a range of observations (Kleinfeld et al., 1998; Lipowsky, 2005; Villringer et al., 1994; Vovenko, 1999). The diameters of upstream arterioles and downstream venules were increased to decrease the integrated cross-sectional area at each branch level to result in a greater velocity in the larger vessels, consistent with the review of Lipowsky (2005). Furthermore, venous diameters were made larger to be more consistent with the observations of venous versus arteriole diameters (Vovenko, 1999) and to provide a greater volume fraction of blood in the veins versus the arteries, roughly 40% versus 30% with the remaining 30% in the capillaries (consistent with Duong and Kim, 2000). The capillary length was set to 250 μm, consistent with that of Vovenko (1999). Lacking literature support for the lengths of different arterioles and venules, these lengths were set to 100 μm. Blood viscosity depends on plasma viscosity and hematocrit as reviewed in Pries et al. (1992). Hematocrit varies with vessel diameter as reviewed in Lipowsky et al. (1980). These values taken together give the resistance of each vascular segment, and with the vascular topology and pressure drop allow the circuit equation calculation of flow, velocity, and transit time within each segment. While our balanced network of branching arterioles, capillaries, and converging venules is a simplistic representation of a true vascular network, the calculated velocity and blood pressure at each branch level agreed well with literature values (Fig. 3). Increasing the lengths of the arterioles and venules by 50% produces a moderate increase in resistance and only a small decrease in the velocity at each branch level, while changing the diameters at each branch level has a dramatic effect owing to the fourth power dependence on diameter.
Oxygenated hemoglobin is delivered by the feeding arteriole and flows through the vascular network. Oxygen bound to hemoglobin is in equilibrium with the plasma blood concentration and diffuses by the plasma, through the vascular wall, down an oxygen concentration gradient to the mitochondria, where it is metabolized. Increasing blood transit time or decreasing blood flow results in a greater fraction of oxygen extracted from the blood, thus decreasing the hemoglobin oxygen saturation. Increasing vascular length increases transit time and results in a modest decrease in oxygen saturation. Decreasing vessel diameters significantly reduces blood flow and oxygen saturation.
Comparing our VAN model predictions of velocity and oxygen saturation with literature values indicates strong model sensitivity to vascular diameters and weaker sensitivity to vascular lengths. Closer inspection of the data reveals differences in the venous velocity, arteriole pressure drop, and arteriole oxygen saturation. The model venous velocity is greater than literature values. This can be brought into better agreement by further increasing the diameter of the venules with little consequence on the oxygen saturation (data not shown) or considering venous topology giving rise to more ascending venules per descending arteriole. The literature arteriole pressure drops more rapidly than the model, suggesting greater resistance in the arterioles than is captured by our model. Obtaining better model agreement on the pressure will require exploration of more realistic vascular topologies and consideration of the resistance arising from deformation of red blood cells entering the capillaries (Lipowsky, 2005). Our model arteriole oxygen saturation drops more rapidly than literature values. Given our good agreement on the arteriole velocities, this suggests that our model of oxygen flux from arterioles to tissue needs modification. Modifications include properly accounting for the vascular wall thickness in Eqs. (11) and (13), and considering that tissue oxygen content may be greater near arterioles than near capillaries and venules. Our model currently uses the same oxygen extraction coefficient (Eq. (13)) for all vascular segments and assumes that tissue oxygen content is the same around every vascular segment. The former is an important issue considered by Herman et al. (2006). The latter is known to be false as Vovenko (1999) has measured increased oxygen content around arterioles than around capillaries and venules.
Our VAN model response to a transient arteriole dilation and oxygen consumption produced results consistent with literature in terms of relative times-to-peak and relative amplitudes for the hemodynamic variables. Blood volume (i.e., total hemoglobin concentration) and blood flow peak was followed first by a delayed peak in oxygenated hemoglobin and then by a delayed peak in deoxygenated hemoglobin. Blood volume and flow peak closely following the peak in arteriole dilation. Oxy- and deoxy-hemoglobin exhibit delayed peaks due to the venous washout of deoxy-hemoglobin. The peak for oxygenated hemoglobin is not as delayed as for deoxy-hemoglobin because of the differing relative contributions of arteriolar and venous contributions to each.
As shown in Fig. 5, the spatial extent of the blood volume response is better localized than for the flow and oxy- and deoxy-hemoglobin responses. This is consistent with the findings of Culver et al. (2005), Gautama et al. (2003), Leite et al. (2002), Mandeville et al. (2001), Sheth et al. (2004a), and Vanzetta et al. (2004). Oxy- and deoxy- hemoglobin are more spatially extended because of washout down the veins. Flow might be expected to have spatial localization comparable to that of blood volume, depending on the compliant blood volume response of the capillaries and veins, but our results indicate that the arteriolar blood volume change comprises the major portion of the blood volume change. This results in a strong localization of the blood volume response in the arterioles, whereas the blood flow response is uniform across arterioles, capillaries, and venules. The significant contribution of arteriolar dilation to the blood volume response and relatively insignificant venous contribution has been recently observed experimentally (Hillman et al., 2007; Kim et al., 2007; Lee et al., 2001). This is an important issue to be explored further as interpretation of the BOLD signal in functional MRI, particularly post-stimulus behavior, has relied on venous volume changes (Buxton et al., 1998; Kong et al., 2004; Mandeville et al., 1999).
Our VAN model showed that localized arterial dilation induced a chain of passive events in both a center ROI and in the surrounding vessels such that, concurrent with increases in flow and oxy-hemoglobin in the center, there was a smaller decrease in flow and oxy-hemoglobin in the surround (Figs. 5 and and6).6). Center-surround hemodynamic activity patterns have been measured previously in response to stimuli that activate localized areas of cortex (Devor et al., 2007; Devor et al., 2005; Sheth et al., 2005; Takashima et al., 2001), as have center-surround patterns of blood flow (Cox et al., 1993; Woolsey et al., 1996) and arterial dilation (Devor et al., 2007; Liang et al., 1995). Our model predicted that the observed surround decreases are at least partially due to passive vessel properties that create a local redistribution of blood. However, a comparison with experimental data suggested that the incorporated vessel properties are not sufficient to account for the magnitudes of decrease in the surround (Fig. 10). Several factors may contribute to this discrepancy. The first of these is the active vasoconstriction of arteries and arterioles through activation of the smooth muscles that line them (Devor et al., 2007; Hamel, 2006; Iadecola, 2004; Lauritzen, 2005), or the active changes in the diameter of capillaries (Peppiatt et al., 2006). Active vasoconstriction could be induced by the release of neurotransmitters or neuropeptides, and may be associated with center-surround electrophysiological patterns of activity consisting of center excitation and surround inhibition (Cauli et al., 2004; Derdikman et al., 2003; Devor et al., 2007; Faraci and Breese, 1993; Faraci and Sobey, 1998; Hamel, 2006; Takashima et al., 2001). The second possible factor is the simplicity in our VAN model architecture, which does not fully represent the three-dimensional (3D) asymmetric structure of a branched vascular network. In particular, inclusion of additional draining vessels for each feeding arteriole could increase the magnitude of the surround decreases in oxy-hemoglobin and blood flow. Third is arterial compliance, which is another passive mechanism that has been reported to be coupled to the regulation of blood flow (Behzadi and Liu, 2005). The relative contributions of each of these mechanisms to the timing and magnitude of the observed center-surround activity patterns are likely to be tied to the degree of center activation. A future investigation of the influence of these mechanisms with our VAN model will likely prove its further value as a non-invasive tool to study spatial–temporal neural–vascular coupling.
Despite being a simplified branching network of diverging arterioles and converging venules, the model produced spatio-temporal hemodynamic responses that agree well with experimental observations. Future improvements are needed to better account for the true 3D structure of the vessels from the pial arteries, through the descending arterioles, capillary networks, and ascending veins, up to the pial veins. This will provide information for the lengths of each vascular branch and branching details, including a number of branches and diameter change at each branch. Detailed statistics of these parameters are not provided in the literature, and are likely to vary between brain areas and even within an area (Harrison et al., 2002).
Calculation of the oxygen diffusion into the tissue will be improved by considering the true 3D structure of the vessels, allowing oxygen to diffuse to a given point in the tissue from multiple vascular segments. Our VAN model presently considers an independent tissue compartment for each vascular segment. As such, oxygen is delivered to each tissue compartment from only a single vascular segment. In addition, we fixed the initial tissue PO2 in these compartments with the result that the baseline oxygen delivery (and thus oxygen consumption) of each compartment is different. The 3D structure will allow us to specify baseline oxygen consumption as a function of spatial coordinate to, for instance, account for laminar differences. Measures of the vessel wall thickness will also be important to better model the oxygen flux from vessels to tissue (Eq. (11)). Further insight into the delivery of oxygen to the tissue will come from the comparison of these model predictions with direct experimental observation of vascular and tissue oxygen content as afforded by oxygen-quenched phosphorescent dyes (Dunphy et al., 2002), as utilized by Ances et al. (2001) and Vanzetta and Grinvald (1999).
Finally, Piechnik et al. (2008) have recently published a model of 19 vascular compartments from feeding large arteries to draining large veins based on the statistics of the vascular tree topology. Their model was developed to explore the relationship between cerebral blood flow and volume during CO2 manipulation and thus considers the CO2 reactivity of individual vascular compartments as derived from a review of the literature. The details of this flow–volume relationship is important for calibration of the fMRI BOLD signal to enable quantitative estimate of evoked changes in CMRO2 (Buxton et al., 2004; Davis et al., 1998; Hoge et al., 1999a). While their model is topologically similar to ours, it does not explore the spatio-temporal response of the vessels and oxygenation to the localized modulation of arteriolar diameter and oxygen consumption. Further vascular modeling improvements will follow from combining our spatio-temporal modeling of individual vascular segments with details of the vascular reactivity at different branch levels from Piechnik et al. (2008), and with explicit modeling of the neuronal modulation of arteriole compliance to induce arteriolar diameter changes (Behzadi and Liu, 2005).
We have developed a vascular anatomical network (VAN) model of blood flow and oxygen delivery based on physically measurable parameters. This model accurately predicts the blood pressure and hemoglobin oxygen saturation at different branch levels of the arterioles, capillaries, and venules. Further, localized arteriole dilation produces a dynamic evolution of blood flow, volume, and oxygenation consistent with experimental studies. The model-predicted spatial characteristics of the blood flow, volume, and oxygenation response reveal a passive surround decrease in flow and oxygenation. This model observation underestimates recent experimental observations of surround negativity suggesting the active role of vasoconstrictive mediators. Interpretation of experimental studies of the neurovascular unit and its response to stimuli will be improved when analyzed within the framework of a VAN model. In particular, we expect to derive a more comprehensive understanding of oxygen delivery to the tissue. Finally, simplistic lumped parameter vascular models, like the balloon and Windkessel models, are providing precise interpretation of the physiological factors driving the hemodynamic response to brain activation, i.e., arteriole dilation and oxygen consumption. This VAN model can be used to help determine the quantitative accuracy of the inferences from these lumped parameter models.
This work was supported by NIH R01EB002482, R01NS051188, R01NS057476, and R01EB000790.