Abstract
Brief neural stimulation results in a stereotypical pattern of vascular and metabolic response that is the basis for popular brain-imaging methods such as functional magnetic resonance imagine. However, the mechanisms of transient oxygen transport and its coupling to cerebral blood flow (CBF) and oxygen metabolism (CMRO2) are poorly understood. Recent experiments show that brief stimulation produces prompt arterial vasodilation rather than venous vasodilation. This work provides a neurovascular response model for brief stimulation based on transient arterial effects using one-dimensional convection-diffusion transport. Hemoglobin oxygen dissociation is included to enable predictions of absolute oxygen concentrations. Arterial CBF response is modeled using a lumped linear flow model, and CMRO2 response is modeled using a gamma function. Using six parameters, the model successfully fit 161/166 measured extravascular oxygen time courses obtained for brief visual stimulation in cat cerebral cortex. Results show how CBF and CMRO2 responses compete to produce the observed features of the hemodynamic response: initial dip, hyperoxic peak, undershoot, and ringing. Predicted CBF and CMRO2 response amplitudes are consistent with experimental measurements. This model provides a powerful framework to quantitatively interpret oxygen transport in the brain; in particular, its intravascular oxygen concentration predictions provide a new model for fMRI responses.
Keywords
INTRODUCTION
Neural activity causes local changes in cerebral oxygen consumption (CMRO2), blood flow (CBF), and deoxyhemoglobin concentration.1,2 This combination of neurometabolic and neurovascular coupling enables the use of hemodynamic imaging methods to infer brain activity. Of particular interest is the response to brief periods (few seconds) of neural activity, the hemodynamic response function (HRF). The stereotypical HRF permits the use of powerful linear analysis methods.
Despite the utility of the HRF for quantifying brain activity, the mechanisms of transient oxygen transport remain poorly understood. There are at least three controversial issues. First, the role of venous volume changes is not clear, particularly in response to brief neural stimulation.3–5 Second, there is disagreement about the relative role of capillaries and arterioles in oxygen transport to tissue.6–8 Third, experiments in oxygen mass balance have sometimes shown that metabolic consumption is insufficient to account for global oxygen loss from the vasculature. It was therefore suggested that tissues at or adjacent to the vascular wall may consume the missing oxygen.9,10 Thus, questions remain about oxygen transport both in steady state and as a consequence of brief stimulation.
Various computational models have attempted to explain the dynamics of the HRF. The most popular model is the ‘balloon model,’ which postulates non-linear venous inflation with increased CBF;3,11 such models have become increasingly sophisticated and complex.12,13 These models provide accurate fits for the HRF. However, recent brief stimulation experiments demonstrate significant arterial volume changes with minimal venous volume changes.4,5,14 Accordingly, a ‘proximal-integration’ hypothesis was suggested, in which a vasodilation signal proximal to neural activity rapidly propagates upstream to pial arteries. 15 These arterial transients, rather than venous compliance, must have a key role in the dynamics of the HRF. Accordingly, for brief neural activity, the current theoretical framework around non-linear venous dilation needs to be revised.
Previously, we proposed a transient oxygen-transport model driven by upstream arterial flow changes. 16 We asserted that such dilation could generate pressure fluctuations in pial arteries that couple to microvascular flow changes. We described this flow response using a lumped linear model. The earlier model was able to fit a heterogeneous set of transient tissue oxygen time-series responses measured in cat visual cortex. 17 However, this model had some limitations. Specifically, it did not include the dynamics of oxygen dissociation, used only a rectangular function for the CMRO2 response, and did not operate from a fixed steady-state condition. Moreover, without considering hemoglobin oxygen dissociation, the model was unable to predict absolute oxygen concentration.
Here, we substantially extend the previous model. First, the transport model now includes realistic treatment of oxygen dissociation, which enables calculation of oxygen concentration in absolute units and a full treatment of steady-state oxygen mass balance. Second, the model can thereby start from a realistic steady state, decoupled from the transient dynamics. Previously, the steady-state conditions were linked to the transient modeling by modifying the transport geometry. Third, we include the possibility of wall consumption (or similar loss mechanism), which can be necessary to obtain oxygen mass balance depending upon the chosen set of the steady-state parameters. Finally, instead of using a rectangular function for transient CMRO2 response, the model uses a gamma function to match experimentally observed temporal characteristics of the CMRO2 response. This feature provides physiologic insight into how the CMRO2 response affects the HRF.
All in all, we provide a tractable oxygen-transport model based on prompt arterial dilation evoked by brief neural activity capable of predicting absolute oxygen concentration. Our model provides accurate and detailed predictions of transient extravascular oxygen time courses with heterogeneous combinations of initial dips, hyperoxic peaks, undershoots, and ringing in cat visual cortex using only six optimization parameters. It also provides CBF and CMRO2 responses that explain the physiology of the HRF in detail.
MATERIALS AND METHODS
Transport Model
A one-dimensional (1D) cylindrical geometry is employed to represent a single, stereotypical element of the microvascular network in cerebral cortex, a ‘vascular unit’ for oxygen transport starting with small arterioles through its associated capillary beds (Figure 1A). Experimental findings on oxygen transport in rat brain demonstrated that oxygen transport in the brain occurs principally from very small vessels (eg., small arterioles and capillaries).8,18 Accordingly, we assumed that the oxygen transport begins principally from third-order arterioles. These studies also motivated our assumption of a uniform cylindrical geometry.

