Microscopic Quantification of Oxygen Consumption across Cortical Layers

The cerebral cortex is organized in cortical layers that differ in their cellular density, composition, and wiring. Cortical laminar architecture is also readily revealed by staining for cytochrome oxidase – the last enzyme in the respiratory electron transport chain located in the inner mitochondrial membrane. It has been hypothesized that a high-density band of cytochrome oxidase in cortical layer IV reflects higher oxygen consumption under baseline (unstimulated) conditions. Here, we tested the above hypothesis using direct measurements of the partial pressure of O2 (pO2) in cortical tissue by means of 2-photon phosphorescence lifetime microscopy (2PLM). We revisited our previously developed method for extraction of the cerebral metabolic rate of O2 (CMRO2) based on 2-photon pO2 measurements around diving arterioles and applied this method to estimate baseline CMRO2 in awake mice across cortical layers. To our surprise, our results revealed a decrease in baseline CMRO2 from layer I to layer IV. This decrease of CMRO2 with cortical depth was paralleled by an increase in tissue oxygenation. Higher baseline oxygenation and cytochrome density in layer IV may serve as an O2 reserve during surges of neuronal activity or certain metabolically active brain states rather than baseline energy needs. Our study provides the first quantification of microscopically resolved CMRO2 across cortical layers as a step towards informed interpretation and modeling of cortical-layer-specific Blood Oxygen Level Dependent (BOLD) functional Magnetic Resonance Imaging (fMRI) signals.


