Real-time tracking of brain oxygen gradients and blood flow during functional activation

Abstract. Significance Cerebral metabolic rate of oxygen (CMRO2) consumption is a key physiological variable that characterizes brain metabolism in a steady state and during functional activation. Aim We aim to develop a minimally invasive optical technique for real-time measurement of CMRO2 concurrently with cerebral blood flow (CBF). Approach We used a pair of macromolecular phosphorescent probes with nonoverlapping optical spectra, which were localized in the intra- and extravascular compartments of the brain tissue, thus providing a readout of oxygen gradients between these two compartments. In parallel, we measured CBF using laser speckle contrast imaging. Results The method enables computation and tracking of CMRO2 during functional activation with high temporal resolution (∼7  Hz). In contrast to other approaches, our assessment of CMRO2 does not require measurements of CBF or hemoglobin oxygen saturation. Conclusions The independent records of intravascular and extravascular partial pressures of oxygen, CBF, and CMRO2 provide information about the physiological events that accompany neuronal activation, creating opportunities for dynamic quantification of brain metabolism.


General information
All solvents and reagents were purchased from standard commercial sources and used as received.

Oxygen quenching properties of Oxyphor PtR4
The oxygen quenching plots of Oxyphor PtR4 were measured in buffered aqueous solutions (20 mM phosphate, pH 7.2) in the absence and presence of bovine serum albumin (BSA, 5%) and found to be nearly identical (Fig. S2). The data were fitted to an empirical bi-exponential form (shown in Fig. S2), and its coefficients were used to convert the phosphorescence lifetime measured in vivo to pO2. The phosphorescence lifetime of Oxyphor PtR4 changes from ~8 μs at air saturation to ~43 μs at zero-oxygen at (36.8°C).

Two-channel phosphorometer
The instrument for two-color phosphorimetric measurements was constructed in house. The excitation sources were modulated laser diodes (Power Technology) operating at λmax=630 nm (15 mW) and λmax=517 nm (10 mW) for excitation of PtG4 and PtR4, respectively. Both diodes have the rise time of ~50 ns and can be controlled by standard TTL signals. The detectors were avalanche photodiodes (C12703-01, Hamamatsu; rise time 3.3 µs). The control of the data acquisition was performed using a digital-to-analog (DA)/analog-to-digital (AD) board (NI USB-6351, National Instruments; 1 MHz digitization frequency). Figure S3. Optical absorption and emission spectra of the probes and associated laser lines and optical filters in Channels 1 and 2 of the phosphorometer. The spectral ranges seen by the detectors (APDs) are shown by shaded areas. 'lp' is the abbreviation for 'long-pass filter'.

S9
The detectors, the digital board and the power supplies were located inside the instrument box, which was connected to a control computer by a USB cable. The front panel of the box contained two terminals for the cables conducting power and signals to the laser diodes and two optical ports for the optical fibers conducting the emitted light (phosphorescence) to the APDs, which were located inside the instrument box. The acquisition control and data analysis were performed by a program written in C/C++ (Qt, Nokia).
Each laser diode was placed inside a plastic holder tube, equipped with a lens (f=60 mm, ∅=12 mm), such that the beam could be focused into a spot (∅~0.2 mm) ~6 cm away from the tube's end. A single core plastic fiber (∅=4 mm, Fiberoptic Technology) ran on the side of the tube and held ~2 mm away from the excitation focus. The emitted light that entered the fiber was forwarded to an optical port in the instrument box. Inside the port the fiber tip was positioned next to a spherical lens (∅=10 mm) for light collimation, after which the light passed through a set of optical filters (∅=2.54 cm) and focused by another spherical lens (d=10 mm) on the entrance aperture of the APD (∅=3 mm). The optical filters configuration in the channels is shown in Fig. S3.