The oxygen-transport and flow models. (
Our transport model is based on the mechanisms of flow and diffusion in multiple concentric cylindrical compartments. We illustrate the relevant dynamics with a simplified example with only two compartments: intravascular and extravascular. Temporally variable CBF occurs in the intravascular compartment (subscript
Thus, the rate of change of oxygen is the sum of oxygen delivered by flow balanced by losses from transmural diffusion. These concepts form the core of our analysis.
In the current model, we use a concentric three-compartment geometry to include the dynamics of oxygen dissociation (Figure 1B). This represents the small arterioles and capillaries where oxygen transport mainly occurs. The intravascular compartment is separated into erythrocyte (subscript
where is volume per unit length,
Equation (2) states that the rate of change in erythrocyte oxygen concentration is balanced with oxygen-dissociation efflux to the plasma compartment (first term on the right side) and oxygen delivery in
Analytic solution of this set of equations is not possible in general, so we utilized numerical methods. Hemoglobin oxygen dissociation is often approximated using the Hill equation,
19
which has a sigmoidal form that complicates numerical solution of the mass-balance equations. We utilized a quadratic fit to the Hill equation to describe oxygen dissociation within the range of oxygen concentrations relevant to
Setting all the time derivatives to zero in steady state yields a pair of ordinary differential equations that can be solved analytically (Appendix 2). The results relate the blood vessel length (L) and wall consumption (W) to a set of experimental parameters obtained from the literature (Table 1). The steady-state solution provides the starting point for the evaluation of the dynamic changes evoked by neural stimulation.
Steady-state parameters
CBF, cerebral blood flow; CBV, cerebral blood volume; CMRO2, cerebral oxygen consumption.
The choice of parameters strongly affects the steady-state solution, so we must consider our choices carefully. In particular, there is uncertainty or wide variability in measured values for steady-state CMRO2 (Γ
CBF and CMRO2 Response
We model the transient response in blood flow based upon the experimental findings that suggested the proximal-integration hypothesis. 15 During transient neural activity, vasodilation starts in the microvasculature local to the neural activity (Figure 2). This dilation causes a local increase in blood flow that generates mechanical shear stress on the vascular wall. This shear stress, in turn, produces a back-propagating wave of vasodilation into the upstream arterial vasculature, which then creates an upstream pressure fluctuation. The time scale of this process is fast (~ 1 second), 14 so that it can be treated as an impulse on the time scale of the blood-flow dynamics.

Sequence of the flow response caused by proximal integration. (1) A brief period of neural activity triggers increased microvascular flow, leading to (2) a wave of vasodilation that propagates up to (3) the pial arteries. This transient dilation, in turn, drives (4) a pressure fluctuation that couples to (5) a downstream viscous-flow response.
We use a lumped linear flow model to describe the effects of this pressure fluctuation in the pial arteries (Figure 1C). It has been shown that an inertial element is necessary when considering hemodynamics at the scale of the pial arteries. 27 We, therefore, model the upstream vasculature by including an inertive element to represent the kinetics of blood flow (red outline). The effects of the downstream venous vasculature are modeled by including resistive and compliant elements to represent its aggregate dynamics (blue outline). This lumped linear flow model is also known as a 4-element windkessel (FEW), which has been widely used for cardiovascular flow modeling. 27 The flow response in this linear system is described by an ordinary differential equation:
where
Solution of Equation (5) thus gives the microvascular flow response, which has two forms depending upon the time constants of the system:
Three parameters characterize the flow impulse response: frequency
Spatially, uniform oxygen uptake from the surrounding tissue is assumed along the whole length of the modeled cylindrical blood vessel. In several experiments, the CMRO2 response to brief neural activity consists of a prompt fast increase and relatively slow decay.30–32 A gamma function is therefore used to model the transient CMRO2 response:
where
We assume that these two transient impulse responses in CBF and CMRO2 are driven by the visual stimulus, which we model as a rectangular pulse, rect(
Fitting Tissue-Oxygen Measurements
We applied the model to experimental tissue-oxygen measurements obtained in area 17 of cat visual cortex. 17 Briefly, in those experiments, a combined microelectrode and polarographic sensor was used to measure the dynamics of neural activity and tissue oxygenation, respectively. Spiking activity recorded with the electrode from a particular local cell was used to guide the visual stimulus selection. Stimuli were drifting full-contrast gratings applied for 4-second periods to the receptive field of each local cell. There were eight stimulus conditions: seven grating orientations that formed a rough orientation-tuning curve, including the ‘preferred’ (largest evoked spiking) and orthogonal orientations presented to the dominant eye; only the preferred orientation was presented to the non-dominant eye. The field of sensitivity of the oxygen sensor was a sphere ~60 μm in diameter. Both neuronal and oxygen measurements were collected for 40-second periods starting 10-seconds before stimulus onset. The measurements gave only relative rather than absolute oxygen concentration dynamics.
We obtain the full transient evolution of the system starting from a particular steady state. Heun’s modified Euler’s method
33
is used to accurately evolve numerical solutions of the transient concentrations in each compartment. The system is spatially gridded with at least 80 steps between 0 and
Optimization and Data Analysis
We use the model to fit the 166 individual oxygen-concentration time series
17
described above. Model fits start from a particular choice of steady-state, and best-fit transient response predictions are then computed using six optimization parameters. Three parameters specify a CBF response: frequency
For successful fits, we also investigate correlations between model optimization parameters, and the observed neuronal activity. Neuronal activity is calculated as the difference between the mean of the steady-state neuronal firing rate and the mean firing rating during stimulus presentation. Because each electrode sampled individual units, there was substantial cell-to-cell variability in the neuronal firing rates (a population average across columns would show much less variability). To deal with this variability, we normalize the measured neuronal firing rate and each parameter by a ‘cocktail mean’ for each cell, that is, the mean value of each parameter across all stimulus orientations.
RESULTS
Steady-State Conditions
We use our analytic result for the steady-state intra- and extravascular spatial oxygen concentrations (Appendix 2). Both concentrations decrease from proximal to distal, typically starting with a non-linear exponential decrease at high concentration, followed by a comparatively linear decrease at lower concentrations (Figure 3). The steady-state plasma and extravascular concentration in the spatial sampling region (blue shading) are 0.053 and 0.044 mmol/L, respectively (38.13 and 31.65 mm Hg).

Normalized spatial oxygen concentrations in plasma and extravascular compartments for our nominal steady state. Proximal plasma concentration, 0.093 μmol/L (66.9 μm Hg) at z = 0 was calculated from hemoglobin dissociation assuming an erythrocyte concentration of 8.56 μmol/L. From proximal to distal, there is a non-linear decrease at high concentrations, followed by a relatively linear decrease at lower concentrations. Blue shaded region represents the spatial sampling used to compare modeling results with the measured transient extravascular oxygen concentrations.
We investigate the effect of the initial transmural gradient (σ), steady-state CMRO2 (Γ

Effects of three parameters on wall consumption (W) and cylinder length (L) in steady state. Panels
Figures 4C and 4D show the effects of varying Γ
Transient Oxygen Response Fitting
The temporal evolution of the oxygen responses are simulated corresponding to various choices of CMRO2 and CBF responses. The cylinder length,
Most fits exhibit a stereotypical response that consists of three temporal regimes: initial oxygen decrease (initial dip), hyperoxic peak, and then undershoot and late-time oscillation (ringing). However, large variations of the relative magnitudes and timings of these phases are observed (Figure 5). The model is capable of fitting this wide range of variations in the measured oxygen response. In all cases, the transient CMRO2 response, Γ

Representative examples of model fits to measured extravascular oxygen time courses (
Some measurements exhibit a brief initial dip and strong hyperoxic peak (Figures 5A and 5B). In such cases, the effects of the large flow responses (Figures 5E and 5F, red lines) quickly rise to overwhelm the relatively weak extravascular oxygen demand (green lines). In contrast, there are cases of strong initial drop followed by slow recovery (Figures 5C and 5D). Here, the model fit consists of relatively large and slow oxygen responses (green lines, Figures 5G and 5H), while flow responses are comparatively weak. The model also fits either strong (Figures 5A and 5C) or weak ringing (Figures 5B and 5D) by adjusting the damping time of the flow response.
The model successfully fit ~97% of the measured oxygen time courses. Variations in the parameters follow the heterogeneity of the hemodynamic response behavior. Flow response parameters have the following mean and standard deviation values: amplitude, 25.4%±11.9%; damping time, 11.1 ± 10.9 seconds; frequency, 0.064 ±0.016 Hz. CMRO2 response parameters are: amplitude, 19.7%±9.0%; damping time 1.6 ± 1.2 seconds; and shape coefficient, 2.0 ± 1.0.
The flow-response amplitudes show a complex distribution with peaks at 13.6%, 23%, and 41.8%; most flow responses are < 45% (Figure 6A). The flow damping times show a Rayleigh-like distribution with a peak at 4.6 seconds (Figure 6B). Flow frequencies have a bimodal distribution with modes at 0.048 and 0.068 Hz (Figure 6C). The CMRO2 response amplitudes have a complex distribution with peaks at 8.6%, 18.7%, and 25.4%; most CMRO2 responses are <30% (Figure 6D). The CMRO2 damping-times include two populations, with a cluster of small values (peak at B0.63 seconds) corresponding to sharp decay of metabolic responses, while a distribution of larger values (peak at ~2.5 seconds) corresponds to slower decreases in CMRO2 (Figure 6E). The CMRO2 shape coefficients also exhibit two populations: small values (peak at ~1) corresponding to prompt CMRO2 increase, while larger values (>2) correspond to slower onsets to the CMRO2 response (Figure 6F).

Histograms for six optimization parameters. Three cerebral blood flow parameters are shown in upper panels. (
We repeated the fitting procedure using six different sets of steady-state parameters, adjusting the initial transmural gradient (à), steady-state CMRO2 (Γ
Means and standard deviations of optimization parameters for fits to the same experimental data using 7 different sets of steady-state parameters
For successful fits, we examine correlations between the observed neuronal activity and the optimization parameters. We find strong and significant correlation (

Correlations between optimization parameters and neuronal activity. (
DISCUSSION
We provide a novel model based upon 1D convection-diffusion oxygen transport in a cylindrical three-compartment geometry. It accurately accounts for hemoglobin dissociation (neglected in previous efforts3,11,12,16), which allows predictions of absolute oxygen concentrations both in steady state and during the response to brief neural stimulation. Neural activity is assumed to produce responses in CBF and CMRO2 based on linear models, and these responses compete with each other to predict measured extravascular oxygen concentrations. Thus, the model provides a detailed physiologic interpretation of neurovascular and neurometabolic coupling in the brain. Accurate fits to experimental data are possible starting from a wide range of plausible steady states, and the model permits testing of a wall-consumption hypothesis.
The proximal-integration hypothesis motivates use of a lumped linear flow model, a FEW. This flow impulse response has the form of a damped harmonic oscillator. Previous models used oscillatory descriptions to model the flow response.28,29 Here, we provide a physiologic basis for this oscillatory character in terms of energy exchange between the kinetic energy of moving arterial blood, and energy stored in the compliance of the downstream venous vasculature.
The oscillatory character of this flow model permits fits to the undershoot and ringing observed in our measured extravascular oxygen time courses. The flow frequency has a bimodal distribution at ~0.05 and ~0.07 Hz. The modeled frequencies are ~0.05 Hz lower than the spectral centroids of the raw experimental data suggesting non-linear coupling between the CBF and CMRO2 responses. The observed ringing is also consistent with underdamped flow responses measured using optical imaging spectroscopy, 36 two-photon imaging, 4 and laser speckle. 37 Moreover, we find a mild correlation between neural activity and flow-damping time (Figure 7C). Greater neuronal activity generates greater oxygen demand in the neuropil, and the observed increase in ringing may reflect a decrease in damping time that speeds the delivery of oxygenated blood to the activated tissue.
However, most experiments in rodents do not show strong ringing. This could be a consequence of anesthesia, which strongly dampens the HRF.36,38 Also, we show here that stimulus-evoked ringing has substantial site-to-site variations in frequency, so the ringing would be masked to spatially averaged measurements such as laser speckle or Doppler flowmetry.31,39 Regardless, our model is also capable of fitting non-oscillatory flow responses using the overdamped form of the FEW impulse response.
Low-frequency fluctuations have been observed extensively in the cerebrovasculature and are sometimes termed Mayer waves, which exhibit several frequency bands.40,41 Our predicted flow fluctuations are in good correspondence to the LF (low frequency) and VLF (very low frequency) bands. Mayer waves in these bands are usually associated with vasomotion, rhythmic fluctuations in vascular tone caused by intravascular calcium variations. 42 Our results suggest that such slow hemodynamic oscillations can be driven indirectly by neural activity in the parenchyma by way of the proximal-integration response. In this view, neural activity would generate local vasodilation, leading to mechanical fluid shear that drives waves of vascular tone changes, a phenomenon similar to vasomotion in the brain.
The stimulus-evoked CBF response amplitudes are consistent with experimental measurements. Mean CBF response amplitude is 25.4%, which is similar to experimental measurements (~30%) by PET and functional magnetic resonance imaging (fMRI). 43 However, there is no positive correlation between neuronal activity and CBF response amplitude (Figure 7B). This suggests that the flow response is not tightly localized, but occurs on a coarser spatial scale than the localized neuronal activity, which in these experiments varies on the 0.7-mm scale of the feline ocular dominance columns. A less localized flow response is consistent with the proximal-integration hypothesis, which suggests that local neural activity gives rise to an upstream vascular response that spreads in spatial extent.
The predicted activity-induced increase in CMRO2 has a mean value, ~19%, which is similar to values (~10%) measured using PET 44 and fMRI. 45 We find a strong correlation between neuronal activity and CMRO2-response amplitude (Figure 7A), consistent with the assertion that neuronal activity gives rise to localized oxygen demand in the tissue. 17
We simulate CMRO2 responses using a gamma function. Measurements of oxygen consumption by brief neuronal activation have demonstrated relatively fast increases and slow decays, 32 which we can model using appropriate choices for shape coefficient and time constant. The prompt increase in CMRO2 directly causes the initial dip in extravascular oxygen measurements. In all cases, the CMRO2 response begins first, followed by the CBF response, consistent with immediate consumption of oxygen by the neuropil due to the transient neural activity. This initially leads to a drop of oxygen in the tissue that is only slowly remediated by the sluggish CBF response. Later, the continued competition between the CMRO2 response and the oscillatory CBF response determines how the extravascular oxygen concentration returns to steady state, thus describing the observed heterogeneous late-time responses.
The model provides an analytic solution for the cerebrovascular steady state, which allows interesting predictions. We are particularly interested in the effects of oxygen consumption by the blood-vessel wall (or apposed tissues) on steady-state mass balance. We chose a small transmural gradient consistent with previous polarographic measurement in cerebral cortex,8,18 which corresponds to a relatively small wall consumption (19% of total oxygen consumption). However, wall consumption can contribute up to ~90% of oxygen consumption, depending upon the choice of transmural gradient and steady-state CMRO2 (Figures 4A and 4D). Additional experiments that measure absolute oxygen concentrations will be needed to determine how much wall consumption is actually present.
The steady-state solution also determines the length of the simulated blood vessel, which provides a rough estimate of where oxygen transport occurs in the parenchyma. Our steady-state model predicts a range of cylinder lengths, 26 to 1,017 μm, depending upon parameter choices. Shorter lengths are consistent with oxygen delivery from capillaries (<200 μm), while longer lengths suggest delivery from arterioles and capillaries (~900 μm). The model is thus consistent with previous experimental oxygen-transport findings.8,10,19
We investigated the effect of transmural gradient, steady-state CBF, and CMRO2 on the optimization parameters. The changes modify the balance between CBF and CMRO2 responses, but both remain in a physiologically reasonable range, and excellent fits to our experimental measurements are still obtained. Further experiments, therefore, will be necessary to elucidate the appropriate choices for the steady state.
Our model can be broadly used to interpret and predict oxygen transport in the brain. For example, the model can be applied to localized measurements of CBF as well as intra- and extravascular oxygen concentrations that are now becoming available using two-photon imaging methods.46–48 From these measurements, the model can provide estimates of temporal CMRO2 responses. Moreover, the model can also predict blood oxygen-level dependent (BOLD) fMRI responses. Recent experiments suggest that there is no significant venous blood volume change for vascular and metabolic responses evoked by brief neural stimulation.4,5,32 Therefore, the BOLD response is primarily a function of the intravascular oxygen concentration and can be estimated by our model (Appendix 3). Further brief-stimulation experiments with high-resolution fMRI that are fully localized to the gray matter will be needed to test this new BOLD response model. For such measurements, our model will be able to separate the effects of CBF and CMRO2 on the BOLD HRF.
The model has some limitations. Although we use weighted parameters to consider the volume effect of different vessel sizes, we underestimate the surface-area-to-volume ratio by representing the vastly branching capillary bed by a single cylinder. We can remediate this flaw by using a parallel network that breaks the effective blood volume into many identical branches. Including a longitudinal series of compartments to represent multiple orders of arterioles and the capillary bed would also increase the accuracy of the model.
This model for the HRF does not attempt to explain the mechanisms of longer neural stimulation that do lead to microvascular venous volume changes. Previous models that postulate venous inflation are motivated by the slow hemodynamic responses of oxygen autoregulation.3,11,12,49,50 Such mechanisms are likely essential to describe the hemodynamic response to extended periods of neural activity (tens of seconds). We believe that both the transient mechanisms described by our model, and the slower mechanisms that yield venous dilation are likely to operate in parallel, and both hemodynamic response mechanisms must be included depending upon time scale.
To interpret the experimental data from a polarographic oxygen sensor, we had to make a spatial sampling assumption. Because oxygen transport mainly occurs in capillaries (~52% of total oxygen transport),8,18 we chose to use the predicted extravascular oxygen concentration over a 60-mm region at the distal end of the model blood vessel, which corresponds to the mean length of a capillary segment. 34 However, we also tested two other sampling approaches: permitting the sampling to vary along the length of the cylinder as an optimization parameter, and averaging over the entire blood-vessel length. These approaches also provide good fits to the measurements that were qualitatively identical to those presented in the Results.
In conclusion, our computational model uses a simple geometry and tractable computations to provide accurate and detailed predictions of transient oxygen responses with heterogeneous combinations of initial dips, hyperoxic peaks, undershoots, and ringing in cat visual cortex using only six optimization parameters. The FEW flow model provides a mechanism by which neural activity couples to parenchymal CBF in a substantially linear fashion. The model is flexible, as overdamped flow responses could also fit experimental data that do not show ringing. The realistic steady-state formulation is adjustable, permitting tests of hypotheses for oxygen transport mechanisms. Altogether, the model provides a theoretical framework to explain the full panoply of vascular and metabolic responses evoked by brief neural stimulation in the gray matter of the brain.
DISCLOSURE/CONFLICT OF INTEREST
The authors declare no conflict of interest.
Footnotes
ACKNOWLEDGEMENTS
The authors thank Ralph Freeman for providing the experimental data, and Rick Buxton for providing thoughtful comments and advice.
APPENDIX 1
APPENDIX 2
APPENDIX 3
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