Introduction
Dramatic improvements in functional Magnetic Resonance Imaging (fMRI) technology in recent years, including significant advances in hardware and image reconstruction (Polimeni & Wald 2018, Ugurbil 2018, Wald 2012, Yu et al 2014, have enabled resolution of cortical layers bringing us one step closer to the spatial scale of local neuronal circuits (Dumoulin et al 2018, Goense et al 2016, Polimeni et al 2010, Turner 2016. Indeed, neuronal circuits in cerebral cortex are organized in layers that differ in their cellular density, composition and wiring (Helmstaedter et al 2007, Woolsey & Van der Loos 1970, as well as the density of mitochondrial cytochrome oxidase, a marker of O 2 metabolism (Wong-Riley 1989). In primary cortices, the highest density of cytochrome oxidase is found in layer IV, and the lowest -in layer I (Land & Simons 1985, Weber et al 2008. Therefore, it is commonly believed that layer IV has higher cerebral metabolic rate of O 2 (CMRO 2 ) compared to other layers. These differences have been specifically hypothesized to reflect laminar variation in metabolic costs under the baseline (unstimulated) conditions (Weber et al 2008). Experimentally addressing this hypothesis is important for physiological interpretation of the Blood Oxygenation Level Dependent (BOLD) contrast used in cortical layer-resolved fMRI studies (Dumoulin et al 2018, Goense et al 2012, Hirano et al 2011, Siero et al 2015, Yu et al 2014. Experimental measurement of layer-specific (i.e. laminar) CMRO 2 has been challenging, because in common practice it requires information about both blood oxygenation and blood flow (Devor et al 2012a, Devor et al 2012b, Dunn et al 2005, Royl et al 2008. Previously, we introduced a method for extraction of CMRO 2 (Sakadzic et al 2010) based on a single imaging modality: 2-photon phosphorescence lifetime microscopy (2PLM) (Finikova et al 2008, Lecoq et al 2011, Sakadzic et al 2010 providing measurements of the partial pressure of O 2 (pO 2 ). This method relies on the Krogh-Erlang model of O 2 diffusion from a cylinder (Krogh 1919), which assumes that a volume of tissue within a certain radius around a blood vessel gets all its O 2 from that blood vessel. Previously, we applied our method to estimate CMRO 2 in the rat cerebral cortex, where capillaries were absent within an ~100-µm radius around penetrating arterioles satisfying the model assumption. The laminar profile of CMRO 2 , however, has not been addressed due to limited penetration of our imaging technology at that time.
In the present study, we increased our penetration depth by (i) utilizing a new pO 2 probe Oxyphor 2P with red-shifted absorption and emission spectra , (ii) switching from rats to mice since mice have a thinner cortex and smaller diameter of surface blood vessels attenuating light, and (iii) optimizing our procedure of probe delivery to avoid spilling of the probe on the cortical surface exacerbating out-of-focus excitation. With these improvements, we performed pO 2 measurements in fully awake mice in cortical layers I-IV.
In the mouse cortex, periarteriolar spaces void of capillaries are narrower than those in the rat (Kasischke et al 2011, Uhlirova et al 2016b. This challenges the basic Krogh-Erlang assumption of the center arteriole serving as the sole O 2 source. Therefore, we revisited the model to account for the contribution of the capillary bed and then applied the revised model to quantify the baseline laminar CMRO 2 profile. Our results show that, contrary to the common notion (Weber et al 2008), baseline CMRO 2 in layer IV does not exceed that in the upper cortical layers. We speculate that higher cytochrome oxidase density in layer IV may reflect higher O 2 metabolism during transient surges of neuronal activity or certain metabolically active brain states rather than baseline energy needs.

Depth-resolved measurements of tissue pO 2 in awake mice with chronic optical windows
We used 2PLM (Fig. 1A) (Finikova et al 2008, Lecoq et al 2011, Sakadzic et al 2010 in combination with a new pO 2 probe Oxyphor 2P  to image baseline interstitial (tissue) pO 2 in the primary cerebral cortex (SI) around diving arterioles at different depths: from the cortical surface down to ~500 µm deep. All measurements were performed in awake mice with chronically implanted optical windows ( Fig. 1B-C and Suppl Fig. 1). Oxyphor 2P was pressure-microinjected into the cortical tissue ~400 µm below the surface and was allowed to diffuse resulting in sufficient signal intensity in cortical layers I-IV within ~40 min after the injection. Sulforhodamine 101 (SR101) was co-injected with Oxyphor 2P to label astrocytes as a means for real-time monitoring of tissue integrity. In addition, fluorescein isothiocyanate (FITC)labeled dextran was injected intravenously to visualize the vasculature. We imaged pO 2 using either square or radial grids of 400 points in the XY plane covering an area around a diving arteriole up to ~200 µm in the radial distance. Overall, pO 2 measurements were acquired at 24 planes along 6 diving arterioles in 4 subjects.
We traced individual diving arterioles from the cortical surface and acquired grids of pO 2 points with an arteriole at the center ( Fig. 1D and Suppl Fig. 2). Figure 1E shows an example imaging plane crossing a diving arteriole. For each plane, a corresponding "reference" vascular image of intravascular FITC fluorescence was acquired immediately before and immediately after the pO 2 measurements for coregistration of the measurement points in the coordinate system of the vascular network. The measurements are color-coded according to the pO 2 level (in mmHg) and superimposed on the reference vascular image (Fig. 1F).
We grouped pO 2 measurements according to the following depth categories: layer I (50-100 µm), layer II/III (150-300 µm) and layer IV (320-500 µm). Figure 1G shows distributions of 7 tissue pO 2 values for the specific example shown in Figure 1E-F, and Figure 1H shows similar histograms for the overall dataset pooled across subjects. The high tail of the overall pO 2 distribution (the upper 35%), reflecting oxygenation in tissue in the immediate vicinity of arterioles, did not change significantly as a function of depth. This lack of dependence of the peri-arteriolar pO 2 on depth is in agreement with recent depth-resolved intravascular pO 2 studies in awake mice by us and others showing only a small decrease in the intravascular pO 2 along diving arteriolar trunks (Li et al 2019, Lyons et al 2016, Sencan et al 2020. In contrast, there was a significant increase in the low tail (the bottom 35%) of the pO 2 distribution with depth ( Fig. 1I, p<0.05) presumably indicating an increase in oxygenation of tissue in between capillaries in layer IV compared to layer I. This depth-dependent increase questions the common believe that layer IV has higher baseline CMRO 2 compared to other layers. In addition, the median of the pO 2 distribution shifted to the higher pO 2 values below 200 µm (Fig. 1H). The median increase can be explained, at least in part, by branching: at certain depths, diving arterioles branched leading to the presence of highly oxygenated vessels (arterioles and capillaries) in the vicinity of diving arteriolar trunks.

ODACITI model for CMRO 2 extraction
To quantify the laminar profile of CMRO 2 in mouse cerebral cortex, we revised our model for estimation of CMRO 2 based on peri-arteriolar pO 2 measurements. The vascular architecture in the cerebral cortex features the absence of capillaries around diving arterioles (Kasischke et al 2011, Sakadzic et al 2014, Sakadzic et al 2016. Previously, we and others have argued that this organization agrees with the Krogh-Erlang model of O 2 diffusion from a cylinder (Krogh 1919), where a diving arteriole can be modeled as a single O 2 source for the tissue in the immediate vicinity of that arteriole where no capillaries are present. Indeed, radial pO 2 gradients around cortical diving arterioles in the rat SI approached zero at the tissue radii corresponding 8 to the first appearance of capillaries, and application of this model to estimate CMRO 2 yielded physiologically plausible results (Sakadzic et al 2016). In the mouse, however, capillary-free peri-arteriolar spaces are smaller compared to those in the rat, and O 2 delivery by the capillary bed cannot be ignored (Moeini et al 2019). Therefore, we revised the model to explicitly include both arteriolar and capillary delivery as we describe below.
In the Krogh-Erlang model, a blood vessel -approximated as an infinitely long cylinderhas a feeding tissue territory with a radius . Outside of this cylinder, we enter a feeding territory of other vessels. Therefore, according to the Krogh-Erlang model, tissue pO 2 will monotonically decrease with r while moving away from the vessel until we reach . The minimum occurs at , where the first spatial derivative along the radial direction is zero ( . Outside of this cylinder, the predicted pO 2 for has no physiological meaning as the model is defined for . However, Krogh-Erlang formula for predicts rapid pO 2 increase, which is not compatible with the experimental observations of small tissue pO 2 changes in the capillary bed and it complicates the fitting procedure (Sakadzic et al 2016).
Importantly, in the cerebral cortex, not only arterioles but also capillaries contribute to tissue oxygenation, and this contribution is larger in awake animals (Li et al 2019, Sencan et al 2020 compared to those under anesthesia (Sakadzic et al 2014 Capillaries Into TIssue), is schematically illustrated in Figure 2A. As for the Krogh-Erlang model, the central arteriole (red) with the radius is feeding the peri-arterial tissue ( Fig.   2A, pink) with the radius . Beyond that, O 2 is supplied by the capillary bed so that the second spatial derivative is zero for ( Fig. 2A, blue). With this set of assumptions, we derived a new specific solution to the forward Poisson problem in cylindrical coordinates (see

Validation of the ODACITI model using synthetic data
To validate the ODACITI model, we used synthetic data where the ground truth (spaceinvariant) CMRO 2 was preset and thus known. First, we compared the functional form of the radial pO 2 gradient obtained with ODACITI to that obtained with the Krogh-Erlang equation. To that end, we generated synthetic data by analytically solving the Krogh-Erlang and ODACITI equation for a given CMRO 2 and . As expected, the models agree for producing the same monotonic decrease in pO 2 with increasing distance from the center arteriole. Beyond , the Krogh-Erlang formula produces an unrealistic increase in pO 2 , while ODACITI is able to take into account measured tissue pO 2 in the capillary bed and it produces more realistic, nearly constant level of pO 2 (Fig. 2B). Next, we tested ODACITI on synthetic data with added noise to mimic noise of experimental measurements. To that end, we (1) amended the Krogh-Erlang solution by imposing a constant pO 2 equal to pO 2 (R t ) for to represent the capillary bed, and (2) added Gaussian noise (σ=2). With these synthetic data, we were able to sufficiently recover the true CMRO 2 with ODACITI ( Fig. 2C). In another test, we performed fitting while selecting a constant between 60 -80 µm. This resulted in an 18% error in the estimated CMRO 2 (Suppl. Fig. 3). In comparison, the same variation in resulted in a 45% error while using a fitting procedure developed in Sakadzic et al. 2014 to fit for CMRO 2 based on the Krogh-Erlang model (Sakadzic et al 2014). Thus, ODACITI is less sensitive to the error in compared to the Krogh-Erlang model.
ODACITI assumes that the first spatial derivative monotonically decreases in the proximity of the arteriole and remains near zero in the capillary bed. In reality, however, this derivative varies in different directions from the center arteriole due to asymmetry of the vascular organization. Moreover, for some radial directions, pO 2 may remain high due to the presence of another highly oxygenated vessel (a small arteriolar branch or highly oxygenated capillary (Li et al 2019, Sakadzic et al 2014. To model this situation, we placed additional O 2 sources representing highly oxygenated blood vessels around the center arteriole (Fig. 2D).
When multiple nearby vessels serve as O 2 sources, no analytical solution for the tissue pO 2 map is available. Therefore, we solved the Poisson equation, which relates CMRO 2 and tissue pO 2 , by means of finite-element numerical modeling implemented in the software package  Fig. 2D illustrates an example FEniCS output for a center 15-µm-diameter "arteriole" surrounded by seven 5-µm-diameter "capillaries" located 80-130 µm away from the diving arteriole (Fig. 2D), which is a typical size of the region around diving arterioles void of capillaries in mouse cerebral cortex. The intravascular pO 2 was set to 70 mmHg and 35 mmHg for the arteriole and capillaries, respectively. In addition, we placed another "arteriole" with intravascular pO 2 was set to 60 mmHg 110 µm away from the center arteriole to simulate the occasional presence of highly oxygenated vessels within our measurement grid.
These simulated data were used to devise a procedure for segmenting a region of interest (ROI) consistent with the assumption that pO 2 should decrease monotonically from the center arteriole until reaching a certain steady-state value within the capillary bed (see Methods). We then used the data included in the ROI to extract the pO 2 profile as a function of the radial distance from the center arteriole. Finally, we applied the ODACITI model to these radial pO 2 profiles to calculate CMRO 2 . We calculated CMRO 2 in a noise-free case as well as with additive Gaussian noise (see Methods). In both cases, the model was able to recover the true CMRO 2 ( Fig. 2E-F).

Estimation of CMRO 2 across cortical layers
Following validation with synthetic data (Fig. 2), we used ODACITI to calculate CMRO 2 across the cortical depth keeping the same depth categories as in Figure 1. For each imaging plane, we registered a pO 2 measurement grid with the corresponding vascular reference image.
As with synthetic data, we segmented an ROI for each measurement grid and collapsed the data within each ROI into a corresponding radial pO 2 gradient ( Figure 3C shows laminar (depth-resolved) CMRO 2 profiles along 6 diving arterioles in 4 subjects. We used a mixed effects model implemented in R to quantify CMRO 2 across the layers. This analysis revealed a significant reduction of CMRO 2 with depth (p=0.001, Fig. 3D). For this calculation, was fixed at 80 µm (see Methods); allowing to vary while fitting for CMRO 2 did not reveal dependence of on cortical depth (not shown).
For each arteriole, we also acquired pO 2 measurements at the cortical surface ( Suppl   Fig. 4). The surface measurements, however, were not used for estimation of CMRO 2 because of violation of the axial symmetry assumption (i.e., absence of cerebral tissue above the imaging plane).
Graphically, CMRO 2 scales with the steepness of the descent of peri-arteriolar radial pO 2 profile (Kasischke et al 2011, Krogh 1919, Saetra et al 2020, Sakadzic et al 2016. In Figure   3E, we overlay grand averaged radial profiles for each layer. Each curve was obtained by averaging all data for that layer (Suppl Fig.5). As can be appreciated by visual inspection, the gradient corresponding to layer I has the steepest descent, while the descent of gradient in layer IV is more relaxed. This indicates that CMRO 2 in layer IV is lower compared to the upper layers, which is consistent with the model-based estimation shown in Figure 3D.
Taken together, our results indicate that an increase in the low tail of the pO 2 distribution with depth, which likely reflects an increase in oxygenation of tissue in between capillaries in layer IV compared to upper layers ( Fig. 1H-I), can be explained at least in part by a depthdependent decrease in CMRO 2 from layer I to layer IV ( Fig. 3D-F). This is in contrast to a wellestablished fact that capillary, cytochrome oxidase, and mitochondrial densities in the mouse SI all peak in layer IV (Blinder et al 2013, Santuy et al 2018 implying that neither of these densities can serve as a proxy for oxidative energy metabolism under baseline (unstimulated) conditions.

Discussion
In this work, we have achieved depth-resolved tissue pO 2 measurements in fully awake mice using 2PLM in combination with the O 2 -sensitive probe Oxyphor-2P. We devised the ODACITI model for estimation of CMRO 2 from peri-arteriolar pO 2 gradients accounting for both arteriolar and capillary O 2 supply. We validated the model using synthetic data and then applied it to estimate the laminar CMRO 2 profile during the baseline level of neuronal activity (i.e., in the absence of external stimulation). Our results demonstrate a significant reduction of CMRO 2 with depth from layer I to layer IV. While these results are based on a model and its specific assumptions, they strongly suggest that, at the very least, the baseline CMRO 2 in layer IV does not exceed that in the upper cortical layers.
In cerebral cortex, neuronal cell type distribution as well as cellular, synaptic and microvascular densities vary between cortical layers (Blinder et  in the respiratory electron transport chain located in the mitochondrial membrane (Weber et al 2008, Wong-Riley 1989. In SI, cytochrome oxidase distribution features a high-density band in cortical layer IV, which has been hypothesized to reflect higher CMRO 2 compared to other layers. Further, cytochrome oxidase was found to have a strong correlation with microvascular but not synaptic density -a finding that was interpreted to reflect high layer IV CMRO 2 in the "idling state" (Weber et al 2008), where about half of the energy expenditure reflects other processes than synaptic and spiking activity (Attwell & Laughlin 2001, Engl & Attwell 2015. The present measurements and calculations do not support the idea of high baseline CMRO 2 in layer IV. If so, what is the purpose of high cytochrome oxidase density in layer IV? In our recent study, we found that intravascular pO 2 changes (∆pO 2 ) in response to a sensory stimulus across layers were conserved (Sencan et al 2020) indicating a conserved ratio between the demand and supply. Taken together with existing data on larger stimulus-induced increases in CBF (∆CBF) in layer IV compared to other layers (Srinivasan & Radhakrishnan 2014), this finding implies that the laminar ∆CMRO 2 profile should also peak in layer IV in order for ∆pO 2 to remain invariant. Thus, cytochrome oxidase may reflect layer-specific differences in peak energetic demands during transient neuronal signaling events. These events may include not only taskinduced computation performed by local neuronal circuits (e.g., Barrel cortex response to a whisker touch (Feldmeyer et al 2013)) but also other dynamic neuronal processes occurring on larger spatiotemporal scales such as neuromodulation and sleep (Madsen et al 1991).
In this study, we observed a shift towards higher oxygenation of tissue within the capillary bed -quantified as the mean of the low tail of pO 2 distribution -with increasing depth.
Because our pO 2 grids were always placed around diving arterioles, the low tail of the pO 2 distribution may not reflect oxygenation of the capillary bed remote from arteriolar O 2 sources.
This observation, however, is potentially in line with our recent intravascular pO 2 study in awake mice, where we showed that mean capillary pO 2 in layer IV was ~15% higher than that in the upper layers (Li et al 2019). It is also consistent with a recent study where tissue pO 2 measurements were performed at different depths while a Clark-type polarographic electrode was advanced throughout cortical layers . We speculate that higher baseline oxygenation in layer IV may serve as an O 2 reserve for transient increases in neuronal activity.
In principle, higher capillary density in layer IV (Blinder et al 2013) may effectively shrink the peri-arteriolar cylinder where all O 2 is provided by the diving arteriole. This hypothetical scenario, however, would not affect the pO 2 gradient around the arteriole (which reflects CMRO 2 ). Rather, it would decrease the radius of this cylinder, . Although in our calculations of CMRO 2 was a fixed parameter, we show that sensitivity of ODACITI to varying was relatively low. In addition, allowing to vary while fitting for CMRO 2 did not reveal depthdependence.
The current study is part of our ongoing effort to improve 2PLM ( In principle, if we knew space-resolved tissue pO 2 as well as intravascular pO 2 for each blood vessel within that tissue, we could solve for CMRO 2 (that may vary in space and time) by finding an inverse solution for the Poisson diffusion equation (Saetra et al 2020). In practice, however, the signal to noise ratio (SNR) of our measurements may not be sufficient to accurately resolve pO 2 gradients around capillaries. In addition, our current spatial resolution limits simultaneous measurements of capillary and tissue pO 2 , and sampling pO 2 in 3D also remains a challenge. To mitigate these limitations, similar to the Krogh-Erlang model, the ODACITI model relies on radial pO 2 gradients around diving arterioles assuming cylindrical symmetry. In reality, highly oxygenated arteriolar branches and/or low branching order capillaries can often be present on one side of a diving arteriole requiring ROI segmentation.
Applied to synthetic data, our data analysis stream -including the segmentation algorithm followed by ODACITI estimation -was able to recover the "ground truth" CMRO 2 validating the method. In the future, further improvements in pO 2 probes and 2PLM technology may enable estimation of CMRO 2 in 3D using numeric methods without the need to segment the data (Saetra et al 2020).
Our measurements were limited to cortical layers I-IV. In 2-photon microscopy, depth penetration is fundamentally limited by out-of-focus excitation (Theer & Denk 2006). Because light is scattered and absorbed by tissue, the laser power is ought to increase with depth to overcome these effects and deliver a sufficient number of photons to the focal volume. Knowing layer-specific CMRO 2 in cerebral cortex is important for better understanding of the normal brain physiology as well as pathophysiology in diseases that affect cerebral microcirculation (Berthiaume et al., 2018;Iadecola, 2016Iadecola, , 2017Pantoni, 2010;Zlokovic, 2011).
It also is important for interpretation and modeling of the BOLD fMRI signals that are affected by both CBF and CMRO 2 (Buxton 2010, Gagnon et al 2015, Uhlirova et al 2016a. Of particular relevance to the current results, the BOLD response is sensitive not only to the balance of the ∆CBF and ∆CMRO 2 but also to their baseline state that can vary with age, disease, or even after consuming a cup of coffee (Griffeth et al 2011, Perthen et al 2008 altering the hemodynamic response for the same neuronal reality. Our study provides the first quantification of baseline CMRO 2 across cortical layers as a step towards informed interpretation and modeling of high resolution BOLD signals.

O 2 probe
The synthesis and calibration of the O 2 probe Oxyphor 2P were performed as previously

Two-photon imaging
Images were obtained using an Ultima 2-photon laser scanning microscopy system (Fig.   1A) (Bruker Fluorescence Microscopy). Two-photon excitation was provided by a Chameleon Ultra femtosecond Ti:Sapphire laser (Coherent) tuned to 950 nm.
For 2-photon phosphorescence lifetime microscopy (2PLM), laser power and excitation gate duration were controlled by two electro-optic modulators (EOMs) (350-160BK, Conoptics) in series with an effective combined extinction ratio > 50,000. This was done to mitigate bleedthrough of the laser light and reduce the baseline photon count (i.e., during gate-OFF periods), which was critical for accurately fitting phosphorescence decays.
We used a combination of Zeiss 5x objective (Plan-NEOFLUAR, NA=0.16) for a coarse At each imaging plane up to 400 µm below the cortical surface, a ~300 x 300 µm field of view (FOV) was selected for pO2 measurements. At each plane, pO 2 measurements were performed serially arranged in a square or radial grid of 400 points.
The phosphorescence was excited using a 13-µs-long excitation gate, and the emission decay was acquired during 287 µs. The photon counts were binned into 2-µs-long bins.
Typically, 50 decays were accumulated at each point with the total acquisition time of 15 ms per point. With 400 points per grid, one grid was acquired within 6 s. All points were revisited 20 times yielding a total of 1000 excitation cycles per point (50 cycles x 20 repetitions) acquired within 120 s. Dividing acquisition of 1000 cycles per point in blocks of 20 repetitions increased tolerance against motion (see Motion correction below) at a price of a slightly reduced sampling rate due to settling time of the galvonometer mirrors while moving from point to point.

Implantation of the cranial window and headpost
All animal procedures were performed in accordance with the guidelines established by the Institutional Animal Care and Use Committee (IACUC) at University of California, San Diego. We used 13 adult C57BL/6J mice of either sex (age: 4-7 months). Seven of them were rejected due to imperfect healing of the cranial window, and 6 were used for pO 2 measurements. Mice had free access to food and water and were held in a 12 h light/12 h dark cycle.
The surgical procedure for implantation of a chronic optical window was performed as  (Fig.1B and Suppl. Fig. 1). The window and the headpost were fixed to the skull in a predetermined orientation such that, when the mouse head was immobilized in the mouse holder ("hammock"), the window plane would be horizontal. Dextrose saline (5% dextrose, 0.05 ml) was injected subcutaneously before discontinuing anesthesia.
Post-operative analgesia was provided with buprenorphine (0.05 mg/kg subcutaneously) injected ~20 min before discontinuing anesthesia. A combination of sulfamethoxazole/trimethoprim (Sulfatrim) (0.53 mg/mL sulfamethoxazole and 0.11 mg/mL trimethoprim), and ibuprofen (0.05 mg/ml) was provided in drinking water starting on the day of surgery and for 5 days after surgery. Generally, full recovery and return to normal behavior were observed within 48 h post-op.

Delivery of optical probes
Oxyphor 2P and astrocytic marker SR101 were delivered by intracortical microinjection through the silicon port in the glass window ( Fig. 1B and Suppl. Fig. 1). (thereafter referred to as the "dye") and fixed in a micromanipulator at an angle of ~35°. The pipette was guided under visual control through the silicon port to its final location in tissue while viewing the SR101 fluorescence through the microscope eyepiece. First, to target the tissue surrounding a diving arteriole, the pipette was oriented above the window such that its projection onto the window coincided with the line connecting the diving point with the silicon port. The pipette was then retracted and lowered to touch the surface of the port. There, a drop of the dye solution was ejected to ensure that the pipette was not clogged. A holding positive pressure was maintained (using PV830 Picopump, World Precision Instruments) to avoid clogging of the pipette while advancing through the port. An experimenter could recognize the pipette emerging below the port by the dye streaming from the tip. At this point, the holding 21 pressure was quickly set to zero in order to avoid spilling of the dye on the cortical surface.
Next, the pipette was advanced for ~600 µm along its axis at 35°. Below ~100 µm, some holding pressure was applied to allow leakage of the dye that was used to visualize the pipette advancement. The pressure was manually adjusted to ensure visible spread of the dye without displacing cortical tissue. When the target artery was reached, the pipette was withdrawn by ~50 µm, and holding pressure was maintained for ~20 min to allow slow diffusion of the dye into the tissue. A successful injection was recognized by the appearance of SR101 in the perivascular space around the targeted arteriole but not around its neighbors. The contrast between the targeted and neighboring arterioles was lacking when the dye was injected too shallow or too vigorously, in which case the experiment was aborted. After the loading, the pressure was set to zero and the pipette was withdrawn. The mouse was placed back in its home cage to recover for ~40 minutes until the start of the imaging session. At that time, the dye usually diffused within a 300-600-µm radius around the targeted arteriole.
FITC-dextran (FD2000S, 2 MDa, Sigma Aldrich) was used to visualize cortical vasculature. Mice were briefly anesthetized with isoflurane, and 50 µl of the FITC solution (5% in normal saline) was injected retro-orbitally prior to intracortical delivery of Oxyphor 2P for each imaging session.

Habituation to head fixation
Starting at least 7 days after the surgical procedure, mice were habituated in one session per day to accept increasingly longer periods of head restraint under the microscope objective (up to 1 h/day). During the head restraint, the mouse was placed on a hammock. A drop of diluted sweetened condensed milk was offered every 20-30 min during the restraint as a reward. Mice were free to readjust their body position and from time to time displayed natural grooming behavior.

Estimation of pO 2
All image processing was performed using custom-designed software in MATLAB (MathWorks). For each point, cumulative data from all excitation cycles available for that point were used for estimation of the phosphorescent decay. Starting 5 µs after the excitation gate to minimize the influence of the instrument response function, the decay was fitted to a singleexponential function: , where N(t) is the number of photons at time t, N 0 is the number of photons at t=0 (i.e., 5 µs after closing the excitation gate), τ is the phosphorescence lifetime and x is the offset due to non-zero photon count at baseline. The fitting routine was based on the nonlinear least-squares method using the MATLAB function lsqnonlin. The phosphorescence lifetime τ was then converted into absolute pO 2 using an empirical biexponential form: (2) Where parameters A 1 , t 1 , A 2 , t 2 and y 0 were obtained during independent calibrations .
The locations of pO 2 measurements were co-registered with FITC-labeled vasculature.
Color-coded pO 2 values are overlaid on vascular images (grayscale) in all figures displaying pO 2 maps.

Identification of data corrupted by motion
To exclude data with excessive motion, an accelerometer (ADXL335, Sainsmart, Analog Devices) was attached below the mouse hammock. The accelerometer readout was synchronized with 2-photon imaging and recorded using a dedicated data acquisition system (National Instruments). During periods with extensive body movement (e.g., grooming behavior), the accelerometer signal crossed a pre-defined threshold, above which data were rejected. In pilot experiments, in addition to the accelerometer, a webcam (Lifecam Studio, Microsoft; IR filter removed) with IR illumination (M940L3-IR (940 nm) LED, Thorlabs) was used for video recording of the mouse during imaging. The video and accelerometer reading were in general agreement with each other. Therefore, accelerometer data alone were used to calculate the rejection threshold. Because every point was revisited 20 times, typically at least 10 repetitions (or 500 excitation cycles) were unaffected by motion and were used to estimate τ.
The MATLAB function lsqnonlin used to fit phosphorescent decays returned the residual error for each fit, that we plotted against the number of cycles (Suppl. Fig. 6). At around 500 cycles, the error stabilized at a low level. Therefore, we quantified points where at least 500 cycles were available.

Estimation of CMRO2
The O 2 transport to tissue is thought to be dominated by diffusion. In the general case, the relationship between pO 2 denoted as and the O 2 consumption rate denoted can be described by (3) where is the Laplace operator in three spatial dimensions, and D and α are diffusion coefficient and solubility of O 2 in tissue, respectively. For all our calculations, we assumed α = 1.39 µM mmHg −1 and D = 4x10 −5 cm 2 s −1 (Goldman 2008). This equation is only applicable outside the blood vessels supplying O 2 to the tissue. Oxygen supplied by a blood vessel is represented by a boundary condition of pO 2 imposed at the vessel wall. When the system is in steady state, the term can be neglected. If we also assume that there is no local variation of pO 2 in the vertical z-direction, that is, the direction along the cortical axis parallel to penetrating arterioles, Eq. 3 simplifies to (4) where P(r) represents pO 2 measured at the radial location r and refers to the twodimensional Laplace operator.
Here, we derive a specific solution to the forward problem of this partial differential equation using the following assumptions: First, as with the Krogh-Erlang model (Krogh 1919), we assume that tissue CMRO 2 in the vicinity of a diving arteriole is space-invariant. Second, we assume that within a capillary-free region around the arteriole, all O 2 is provided by that arteriole. That means that if we have an arteriole with the radius r=R ves and a region around the arteriole void of capillaries within a radius r=R t , for the radial distance R ves <r<R t , the rate of O 2 consumption is equal to CMRO 2 (that we are trying to find). This is because no local O 2 sources are present within this region except for the arteriole in the center. Third, we assume that for r>R t the rate of O 2 delivery by the capillary bed is uniformly distributed and equal to the rate of tissue O 2 consumption. With these assumptions, we can derive an analytic solution to Eq. 4 that we describe below.
We generally denote the partial differential equation in cylindrical coordinates as , where pO 2 (r) is the pO 2 as a function of the radial distance from the center of a diving arteriole, and H(r) is a Heaviside step function . Therefore, Eq. 5 considers three spatial regions , , and , where only in the middle region (e.g., for ) the right hand side of Eq. 5 is non zero. Additionally, we required that is a continuous function at (e.g., and at (e.g., and a smooth function at (e.g., ). Importantly, to better account for the O 2 delivery by the capillary bed that surrounds the arteriole, we assumed that for some , where . Finally, pO 2 values that we measured at did not show any diffusion gradients because they originated from extracellular dye in the tissue next to the vessel. Therefore, in our model we assumed that . We considered the following solution of Eq. 5 which satisfies the above conditions: and total O 2 flux through any cylindrical surface with length and radius r from the vessel is equal to the product of metabolic rate of O 2 and volume of the tissue that remains outside the cylinder with radius r: If we denote the O 2 flux in ODACITI as , total O 2 flux through the wall of the central arteriole, for the arteriolar segment of length , is given by (12) For any cylindrical surface with radius r between and , the total O 2 flux in ODACITI is given by (13) which represents a difference between the oxygen diffused from the vessel out ( ) and the oxygen consumed in the inner tissue cylinder ( ). Finally, total O 2 flux through a cylindrical surface with is constant: , which is in agreement with the model assumption that for r>R t , CMRO 2 is exactly equal to O 2 delivery rate of capillary bed. In ODACITI, if oxygen flux at r=R t is zero ( ), then , all oxygen delivered by the arteriole is consumed inside the r<R t tissue volume (e.g., ) and ODACITI scales back exactly to the Krogh-Erlang model with additions that pO 2 = const for r<R ves and r>R t . However, in ODACITI it is possible that the radius where can be smaller than , which can better account for the influence of the capillary bed. In addition, this model allows to include more measured data points from small (r<R ves ) and large radii (r>R t ) into the fitting procedures as compared to using the Krogh-Erlang model.

Segmentation of regions of interest
ODACITI describes pO 2 gradients around diving arterioles for an ideal case of radial symmetry. In practice, however, capillary geometry around diving cortical arterioles is not perfectly radially symmetric. The model also assumes that pO 2 decreases monotonically with increasing radial distance from the arteriole until reaching a certain stable level within the capillary bed. In reality, however, at certain distance pO 2 may increase again due to presence of another highly oxygenated blood vessel that can be (a branch of) a diving arteriole or a postarteriolar capillary. To mitigate this issue, we devised an algorithm for automated segmentation of a region of interest (ROI) for each acquired pO 2 grid. First, tissue pO 2 measurements were co-registered with the vascular anatomical images to localize the center of the arteriole. Next, 20 equally spaced radii were drawn in the territory around a center arteriole and a pO 2 vector was interpolated as a function of the radial distance from the arteriole. These vectors were spatially filtered using the cubic smoothing spline csaps MATLAB function with smoothing parameter p=0.001 to reduce the effect of noise in the measurements. For each vector, we calculated the first spatial derivative and included all pO 2 data points with radii where . Finally, we added all pO 2 measurements within the low 35% tail of each pO 2 map to include capillary bed data points. Taken together, these steps resulted in exclusion of data points around additional oxygen sources within the capillary bed deviant from the model assumptions. Included pO 2 data were then collapsed to generate a radial gradient for estimation of CMRO 2 with the ODACITI model for each peri-arteriolar grid. When multiple acquisitions of the same plane were available, CMRO 2 estimates were averaged.

Synthetic data
To validate estimation of CMRO 2 using ODACITI, we generated synthetic pO 2 data by solving the Poisson equation for a given (space-invariant) value of CMRO 2 and chosen geometry of vascular sources and measurements points. This was done numerically using the finite element software package FEniCS (Logg et al 2012) as described in our recent publication (Saetra et al 2020 Experimental data, in contrast, are measured on a Cartesian grid. Therefore, we transferred the synthetic data generated by FEniCS to a Cartesian grid similar to that used in our experiments. Additive Gaussian noise was implemented using the normrnd MATLAB function. For each synthetic data point in the grid, we drew a random number from a Gaussian distribution with the pO 2 value at this point as the mean (µ) and standard deviation of σ=2. Afterwards, we replaced the pO 2 value by (µ + σ).

Statistics
Statistical analysis was performed using The R Project for Statistical Computing (www.rproject.org) using a linear mixed-effect model implemented in the lme4 package, where trends observed within a subject (e.g., the dependence of CMRO 2 on depth) were treated as fixed effects, while the variability between subjects was considered as a random effect. Observations within an animal subject were considered dependent; observations between subjects, independent. 900SP and 945SP -short pass optical filters with a cutoff at 900 and 945 nm, respectively;

B.
Schematics of the chronic cranial window with a silicon port for intracortical injection of Oxyphor 2P.
C. An image of surface vasculature calculated as a maximum intensity projection (MIP) of a 2photon image stack 0-300 μm in depth using a 5x objective. Individual images were acquired every 10 μm. Fluorescence is due to intravascular FITC.
D. An example set of images tracking a diving arteriole throughout the top 500 μm of cortex (red arrows). Fluorescence is due to intravascular FITC. C. Application of ODACITI to synthetic data with added Gaussian noise (σ=2). In these data, pO 2 = P ves for r< R ves and pO 2 = pO 2 (R t ) for r>R t . For R ves <r<R t , we solved for pO 2 using the Krogh-Erlang equation. Three cases with CMRO 2 of 1, 2 and 3 μmol cm-3 min-1 (color-coded) are superimposed. The ODACITI fit (lines) are overlaid on the data (points).

D.
Simulated pO 2 data generated by solving the Poisson equation in 2D for a given geometry of vascular O 2 sources, including one highly oxygenated vessel (on the right) and given CMRO 2 =2 μmol cm-3 min-1.
E. The pO 2 gradient as a function of the distance from the center arteriole. The green line shows ODACITI fit.

Figure 3. CMRO 2 estimation across cortical layers
A. An imaging plane 100 μm below the surface; a grid of pO 2 points is overlaid on the vascular image. pO 2 values were interpolated along the radial yellow lines. The yellow shaded area extends along each direction until the first derivative becomes zero. Segmented ROI used for extraction of the radial pO 2 profile is shown on the right. This ROI includes in addition values within the low tail of the pO 2 distribution.
B. The radial pO 2 profile extracted from (A). ODACITI fit is overlaid on the data. The fitted CMRO 2 value is 1.5 μmol cm -3 min -1 .
C. Estimated CMRO 2 for the entire dataset of 24 planes along 6 diving arterioles in 4 subjects.
Measurements along the same arteriole are connected in a line. Subjects are color-coded.

D.
Quantification of CMRO 2 across layers using the data in (C). Error bars show standard error calculated using a linear mixed effect model implemented in R.
E. Superimposed color-coded grand averaged radial pO 2 profiles for layers I, II/III and IV (mean±SEM for 5-µm bins). For each layer, the profile was obtained by averaging all measurements within that layer across arterioles and subjects.

Supplementary Figure 1. Optical window with a silicon access port
A. Top view of the window before filling the injection port with silicon. There are three 3-mm glass coverslips stacked together and glued to a 5-mm glass coverslip. The 3-mm stack is beveled to guide the pipette at an angle through the port.

B.
A side view of the window; the port is filled with silicon.

C.
A top view of an implanted window; surface blood vessels are visible under the glass.
D. An image of surface vasculature within the same window calculated as a maximum intensity projection (MIP) of a 2-photon image stack 0-300 μm in depth using a 5x objective. Individual images were acquired every 10 μm. Fluorescence is due to intravascular FITC. C. Phosphorescence decays for two points labeled in (B). Data (points) and fit (lines) are overlaid. Lower pO 2 corresponds to slower decay (blue).

D.
Residual error over time.
E. The squared residual norm (output of lsqnonlin) plotted against the number of excitation cycles that were used to sum the photon counts. .