Laser speckle imaging optics
Laser speckle images were collected with an infinity-corrected optical system composed of two lenses with focal lengths of f=135 mm (Mitakon Zhongyi Mark II Lens for Nikon F, Zhong Yi Optics, 1.8 times.

Stimulus protocol
Each stimulus cycle was 1 min-long and it consisted of 4s-long of baseline collection period, 4slong electric forepaw stimulation period and 52s-long post-stimulus resting period (Fig. S4). Figure S4. Timeline of the forepaw stimulation cycle.

Modeling of oxygen diffusion and calculation of CMRO2
Oxygen  O2 is delivered to the activation region by the blood from the descending arteriole. From the arteriole, the blood is dispersed into a capillary network and is ultimately collected by an ascending venule.
O2 dissolved in the blood plasma is in equilibrium with O2 bound to hemoglobin, and the equilibrium is described by the Hill equation. 52 Oxygen transport in and out of the capillary network is modeled by Eq. S1 (also Eq. 1 in the main text).
The rate of change in the concentration of oxygen molecules in the capillaries (Cc), including bound and unbound O2, is given on the left-hand-side (LHS) of Eq. S1. Here V i denotes the intracapillary volume, which is typically given per 100g of tissue, i.e. mL/100g tissue.
The first term on the right-hand-side of Eq. S1 describes the influx of oxygen into the capillaries, where CBF(t) is the blood flow at time t. C a and C v are the arterial and venous concentrations of oxygen, respectively, and they include both bound and unbound O2. The second term represents the rate of extraction of oxygen from the capillary compartment into extravascular tissue. PS c is the O2 mass transfer proportionality parameter that is often written as the product of the capillary O2 permeability, P [cm/s], and capillary surface area, S c [cm 2 /100g tissue]. Note, however, that the precise microscopic assignment of P and S c is still debated, and some researchers 49 prefer to define and use another proportionality constant (KTO2) which does not impose such a strict microscopic interpretation. C i [µM] and C e [µM] represent the S12 volume-averaged concentrations of dissolved (unbound) oxygen in the intravascular and extravascular compartments, respectively.
The extraction rate of oxygen from the capillary is proportional to the concentration gradient between the intravascular and extravascular compartments. Thus, on general grounds, the capillary permeability ( P ) is expected to be proportional to the oxygen diffusion coefficient and inversely proportional to a characteristic length scale on the order of the capillary radius (but not necessarily equal to it).
A similar equation describes the rate of change in the extravascular oxygen concentration dCe/dt (also see Eq. 2 in the main text): where V e (also defined per 100g of tissue] denotes the extravascular volume, so that V e = V t − V i , where V t is the total tissue volume. 53 In Eq. S2, the net-flow of oxygen includes its diffusion from the intravascular compartment (first term on RHS) and its consumption in tissue (second term on RHS). The tissue oxygen consumption rate is called the Cerebral Metabolic Rate of Oxygen (CMRO 2 ), and it is typically measured in µmole/s/100g tissue. The variables and parameters in Eqs. S1 and S2 are summarized in Table S1. S13  S6) specifies volume-averaged concentrations C i and C e and provides a form for the dependence of oxygen concentration on the the radial distance (ρ), C(ρ). For a single cylinder, the capillary radius is denoted as S14 Figure S5.
Here, M is the CMRO 2 density in the extravascular volume: M = CMRO 2 /V e . The model assumes that M and D are spatially uniform in the extravascular volume. The model also assumes that the intracapillary oxygen concentration is also uniform, which effectively implies that oxygen inside the capillary (including the capillary wall) equilibrates rapidly throughout the volume compared to the equilibration throughout the S15 extravascular compartment, and that the oxygen transport properties of the capillary wall are the same as those of the capillary interior.
As a boundary condition, we assume that the diffusive O2 flux at the extravascular outer cylinder boundary is zero: dC dρ � ρ = ρ 2 = 0. We also assume that the concentration on the capillary wall edge, C(ρ 1 ) = C i . The steady-state solution of this problem is given by Krogh 18 and, in more detail, by Popel: 54 The volume-averaged extravascular, C e , is then found by computing the spatial average of C(ρ) over the extravascular compartment volume.
An exemplary profile of the partial pressure of oxygen in the extravascular cylinder is shown in Fig. S7. The volume-averaged pO2 is also indicated. Figure S7. Oxygen partial pressure (pO2) profile along the radial distance from the wall (ρ1=3 µm) to the outer surface of the tissue cylinder (ρ2=30 µm) in the Krogh cylinder model. CMRO2 and D are set to be 3.4 µmol/s/100g tissue and 1.7x10 -5 cm 2 /s respectively. The volume-averaged pO2 is 22 mmHg.

