A Model for Transient Oxygen Delivery in Cerebral Cortex

Popular hemodynamic brain imaging methods, such as blood oxygen-level dependent functional magnetic resonance imaging (BOLD fMRI), would benefit from a detailed understanding of the mechanisms by which oxygen is delivered to the cortex in response to brief periods of neural activity. Tissue oxygen responses in visual cortex following brief visual stimulation exhibit rich dynamics, including an early decrease in oxygen concentration, a subsequent large increase in concentration, and substantial late-time oscillations (“ringing”). We introduce a model that explains the full time-course of these observations made by Thompson et al. (2003). The model treats oxygen transport with a set of differential equations that include a combination of flow and diffusion in a three-compartment (intravascular, extravascular, and intracellular) system. Blood flow in this system is modeled using the impulse response of a lumped linear system that includes an inertive element; this provides a simple biophysical mechanism for the ringing. The model system is solved numerically to produce excellent fits to measurements of tissue oxygen. The results give insight into the dynamics of cerebral oxygen transfer, and can serve as the starting point to understand BOLD fMRI measurements.


INTRODUCTION
Detailed understanding of cerebrovascular oxygen transport mechanisms is particularly relevant to popular brain imaging methods that rely on hemodynamics to infer brain activity, such as functional magnetic resonance imaging (fMRI). Generally, such methods use the positive portion of the blood oxygen-level dependent (BOLD) response, which reaches its peak several seconds after a brief period of neural activation. However, there is also interest in the earlier phase of the hemodynamic response. In this period, during the fi rst few seconds after brief activation, there is evidence for a transient drop in intravascular oxygen concentration (Malonek et al., 1997;Vanzetta and Grinvald, 1999). Moreover, measurements of tissue oxygen concentration clearly show the early deoxygenation, subsequent overshoot, and oscillatory late-time behavior (Thompson et al., 2003(Thompson et al., , 2005. Here, we provide a model that explains the full time course of these transient changes in tissue oxygen concentration.
These detailed mechanisms of microvascular oxygen transport underlie the BOLD response observed using T 2 *-weighted magnetic resonance imaging. The time course of the detected contrast changes are, in fact, critically related to changes in intravascular oxygen concentration. An understanding of oxygen transport is a critical requirement, although a complete treatment of the BOLD response naturally involves other factors, such as the signal component associated with extravascular volume changes.
Most attempts to model cerebrovascular oxygen transport have taken rate-equation approaches, e.g., (Buxton and Frank, 1997;Buxton et al., 1998;Hayashi et al., 2003;Herman et al., 2006;Hudetz, 1999;Hyder et al., 1998;Mintun et al., 2001;Zheng et al., 2005). These models generally combine a non-linear fl ow model The model's most novel feature is the inclusion of the inertial effects of moving arterial blood. As the vascular system transitions from small arteries to penetrating arterioles at the pial surface, blood fl ow makes a transition from inertial to viscous fl ow as vessel diameters and fl ow velocities decrease. The Reynolds number, R = UD/n, describes this transition, where U is the fl ow velocity, D is the vessel diameter, and n is the kinematic viscosity. R is the ratio of dynamic pressure, a measure of inertial energy content, to shearing stress, a measure of frictional energy loss. Generally, larger cerebral arteries have R > 1, while arterioles have R < 1. In fact, small arteries at the pial surface have R ∼1. For example, if we consider a vessel with D = 0.2 mm transporting blood with U = 20 mm/s (Kobari et al., 1987), and take n = 4 centistokes (mm 2 /s), we obtain R = 1. If the transients of the hemodynamic response extend into these relatively upstream regions with similar fl ow velocities and diameters, they will be accompanied by inertial effects that have been neglected in previous models.
A critical feature of the new work is a simple linear fl ow model that includes inertial effects, the four-element windkessel (FEW, Figure 1). The FEW model includes an inertive element, which represents the dynamics of fl owing blood, together with the usual compliant and resistive elements, which represent the elasticity of the venous system and frictional losses, respectively. In response to brief neural stimulation in cerebral cortex, we postulate that a correspondingly brief vasodilation in the larger pial vessels gives rise to pressure fl uctuations that couple directly to the smaller arterioles, producing corresponding fl ow velocity perturbations. The FEW model allows blood fl ow to increase and decrease multiple times in response to a single brief perturbation, an effect analogous to the "sloshing" observed when an elastic container of fl uid is sharply disturbed.
We note that the FEW model is a substantial simplifi cation of the complex fl ow dynamics associated with the whole of the cerebrovasculature. However, within the context of linear network models, this kind of analogy has been utilized extensively to solve many problems. Complex networks of linear mechanical or fl uid elements can always be converted to a Norton or Thevenin equivalent with a multi-pole admittance or impedance (Olson, 1943), and this approach has been utilized extensively in cardiovascular applications (Burattini and Di Salvia, 2007;Stergiopulos et al., 1999;Thiry and Roberge, 1976). In this case, we have utilized a second-degree, two-pole model for the impedance transfer function.
As a fi rst step in examining the effects of such inertial dynamics on cerebral oxygen transport, we combine the FEW fl ow model with a simple one-dimensional, three-compartment, convectiondiffusion model for oxygen transport in order to explain measurements of tissue oxygen dynamics during visual cortex stimulation experiments. This combination of the FEW fl ow perturbations with convection-diffusion dynamics produces a simple but powerful model capable of explaining a wide range of experimental tissue oxygen time-course measurements in mechanistic terms.

MODEL
We utilize a simple model geometry, a pair of uniform cylinders of length L that divides space into intravascular, extravascular, and intracellular compartments (Figure 2). Although this geometry cannot fully approximate the complex cerebrovascular architecture, it serves as an elementary starting point, roughly representing the aggregate dynamics of the vascular tree from arterioles through the capillaries.
The walls of the inner cylinder permit diffusive transfer of oxygen from the intravascular to the extravascular compartment. In steady-state equilibrium, we assume constant fl ow of blood inside the cylinder (U 0 ), creating a linear oxygen concentration gradient along its length; the oxygen concentration is higher in the intravascular compartment so that oxygen diffusion matches the uptake of oxygen by the extravascular tissue (Γ 0 , steady-state CMRO 2 ); the intracellular compartment serves purely as an oxygen sink. Note that this one-dimensional model has no explicit radial variation. Both inner compartments have radially uniform oxygen concentrations with a FIGURE 1 | The FEW fl ow model. (A) Lumped-element equivalent circuit schematic. Coupling to the "test" arteriole at T is represented by the load R A >> R 2 , R 1 ; the steady state pressure at T gives rise to fl ow U 0 . At t = 0, a brief upstream vasodilation is represented by the closing of the switch around R 1b . The inertive element maintains fl ow continuity by introducing a pressure transient shown as an impulsive source in (B), which then couples to the load to create fl ow perturbation U 1 . Transient oxygen-delivery model jump discontinuity between them that drives the diffusive oxygen transfer. In the interest of simplicity we will neglect the non-linear relationship between oxygen concentration in the blood plasma and the erythrocytes. This assumption may distort details of the intravascular modeling, but should permit investigation of the convectiondiffusion dynamics on the extravascular oxygen concentration. At the proximal end of this cylinder we postulate a connection to the larger vascular network in which vasodilation (or other regulatory transient) produces a pressure fl uctuation that we can model using the FEW. This linear model allows exchange between the kinetic energy of the blood fl ow, and the stored energy corresponding to small elastic distortions of vessel walls within the downstream venous compartment, damped oscillatory behavior often called "ringing." These upstream pressure fl uctuations couple directly to fl ow perturbations, U 1 (t), within the model intravascular compartment.
In addition to the fl ow perturbation, we also assume that neural stimulation produces a prompt rectangle-function increase in cerebral oxygen consumption (CMRO 2 ) throughout the extravascular compartment (Γ 1 ). So long as the stimulation duration is short compared to the time scale of the hemodynamic response, this step-function serves as a good approximation, even though the precise time course of the oxygen consumption may be more complex. In modeling results for a few cells, we confi rmed that there was no signifi cant difference between the rectangle-function assumption and a time-course that more precisely follows the local neuronal response.
Finally, we assume that a sensor at the distal end of the cylinder measures the effects of these perturbations in the extravascular compartment. Because this sensor measures the oxygen concentration produced by many local capillaries within the gray matter parenchyma, we average the distal portion of model oxygen concentration spatial profi le over a length of a capillary segment, ∼60 µm (Cassot et al., 2006). Altogether, the model utilizes a total of ten parameters, fi ve of which are generally fi xed when the model is applied to fi t experimental measurements. Thus, we present a fi veparameter model capable of fi tting a large experimental dataset.
In the interests of simplicity, we assume that oxygen concentrations are radially uniform and variation only occurs along the longitudinal (z) axis of the cylinder. Thus, the fl uid has oxygen concentration Q i (z, t) in the intravascular compartment (r < r 1 ), and Q e (z, t) in the extravascular compartment (r 2 > r > r 1 ); the jump discontinuity between these compartments drives the diffusive coupling. We will also neglect longitudinal diffusion in both compartments, a reasonable assumption because longitudinal oxygen concentration gradients will be small compared to the near-discontinuity across the vessel wall. The intravascular fl uid has a spatially uniform but time-varying fl ow velocity, U(t).
Consider a small (differential) length element of the cylinder. In each length element, the change in intravascular oxygen concentration is a balance between the downstream fl ow of oxygenated blood, and the radial diffusion of oxygen across the vessel wall. The diffusion rate is specifi ed by coeffi cient A, which is the diffusion conductance weighted by the surface-to-volume ratio of the intravascular compartment (see Appendix for details). A coupled pair of convection-diffusion equations then specifi es the relationship between Q i and Q e in this simple system: where Γ(t) represents oxygen uptake from the tissue (CMRO 2 ). We will use steady-state equilibrium conditions to specify A, so we need only specify the differential volume ratios between the compartments, α = − r r r We use the steady-state solution to obtain relationships between certain parameters of the model. The model is normalized to the intravascular oxygen concentration at the control point, so that Q i0 (0) = 1, and oxygen concentration becomes a dimensionless quantity; this choice gives Γ and A units of inverse time. The oxygenextraction fraction, E, is specifi ed as the relative oxygen concentration at the downstream end of the intravascular compartment, z = L, so that Q i0 (L) = 1 − E, and Γ 0 = αEU 0 /βL. We can then specify the radial oxygen diffusion by choosing the initial relative tissue-oxygen concentration, Q e0 (0) = Q e00 , so that A = EU 0 /L(1−Q e00 ). Thus, specifying α, β, E, U 0 , L, and Q e00 sets the initial state of the system.
We simulate the experimental response in this system by assuming an initial constant-fl ow equilibrium, and introduce perturbations in fl ow, U(t) = U 0 + U 1 (t), and oxygen uptake rate, Γ(t) = Γ 0 + Γ 1 (t). Equation 1 could then be linearized to fi rst order in these perturbations, but we will use the full dynamics to properly model its non-linearities.
We model the fl ow perturbation using a four-element windkessel (FEW), a lumped linear equivalent-circuit model for the blood fl ow dynamics (Figures 1A,B). The inertive element, I, models Newton's law regarding the fl ow of blood in the pial arteries. The series resistance, R 1 = R 1a + R 1b , represents frictional losses associated with the parallel combination of all downstream fl ow through the fractal vascular tree. The combination of the compliant element, C, and its parallel resistance R 2 , represent the lossy compliance of the downstream venous vasculature. In steady state, and neglecting cardiac pulsatility, there is a constant pressure at point T, which represents the input to our model three-compartment geometry ( Figure 1A). The model cylinder couples resistively to the larger network, but the load resistance R A i s much greater than R 1 and R 2 , so it does not perturb the global dynamics. The transient behavior of this system is modeled by the switch across R 1b , which briefl y closes to simulate a vasodilation that decreases the upstream vascular resistance. Because there was steady-state fl ow through I, the vasodilation creates a brief positive pressure fl uctuation at T, which corresponds to the impulse response of the network in Figure 1B.
Because this system is linear, we need only derive its impulse response function to fully characterize its behavior. The fl ow perturbation impulse response has three forms depending upon the time constants of the system: where the fi rst form is under-damped, the second critically damped, and the third over-damped ( Figure 1C). Thus, the fl ow perturbation is specifi ed by the three parameters: f, τ, and amplitude U 10 . All of the experimental fi ts utilized fl ow perturbations that varied from under-damped to critically-damped; overdamped response fi ts were not observed.

SIMULATIONS OF EXPERIMENTAL MEASUREMENTS
We applied the model to a large dataset of oxygen measurements obtained in cat brain area 17 (Thompson et al., 2003). In those experiments, visual stimuli were full-fi eld drifting gratings presented monocularly for 4-s periods. During a series of penetrations in three cat brains, combined neuronal and oxygen measurements were collected for a 40-s period starting 10 s before the stimulus onset. Recordings were collected for eight stimulus conditions: seven grating orientations, including the "preferred" (largest evoked spiking orientation, presented to the dominant eye), and only the preferred orientation presented to the non-dominant eye. We modeled this experiment with two perturbations. First, to model a step-like increase in CMRO 2 during the stimulus period, we set Γ 1 (t) = Γ 10 H(t/T − 0.5), where H(t) is the rectangle function and T = 4 s is the stimulus duration. Second, to model the fl ow perturbation produced by a step-like dilation of a sphincter control point, we set U 1 (t) = U 1,FEW *H(t/T − 0.5), where the "*" symbol denotes convolution. During the fi tting procedure, the fl ow was limited to be ≥0 to avoid fl ow reversals. Although fl ow reversals have been observed in individual capillaries, they were not a dominant feature, nor were they observed in larger vessels (Kleinfeld et al., 1998). Because our model lumps together the aggregate effects of a portion of the arteriolar tree and capillaries, we consider it reasonable to neglect the effects of fl ow reversals on the oxygen transport.
We obtained the temporal evolution of the system of Eq. 1 from its initial linear equilibrium by numerical integration using a standard leapfrog technique. The absolute O 2 concentration was limited to be ≥0 in both compartments. The two-compartment system was gridded longitudinally into at least 40 steps between 0 and L; more steps were used if needed to enforce a minimum step size of 20 µm. Stability was ensured by making the temporal increment small compared to the steady-state convection, Δt = 0.2 Δz/U 0 . Stability was assured by running lengthy simulations (e.g., 1000 s) and noting that the equilibrium state persisted with negligible numerical noise when unperturbed, and returned to the equilibrium state after perturbation effects subsided. To improve speed, the numerical integration was coded in C, and then linked to Matlab (Mathworks Corp., Natick, MA, USA). To compare the model to the experimental measurements, the temporal evolution of the relative perturbation of Q e was followed from t = 0-25 s, and the results were sampled over a 60-µm region about z = L.
The full model presented above involves a total of 10 parameters: E, U 0 , L, Q e00 , Γ 10 , U 10 , f, τ, α, and β. We simplifi ed matters by fi xing the fi rst three parameters based on experimental observations in the literature. We set the oxygen-extraction fraction E to 0.44 (Ito et al., 2004;Tsai et al., 2003;Vovenko, 1999), the steady-state fl ow velocity U 0 to 0.77 mm/s (Kleinfeld et al., 1998), and Q e00 to 0.66 (Tsai et al., 1998(Tsai et al., , 2003. The choice of an appropriate value for α requires consideration of gray matter cytoarchitecture. The blood volume fraction in the cortex (CBV) is close to 0.06 (Ladurner et al., 1976;Phelps et al., 1979). If the oxygen probe uniformly samples the whole of the extravascular and intracellular volumes, we would choose α = CBV = 0.06. However, the probe is mechanically restricted to the extracellular volume. The extracellular volume fraction in gray matter has been measured by tracer diffusion methods to be 0.2 (Nicholson and Sykova, 1998;Rusakov and Kullmann, 1998). So, if the oxygen probe only samples the portion of the volume that is both extracellular and extravascular, we would choose α = 0.06/0.2 = 0.3. We cannot resolve this ambiguity, so we chose the mean between these two extremes, α = 0.18. This then requires β = α/CBV = 3.
The remaining fi ve parameters were then initially adjusted manually to fi nd an approximate match to the extravascular oxygen measurements associated with maximal neuronal activity in each cell. Adjustment began with a standard set of parameters values, similar to the median parameter values reported in the Results. These parameters were adjusted so that the model time-course roughly captured the major features of the measurement. For example, the CMRO 2 perturbation Γ 1 was adjusted to roughly match the magnitude of the initial dip in tissue oxygen concentration, the fl ow perturbation amplitude U 1 was adjusted to match the hyperoxic peak, the frequency f adjusted to match the observed oscillatory period, and so forth. Using this parameter set as a starting point, we then used non-linear optimization (Lagarias et al., 1998) to improve the fi ts to oxygen measurement time series across all stimulus conditions of that cell.
We repeated this adjustment and optimization process using three individuals, and found no signifi cant differences in the results. Specifi cally, median parameter values varied by <10% of the standard deviation of each parameter value across the ensemble of measurements. Thus, individual biases in the initial fi tting strategy had a negligible effect on the results.
We arrived at the best-fi t parameters in the maximum likelihood sense by minimizing the negative sum of the log likelihood of the model at all time points, given the data. The likelihood was computed assuming a Gaussian error distribution defi ned by the mean and standard error of the mean based on the data provided by Thompson and colleagues. The quality of a fi t was assessed by computing the mean deviation between the fi t of the measurement as a fraction of the experimental standard error. Fits where the mean deviation was <1 were considered successful. We also performed the optimization by minimizing a simple unweighted root-mean-square (RMS) error metric and arrived at similar results. Quality-of-fi t was then assessed by comparing the ratio of the RMS fi t error to the RMS standard error; once again, fi ts where this ratio was <1 were considered successful.
The tissue-oxygen measurements included a 10-s period before stimulation. For most of the modeling, we subtracted the mean value of this period to establish a baseline for the subsequent measurements. However, for a small fraction (∼20%) of the measurements, there appeared to be substantial linear baseline trends and these were removed before fi tting. Trends were estimated by fi tting a straight line to the data points near the beginning (−5 to 0 s) and the end (25-30 s) of the measurement period.

CHARACTERISTICS OF THE MODEL
The model produced plausible results for the extravascular oxygen concentration corresponding to the two perturbations (Figure 3, solid lines). The temporal responses for the extravascular oxygen concentration perturbation are shown in Figure 3A for the median parameter set corresponding to the experimental fi ts to the entire dataset of 21 cells described below (L = 0.9 mm, Γ 10 /Γ 0 = 0.21%, U 10 /U 0 = 29%, f = 0.095 Hz, and τ = 8.3 s). The fi ts show an initial oxygen decrease caused by the 4-s-long pulse of CMRO 2 increase, followed by a positive phase of oxygen increase produced by the fl ow response together with subsequent under-damped oscillations. The intravascular fl ow response, shown in red, is fairly prompt, corresponding to the damped sinusoidal impulse response convolved with the 4-s-duration rect-function stimulus duration. We also illustrate the corresponding intravascular oxygen concentration  (Figure 3, dashed lines). The oxygen concentration results are distinctly different in the intravascular compartment, in that the deoxygenation phase is smaller or absent.
We investigated the effects of the fi ve parameters by plotting the model output as each parameter was increased and decreased from its median value by its standard deviation across the experimental data (Figures 3B-F). The length parameter (L, Figure 3B) has a global impact, with greater lengths increasing the effects of the extravascular CMRO 2 perturbation as well increasing transit-delay and smoothing. The magnitude of the CMRO 2 perturbation (Γ 1 , Figure 3C) has a much more temporally specifi c effect, directly modulating the magnitude of the early deoxygenation portion of the response. Note how the deoxygenation effects oppose the fl ow effects, so that an increased CMRO 2 perturbation decreases the magnitude of the initial fl ow-generated hyperoxic peak and viceversa. The long delays to the hyperoxic peak are also a consequence of this competition. The fl ow frequency parameter (f, Figure 3D) modulates the frequency of the underdamped ringing, but also affects the width of the initial oxygen overshoot. The fl ow damping parameter (τ, Figure 3E) directly affects the amount of ringing evident in the oxygen response; most measurements showed substantial ringing. The fl ow perturbation magnitude parameter (U 10 , Figure 3F) globally modulates the fl ow effects over the whole of the time course. Note that only when the fl ow is very small (Figure 3F, lower panel) does the intravascular oxygen concentration drop below the equilibrium concentration. Generally, our model does not predict an "initial dip" in the intravascular compartment.

FITS TO EXPERIMENTAL DATA
The model was capable of fi tting the diverse dynamics of the experimentally observed tissue oxygen variations, as illustrated by examples from four different cell recordings (Figure 4). Each row corresponds to a different cell. The leftmost column is the measurement and fi t corresponding to the grating orientation that evoked the maximal neuronal response ("preferred" stimulus), the center column for the orthogonally oriented grating, and right column for the preferred-orientation stimulus presented to the non-dominant eye.
Fits to the measurements highlighted in Thompson et al. (2003) are shown in Figure 4A, and corresponding intravascular fl ow and oxygen perturbations are shown in Figure 5A. During the stimulus period, the CMRO 2 increases promptly, competing with the comparatively sluggish infl ow of oxygenated blood. Consequently, there is an initial decrease in extravascular oxygen concentration that is subsequently overwhelmed by infl ow to produce an increase in oxygen concentration. Later, there is an undershoot caused directly by "ringing" in the fl ow perturbation: blood fl ow periodically decreases below baseline levels, and these changes cause similar decreases in oxygen concentration. As the fl ow perturbation damps out, the oxygen concentration returns to baseline.
Measurements in the vicinity of three other cells are also shown. Most exhibited stronger perturbation responses than those from Cell 20. While a few other cells also showed relatively strongly damped behavior (Figures 4B and 5B), strong "ringing" was often evident (Figures 4C,D). The accurate fi ts to this behavior are a direct consequence of the FEW fl ow model, which allows underdamped fl ow perturbations (Figures 5C,D).
The corresponding fi t parameters for these four cells are shown in Table 1. The parameters illustrate the diversity of the behavior and the fi ts. CMRO 2 perturbations vary from 0.07 to 0.80%; length parameter vary from 0.1 to 3.5 mm, damping times from 4.3 to 50 s, and fl ow perturbations from 5.1 to 60.1%. The frequency is comparatively stable, varying from 0.063 to 0.106. Despite the substantial variability of the fi t parameters, most of the fi t parameters show signifi cant correlations with neuronal activity (see below).
Our model successfully fi t ∼93% (156 out of 166) of the experimental measurements using the maximum-likelihood metric, and 92% (154/166) using the least-squares metric. Consideration of the fi ts to the entire dataset (Figures 1-21 in Supplementary Material), show an interesting variety of dynamics in the time courses. For example, fi ts to the oxygen measurements near Cell 3 (Figure 3 in Supplementary Material) exhibited the combination of relatively modest ringing and long length so that the deoxygenation phase decayed slowly and dominated the dynamics when the neuronal activity was strong. For this same cell, when the stimulus was switched to the non-dominant eye, the ringing became much more pronounced, suggesting that a different microvascular "circuit" had been engaged. Cell 10 ( Figure 10 in Supplementary Material) exhibits a faster response (modeled as a shorter length), and a weak fl ow response. The ability of the model to fi t such a variety of response waveforms gave us additional confi dence in its utility. This simple model, a combination of early oxygen diffusion in competition with underdamped fl ow oscillation, offers a quantitative mechanistic description for over a hundred detailed oxygen time series.
In fact, even the unsuccessful fi ts were still fairly good, with errors never exceeding the success threshold by more than 39%. Measurements near Cell 14 (Figure 14 in Supplementary Material) were least well fi t by the model, with fi ve out of the eight stimuli successfully fi t by the stricter (least-squared-error) metric. This cell exhibited a different dynamic: several early deoxygenation periods with the largest occurring later in time after the visual stimulus had ceased. It is notable that this cell also exhibited a high baseline-fi ring rate (5-7 spikes/s) compared to most other cells. Most of the other measurements that were not fi t satisfactorily did not return to zero toward the end of the measurement period; these measurements were also relatively noisy.
The parameter values of the successful fi ts varied substantially across the many measurements of the full dataset. The variations were substantial both cell-to-cell and within stimulus type. We examined the cell-to-cell variations by forming a "cocktail mean" for each cell, that is, the mean value of each parameter across all stimulus orientations, and then calculated the standard deviation of these cocktail means across the 21 cells. The standard deviations were: length, 45%; O 2 rate, 42%; fl ow frequency, 19%; damping time, 77%; fl ow magnitude, 61%. Within-cell variations across stimuli were obtained by normalizing each cell's set of parameter values by its cocktail mean, and calculating the standard deviation; we then calculated the mean within-cell standard deviation across the ensemble of cells. The values were as follows: length, 56%; O 2 rate, 33%; fl ow frequency, 15%; damping time, 51%; fl ow magnitude, 58%. With the exception of the length parameter, cellto-cell variations were somewhat larger than within-cell (stimulusto-stimulus) variations.

Ress et al. Transient oxygen-delivery model
We also examined the effects of varying four other parameters: volume ratio (α), equilibrium fl ow velocity (U 0 ), oxygen-extraction fraction (E), and equilibrium extravascular oxygen concentration (Q e00 ). In all cases, the results were similar for the quality of the fi ts to the extravascular oxygen measurements. However, adjusting these parameters had substantial effects on some of the fi ve "free" model parameters. Because of the uncertainty in the sampling volume for oxygen probe, we refi t the entire data set for α = 0.1 and 0.3, holding CBV constant at 0.06. Increasing the volume ratio decreases the "loading" of the extravascular volume on the diffusive losses from the intravascular volume. Consequently, the fi ts exhibit much greater length parameters (L) in order to match the temporal dynamics. Increasing α also causes a less pronounced decrease in the fi tted fl ow perturbation magnitude (U 10 ). Decreasing α has the reverse effect on both parameters. Each of the other three parameters was varied both above and below its nominal value by ∼15%, and the fi tting procedure repeated. The main effect of changing Q e00 was to scale the magnitude of the perturbations in fl ow (U 10 ) and CMRO 2 (Γ 10 ) required by the fi ts. Increasing Q e00 increases the equilibrium extravascular oxygen concentration near the end of the model cylinder. Consequently, the model fi ts require larger perturbations to fi t the measurements. Decreasing E has much the same effect as increasing Q e00 , but also reduces the longitudinal oxygen-concentration gradient. This, in turn, reduces the effect of convection relative to diffusion, changing the scaling between the two perturbations required to fi t the measurements. Increasing Transient oxygen-delivery model U 0 also increases the effect of convection relative to diffusion, and modifi es the scaling between the two perturbations. All of these changes slightly modify the length parameter (L), but do not much affect the fl ow frequency or damping parameters. Thus, although these parameters had little effect on the quality of the fi ts to the extravascular oxygen measurements, they can substantially affect the intravascular dynamics.

CORRELATIONS WITH NEURONAL ACTIVITY
We examined the correlation of the parameters for successful fi ts with neuronal activity in the following fashion. Spike-rate modulations were calculated as the difference between the mean baseline fi ring rate and the mean fi ring rate during the stimulus presentation period. We divided each cell's collection of spike-rate modulations and corresponding parameter values by their cocktail mean to compensate for cell-to-cell variations. We then plotted each normalized parameter value against the corresponding normalized spike-rate modulation (Figure 6). All of the fi ve parameters showed signifi cant correlations with neuronal activity. The length parameter, L, had the strongest correlation (R = 0.38, infi nitesimal p). indicating that stronger fi ring is associated with slower oxygen dynamics. The oxygen-transfer rate, Γ 10 , showed a weaker but still signifi cant correlation (R = 0.28, infi nitesimal p), consistent with the correlations reported by Thompson et al. (2003). A positive correlation (R = 0.27, p = 0.031) with neuronal fi ring was also evident for the decay time constant, τ, indicating that the "ringing" increases with fi ring-rate modulation. A weak correlation (R = 0.19, p = 0.026) was also observed for the normalized Each graph shows the variations for same three stimulus confi gurations presented in Figure 4: preferred, orthogonal, and non-dominant. Figure 4. Ress et al. Transient oxygen-delivery model fl ow magnitude, U 10 , suggesting that neuronal activity also drives the fl ow perturbation. Finally, the fl ow frequency shows a positive correlation (R = 0.18, infi nitesimal p), indicating that the ringing frequency increases with neuronal activity. However, the slope of this correlation is very small, so this is a weak effect.

DISCUSSION
Our model for transient tissue oxygen delivery in cerebral cortex successfully fi ts a large dataset of extravascular oxygen measurements corresponding to brief visual stimulations. The model offers several useful features. First, it provides a very concise description of oxygen transport as a competition between prompt local oxygen diffusion and a spatially distributed underdamped fl ow response.
Although elements of this model are similar to some previous models, the incorporation of the FEW fl ow model tremendously improves its explanatory power. Because the model offers accurate fi ts in a concise, fi ve-parameter framework, certain aspects of the hemodynamic response can be understood by examining variations of the model parameters with neuronal fi ring-rate. Second, the model offers a mechanism that unifi es the hemodynamic response with some components of the frequently observed phenomena of "spontaneous fl uctuations" or "vasomotion." This FEW fl ow model naturally incorporates underdamped fl ow perturbations, providing a simple mechanism for the observed late-time oscillations at ∼0.1 Hz that have not been well understood and modeled previously. Finally, the model also provides a description of the intravascular oxygen response that could yield some new insights into the hemodynamic response measured by BOLD functional magnetic resonance imaging (fMRI).
Our model provides mechanisms for the four temporal phases of the hemodyamic response: (1) The initial period of latency or dip observed in oxygen-sensitive hemodynamic measurements. The model gives a quantitative mechanism for this period as a competition between infl owing oxygenated blood, driven by the fl ow transient and delayed by propagation, and oxygen uptake. (2) The hyperoxic peak, also driven by the fl ow transient, which competes with and eventually overwhelms the oxygen consumption. (3) The undershoot, which is explained here as a brief fl ow decrease caused by the "ringing" dynamics of the FEW. (4) Subsequent oscillatory phenomena, which is explained as continued ringing in the fl ow model. Valabregue et al. (2003) suggested a model with similar geometry and oxygen-transport dynamics. Their model included the non-linear relationship between oxygen and oxyhemoglobin concentration in the intravascular compartment, as compared to the linear assumption of our own model. However, they did not include a fl ow model corresponding to a brief neural stimulation, nor did they attempt to fi t experimental measurements of the transient hemodynamic oxygen response. To our knowledge, we offer the fi rst attempt to fi t such a model to the slow, event-related measurements of the sort obtained by Thompson et al. (2003).
A model combining convection-diffusion oxygen transport with windkessel fl ow was fi t to experimental optical imaging data by Huppert et al. (2007) with excellent results. However, there were some notable differences between that effort and our own. First, they used a delayed-compliance windkessel model rather than the FEW. Second, they used a more complex model that involved many parameters, so a much more intricate fi tting procedure was required. Third, and most important, they applied their model to data from a fast event-related stimulus design and only fi t the fi rst 7 s of the experimental response. This approach has the effect of fi ltering out low-frequency variations that were so prevalent in the experimental results of Thompson et al. (2003). It is therefore diffi cult to compare the effi cacy of these two models across such disparate experimental data sets.
Our model does not utilize a non-linear fl ow-volume scaling relationship, such as Grubb's power-law relationship (Grubb et al., 1974), used in many other models (Buxton and Frank, 1997;Buxton et al., 1998Buxton et al., , 2004Hyder et al., 1998;Zheng et al., 2002Zheng et al., , 2005. Such scaling relations were obtained during experiments that took place on much slower time scales than the present experimental dataset. We suggest that the hemodynamics modeled here are distinct from cerebral auto-regulation, with the former mediated by microvascular control mechanisms on a relatively fast time scale (s), while the latter is mediated by smooth-muscle fi bers that surround larger arteries on a slower time scale (tens of seconds). If so, it is likely that the mechanisms work in parallel, and both must be considered during experiments that manipulate the hemodynamic response.
The mean parameter values for our model present an interesting picture of microvascular hemodynamics. The median length parameter is 0.9 mm, which is consistent with oxygen transport from arterioles as well as capillaries, as has been previously suggested (Sharan et al., 1989;Vovenko, 1999;Zheng et al., 2005), with upstream control points located along the arteriolar tree. This result is also consistent with tissue oxygen measurements from the lateral geniculate nucleus suggest that the point spread function of the fl ow response has a width of approximately 1.4 mm (Thompson et al., 2005). However, our length parameters are only rough estimates, because the simple geometry of our model does not take into account the variations in diameter and fl ow velocity that occur within the vascular tree. The median frequency parameter in the fl ow model, 0.10 Hz, was quite stable, suggesting that the pial vascular network behaves as a second-degree, roughly linear fi lter that is mechanically tuned to a particular frequency. Moreover, the median damping time of ∼8 s suggests that this network has an underdamped, bandpass character, with several cycles of "ringing" a normal aspect of the hemodynamic response.
The variation of our model parameters with fi ring-rate modulation provides some insight into the character of the hemodynamic response to brief neural stimulation. Our model indicates signifi cant correlations between both oxygen demand and the fl ow perturbation with local neuronal fi ring (Figures 6A,D). Both correlations are consistent with the usual notions of neurovascular coupling. The correlation of our length parameter L on fi ring-rate modulation ( Figure 6A) enriches this picture. The tendency toward larger length parameters with increasing fi ring-rate modulation suggests that stronger neuronal activity tends to attract blood fl ow from a broader region of vascular control. Finally, the signifi cant correlation of damping time with spike rate modulation ( Figure 6B) suggests that the observed underdamped ringing effects are an integral part of the hemodynamic response. Strong spike-rate modulation appears to drive strong ringing. Our fl ow model suggests an alternative explanation for some of the low-frequency oscillations or "vasomotion" observed with various techniques, e.g. (Katura et al., 2006;Kleinfeld et al., 1998;Koepchen, 1984;Mayhew et al., 1996;Spitzer et al., 2001), including the present tissue oxygen measurements (Thompson et al., 2003). Generally, these slow hemodynamic fl uctuations have been attributed to other mechanisms such as blood pressure feedback from baroreceptors or neurogenic mechanisms. These fl uctuations, sometimes called Mayer waves, are observed in three frequency bands. The so-called "low-frequency" (LF) band of such fl uctuations occurs at ∼0.1 Hz, making it particularly relevant to the measurements and model presented here. The FEW has a bandpass response to transient changes in cerebrovascular control. Thus, unstimulated or spontaneous neural activity would be fi ltered by the hemodynamic network to appear as low-frequency oscillations in fl ow, oxygen transport, and other aspects of the hemodynamics. This hypothesis is consistent with the use of low-frequency fl uctuations as a "resting-state" measurement of neural activity using fMRI (Biswal et al., 1997;Fransson, 2005;Peltier and Noll, 2002). The purpose of these oscillations is not yet understood. It has been suggested that low-frequency oscillations could increase oxygen transport effi cacy (Goldman and Popel, 2001;Kislukhin, 2004). Alternatively, the advantage of underdamping may simply be that of speed. Given our simple second-degree fl ow model, an underdamped fl ow response provides more immediate and copious delivery of oxygenated blood than a critically damped response.
It could also be that several mechanisms are responsible for these fl uctuations, with the FEW responsible only for fl uctuations in hemodynamics associated with transient neural activation, while the usual neural feedback mechanisms are responsible for oxygen auto-regulatory functions. It is also possible that the two mechanisms are related, with the tuned bandpass character of the cerebral vasculature actually a consequence of the putative neurogenic feedback mechanisms. Although in global scope, the mechanism of the LF fl uctuations may be complex and non-linear (Seydnejad and Kitney, 2001), at the level of the present cerebrovascular oxygen measurements they may be satisfactorily described by the FEW linear model. In this view, the dynamics of the FEW are an abstract model of more complex feedback mechanisms rather than a literal model of simple mechanical processes. Further experiments will be necessary to elucidate the relationship between the FEW and LF fl uctuations.
The CMRO 2 changes predicted by the model are much smaller than measured values [e.g., ∼10% observed using PET (Ito et al., 2005)]. This underestimate could refl ect the model's assumption that oxygen demand increases uniformly along the entire length of the model. In reality the demand could increase in a more localized fashion corresponding only to the active portions of the parenchyma. Moreover, this simple model represents the dynamics of the entire arteriolar tree with a single fl ow velocity, so it may not be possible to consistently match predictions for both velocity and CMRO 2 perturbations to measurements. Finally, this distortion could also be caused by our assumption of linear hemoglobin oxygen dissociation.
Our model also makes interesting but tentative predictions about the character of the intravascular fl ow and oxygen responses to brief neural stimulation. With the baseline set of fi xed parameters (α, Q e00 , E, and U 0 ), the predicted fl ow changes have a mean value of 26%, which is similar to measurements, e.g., the ∼30% fl ow changes obtained by positron emission tomography and MRI (Feng et al., 2004;Fox and Raichle, 1984). However, neglecting the non-linearity in oxygen dissociation could signifi cantly distort the model's detailed predictions of the intravascular dynamics. This could bear on the model's prediction of little or no "initial dip" in the intravascular oxygen concentration. Because the non-linearity "stiffens" the oxygen transfer from the erythrocytes to the vessel wall, the response to the initial period of oxygen demand may be larger than predicted.
The intravascular aspects of our model will form the basis for a model of the BOLD fMRI hemodynamic impulse response. To do so, the current model's predictions of oxygen concentration would be coupled to the additional dynamics of intravascular oxygen dissociation to obtain estimates of oxygenated and deoxygenated hemoglobin concentrations. These estimates, in turn, would be linked to a model of MR sampling to quantify the intravascular and extravascular components of the BOLD response. Finally, the current model's linear fl ow response system would form the basis for an estimate of the fMRI volume response component.

AUTHOR CONTRIBUTIONS
D. Ress conceived this model, performed the mathematical analysis, wrote the initial numerical solution and fi tting procedures, and performed the fi ts to the experimental data. J. Thompson collected the experimental data, organized it into a convenient form, and provided critical assistance on its use and interpretation. B. Rokers converted the initial Matlab code into C, tremendously increasing the speed of the fi tting process. R. Khan performed an