2.1.

Mathematical Modeling of CMRO₂ and pO₂

The blood-tissue ${\mathrm{O}}_2$ transport is thought to be dominated by diffusion.¹² The relationship between ${\mathrm{pO}}_2$ values denoted as $P(\mathbf{r},t)$ and the net rate of oxygen consumption $s(\mathbf{r},t)$ in the tissue can then, in the general case, be described by^6,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 ${\mathrm{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(\mathbf{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 $\partial P(\mathbf{r},t)/\partial t$ can be neglected and (ii) there is no variation of ${\mathrm{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 ${\mathrm{CMRO}}_2$ in Ref. 4. In this case, the diffusion equation [Eq. (2)] simplifies to

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:

2.2.

Inverse Problem of Estimating CMRO₂ from pO₂ Measurements

We estimate ${\mathrm{CMRO}}_2$ by solving the inverse diffusion problem, that is, the problem where $P(\mathbf{r})$ is known from experiments, and the net rate of oxygen consumption $s(\mathbf{r},t)$ is the unknown function of interest. It follows from Eq. (2) that $s(\mathbf{r},t)$ based on ${\mathrm{pO}}_2$ measurements $P_{\text{data}}(\mathbf{r},t)$ is given by

Equation 9 says that given a data set of oxygen partial pressure $P_{\text{data}}$ measured on a 2D spatial grid, $M$ can be estimated by taking the Laplacian of $P_{\text{data}}$ (or in practice, a smoothed version of $P_{\text{data}}$). We refer to this approach as the diffusion-operator method for ${\mathrm{CMRO}}_2$ estimation. If one wants estimates for $s_{\mathrm{est}}$, values of the diffusion coefficient $D$ and the solubility $\alpha$ 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 ${\mathrm{pO}}_2$ profiles. Thus to reduce adverse effects of noise in the ${\mathrm{pO}}_2$ measurements, we pursue a method that spatially smooths $P_{\text{data}}$ before application of the Laplace operator.

2.2.3.

Choice of smoothing parameter

The effect of the csaps smoothing function can be characterized by a smoothing length ${\widehat{d}}_q$, which describes how much a spatial $\delta$-function is smeared out in space. By numerical exploration, we found that this characteristic smoothing length depends on $q$ and ${\widehat{d}}_{\text{data}}$ through the relationship

Eq. (14)

$${\widehat{d}}_q=k{(q{\widehat{d}}_{\text{data}})}^{1/4},$$

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 $\delta$-function was then plotted, for a fixed value of ${\widehat{d}}_{\text{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 length ${\widehat{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 estimated ${\widehat{d}}_q$ on $q$ (for a fixed ${\widehat{d}}_{\text{data}}$ of 0.005). We observe that ${\widehat{d}}_q$ increases slowly with $q$, that is, when $q$ is increased by a factor ${10}^4$, ${\widehat{d}}_q$ increases only by a factor 10. Figure 1(c) shows the smoothed $\delta$-function when instead the value of $q$ is fixed, while ${\widehat{d}}_{\text{data}}$ has different values. Again, when ${\widehat{d}}_q$ is read out from the curve and plotted as a function of ${\widehat{d}}_{\text{data}}$ [Fig. 1(d)], we see that ${\widehat{d}}_q$ increases slowly with ${\widehat{d}}_{\text{data}}$, that is, when ${\widehat{d}}_{\text{data}}$ is increased by a factor ${10}^4$, ${\widehat{d}}_q$ increases only by a factor 10.

Fig. 1

Choice of smoothing parameter in csaps. The effect of the smoothing function csaps is characterized by a smoothing length ${\widehat{d}}_q$ that is related to the smoothing factor $q$ and the spatial spacing ${\widehat{d}}_{\text{data}}$ through Eq. (14). We found this relationship by smoothing a 2D spatial $\delta$-function using different values of $q$ and ${\widehat{d}}_{\text{data}}$, and plotting the result as a function of the distance $\widehat{r}$ from the position of the $\delta$-function. (a) and (c) The normalized smoothed $\delta$-function [${\delta }_{\text{smooth}}(\widehat{r})$] for different values of $q$ (${\widehat{d}}_{\text{data}}$ fixed) and ${\widehat{d}}_{\text{data}}$ ($q$ fixed), respectively. The characteristic smoothing length ${\widehat{d}}_q$ is defined as the distance corresponding to ${\delta }_{\text{smooth}}=0.5$ (dotted lines) and is plotted as a function of $q$ and ${\widehat{d}}_{\text{data}}$ in (b) and (d), respectively. In (e), we demonstrate how different sets of $q$ and ${\widehat{d}}_{\text{data}}$-values correspond to the same ${\widehat{d}}_q$, that is, the same smoothing effect.

The detailed value of the constant $k$ in Eq. (14) is not critical for our purpose. We set it by reading out the value for ${\widehat{d}}_q$ from the graph for the case with ${\widehat{d}}_{\text{data}}=5\times {10}^{-3}$ and $q=5\times {10}^{-4}$ as shown with a blue line in Fig. 1(e). The readout value, ${\widehat{d}}_q\approx 5.6\times {10}^{-2}$, was then used to calculate $k$ from Eq. (14). After rounding to one decimal, this gave $k=1.4$.

Thus given ${\widehat{d}}_{\text{data}}$ and a chosen value of ${\widehat{d}}_q$, we can find which $q$ to use in csaps in Eq. (11) through the following equation:

Eq. (15)

$$q={\left(\frac{{\widehat{d}}_q}{1.4}\right)}^4\frac{1}{{\widehat{d}}_{\text{data}}}.$$

This equation tells us that if, say, ${\widehat{d}}_{\text{data}}$ increases from $5\times {10}^{-3}$ to $1\times {10}^{-2}$, then $q$ must decrease from $q=5\times {10}^{-4}$ to about $q=2.6\times {10}^{-4}$ to keep the same smoothing effect, that is, give the same value of ${\widehat{d}}_q$. The dotted orange line in Fig. 1(e) illustrates that this is indeed the case.

2.3.

Forward Modeling of Ground Truth pO₂ Data

To validate the ${\mathrm{CMRO}}_2$ estimation method, we generate synthetic data of oxygen partial pressure $\widehat{P}(\widehat{\mathbf{r}})$ by solving Eq. (4) for chosen values of $\widehat{M}(\widehat{\mathbf{r}})$ and chosen geometries of vascular sources and measurements points. The synthetic data work as a “ground truth.” Since we know its true value of $\widehat{M}(\widehat{\mathbf{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.

2.4.

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

The RMSE combines both bias and precision as its squared value MSE is equal to the SD squared plus the bias squared: $\mathrm{MSE}={\mathrm{SD}}^2+{\text{bias}}^2$.¹⁵

3.1.

Illustration of the Diffusion-Operator Method

The principle of the diffusion-operator method for estimation of the net oxygen consumption $M(\mathbf{r})$ from ${\mathrm{pO}}_2$ measurements $P(\mathbf{r})$ is illustrated in Fig. 2. In this example, we assume the spatial map of ${\mathrm{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.

Fig. 2

Illustration of diffusion-operator estimation method. (a) An example of a synthetic ${\mathrm{pO}}_2$ profile calculated using the Krogh–Erlang formula in Eq. (17) without noise. (b) The corresponding 2D representation of this ${\mathrm{pO}}_2$ data set with use of dimensionless parameters. (d) A map of additive Gaussian noise ${\widehat{P}}_{\sigma }$ and (e) the corresponding ${\mathrm{pO}}_2$ map where this noise has been added. (g) A data set where smoothing has been applied. (c), (f), and (h) Estimated $\widehat{M}\mathrm{s}$ calculated from the ${\mathrm{pO}}_2$ data in (b), (e), and (g), respectively. Parameter values: all panels: $P_{\mathrm{ves}}=80\mathrm{mmHg}$, $R_{\mathrm{ves}}=6\mu \mathrm{m}$, $R_{\mathrm{t}}=200\mu \mathrm{m}$, $M=0.001\mathrm{mm}\mathrm{Hg}\mu {\mathrm{m}}^{-2}$. (a) $d_{\text{data}}=1\mu \mathrm{m}$; (b)–(h) ${\widehat{d}}_{\text{data}}=0.007$; $r^{*}=141\mu \mathrm{m}$, $M^{*}=M$; (d)–(h) ${\widehat{\sigma }}_P=5\times {10}^{-4}$; and (g), (h) ${\widehat{d}}_{\mathrm{est}}=0.001$, ${\widehat{d}}_q=0.04$.

Figure 2(a) shows the pressure profile in the radial directions as described by this formula with example parameters $P_{\mathrm{ves}}$, $M$, and $R_{\mathrm{ves}}$ chosen to be in qualitative agreement with example data from Ref. 4, that is, $P_{\mathrm{ves}}=80\mathrm{mmHg}$, $M=0.001\mathrm{mmHg}\mu {\mathrm{m}}^{-2}$, and $R_{\mathrm{ves}}=6\mu \mathrm{m}$, and $R_t$ set to $200\mu \mathrm{m}$. Figure 2(b) shows a contour plot of this synthetic ${\mathrm{pO}}_2$ map in 2D. Here dimensionless parameters (cf., Sec. 2) are used with $r^{*}=141\mu \mathrm{m}$ and the convenient choice $M^{*}=M$ so that the maximal pressure $P_{\mathrm{ves}}$ corresponds to ${\widehat{P}}_{\mathrm{ves}}\approx 4.1$ and $\widehat{M}=1$. We show the ${\mathrm{pO}}_2$ map in a square window with side lengths of $282\mu \mathrm{m}$ so that the dimensionless position coordinates extends from $-1$ to 1 along the $\widehat{x}$ and $\widehat{y}$ axes. With this choice, the corners of the square correspond to a radial distance equal to ${\widehat{R}}_{\mathrm{t}}$, the radius of the tissue cylinder.

The problem of ${\mathrm{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 ${\mathrm{pO}}_2$ is “measured” is set to $\widehat{d}=0.007$, corresponding to a physical grid-point distance of about $1\mu \mathrm{m}$. It is seen that some distance away from the vessel, the estimator predicts $\widehat{M}$ very close to 1, that is, $M\simeq M^{*}$, as it should.

However, close to the vessel, that is, for $\widehat{r}\gtrsim {\widehat{R}}_{\mathrm{ves}}$, clearly incorrect values of $\widehat{M}$ 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 ${\nabla }^{2}P\ne 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 distance $\widehat{d}$ 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 of $\widehat{d}$, 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 ${\mathrm{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 estimator ${\widehat{M}}_{\mathrm{est}}$ in Eq. (13) is used sufficiently far away from the vessel wall. Experimental ${\mathrm{pO}}_2$ data will always contain noise, however, and Fig. 2(d) shows a map of additive Gaussian noise ${\widehat{P}}_{\sigma }$ with zero mean and SD ${\widehat{\sigma }}_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). When ${\widehat{M}}_{\mathrm{est}}$ in Eq. (13) is applied on these synthetic data, the estimated values of $\widehat{M}$ are wildly inaccurate [Fig. 2(f)]. Not only does the estimated values of $\widehat{M}$ have much larger magnitudes than the true value of $\widehat{M}=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 estimated $\widehat{M}$ can be reduced by the use of spatial smoothing, that is, low-pass filtering, of the data $\widehat{P}$ before application of ${\widehat{M}}_{\mathrm{est}}$. While the smoothed map ${\widehat{P}}_{\text{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 of $\widehat{M}$ 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 ${\mathrm{pO}}_2$ data before using the Laplace estimator ${\widehat{M}}_{\mathrm{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.

3.2.

Noise Removal versus Bias

Figure 3 illustrates the dilemma when choosing the right level of low-pass smoothing of the ${\mathrm{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_{\text{smooth}}$ while at the same time fitting the synthetic data $P_{\text{data}}$. The level of smoothing is set by the smoothing length $d_q$ (or ${\widehat{d}}_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 $\delta$-function) will be smeared out in space. Thus the larger $d_q$ is, the more the ${\mathrm{pO}}_2$ map will be smeared out or smoothed.

Fig. 3

Illustration of noise removal versus bias. (a), (d), (g), and (j) Bias of $M_{\mathrm{est}}$ for different values of $d_q$. (b), (e), (h), and (k) SD of $M_{\mathrm{est}}$ for different values of $d_q$. (c), (f), (i), and (l) RMSE of $M_{\mathrm{est}}$ for different values of $d_q$. Bias is computed from Eq. (19) for the case without noise ${\widehat{\sigma }}_P=0$ so that a single estimate of ${\widehat{M}}_{\mathrm{est}}$ is sufficient, that is $N=1$ in Eq. (19). SD is computed from Eq. (20) with ${10}^4$ estimates of ${\widehat{M}}_{\mathrm{est}}$, that is, $N={10}^4$. In the computation of SD and RMSE, ${\widehat{\sigma }}_P=5\times {10}^{-4}$. All performance measures are given as the percentage of the ground truth value $\widehat{M}=1$. Note also that the MATLAB routine csaps is used also for the case without smoothing (${\widehat{d}}_q=0$) with $q=0$ inserted in Eq. (11). Other parameter values: ${\widehat{d}}_{\text{data}}=0.007$, $P_{\mathrm{ves}}=80\mathrm{mmHg}$, $R_{\mathrm{ves}}=6\mu \mathrm{m}$, $R_{\mathrm{t}}=200\mu \mathrm{m}$, $M=0.001\mathrm{mmHg}\mu {\mathrm{m}}^{-2}$, $r^{*}=141\mu \mathrm{m}$, and $M^{*}=M$.

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_{\mathrm{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_{\sigma }=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 [${\widehat{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 [${\widehat{d}}_q=0.02$, Fig. 3(d)], the bias around the vessel is increased. For ${\widehat{d}}_q=0.04$ [Fig. 3(g)] and ${\widehat{d}}_q=0.08$ [Fig. 3(j)], this tendency of increased bias with increasing ${\widehat{d}}_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_{\sigma }$. In the present example in Fig. 3, a Gaussian noise with a SD of ${\widehat{\sigma }}_P=5\times {10}^{-4}$ is used. With $r^{*}=141\mu \mathrm{m}$ and $M^{*}=0.001\mathrm{mmHg}\mu {\mathrm{m}}^{-2}$ as in Fig. 2 this corresponds to a physical noise level of ${\sigma }_P\approx 0.01\mathrm{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 (${\widehat{d}}_q=0$), (ii) the SD decreases with increasing ${\widehat{d}}_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 pressure ${\widehat{\sigma }}_P$. Thus if ${\widehat{\sigma }}_P$ was doubled to 0.001, the SDs in Figs. 3(b), (e), (h), and (k) would be doubled as well.

The accuracy of the estimator $M_{\mathrm{est}}$ is measured by the RMSE [Eq. (21)], which incorporates both the bias and precision (SD) through the relation

The large RMSE close to the blood vessel even for the “best” choice of ${\widehat{d}}_q$ in Fig. 3(I) reflects the large bias at these locations [Fig. 3(g)].

3.3.

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_{\text{data}}$. Since the bias is independent of the noise level, and the SD is linearly proportional to the SD ${\sigma }_P$ of the noise, it is convenient to first explore the interplay between $d_q$ and $d_{\text{data}}$ for the bias and SD separately.

In Fig. 4, we show how the bias varies with $d_{\text{data}}$ and $d_q$ for three choices of parameter values of each: ${\widehat{d}}_{\text{data}}=0.0035$, 0.007, 0.014 (here corresponding to physical grid resolutions of approximately 0.5, 1, and $2\mu \mathrm{m}$, respectively), ${\widehat{d}}_q=0$, 0.02, and 0.04 (corresponding to physical smoothing lengths of approximately 0, 3, and $6\mu \mathrm{m}$, respectively). For the case with no smoothing [Figs. 4(a), (d), and (g)], we observe that the bias increases with increasing ${\widehat{d}}_{\text{data}}$. This illustrates that the error due to the discreteness of the Laplace estimator is sensitive to $d_{\text{data}}$ even when $d_{\mathrm{est}}$ is set to a very small number (${\widehat{d}}_{\mathrm{est}}=0.001$, cf., Sec. 2). This is not surprising because decreasing the grid resolution from ${\widehat{d}}_{\text{data}}$ to ${\widehat{d}}_{\mathrm{est}}$ means that we estimate $\widehat{P}$ 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 of ${\widehat{d}}_q$ is, the larger the bias is. (Note the difference in color scales in this figure.)

Fig. 4

Bias for different smoothing. (a)–(i) Bias of $M_{\mathrm{est}}$ for different values of $d_q$ (increasing from left to right) and ${\widehat{d}}_{\mathrm{est}}=0.001$ (increasing from top to bottom) computed from Eq. (19) and given as the percentage of the ground truth value $\widehat{M}=1$. There was no noise added to the pressure data so that a single estimate of ${\widehat{M}}_{\mathrm{est}}$ is sufficient, that is $N=1$ in Eq. (19). Parameter values: $P_{\mathrm{ves}}=80\mathrm{mmHg}$, $R_{\mathrm{ves}}=6\mu \mathrm{m}$, $R_{\mathrm{t}}=200\mu \mathrm{m}$, $M=0.001\mathrm{mmHg}\mu {\mathrm{m}}^{-2}$, $r^{*}=141\mu \mathrm{m}$, and $M^{*}=M$.

In Fig. 5, we correspondingly show how the SD varies with $d_{\text{data}}$ and $d_q$ for the same set of parameters as in Fig. 4 for a fixed level of noise in the data, ${\widehat{\sigma }}_P=5\times {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_{\text{data}}$ (i.e., making the grid of measurements more sparse).

Fig. 5

SD for different smoothing—fixed noise level. (a)–(i) SD of $M_{\mathrm{est}}$ for different values of $d_q$ (increasing from left to right) and $d_{\text{data}}$ (increasing from top to bottom) computed from Eq. (20) with $N={10}^4$. Values are given as the percentage of the ground truth value $\widehat{M}=1$. (Note that the grid-like pattern visible in some of the panels is a numerical artifact resulting from the application of the MATLAB routine csaps.) Parameter values: ${\widehat{\sigma }}_P=5\times {10}^{-4}$, $P_{\mathrm{ves}}=80\mathrm{mmHg}$, $R_{\mathrm{ves}}=6\mu \mathrm{m}$, $R_{\mathrm{t}}=200\mu \mathrm{m}$, $M=0.001\mathrm{mmHg}\mu {\mathrm{m}}^{-2}$, $r^{*}=141\mu \mathrm{m}$, and $M^{*}=M$.

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 the ${\widehat{d}}_{\text{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 choice ${\widehat{d}}_q=0.02$ (second column) gives a good estimate for ${\widehat{d}}_{\text{data}}=0.0035$, that is, low RMSE, for large parts of the map. For ${\widehat{d}}_{\text{data}}=0.007$ and especially ${\widehat{d}}_{\text{data}}=0.014$ the SD is not sufficiently reduced, and the RMSE is overall high. For the case with a larger smoothing (${\widehat{d}}_q=0.04$, third column) the SD is much reduced for all values of ${\widehat{d}}_{\text{data}}$. However, the region with large bias around the vessel is increased, and the spatial region in which RMSE values are small is shrunken.

Fig. 6

Root-mean-square error for different smoothing—fixed noise level. (a)–(i) RMSE computed from Eq. (21) for the bias and SDs shown in Figs. 4 and 5, respectively. Values are given as the percentage of the ground truth value $\widehat{M}=1$.

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, ${\widehat{\sigma }}_P=5\times {10}^{-4}$. However, the SD is proportional to the noise level, so a doubling of ${\widehat{\sigma }}_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).

3.4.

Estimation of CMRO₂ 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 ${\mathrm{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.⁴ 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.

3.4.1.

Spatially varying CMRO₂

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 ${\mathrm{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 ${\mathrm{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 ${\mathrm{pO}}_2$ maps found using the FEniCS numerical solver (see Sec. 2). Figure 7(a) shows a 1-D representation of this ${\mathrm{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 ($|\widehat{r}|=0.7$) where the value of $\widehat{M}$ changes. With the characteristic length $r^{*}$ used throughout this paper, this corresponds to a physical distance of $\sim 100\mu \mathrm{m}$, which is a typical size of the region around diving arterioles void of capillaries in the rat cortex.⁴ We see in Fig. 7(a) that beyond this distance, there is almost no decay in the ${\mathrm{pO}}_2$ compared to that within the capillary-free region.

Fig. 7

Estimation of spatially varying $M$. Diffusion-operator estimation of $\widehat{M}$ for the case with a single oxygen-releasing vessel in the center with a larger $\widehat{M}$ close to the vessel ($\widehat{r}<0.7$ and $\widehat{M}=2$) than far away ($\widehat{r}>0.7$ and $\widehat{M}=0.5$). The synthetic ${\mathrm{pO}}_2$ maps were calculated using the FEniCS numerical solver (see Sec. 2). (a) 1-D illustration of ${\mathrm{pO}}_2$ map for the case without noise (${\widehat{\sigma }}_P=0$). The dotted lines mark the boundary between different levels of $\widehat{M}$. (b) 2D illustration of the case with noise added (${\widehat{\sigma }}_P=0.0005$). (c) Estimated $\widehat{M}$ from the noise-less data without use of smoothing. (d) Estimated $\widehat{M}$ from the data in (b) (where noise is present) with use of smoothing (${\widehat{d}}_q=0.04$). Other parameter values: ${\widehat{d}}_{\text{data}}=0.007$, $P_{\mathrm{ves}}=80\mathrm{mm}\mathrm{Hg}$, $R_{\mathrm{ves}}=6\mu \mathrm{m}$, $r^{*}=141\mu \mathrm{m}$, and $M^{*}={10}^{-3}$.

When using the Laplace estimator on the noise-free data, we obtain excellent estimates of $M$, that is, ${\widehat{M}}_{\mathrm{est}}\approx 2$ within the capillary-free region and ${\widehat{M}}_{\mathrm{est}}\approx 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 ${\mathrm{pO}}_2$ data used in the estimation (${\widehat{d}}_{\text{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 of $\widehat{M}$ 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 constant $\widehat{M}$ in Fig. 6.

3.4.2.

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 ${\mathrm{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 ${\mathrm{pO}}_2$ map seen in the earlier examples is broken around the vessels, but the Laplace estimator is still able to accurately estimate $\widehat{M}$ except in locations close to the vessels [Fig. 8(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 ${\mathrm{pO}}_2$ map was calculated using the FEniCS numerical solver (see Sec. 2). Here $P_{\mathrm{ves}}$ is set to 80, 70, and 50 mmHg for the vessel on the left, lower right, and upper right, respectively, whereas $R_{\mathrm{ves}}$ is set to $6\mu \mathrm{m}$ for all vessels. Noise is added to the synthetic data in (a) (${\widehat{\sigma }}_P=0.0005$), and ${\widehat{d}}_q=0.04$ is used in the smoothing to provide the estimates of $\widehat{M}$ in (b). Other parameter values: ${\widehat{d}}_{\text{data}}=0.007$, $M=0.001\mathrm{mm}\mathrm{Hg}\mu {\mathrm{m}}^{-2}$, $r^{*}=141\mu \mathrm{m}$, and $M^{*}=M$.

3.5.

Estimation of Spatially-Averaged $M$

So far, we have used the Laplace estimator to estimate spatial maps of $M$. The Laplace estimator can give accurate estimates as long as the noise level is not too large, but the estimates of $M$ in the immediate vicinity of the oxygen-releasing blood vessels are typically inaccurate due to the bias introduced by the smoothing procedure.

In situations where the ${\mathrm{pO}}_2$ data are too noisy to give reliable spatially resolved maps of estimated $M$, one can still obtain estimates of spatially averaged values of $M$ (as when estimating ${\mathrm{CMRO}}_2$ based on fitting the Krogh–Erlang model in Eq. (16) to experimental data⁴). The obvious procedure for estimating such average values $M_{\mathrm{est},\mathrm{av}}$ is to take the average over spatially resolved values of $M_{\mathrm{est}}$, that is

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_{\mathrm{est},\mathrm{av}}$ substantially. Larger values of the smoothing length ${\widehat{d}}_q$ give larger regions of large bias around the vessel (Fig. 4). Thus larger areas around the vessel, parameterized by the diameter ${\widehat{d}}_{\text{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) and thus a larger value of SD of $M_{\mathrm{est},\mathrm{av}}$. Again, a compromise between the bias and the SD must be found to get the most accurate estimate.

Fig. 9

Estimation of spatially averaged $M$. Illustration of accuracy of the estimation of spatially averaged $M$ for different values of the diameter ${\widehat{d}}_{\text{cut}}$ of the circular disc removed from the average in Eq. (23). $N=1000$ has been used in the estimation of the SD [Eq. (20)]. Other parameter values: All panels: ${\widehat{d}}_{\text{data}}=0.035$, $P_{\mathrm{ves}}=80\mathrm{mmHg}$, $R_{\mathrm{ves}}=6\mu \mathrm{m}$, $R_{\mathrm{t}}=200\mu \mathrm{m}$, $M=0.001\mathrm{mmHg}\mu {\mathrm{m}}^{-2}$, $r^{*}=141\mu \mathrm{m}$, and $M^{*}=M$. (a) ${\widehat{\sigma }}_P=0$; (b)–(d) ${\widehat{\sigma }}_P=5\times {10}^{-4}$; and (e)–(g) ${\widehat{\sigma }}_P=5\times {10}^{-2}$, ${\widehat{d}}_q=0.1$. Note that for figure clarity, only the circles corresponding to ${\widehat{d}}_{\text{cut}}=0.1$, 0.2, and 0.5 are shown in (b) and (e).

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_{\mathrm{est},\mathrm{av}}$ becomes $<1\%$ for all averaging areas considered, that is, all choices of ${\widehat{d}}_{\text{cut}}$ [cf., ${\widehat{d}}_q=0$ in Fig. 9(c)]. With smoothing applied, the SD of $M_{\mathrm{est},\mathrm{av}}$ becomes even smaller, much $<0.1\%$ [Fig. 9(c)]. We also note that the SD is largest for the largest value of ${\widehat{d}}_{\text{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_{\mathrm{est},\mathrm{av}}$. The lowest RMSEs are obtained for ${\widehat{d}}_q\approx 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 (${\widehat{\sigma }}_P$ a factor 100 larger, that is, ${\widehat{\sigma }}_P=5\times {10}^{-2}$) is shown in Figs. 9(e)–(g). The spatially resolved RMSE using a smoothing factor of ${\widehat{d}}_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_{\mathrm{est},\mathrm{av}}$ [Fig. 9(f)] is seen to be on the order of 50% for the case without smoothing (${\widehat{d}}_q=0$), and a smaller RMSE can thus be obtained with smoothing applied [Fig. 9(g)]. The smallest RMSE, less than $\sim 10\%$, is obtained for ${\widehat{d}}_q\approx 0.1$ and ${\widehat{d}}_{\text{cut}}=0.3$.

This high-noise example illustrates how accurate estimates of $M_{\mathrm{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\mathrm{mm}\mathrm{Hg}\mu {\mathrm{m}}^{-2}$ and $r^{*}=141\mu \mathrm{m}$, a ${\widehat{\sigma }}_P$ of $5\times {10}^{-2}$ corresponds to a physical noise level ${\sigma }_P$ of $\approx 1\mathrm{mmHg}$. [Here we have used that ${\sigma }_P={\widehat{\sigma }}_P{M}^{*}{r}^{*2}$, cf., Eq. (6).] For comparison, the corresponding ${\mathrm{pO}}_2$ at the vessel wall in this example would be $P_{\mathrm{ves}}=80\mathrm{mmHg}$.

4.1.

Choice of Inverse-Modeling Method

The present diffusion-operator method is an approach to the inverse diffusion problem in the context of ${\mathrm{CMRO}}_2$ estimation from high-resolution ${\mathrm{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 ${\mathrm{CMRO}}_2$, or more precisely the variable $M(\mathbf{r})$ in principle can be found by applying the Laplace operator on measured ${\mathrm{pO}}_2$ maps $P_{\text{data}}(\mathbf{r})$ and (ii) the use of a smoothing routine on $P_{\text{data}}(\mathbf{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 ${\mathrm{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 ${\mathrm{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 ${\mathrm{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 ${\mathrm{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 diffusion-operator 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.

4.2.

Use of the Diffusion-Operator Method

The noise level and sampling distance in the experimental ${\mathrm{pO}}_2$ data reported in Ref. 4 were too large to allow for reliable estimation of spatially resolved maps of ${\mathrm{CMRO}}_2$ (results not shown). Further advancements in engineering of brighter and more sensitive optical ${\mathrm{pO}}_2$ probes and further development of optical instrumentation will improve the measurement accuracy⁴ and facilitate estimation of such maps. Additionally, other inverse-modeling methods may allow for more accurate spatially resolved ${\mathrm{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.

4.3.

Generalization of the Diffusion-Operator Method

Here the diffusion-operator method has been applied to estimation of ${\mathrm{CMRO}}_2$ for the case with 2D measurements of (assumed) steady state ${\mathrm{pO}}_2$ data. The diffusion-operator method straightforwardly generalizes to the 3D situation and also the nonstationary case where the ${\mathrm{pO}}_2$ varies with time. With time-resolved measurements of ${\mathrm{pO}}_2$ across a 3D volume of brain tissue, spatiotemporally resolved estimates of ${\mathrm{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.