S16
From Eq. S5, we can see that one parameter that we experimentally measure, C e , is predicted to be an explicit function of morphological geometric factors ρ 1 and ρ 2 , as well as of the physiological factors M and D. For fixed capillary volume fraction and capillary radius, the average O2 concentration is thus a function of M and D.
It is useful to rewrite Eq. S5 in other ways. For example, the expression below reveals linear relationship between M and the oxygen concentration difference between the compartments.
If we substitute C(ρ 1 ) for C i , and replace M with CMRO 2 /V e in Eq. S6, we obtain an expression for the CMRO 2 in a steady-state in two-compartment model (Eq. S7).
Equation S7 relates the experimentally measured paramerers, i.e., C i and C e , to CMRO 2 and PS c .
Importantly, PS c depends on microvascular parameters that have been previously measured or estimated, i.e., ρ 1 , ρ 2 , D, and L. (Note, L is a function of V t and ρ 2 ).
Notably, the Krogh-based model gives an explicit formula for PS c : The function G(ρ 2 /ρ 1 ) in Eq. S8 is a function of the ratio of inter-capillary distance to capillary diameter. Therefore, with caveats related to the approximations/simplifications outlined above, PS c represents the oxygen mass transfer coefficient connecting CMRO 2 to the oxygen gradient across intra-and extravascular tissue compartments. Moreover, within this model, P and S c can be separately defined. For example, if we take Sc to be the capillary surface area, i.e., S c = 2πρ 1 L, then we generate an oxygen permeability associated with tissues surrounding the capillary wall, i.e., P = D/[ρ 1 G(ρ 2 /ρ 1 )] . The function, G(ρ 2 /ρ 1 ), is not significantly larger than unity for typical microvascular parameters, e.g., it S17 ranges from ~1.1 to ~1.6 for typical capillary diameter and separation when V i /V t is between 1 and 3 percent.
The Krogh cylinder-based model also provides an explicit formula for the extravascular compartment exponential relaxation time constant, τ = V e /PS c , in terms of ρ 1 , ρ 2 and D. This time constant arises naturally in Eq. S2 and can be important for the description of metabolism dynamics at early times. Interestingly, the relationship between τ and the mean time, Δt, for oxygen to diffuse the distance 2 is well defined in the model; the ratio τ/∆t ranges from approximately 2.1 to 3.2 for typical system parameters, and depends only weakly on ρ 2 /ρ 1 in the physiologically relevant experimental regime (i.e., V i /V t between 1 and 3 percent). Thus, the time for the extravascular compartment to reach a new steadystate following an impulse change in oxygen concentration is 2x-3x longer than the mean time it takes for oxygen to diffuse to the extravascular tissue boundary, ρ 2 .
Limitations of the model. As many biophysical models which give exact results, the present model has important limitations. For example, besides the specific numerical approximations noted above, the model ignores the heterogeneous nature of the capillary network. In practice, the capillary size and especially the inter-capillary separation is not uniform. A more complex theory should average over the distribution of regions with different morphologies. In that case, even if we continue to adopt the core cylinder model, a slightly more complex theory would quantitatively predict a distribution of PS c and τ.
Accounting for these distributions could change the predicted dynamics, i.e., compared to replacing distributions by their average values.
Another limitation derives from the fact that tissue has both neurons and other cells, and they likely have different temporal metabolic responses to stimulation. Still another limitation is the description of the capillary wall as a uniform extension of the inner capillary. For example, the oxygen diffusion coefficient in the capillary wall is likely to be different from that in the surrounding space. All of these factors need to be considered in data interpretation. Notably, our new methodology will make confrontation of these issues possible.

