Spatially resolved estimation of metabolic oxygen consumption from optical measurements in cortex

Abstract. Significance: The cerebral metabolic rate of oxygen (CMRO2) is an important indicator of brain function and pathology. Knowledge about its magnitude is also required for proper interpretation of the blood oxygenation level-dependent (BOLD) signal measured with functional MRI. Despite the need for estimating CMRO2, no gold standard exists. Traditionally, the estimation of CMRO2 has been pursued with somewhat indirect approaches combining several different types of measurements with mathematical modeling of the underlying physiological processes. The recent ability to measure the level of oxygen (pO2) in cortex with two-photon resolution in in vivo conditions has provided a more direct way for estimating CMRO2, but has so far only been used to estimate the average CMRO2 close to cortical penetrating arterioles in rats. Aim: The aim of this study was to propose a method to provide spatial maps of CMRO2 based on two-photon pO2 measurements. Approach: The method has two key steps. First, the pO2 maps are spatially smoothed to reduce the effects of noise in the measurements. Next, the Laplace operator (a double spatial derivative) in two spatial dimensions is applied on the smoothed pO2 maps to obtain spatially resolved CMRO2 estimates. Result: The smoothing introduces a bias, and a balance must be found where the effects of the noise are sufficiently reduced without introducing too much bias. In this model-based study, we explored this balance using synthetic model-based data, that is, data where the spatial maps of CMRO2 were preset and thus known. The corresponding pO2 maps were found by solving the Poisson equation, which relates CMRO2 and pO2. MATLAB code for using the method is provided. Conclusion: Through this model-based study, we propose a method for estimating CMRO2 with high spatial resolution based on measurements of pO2 in cerebral cortex.


Introduction
The level of consumption of oxygen by metabolic processes, that is, the cerebral metabolic rate of oxygen (CMRO 2 ), is an important indicator of brain function and pathology. Further, knowledge about the magnitude of the CMRO 2 is required for a proper interpretation of the blood oxygenation level-dependent (BOLD) signal measured in functional MRI studies. 1 The ability to measure CMRO 2 with high spatial and temporal resolution in cortex is thus crucial. Traditionally, the CMRO 2 has been estimated from several different types of measurements combined with mathematical modeling of the underlying physiological processes. 1 Given the numerous assumptions and experimental limitations typically involved, questions have been raised about the accuracy of the estimates of the CMRO 2 provided by these complex and somewhat indirect approaches. 2 The possibility to optically measure the partial pressure of oxygen (pO 2 ) around cortical diving arterioles with two-photon resolution in vivo 3 has provided a more direct way to estimate the CMRO 2 . Previously, we (Sakadžić, Devor, and collaborators) used measured pO 2 gradients around diving arterioles in rats to estimate the average CMRO 2 in the vessel's vicinity, that is, within a radius of ∼100 μm. 4 We based our estimates on the Krogh-Erlang formula relating the pO 2 to the CMRO 2 in a cylinder section around an arteriole providing the brain tissue with oxygen. 5,6 The Krogh-Erlang formula assumes the pO 2 level to have reached a stationary state, so that the fundamental equation relating the pO 2 and the CMRO 2 in the neural tissue can be described by the Poisson equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 4 8 9 ∇ 2 PðrÞ ¼ MðrÞ; (1) where PðrÞ represents pO 2 measured at the position r, and MðrÞ is a measure that encapsulates the local CMRO 2 . The Krogh-Erlang formula gives a specific solution to the forward problem of this partial differential equation, that is, the radial profile of P, for the case where (i) the CMRO 2 [MðrÞ] is assumed to be a constant and (ii) all the oxygen provided by the center arteriole is assumed to be consumed within a radial distance R t . The problem of estimating MðrÞ based on measured pO 2 profiles PðrÞ is referred to as the inverse problem. In Ref. 4, this inverse problem was solved by fitting the Krogh-Erlang formula to pO 2 data obtained in the close vicinity of a penetrating cerebral arteriole. This approach is global in the sense that it uses all measurements within a radial distance R t to obtain an estimate for an assumed constant value of M.
In this paper, we present a different approach to estimating M based on the same kind of two-photon pO 2 measurements. The solution of inverse source problems for systems described by differential equations is important in many fields of science and technology and has consequently received substantial attention from mathematicians. 7 Equation (1) is known as the Poisson equation, and several approaches have been taken to solve the inverse Poisson problem in different science and engineering contexts. [8][9][10][11] In this study, we develop an approach to the inverse Poisson problem in the context of CMRO 2 estimation. Specifically, we solve the problem by applying the Laplace operator ∇ 2 directly to suitably smoothed pressure maps PðrÞ to obtain a measure of MðrÞ. We will refer to this approach as the diffusion-operator method for CMRO 2 estimation. Unlike the Krogh-Erlang method, the diffusion-operator method provides a spatially resolved map of CMRO 2 estimates around the arterioles and is thus not restricted to estimating an assumed constant value of M. Further, the diffusion-operator method is not restricted to situations with radially symmetric pO 2 maps as when a single arteriole provides all oxygen.
The double spatial derivatives in the Laplace operator make the diffusion-operator method inherently very sensitive to noise in the measured spatial pO 2 maps. In order to have a practical method for CMRO 2 estimation, we smooth the experimental data in two dimensions before application of the Laplace operator to reduce the effects of the noise. Smoothing introduces a bias, that is, a systematic error in the estimates, and a balance must be found where the effects of the noise are sufficiently reduced without introducing too much bias. In the present modelbased study we explore this balance by examining the accuracy of CMRO 2 estimates in situations where the ground truth, that is, spatial maps of MðrÞ are preset and thus known, and the corresponding maps of PðrÞ are found by solving the forward problem of Eq. (1), either numerically or by taking advantage of the Krogh-Erlang formula.
The manuscript is organized as follows. In Sec. 2, we describe the diffusion-operator method, the methods used to provide model-based pO 2 maps used in the method validation, and the metrics used to quantify the accuracy of the resulting estimates. In Sec. 3, we first illustrate the method and the necessary compromise between reducing noise and limiting bias when choosing the level of spatial smoothing. Next, we systematically explore the accuracy of CMRO 2 estimates for a variety of situations with different levels of noise, different grid sizes of the pO 2 measurement, and different levels of smoothing. In these systematic explorations of the efficacy of the method, the simple single-arteriole situation where the Krogh-Erlang formula gives the ground truth, is considered for simplicity. Later, we illustrate the use of the diffusionoperator method in more complicated situations where several arterioles provide the consumed oxygen, or the CMRO 2 varies with position. In Sec. 4, we discuss the diffusion-operator method and its further development and use.