Selection of PSc
Eq. S2 enables computation of CMRO 2 , however the computed profiles are dependent on the choice of parameters V e and PS c . The volume V e is known from experiments. 41,42,53 However, the value of PS c , which affects both CMRO 2 magnitude and dynamics, is not well-known, which reflects the fact that S18 measurements of oxygen diffusion coefficient in tissue us notoriously difficult. PS c values have been measured, 55 computed, 35, 49, 56 and derived by combining measurements and calculations. 57,58,59,60 According to a recent review, 49 the reported value of PS c vary in the range of 35-230 mL/s/100g tissue.
Below we describe how we chose PS c from the available literature data, including teh reported measurements of morphological and physiological parameters.
According to Eq. S8, the relevant morphological parameters are the capillary radius (ρ 1 ), one-half the inter-capillary separation (ρ 2 ) and the volume fraction of intravascular tissue (V i /V t ). These parameters also determine L and V e . Reported values of ρ 1 vary from 2 to 3.5 μm. 60-64 Reported values of ρ 2 vary from 10 to 30 μm, corresponding to average inter-capillary distances ranging from 20 to 60 μm. 39, 61, 65 Reported intravascular volume fraction, V i /V t , in tissue ranges from 1% to 3%. 41,61,62,64,66 Tables S2-S4 summarize the range of possible morphological parameters (ρ 1 , ρ 2 , and L) for V i /V t equal to 0.01, 0.02, and 0.03. Note, all of these reported morphological parameters were measured in rats.

S19
The other relevant physiological parameters are the oxygen diffusion coefficient, D , and the baseline (normal) CMRO 2 . Reported values of D vary from 1.4 to 1.9 [×10 -5 cm 2 /s] at 37°C. 67-70 Note, in the cited papers objects varied from small rodents to a human brain tissue. The reported baseline CMRO 2 of rat brain ranges from 3 to 7.4 [μmol/s/100g tissue]. 44,45,71,72 The cylindrical model for capillary and surrounding tissue enables determination of PS c based on the full range of reported ρ 1 , ρ 2 , V i /V t , and D and using Eq. S8. The resulting "allowed" variation of PS c is large, spanning from 137 to 2470 [mL/s/100g tissue]. Tables S5-S7 show the values of PS c (in bold) as a function of ρ 1 , ρ 2 , V i /V t , and D; note ρ 2 is fixed by choice of ρ 1 and V i /V t . This potential variation of PS c is much larger than reported in the literature.

S21
To further constrain PS c , we next utilize prior measurements of baseline (normal) CMRO 2 and our measurements of intravascular and extravascular pO2. Specifically, our measurements found baseline piO2=33.5 mmHg and peO2=21.3 mmHg. Using the reported range of baseline CMRO 2 and Eq. S7, i.e., CMRO 2 = PS c (C i − C e ), we identify a much smaller subset of PS c 's, which satisfy Eq. S7 subject to constraints of tissue morphology and oxygen diffusion. These results are summarized in Tables S5-S7, which show the predicted CMRO 2 as a function of ρ 1 , ρ 2 , V i /V t , and D for the cases where 3 µmol/s< CMRO 2 < 7.4 µmol/s. The green shaded rectangles in the Tables correspond to the allowed values of CMRO 2 . Note, PS c is given in bold next to CMRO 2 in these green shaded rectangles. Ultimately, to ameliorate this problem it will be valuable to carry out more experiments to further constrain the possible values of PS c .

Timing analysis
To analyze the dynamics of stimulation-induced traces of CBF, piO2, peO2 and CMRO 2 , six timepoints were sought: a "trigger" time (ttrigger), the 20%, 50% and 90% rise-time from baseline (t20%, t50%, t90%), the time at-peak (tpeak), and the 50% fall-time from peak (tfall). To obtain these time points we averaged the 13 single-stimulation traces. Further, to reduce noise in the single-shot stimulation data, we smoothed each single-shot raw data trace utilizing a low-pass-filter de-noising technique. The low-pass filtering kernel was a sliding Gaussian function with a full-width-at-half-maximum of 0.3 s (twice the measurement time-step). After application of the filter, we computed a mean trace from the 13 responses as shown in Fig. 4 in the main text. We also computed a trace of the standard error (standard deviation of the mean) for additional temporal error analysis. Exemplary raw and de-noised data are shown in Fig. S8. Figure S8. De-noised single-stimulation traces. Each trace is low-pass filtered (FWHM=300 ms). In each plot, the raw data is marked in gray. Notice, the peO2 trace is nearly noise-free compared to other traces.

S22
For temporal analysis we adopted an approach published previously. 43 The de-noised average data were normalized, so that the baseline is set to 0, and the peak corresponds to 1 (Fig. 4, main text). The standard error (standard deviation of the mean) is accordingly normalized, i.e., divided by the peak magnitude. Our determination of the trigger time/height sought to identify a peak height that was clearly larger than measurement noise. To this end we first computed the standard deviation ( ) of full-dynamic CMRO 2 during the 4s-long baseline period prior to activation. We chose full-dynamic CMRO 2 for this task because it was our noisiest signal. The resultant baseline standard deviation was = 0.03. We then defined the trigger time/height of full-dynamic CMRO 2 to be the point on the trace where the signal had definitively increased (out of the noise) to 3 = 0.09. Thus, in practice, ttrigger is determined from the 9% rise timepoint with respect to the baseline (t9%) for full-dynamic CMRO 2 and all other traces. Note, the term trigger applies directly to CMRO 2 , CBF, and piO2 traces but not peO2. Determination of t20%, t50%, t90%, tpeak, tfall, and ttrigger= t9% was straightforward. For example, t20% is determined by the intersection of normalized trace with the horizontal line displaced vertically by 0.2 from the horizontal axis. We assigned the timing uncertainties following a standard approach. 73 The results of the timing analysis are summarized in Table   S8. 11. Records of CBF and ABP during stimulation Figure S9. Concurrent records relative CBF (rCBF) and arterial blood pressure (ABP) in a representative experiment.

Potential measurement artifacts: oxygen consumption, probe distribution
Potentially, oxygen measurements by phosphorescence quenching can lead to oxygen consumption.
The consumption is not due to the probe per se, but rather to scavenging of singlet oxygen, a by-product of the triplet quenching reaction, by other organic molecules in the medium. Singlet oxygen is indiscriminate in choosing its substrate, and the probe is certainly one of its potential targets. However, the probe is usually present in micromolar concentrations, while the "concentration" of other organic molecules in tissue (lipids, proteins etc) is in molar range. Hence, virtually all singlet oxygen that gets scavenged, is scavenged by biological matter.
Nevertheless, it is important to realize that the lifetime of singlet oxygen in water is only ~3 µs.
Thus, most of the singlet oxygen generated in blood plasma or interstitial space, decays back to the ground state (triplet O2) faster than it interacts with organic molecules. Only at very high light fluxes, high probe concentrations and in systems with limited supply of oxygen one might observe oxygen depletion upon phosphorescence quenching measurements. It was clearly not the case in the experiments reported in this work. We performed control measurements where a solution of the probe (Oxyphor PtG4, 2 µM) in 1M glycerol, 5mM glucose, 5mM glutathione (GSH), and 10mM HEPES (pH 7.2) and ~1 mM bovine serum albumin (BSA) -a solution mimicking cellular milieu, -in a sealed vial was exposed to the same measurements as in or in vivo experiments, showing that oxygen consumption was not a concern.
Probe concentration does not affect phosphorescence lifetime measurements per se, but it does affect the signal-to-noise ratio. In a highly heterogeneous environment, such as tissue, it is possible that within the volume sampled by light more of the probe would somehow end up in a more oxygenated region vs less oxygenated region or vice versa, i.e. a situation would occur where a probe gradient would affect the average measured pO2 value. However, such a distribution would occur randomly, and since in our S29