Mathematical Modeling of CMRO 2 and pO 2
The blood-tissue O 2 transport is thought to be dominated by diffusion. 12 The relationship between pO 2 values denoted as Pðr; tÞ and the net rate of oxygen consumption sðr; tÞ in the tissue can then, in the general case, be described by 6,12 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 4 8 9 ∂Pðr; tÞ ∂t ¼ D∇ 2 Pðr; tÞ þ sðr; tÞ α : In Eq. (2), ∇ 2 is the Laplace operator in three spatial dimensions (3D). Further, D and α are the diffusion coefficient and solubility, respectively, of oxygen in the tissue. They are assumed to be space-invariant. If warranted, Eq. (2) can be generalized to the case where D depends on position and direction, or when α varies with position. 12 Equation (2) is only applicable outside the arterioles supplying the oxygen to the brain tissue. In the context of this equation, the oxygen supplied to the tissue is represented by a boundary condition of pO 2 imposed at the vessel wall of the arteriole. Note, however, that the effect of an oxygen supply from a bed of small capillary vessels located some distance away from the arteriole may be incorporated in the description. Such an oxygen supply will offset (or even reverse the sign of) the net rate of oxygen consumption sðr; tÞ in this region.
In this paper, we will focus on a special case of this diffusion problem where (i) the system is in a steady-state so that the term ∂Pðr; tÞ∕∂t can be neglected and (ii) there is no variation of pO 2 in the vertical z-direction, that is, the direction along the cortical axis parallel to the penetrating arteriole. These assumptions are also incorporated in the Krogh-Erlang model used to estimate the CMRO 2 in Ref. 4. In this case, the diffusion equation [Eq. (2)] simplifies to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 2 6 6 where ∇ 2 now refers to the 2D Laplace operator (which with Cartesian coordinates reads can be written more compactly as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 1 9 8 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 5 4 MðrÞ ≡ sðrÞ∕Dα: Here MðrÞ is a new position-dependent variable encapsulating the net rate of oxygen consumption in the neural tissue. In principle, Eq. (4) describes the spatial map of pO 2 for any set of oxygen sinks [metabolic consumption, sðrÞ > 0] and sources [i.e., oxygen provided by small capillaries, sðrÞ < 0]. The variable MðrÞ is then proportional to the net rate of oxygen consumption, that is, the difference between oxygen sinks and sources at position r in tissue.
By introducing a characteristic length r Ã and a characteristic oxygen consumption M Ã , we can rewrite Eq. (4) in a dimensionless form which is useful in the further analysis: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 7 1 1∇ 2P ðrÞ ¼MðrÞ: (6) In Eq. (6),r ¼ r∕r Ã ,P ¼ P∕ðM Ã r Ã2 Þ,M ¼ M∕M Ã , and∇ 2 is the Laplace operator in terms of the dimensionless position variables. In this dimensionless form, the number of model parameters is effectively reduced by one, making the further analysis simpler.

Inverse Problem of Estimating CMRO 2 from pO 2 Measurements
We estimate CMRO 2 by solving the inverse diffusion problem, that is, the problem where PðrÞ is known from experiments, and the net rate of oxygen consumption sðr; tÞ is the unknown function of interest. It follows from Eq. (2) that sðr; tÞ based on pO 2 measurements P data ðr; tÞ is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 5 5 5 s est ðr; tÞ ¼ α ∂P data ðr; tÞ ∂t þ αD∇ 2 P data ðr; tÞ: For the stationary 2D case, this reduces to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 5 0 1 s est ðrÞ ¼ αD∇ 2 P data ðrÞ; (8) or E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 4 5 8 where ∇ 2 is the Laplace operator in 2D. Equation 9 says that given a data set of oxygen partial pressure P data measured on a 2D spatial grid, M can be estimated by taking the Laplacian of P data (or in practice, a smoothed version of P data ). We refer to this approach as the diffusion-operator method for CMRO 2 estimation. If one wants estimates for s est , values of the diffusion coefficient D and the solubility α are also required.
The double spatial derivatives in the Laplace operator make the diffusion-operator method inherently very sensitive to noise in the measured spatial pO 2 profiles. Thus to reduce adverse effects of noise in the pO 2 measurements, we pursue a method that spatially smooths P data before application of the Laplace operator.

Smoothing of pO 2 data
To smooth pressure data, we performed cubic smoothing spline interpolation using the csaps function in MATLAB's Curve Fitting Toolbox. The function minimizes the square deviation between the estimated and measured 2D data (so-called L2 norm) while penalizing large double-spatial derivatives in the smoothed data. Other smoothing procedures could have been pursued instead, but a key motivation for this particular choice was the public availability of the tool.
In terms of dimensionless quantities, the csaps function takes a given data setP data ðx;ŷÞ and generates a smoothing splineP smooth ðx;ŷÞ that minimizes E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 1 8 0 Here n and m are the number of entries ofx andŷ, respectively, and q is a smoothing parameter between 0 and 1. q ¼ 0 corresponds to the case with no smoothing, and increasing values of q imply increasing the amount of smoothing. Note that the csaps function takes p ¼ 1 − q as input argument, see MATLAB documentation. This MATLAB function allows for giving more weights to some data points than others in the optimization. We keep the weights identical to 1 for all data points in the present application. The csaps function allows the smoothing splineP smooth to be computed with higher resolution than the spatial resolution of the measurements. This is convenient as it allows for a higher spatial resolution in the maps of estimated M obtained from the discrete Laplace function del2. Assuming that the measurements are taken in a rectangular grid of points, 13 we here refer to the grid spacing between the pressure data points asd data , and the grid spacing of the estimated pressure pointsP smooth asd est . In the smoothing function,d est is set by inserting position vectors for the estimation pointsx est andŷ est with this spacing. Likewise,d data is set by inserting position vectors for the data pointsx andŷ with this spacing. ThenP smooth is estimated from the recorded pressure by the following call of csaps: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 5 8 1P smooth ¼ csapsðfŷ;xg;P; ð1 − qÞ; fŷ est ;x est gÞ: (11) In this paper, we keep a fixed small value ofd est , that is,d est ¼ 0.001, to minimize the error introduced from the discreteness of the Laplace operator used in the estimator presented in the next section. With this choice, the discreteness error is negligible far away from the arteriole and much smaller than other estimation errors close to the arteriole.

Application of Laplace operator
After the smoothing procedure, the net oxygen consumption as described byMðx;ŷÞ can be estimated directly by application of the Laplace operator: WithP smooth given on a square (or rectangular) grid with grid spacingd, we apply the discrete finite difference approximation of the Laplace operator: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 3 7 5M Here the integers i and j represent the grid point positions, that is, In the present application, the MATLAB function del2 is used to compute this discrete finite difference approximation of the Laplace operator. Note that in order to calculate the right-hand side of Eq. (13), one must multiply the output from del2 by 4. Specifically, we use the command 4 Ã del2ðP smooth ;dÞ to calculateM est ðx;ŷÞ.

Choice of smoothing parameter
The effect of the csaps smoothing function can be characterized by a smoothing lengthd q , which describes how much a spatial δ-function is smeared out in space. By numerical exploration, we found that this characteristic smoothing length depends on q andd data through the relationship E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 1 4 0d where k is a constant. This relationship was found numerically by smoothing a square single-entry matrix with one as the center element, and the rest of the elements set to zero. The resulting spatially smoothed δ-function was then plotted, for a fixed value ofd data and different values of q, as a function of the distance r to the center point, as shown in Fig. 1(a). We then defined the characteristic lengtĥ d q to be the distance from the center point, in which the function value had fallen 50% compared to the center value, see dotted lines in Fig. 1(a). Figure 1(b) shows the dependence of the estimatedd q on q (for a fixedd data of 0.005). We observe thatd q increases slowly with q, that is, when q is increased by a factor 10 4 ,d q increases only by a factor 10. Figure 1(c) shows the smoothed δ-function when instead the value of q is fixed, whiled data has different values. Again, whend q is read out from the curve and plotted as a function ofd data [ Fig. 1(d)], we see thatd q increases slowly withd data , that is, whend data is increased by a factor 10 4 ,d q increases only by a factor 10.
The detailed value of the constant k in Eq. (14) is not critical for our purpose. We set it by reading out the value ford q from the graph for the case withd data ¼ 5 × 10 −3 and q ¼ 5 × 10 −4 as shown with a blue line in Fig. 1(e). The readout value,d q ≈ 5.6 × 10 −2 , was then used to calculate k from Eq. (14). After rounding to one decimal, this gave k ¼ 1.4.
Thus givend data and a chosen value ofd q , we can find which q to use in csaps in Eq. (11) through the following equation: Fig. 1 Choice of smoothing parameter in csaps. The effect of the smoothing function csaps is characterized by a smoothing lengthd q that is related to the smoothing factor q and the spatial spacingd data through Eq. (14). We found this relationship by smoothing a 2D spatial δ-function using different values of q andd data , and plotting the result as a function of the distancer from the position of the δ-function. (a) and (c) The normalized smoothed δ-function [δ smooth ðr Þ] for different values of q (d data fixed) andd data (q fixed), respectively. The characteristic smoothing lengthd q is defined as the distance corresponding to δ smooth ¼ 0.5 (dotted lines) and is plotted as a function of q andd data in (b) and (d), respectively. In (e), we demonstrate how different sets of q andd datavalues correspond to the samed q , that is, the same smoothing effect.
This equation tells us that if, say,d data increases from 5 × 10 −3 to 1 × 10 −2 , then q must decrease from q ¼ 5 × 10 −4 to about q ¼ 2.6 × 10 −4 to keep the same smoothing effect, that is, give the same value ofd q . The dotted orange line in Fig. 1(e) illustrates that this is indeed the case.

Forward Modeling of Ground Truth pO 2 Data
To validate the CMRO 2 estimation method, we generate synthetic data of oxygen partial pressurê PðrÞ by solving Eq. (4) for chosen values ofMðrÞ and chosen geometries of vascular sources and measurements points. The synthetic data work as a "ground truth." Since we know its true value ofMðrÞ, we can use it to test our estimation method. In this study, we compute this ground truth data by means of two methods: (i) using the Krogh-Erlang model and (ii) by means of finite-element modeling.

Krogh-Erlang model
In the well-known Krogh-Erlang model, 5 a cylindrical geometry, mimicking a straight segment of a blood vessel, was used to model the metabolic consumption of oxygen provided by capillaries in muscles. In Ref. 4, the same model was used to study metabolic consumption of oxygen provided by penetrating arterioles in brain tissue. The model describes the blood vessel as a small cylinder with radius R ves supplying a tissue cylinder with radius R t with oxygen. The further assumptions are (i) uniform consumption of oxygen in the tissue, that is, constant M outside the vessel, (ii) no axial diffusion of oxygen, (iii) P ¼ P ves at R ves , and (iv) no pressure gradient at the surface of the tissue cylinder, that is, dP∕dr ¼ 0 at R t . With these four assumptions, the solution of Eq. (4) is found to be E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 3 9 6 for R t ≥ r ≥ R ves . This so-called Krogh-Erlang formula predicts the oxygen pressure P in the tissue as a function of the distance r from the vessel's center. For our application, we set PðrÞ ¼ P ves if r < R ves . Equation 16 can be written in dimensionless form as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 3 0 4P ðrÞ ¼ Here we also have introducedP ves Further, the boundary condition dP∕dr ¼ 0 forr ¼R t is assumed.

Finite-element modeling: FEniCS model
The Krogh-Erlang formula relates the oxygen consumption and the partial oxygen pressure under very specific conditions. Another option is to solve Eq. (6) numerically. This allows for the solutions for more general cases, such as a more complicated geometry with, for example, several arterioles providing oxygen, or an inhomogenous oxygen consumption. We implemented Eq. (6) in the finite-element software package FEniCS 14 and verified the implementation by comparing the result to that of the Krogh-Erlang formula. The FEniCS implementation solves the variational formulation of Eq. (6): Let V be a space of test functions fv 1 ; : : : v N g on the computational domain Ω. We aim to findP such that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 6 ; 9 3 where ∂Ω denotes the boundary of the domain, and n is a normal vector pointing out of the domain. This variational form is obtained by multiplying Eq. (6) with the test function v i and integrating over Ω, followed by integration by parts of the Laplacian term. Note that as we apply a fixed value forP by the blood vessel and no pressure gradient at the boundary of the domain, the boundary integral in Eq. (18) vanishes. The solution to Eq. (18) gives usP on an unstructured finite-element mesh. Experimental data are typically measured on a structured Cartesian grid, and to better mimic this we transfer the synthetic data generated by FEniCS to a 2D NumPy array. We do this by first defining a new Cartesian mesh using NumPy with a distanced data between each point. Then in the next step, we pick out values ofP from the FEniCS solution corresponding to these positions and save them to a 2D NumPy array. We set PðrÞ ¼ P ves if r < R ves .

Noise
We add additive Gaussian noise to the synthetic data using the normrnd function in MATLAB. For each valueP of oxygen partial pressure, whether it comes from the Krogh-Erlang equation or the FEniCS solution, we draw a random number from a Gaussian distribution with meanP and standard deviation (SD)σ P , and replaceP by this number.

Performance Measures of the Diffusion-Operator Method
In order to evaluate the performance of the diffusion-operator method, we test it on the synthetic data and calculate its bias, precision, and accuracy. As precision and accuracy measures, we use SD and root-mean-square error (RMSE). The mathematical definitions of these measures are where N is the number of synthetic samples andM est;j is the j'th estimate ofM. The RMSE combines both bias and precision as its squared value MSE is equal to the SD squared plus the bias squared: MSE ¼ SD 2 þ bias 2 . 15 3 Results

Illustration of the Diffusion-Operator Method
The principle of the diffusion-operator method for estimation of the net oxygen consumption MðrÞ from pO 2 measurements PðrÞ is illustrated in Fig. 2. In this example, we assume the spatial map of pO 2 to follow the Krogh-Erlang formula in Eq. (16), mimicking the situation where a single arteriole is the source of the oxygen, and the oxygen consumption M is constant around the arteriole. Figure 2(a) shows the pressure profile in the radial directions as described by this formula with example parameters P ves , M, and R ves chosen to be in qualitative agreement with example data from Ref. 4, that is, P ves ¼ 80 mmHg, M ¼ 0.001 mmHg μm −2 , and R ves ¼ 6 μm, and R t set to 200 μm. Figure 2(b) shows a contour plot of this synthetic pO 2 map in 2D. Here dimensionless parameters (cf., Sec. 2) are used with r Ã ¼ 141 μm and the convenient choice M Ã ¼ M so that the maximal pressure P ves corresponds toP ves ≈ 4.1 andM ¼ 1. We show the pO 2 map in a square window with side lengths of 282 μm so that the dimensionless position coordinates extends from −1 to 1 along thex andŷ axes. With this choice, the corners of the square correspond to a radial distance equal toR t , the radius of the tissue cylinder.
The problem of CMRO 2 estimation now corresponds to estimating M at the different locations inside the square window based on these recordings. Figure 2(c) shows the estimated M (in units of M Ã ) found by applying the Laplace estimator in Eq. (13) on the data in Fig. 2(b). In this example, the dimensionless distance between the grid points, in which pO 2 is "measured" is set tod ¼ 0.007, corresponding to a physical grid-point distance of about 1 μm. It is seen that some distance away from the vessel, the estimator predictsM very close to 1, that is, M ≃ M Ã , as it should.
However, close to the vessel, that is, forr ≳R ves , clearly incorrect values ofM are obtained. One obvious reason is that the discrete Laplace estimator in Eq. (13) will be inaccurate when one or more of the grid points used in the estimation is inside the vessel. Here the pressure P is not described by Eq. (16) and is instead assumed constant so that ∇ 2 P ≠ M, cf. Eq. (4). For the present example, a more important reason is that immediately outside the vessel, the pressure profile drops sharply [due to the last term in the Krogh-Erlang formula in Eq. (16)] so that the discrete Laplace estimator becomes inaccurate when the grid-point distanced is too large. The "flower-like" symmetric pattern of this estimation error in Fig. 2(c) reflects the Cartesian symmetry of the estimator in Eq. (13). This discretization error can be reduced by reducing the value ofd, i.e., using a finer grid. Figure 2(c) illustrates that if the experimental measurements were noiseless, the Laplace estimator in Eq. (13) could be used directly on the pO 2 data, at least when the grid of recordings is finely spaced. This would apply for any distribution of vessels as long as the estimatorM est in Eq. (13) is used sufficiently far away from the vessel wall. Experimental pO 2 data will always (a) contain noise, however, and Fig. 2(d) shows a map of additive Gaussian noiseP σ with zero mean and SDσ P ¼ 0.0005. Figure 2(e) shows the same synthetic data as in Fig. 2(b) where this noise has been added, indistinguishable by eye from the noise-free map in Fig. 2(b). WhenM est in Eq. (13) is applied on these synthetic data, the estimated values ofM are wildly inaccurate [ Fig. 2(f)]. Not only does the estimated values ofM have much larger magnitudes than the true value ofM ¼ 1, they also have both signs and vary strongly between neighboring grid positions (that is, between neighboring pixels in the map). These poor estimates reflect that the double-derivative operation of the Laplacian estimator corresponds to a high-pass spatial filtering that effectively amplifies the effects of the noise in the data. This noise in the estimatedM can be reduced by the use of spatial smoothing, that is, lowpass filtering, of the dataP before application ofM est . While the smoothed mapP smooth in Fig. 2(g) at first glance does not appear to be very different from the unsmoothed version in Fig. 2(e), the effect of the smoothing on the estimated M is dramatic [ (Fig. 2(h)]. With the choice of smoothing used in this example (see figure caption for details), quite accurate estimates ofM are found for a large region of the area around the central vessel [light-colored regions of Fig. 2(h)]. However, the smoothing procedure results in large estimation errors in a sizable region around the blood vessel as well as close to the edges of the square data set.
To summarize, suitable smoothing of the pO 2 data before using the Laplace estimatorM est may dramatically improve the estimation accuracy. However, the choice of smoothing is critical: too little low-pass smoothing will not remove enough of the high-frequency spatial noise; too much smoothing will smooth away spatial information in the data and thus give poor estimates of M. Next, we will investigate this dilemma in more detail. Figure 3 illustrates the dilemma when choosing the right level of low-pass smoothing of the pO 2 data P before using the Laplace estimator in Eq. (13). In the smoothing algorithm, the quantity described in Eq. (10) was minimized to penalize sharp variations in P smooth while at the same time fitting the synthetic data P data . The level of smoothing is set by the smoothing length d q (ord q in dimensionless units) which is related to the smoothing parameter used in the presently used MATLAB function csaps via Eq. (14) (see Sec. 2). This smoothing length describes how much a point (that is, a 2D spatial δ-function) will be smeared out in space. Thus the larger d q is, the more the pO 2 map will be smeared out or smoothed.

Noise Removal versus Bias
To quantify the performance of the estimator, we use the following three performance measures: bias, SD, and RMSE. The bias [Eq. (19)] measures the systematic error in the estimator M est introduced by the smoothing (and discreteness of data points) whether the data is noisy or not. It can be evaluated from noiseless data (that is, with P σ ¼ 0), and the results for different values of smoothing are shown in Figs. 3(a), (d), (g), and (j). In the case of no smoothing [d q ¼ 0, Fig. 3(a)], the only bias comes from the discreteness of the grid of data points and is only observed close to the vessel. With a small amount of smoothing [d q ¼ 0.02, Fig. 3(d)], the bias around the vessel is increased. Ford q ¼ 0.04 [ Fig. 3(g)] andd q ¼ 0.08 [ Fig. 3(j)], this tendency of increased bias with increasingd q is continued, and some bias is also observed close to the edges of the square grid. For the largest smoothing depicted in Fig. 3(j), about one-third or so of the map has a bias with a magnitude larger than 100%.
The SD [Eq. (20)] measures the precision or the error in the estimation due to the presence of noise. This measure obviously depends on the level of noise P σ . In the present example in Fig. 3, a Gaussian noise with a SD ofσ P ¼ 5 × 10 −4 is used. With r Ã ¼ 141 μm and M Ã ¼ 0.001 mmHg μm −2 as in Fig. 2 this corresponds to a physical noise level of σ P ≈ 0.01 mmHg. The SD for different amounts of smoothing is shown in Figs. 3(b), (e), (h), and (k). Three observations of note are that (i) the SD of the estimates is extremely large when no smoothing is applied (d q ¼ 0), (ii) the SD decreases with increasingd q , and (iii) unlike for the bias, the SD has similar values at the different positions.
An essential feature of the SD is that it is proportional to the SD of the noise in the pressurê σ P . Thus ifσ P was doubled to 0.001, the SDs in Figs Fig. 3(I)] offers the best compromise between bias and noise removal and gives the smallest RMSE. For this value ofd q , the RMSE is smaller than 25% for almost all positions except for a region around the blood vessel. The large RMSE close to the blood vessel even for the "best" choice ofd q in Fig. 3(I) reflects the large bias at these locations [ Fig. 3(g)].

Choice of Smoothing Length d q
As illustrated in the previous section, a key question when using the Laplace estimator in Eq. (13) is the choice of the amount of smoothing, or more specifically, the choice of the smoothing length d q . This will not only depend on the noise level, but also the spatial resolution of the data as described by the grid resolution, that is, the distance between adjacent points on the measurement grid, d data . Since the bias is independent of the noise level, and the SD is linearly proportional to the SD σ P of the noise, it is convenient to first explore the interplay between d q and d data for the bias and SD separately. In Fig. 4, we show how the bias varies with d data and d q for three choices of parameter values of each:d data ¼ 0.0035, 0.007, 0.014 (here corresponding to physical grid resolutions of approximately 0.5, 1, and 2 μm, respectively),d q ¼ 0, 0.02, and 0.04 (corresponding to physical smoothing lengths of approximately 0, 3, and 6 μm, respectively). For the case with no smoothing [Figs. 4(a), (d), and (g)], we observe that the bias increases with increasingd data . This illustrates that the error due to the discreteness of the Laplace estimator is sensitive to d data even when d est is set to a very small number (d est ¼ 0.001, cf., Sec. 2). This is not surprising because decreasing the grid resolution fromd data tod est means that we estimateP at a denser grid of points than what is directly available in the data. With smoothing added (two rightmost columns of Fig. 4), the bias increases, and the larger the value ofd q is, the larger the bias is. (Note the difference in color scales in this figure.) In Fig. 5, we correspondingly show how the SD varies with d data and d q for the same set of parameters as in Fig. 4 for a fixed level of noise in the data,σ P ¼ 5 × 10 −4 . Here the most important feature is that the SD is strongly reduced with increased smoothing, that is, increasing d q (from left to right). For the smoothed cases (two rightmost columns), we also observe that SD increases with increasing d data (i.e., making the grid of measurements more sparse). Figure 6 shows the RMSE, computed from Eq. (22), for the example bias and SD shown in Figs. 4 and 5, respectively. For the smoothed cases (two right columns), we observe that the RMSE always increases with thed data . Thus, with the noise level fixed, it is (unsurprisingly) always advantageous to have a dense measurement grid. For the noise level in this example, we see that the choiced q ¼ 0.02 (second column) gives a good estimate ford data ¼ 0.0035, that is, low RMSE, for large parts of the map. Ford data ¼ 0.007 and especiallyd data ¼ 0.014 the SD is not sufficiently reduced, and the RMSE is overall high. For the case with a larger smoothing (d q ¼ 0.04, third column) the SD is much reduced for all values ofd data . However, the region (a) There was no noise added to the pressure data so that a single estimate ofM est is sufficient, that is N ¼ 1 in Eq. (19). Parameter values: with large bias around the vessel is increased, and the spatial region in which RMSE values are small is shrunken.
Note that the SD results in Fig. 5 and the RMSE results in Fig. 6 only pertain to the particular noise level used in the example, that is,σ P ¼ 5 × 10 −4 . However, the SD is proportional to the noise level, so a doubling ofσ P would simply double the SD from what is shown in Fig. 5. RMSE results analogous to Fig. 6 for other noise levels can thus be obtained by appropriate scaling of SD in Eq. (22).

Estimation of CMRO 2 for Other Example Situations
In the examples above, we have applied the diffusion-operator method to the situation with (i) a constant value of M and (ii) a single vessel providing the oxygen so that the pO 2 map is described by the Krogh-Erlang formula in Eq. (16). For these examples, an alternative approach could be to estimate M by fitting the Krogh-Erlang formula directly to measured data. 4 In other situations where, for example, M varies with position or several nearby vessels provide the oxygen so that the circular symmetry assumed in the Krogh-Erlang formula does not hold, this approach would not be applicable. In contrast, the current diffusion-operator method does not assume a constant M and can be applied to cases where multiple arterioles deliver oxygen.

Spatially varying CMRO 2
To illustrate the applicability of the Laplace estimator to the situation with varying M, we consider in Fig. 7 a hypothetical case where a single vessel provides the oxygen, but where the parameter M varies with distance from the vessel. Specifically, the value of M is assumed to be smaller far away from the vessel. This can be due to genuine differences in CMRO 2 . Alternatively, this can mimic the situation where a distant bed of capillaries acts as an oxygen source unaccounted for in the model and leading to an apparent decrease in CMRO 2 . Here the solution of the Poisson equation in Eq. (4) must be found numerically, and in Figs. 7(a) and (b), we illustrate the pO 2 maps found using the FEniCS numerical solver (see Sec. 2). Figure 7(a) shows a 1-D representation of this pO 2 profile in the radial direction for the case without any added noise. Figure 7(b) correspondingly shows a 2D colormap of the same synthetic data when noise has been added. The dotted lines in Fig. 7(a) mark the distance from the vessel (jrj ¼ 0.7) where the value ofM changes. With the characteristic length r Ã used throughout this paper, this corresponds to a physical distance of ∼100 μm, which is a typical size of the region around diving arterioles void of capillaries in the rat cortex. 4 We see in Fig. 7(a) that beyond this distance, there is almost no decay in the pO 2 compared to that within the capillary-free region. When using the Laplace estimator on the noise-free data, we obtain excellent estimates of M, that is,M est ≈ 2 within the capillary-free region andM est ≈ 0.5 outside this region [ Fig. 7(c)]. We only observe sizable errors in the immediate vicinity of the vessel, the errors stemming from the discreteness of the synthetic pO 2 data used in the estimation (d data ¼ 0.007). Further, when using the Laplace estimator on a smoothed version of the data in Fig. 7(b), we still obtain good estimates ofM some distance away from the vessel. This is in agreement with the low values for the RMSE found for suitable smoothing of noisy data for the case with constantM in Fig. 6.

Several vessels providing oxygen
An example of a situation where multiple nearby vessels serve as oxygen sources is shown in Fig. 8. Again, no analytical solution for the pO 2 map is available, and the Poisson equation is instead computed by means of FEniCS. As observed in Fig. 8(a), the circular symmetry of the pO 2 map seen in the earlier examples is broken around the vessels, but the Laplace estimator is still able to accurately estimateM except in locations close to the vessels [ Fig. 8(b)].

Estimation of Spatially-Averaged
The SD of M est;av is then expected to be a factor ffiffiffiffi N p reduced compared to the SD for the spatially resolved estimates M est ðrÞ.
The bias is not reduced by such an averaging procedure, however. To reduce the effects of smoothing-induced bias, one possible procedure is to take the average of M only for positions outside a circular region around the oxygen-delivering vessel. As illustrated in Fig. 9(a), this can reduce the bias in the M est;av substantially. Larger values of the smoothing lengthd q give larger regions of large bias around the vessel (Fig. 4). Thus larger areas around the vessel, parameterized by the diameterd cut , should be removed from the averaging sum in Eq. (23) to keep the bias small. This removal of area from the averaging sum implies a smaller value for N in Eq. (23) (a) (b) Fig. 8 Estimation of M with several vessels providing oxygen. Example of diffusion-operator estimation for a situation where three vessels release oxygen into the tissue. The synthetic pO 2 map was calculated using the FEniCS numerical solver (see Sec. 2). Here P ves is set to 80, 70, and 50 mmHg for the vessel on the left, lower right, and upper right, respectively, whereas R ves is set to 6 μm for all vessels. Noise is added to the synthetic data in (a) (σ P ¼ 0.0005), andd q ¼ 0.04 is used in the smoothing to provide the estimates ofM in (b). Other parameter values:d data ¼ 0.007, and thus a larger value of SD of M est;av . Again, a compromise between the bias and the SD must be found to get the most accurate estimate. This compromise is illustrated in Figs. 9(b)-9(g). Figure 9(b) shows the spatially resolved RMSE for a case with low noise corresponding to no smoothing applied (cf., left column of Fig. 6). Here the noise level is so low that even without smoothing, the SD of M est;av becomes <1% for all averaging areas considered, that is, all choices ofd cut [cf.,d q ¼ 0 in Fig. 9(c)]. With smoothing applied, the SD of M est;av becomes even smaller, much <0.1% [ Fig. 9(c)]. We also note that the SD is largest for the largest value ofd cut , reflecting that here the averaging area [and thus N in Eq. (23)] is the smallest. The corresponding RMSE is shown in Fig. 9(d). For this low-noise situation, there is nothing to gain by doing smoothing when estimating M est;av . The lowest RMSEs are obtained ford q ≈ 0 since smoothing reduces the accuracy of the estimates due to the bias introduced [cf., Fig. 9(a)].
The situation with a much higher noise level (σ P a factor 100 larger, that is,σ P ¼ 5 × 10 −2 ) is shown in Figs. 9(e)-(g). The spatially resolved RMSE using a smoothing factor ofd q ¼ 0.1 is seen to give large lobes with high RMSE values around the vessel [ Fig. 9(e)]. Moreover, the typical RMSE value outside the lobe region is about 120%. The SD of M est;av [ Fig. 9(f)] is seen to be on the order of 50% for the case without smoothing (d q ¼ 0), and a smaller RMSE can thus be obtained with smoothing applied [ Fig. 9(g)]. The smallest RMSE, less than ∼10%, is obtained This high-noise example illustrates how accurate estimates of M av can be obtained even when the spatially resolved estimates for M have a large uncertainty. With the parameter values used here, that is, M Ã ¼ 0.001 mmHg μm −2 and r Ã ¼ 141 μm, aσ P of 5 × 10 −2 corresponds to a physical noise level σ P of ≈1 mmHg. [Here we have used that σ P ¼σ P M Ã r Ã2 , cf., Eq. (6).] For comparison, the corresponding pO 2 at the vessel wall in this example would be P ves ¼ 80 mmHg.

Discussion
In this paper, we have introduced a new method, the diffusion-operator method, to provide spatially resolved maps of CMRO 2 estimates based on two-photon measurements of pO 2 . 3,4 The method has two key steps: (i) spatial smoothing of measured pO 2 maps followed by (ii) application of double spatial derivatives in two spatial dimensions, that is, a Laplace operator. This method is an alternative to the Krogh-Erlang method where a spatially averaged value of CMRO 2 is obtained around arterioles assuming circular symmetry. 4

Choice of Inverse-Modeling Method
The present diffusion-operator method is an approach to the inverse diffusion problem in the context of CMRO 2 estimation from high-resolution pO 2 data obtained with two-photon microscopy. The two key elements of the method are (i) the Poisson equation in Eq. (4) describing how estimates of CMRO 2 , or more precisely the variable MðrÞ in principle can be found by applying the Laplace operator on measured pO 2 maps P data ðrÞ and (ii) the use of a smoothing routine on P data ðrÞ to reduce effects of spatial noise before application of the Laplace operator. The development of the inverse-modeling method was mainly motivated by the need to have a method that is conceptually clear, easy to use, and based on publicly available software.
As the double spatial-derivative operation in the diffusion-operator approach is inherently sensitive to spatial noise, the choice of a suitable smoothing method is thus essential for obtaining accurate CMRO 2 estimates. The ideal smoothing method should reduce the effects of this spatial noise without introducing large biases in the resulting estimates. We performed smoothing using the cubic smoothing spline function csaps from MATLAB's Curve Fitting Toolbox. This method minimizes the square deviation between the estimated and measured data (so-called L2 norm) while penalizing large double-spatial derivatives in the smoothed pO 2 maps [Eq. (10)]. However, other smoothing methods could be used, for example, with norms other than L2 or using different types of splines. Also since CMRO 2 , or more precisely the variable M in Eq. (4), is proportional to double spatial derivatives, the smoothing method inherent in csaps effectively penalizes large magnitudes of M and thus introduces an unwanted bias. An alternative approach could be to penalize instead changes in the spatial derivatives of M, that is, third spatial derivatives of the pO 2 . Finally, while csaps allows for different weighting of different locations within the map, the weighting functions are restricted to be spatially separable in the x and y directions. For the present application, this limitation is not optimal as it would be preferable to exclude only a small region in and around the vessel.
While the exploration of effects of different smoothing methods on estimation accuracy is beyond the current scope, an obvious next step would be to test the accuracy of the diffusionoperator method with other smoothing methods. In particular, it would be interesting to explore to what extent other methods could reduce the size and magnitude of the lobes of large bias seen around the vessel in Fig. 4. The present MATLAB scripts, which can be found online at https:// github.com/CINPLA/CMRO2estimation, are designed to allow for an easy exchange of smoothing methods for such exploration.

Use of the Diffusion-Operator Method
The noise level and sampling distance in the experimental pO 2 data reported in Ref. 4 were too large to allow for reliable estimation of spatially resolved maps of CMRO 2 (results not shown). Further advancements in engineering of brighter and more sensitive optical pO 2 probes and further development of optical instrumentation will improve the measurement accuracy 4 and facilitate estimation of such maps. Additionally, other inverse-modeling methods may allow for more accurate spatially resolved CMRO 2 estimation based on the same set of data.
Pooling of spatially resolved estimates [as described in Eq. (23)] will always improve the accuracy, but this will be at the expense of spatial resolution. This trade-off can be investigated within the present version or future variations of the diffusion-operator method using the scripts accompanying this paper. Estimation accuracy can be studied systematically with model-based ground truth data (either based on the Krogh-Erlang model or based on FEniCS simulations) using the same grid density and noise levels as those in the experimental setting.

Generalization of the Diffusion-Operator Method
Here the diffusion-operator method has been applied to estimation of CMRO 2 for the case with 2D measurements of (assumed) steady state pO 2 data. The diffusion-operator method straightforwardly generalizes to the 3D situation and also the nonstationary case where the pO 2 varies with time. With time-resolved measurements of pO 2 across a 3D volume of brain tissue, spatiotemporally resolved estimates of CMRO 2 can be found by an analogous inverse-modeling problem based on Eq. (7). Also here, model-based validation of the estimation method can easily be pursued with synthetic data generated by finite-element modeling, for example, using FEniCS.

Disclosures
The authors declare that no competing interests exist. Early versions of some of the present work have been published as a Master's thesis 16 and in abstract form. 17