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. |
1.IntroductionThe level of consumption of oxygen by metabolic processes, that is, the cerebral metabolic rate of oxygen (), is an important indicator of brain function and pathology. Further, knowledge about the magnitude of the is required for a proper interpretation of the blood oxygenation level-dependent (BOLD) signal measured in functional MRI studies.1 The ability to measure with high spatial and temporal resolution in cortex is thus crucial. Traditionally, the 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 provided by these complex and somewhat indirect approaches.2 The possibility to optically measure the partial pressure of oxygen () around cortical diving arterioles with two-photon resolution in vivo3 has provided a more direct way to estimate the . Previously, we (Sakadžić, Devor, and collaborators) used measured gradients around diving arterioles in rats to estimate the average in the vessel’s vicinity, that is, within a radius of .4 We based our estimates on the Krogh–Erlang formula relating the to the in a cylinder section around an arteriole providing the brain tissue with oxygen.5,6 The Krogh–Erlang formula assumes the level to have reached a stationary state, so that the fundamental equation relating the and the in the neural tissue can be described by the Poisson equation: where represents measured at the position , and is a measure that encapsulates the local . The Krogh–Erlang formula gives a specific solution to the forward problem of this partial differential equation, that is, the radial profile of , for the case where (i) the [] 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 .The problem of estimating based on measured profiles is referred to as the inverse problem. In Ref. 4, this inverse problem was solved by fitting the Krogh–Erlang formula to 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 to obtain an estimate for an assumed constant value of . In this paper, we present a different approach to estimating based on the same kind of two-photon 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–11 In this study, we develop an approach to the inverse Poisson problem in the context of estimation. Specifically, we solve the problem by applying the Laplace operator directly to suitably smoothed pressure maps to obtain a measure of . We will refer to this approach as the diffusion-operator method for estimation. Unlike the Krogh–Erlang method, the diffusion-operator method provides a spatially resolved map of estimates around the arterioles and is thus not restricted to estimating an assumed constant value of . Further, the diffusion-operator method is not restricted to situations with radially symmetric 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 maps. In order to have a practical method for 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 model-based study we explore this balance by examining the accuracy of estimates in situations where the ground truth, that is, spatial maps of are preset and thus known, and the corresponding maps of 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 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 estimates for a variety of situations with different levels of noise, different grid sizes of the 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 diffusion-operator method in more complicated situations where several arterioles provide the consumed oxygen, or the varies with position. In Sec. 4, we discuss the diffusion-operator method and its further development and use. 2.Methods2.1.Mathematical Modeling of CMRO2 and pO2The blood-tissue transport is thought to be dominated by diffusion.12 The relationship between values denoted as and the net rate of oxygen consumption in the tissue can then, in the general case, be described by6,12 In Eq. (2), is the Laplace operator in three spatial dimensions (3D). Further, 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 depends on position and direction, or when varies with position.12Equation (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 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 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 can be neglected and (ii) there is no variation of in the vertical -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 in Ref. 4. In this case, the diffusion equation [Eq. (2)] simplifies to where now refers to the 2D Laplace operator (which with Cartesian coordinates reads ). Equation 3 can be written more compactly as where Here 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 for any set of oxygen sinks [metabolic consumption, ] and sources [i.e., oxygen provided by small capillaries, ]. The variable is then proportional to the net rate of oxygen consumption, that is, the difference between oxygen sinks and sources at position in tissue.By introducing a characteristic length and a characteristic oxygen consumption , we can rewrite Eq. (4) in a dimensionless form which is useful in the further analysis: In Eq. (6), , , , and 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.2.2.Inverse Problem of Estimating CMRO2 from pO2 MeasurementsWe estimate by solving the inverse diffusion problem, that is, the problem where is known from experiments, and the net rate of oxygen consumption is the unknown function of interest. It follows from Eq. (2) that based on measurements is given by For the stationary 2D case, this reduces to or where is the Laplace operator in 2D.Equation 9 says that given a data set of oxygen partial pressure measured on a 2D spatial grid, can be estimated by taking the Laplacian of (or in practice, a smoothed version of ). We refer to this approach as the diffusion-operator method for estimation. If one wants estimates for , values of the diffusion coefficient 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 profiles. Thus to reduce adverse effects of noise in the measurements, we pursue a method that spatially smooths before application of the Laplace operator. 2.2.1.Smoothing of pO2 dataTo 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 set and generates a smoothing spline that minimizes Here and are the number of entries of and , respectively, and is a smoothing parameter between 0 and 1. corresponds to the case with no smoothing, and increasing values of imply increasing the amount of smoothing. Note that the csaps function takes 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 spline 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 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 as , and the grid spacing of the estimated pressure points as . In the smoothing function, is set by inserting position vectors for the estimation points and with this spacing. Likewise, is set by inserting position vectors for the data points and with this spacing. Then is estimated from the recorded pressure by the following call of csaps: In this paper, we keep a fixed small value of , that is, , 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. 2.2.2.Application of Laplace operatorAfter the smoothing procedure, the net oxygen consumption as described by can be estimated directly by application of the Laplace operator: With given on a square (or rectangular) grid with grid spacing , we apply the discrete finite difference approximation of the Laplace operator: Here the integers and represent the grid point positions, that is, and . 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 to calculate .2.2.3.Choice of smoothing parameterThe effect of the csaps smoothing function can be characterized by a smoothing length , which describes how much a spatial -function is smeared out in space. By numerical exploration, we found that this characteristic smoothing length depends on and through the relationship where 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 of and different values of , as a function of the distance to the center point, as shown in Fig. 1(a). We then defined the characteristic length 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 on (for a fixed of 0.005). We observe that increases slowly with , that is, when is increased by a factor , increases only by a factor 10. Figure 1(c) shows the smoothed -function when instead the value of is fixed, while has different values. Again, when is read out from the curve and plotted as a function of [Fig. 1(d)], we see that increases slowly with , that is, when is increased by a factor , increases only by a factor 10. The detailed value of the constant in Eq. (14) is not critical for our purpose. We set it by reading out the value for from the graph for the case with and as shown with a blue line in Fig. 1(e). The readout value, , was then used to calculate from Eq. (14). After rounding to one decimal, this gave . Thus given and a chosen value of , we can find which to use in csaps in Eq. (11) through the following equation: This equation tells us that if, say, increases from to , then must decrease from to about to keep the same smoothing effect, that is, give the same value of . The dotted orange line in Fig. 1(e) illustrates that this is indeed the case.2.3.Forward Modeling of Ground Truth pO2 DataTo validate the estimation method, we generate synthetic data of oxygen partial pressure by solving Eq. (4) for chosen values of and chosen geometries of vascular sources and measurements points. The synthetic data work as a “ground truth.” Since we know its true value of , 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.3.1.Krogh–Erlang modelIn 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 supplying a tissue cylinder with radius with oxygen. The further assumptions are (i) uniform consumption of oxygen in the tissue, that is, constant outside the vessel, (ii) no axial diffusion of oxygen, (iii) at , and (iv) no pressure gradient at the surface of the tissue cylinder, that is, at . With these four assumptions, the solution of Eq. (4) is found to be for . This so-called Krogh–Erlang formula predicts the oxygen pressure in the tissue as a function of the distance from the vessel’s center. For our application, we set if .Equation 16 can be written in dimensionless form as Here we also have introduced , , , and . Further, the boundary condition for is assumed.2.3.2.Finite-element modeling: FEniCS modelThe 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 FEniCS14 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 be a space of test functions on the computational domain . We aim to find such that where denotes the boundary of the domain, and is a normal vector pointing out of the domain. This variational form is obtained by multiplying Eq. (6) with the test function and integrating over , followed by integration by parts of the Laplacian term. Note that as we apply a fixed value for 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 us 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 distance between each point. Then in the next step, we pick out values of from the FEniCS solution corresponding to these positions and save them to a 2D NumPy array. We set if . 2.3.3.NoiseWe add additive Gaussian noise to the synthetic data using the normrnd function in MATLAB. For each value 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 mean and standard deviation (SD) , and replace by this number. 2.4.Performance Measures of the Diffusion-Operator MethodIn 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 and where is the number of synthetic samples and is the ’th estimate of .The RMSE combines both bias and precision as its squared value MSE is equal to the SD squared plus the bias squared: .15 3.Results3.1.Illustration of the Diffusion-Operator MethodThe principle of the diffusion-operator method for estimation of the net oxygen consumption from measurements is illustrated in Fig. 2. In this example, we assume the spatial map of 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 is constant around the arteriole. Figure 2(a) shows the pressure profile in the radial directions as described by this formula with example parameters , , and chosen to be in qualitative agreement with example data from Ref. 4, that is, , , and , and set to . Figure 2(b) shows a contour plot of this synthetic map in 2D. Here dimensionless parameters (cf., Sec. 2) are used with and the convenient choice so that the maximal pressure corresponds to and . We show the map in a square window with side lengths of so that the dimensionless position coordinates extends from to 1 along the and axes. With this choice, the corners of the square correspond to a radial distance equal to , the radius of the tissue cylinder. The problem of estimation now corresponds to estimating at the different locations inside the square window based on these recordings. Figure 2(c) shows the estimated (in units of ) 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 is “measured” is set to , corresponding to a physical grid-point distance of about . It is seen that some distance away from the vessel, the estimator predicts very close to 1, that is, , as it should. However, close to the vessel, that is, for , clearly incorrect values of 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 is not described by Eq. (16) and is instead assumed constant so that , 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 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 , 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 data, at least when the grid of recordings is finely spaced. This would apply for any distribution of vessels as long as the estimator in Eq. (13) is used sufficiently far away from the vessel wall. Experimental data will always contain noise, however, and Fig. 2(d) shows a map of additive Gaussian noise with zero mean and SD . 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 in Eq. (13) is applied on these synthetic data, the estimated values of are wildly inaccurate [Fig. 2(f)]. Not only does the estimated values of have much larger magnitudes than the true value of , 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 can be reduced by the use of spatial smoothing, that is, low-pass filtering, of the data before application of . While the smoothed map 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 is dramatic [(Fig. 2(h)]. With the choice of smoothing used in this example (see figure caption for details), quite accurate estimates of 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 data before using the Laplace estimator 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 . Next, we will investigate this dilemma in more detail. 3.2.Noise Removal versus BiasFigure 3 illustrates the dilemma when choosing the right level of low-pass smoothing of the data 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 while at the same time fitting the synthetic data . The level of smoothing is set by the smoothing length (or 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 is, the more the map will be smeared out or smoothed. 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 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 ), and the results for different values of smoothing are shown in Figs. 3(a), (d), (g), and (j). In the case of no smoothing [, 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 [, Fig. 3(d)], the bias around the vessel is increased. For [Fig. 3(g)] and [Fig. 3(j)], this tendency of increased bias with increasing 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 . In the present example in Fig. 3, a Gaussian noise with a SD of is used. With and as in Fig. 2 this corresponds to a physical noise level of . 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 (), (ii) the SD decreases with increasing , 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 . Thus if 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 is measured by the RMSE [Eq. (21)], which incorporates both the bias and precision (SD) through the relation This measure describes the total statistical uncertainty of the estimates when is applied on individual data sets. The bias increases with increasing [Figs. 3(a), (d), (g), and (j)], whereas the SD instead decreases with increasing [Figs. 3(b), (e), (h), and (k)]. One would thus expect a suitable intermediate value of to give the smallest RMSE. For the example in Fig. 3, we indeed see that for the values of considered, the intermediate value [Fig. 3(I)] offers the best compromise between bias and noise removal and gives the smallest RMSE. For this value of , 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 of in Fig. 3(I) reflects the large bias at these locations [Fig. 3(g)]. 3.3.Choice of Smoothing LengthAs 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 . 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, . Since the bias is independent of the noise level, and the SD is linearly proportional to the SD of the noise, it is convenient to first explore the interplay between and for the bias and SD separately. In Fig. 4, we show how the bias varies with and for three choices of parameter values of each: , 0.007, 0.014 (here corresponding to physical grid resolutions of approximately 0.5, 1, and , respectively), , 0.02, and 0.04 (corresponding to physical smoothing lengths of approximately 0, 3, and , respectively). For the case with no smoothing [Figs. 4(a), (d), and (g)], we observe that the bias increases with increasing . This illustrates that the error due to the discreteness of the Laplace estimator is sensitive to even when is set to a very small number (, cf., Sec. 2). This is not surprising because decreasing the grid resolution from to means that we estimate 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 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 and for the same set of parameters as in Fig. 4 for a fixed level of noise in the data, . Here the most important feature is that the SD is strongly reduced with increased smoothing, that is, increasing (from left to right). For the smoothed cases (two rightmost columns), we also observe that SD increases with increasing (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 the . 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 (second column) gives a good estimate for , that is, low RMSE, for large parts of the map. For and especially the SD is not sufficiently reduced, and the RMSE is overall high. For the case with a larger smoothing (, third column) the SD is much reduced for all values of . However, the region 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, . However, the SD is proportional to the noise level, so a doubling of 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 CMRO2 for Other Example SituationsIn the examples above, we have applied the diffusion-operator method to the situation with (i) a constant value of and (ii) a single vessel providing the oxygen so that the map is described by the Krogh–Erlang formula in Eq. (16). For these examples, an alternative approach could be to estimate by fitting the Krogh–Erlang formula directly to measured data.4 In other situations where, for example, 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 and can be applied to cases where multiple arterioles deliver oxygen. 3.4.1.Spatially varying CMRO2To illustrate the applicability of the Laplace estimator to the situation with varying , we consider in Fig. 7 a hypothetical case where a single vessel provides the oxygen, but where the parameter varies with distance from the vessel. Specifically, the value of is assumed to be smaller far away from the vessel. This can be due to genuine differences in . 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 . Here the solution of the Poisson equation in Eq. (4) must be found numerically, and in Figs. 7(a) and (b), we illustrate the maps found using the FEniCS numerical solver (see Sec. 2). Figure 7(a) shows a 1-D representation of this 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 () where the value of changes. With the characteristic length used throughout this paper, this corresponds to a physical distance of , 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 compared to that within the capillary-free region. When using the Laplace estimator on the noise-free data, we obtain excellent estimates of , that is, within the capillary-free region and 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 data used in the estimation (). Further, when using the Laplace estimator on a smoothed version of the data in Fig. 7(b), we still obtain good estimates of 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 in Fig. 6. 3.4.2.Several vessels providing oxygenAn example of a situation where multiple nearby vessels serve as oxygen sources is shown in Fig. 8. Again, no analytical solution for the 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 map seen in the earlier examples is broken around the vessels, but the Laplace estimator is still able to accurately estimate except in locations close to the vessels [Fig. 8(b)]. 3.5.Estimation of Spatially-AveragedSo far, we have used the Laplace estimator to estimate spatial maps of . The Laplace estimator can give accurate estimates as long as the noise level is not too large, but the estimates of 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 data are too noisy to give reliable spatially resolved maps of estimated , one can still obtain estimates of spatially averaged values of (as when estimating based on fitting the Krogh–Erlang model in Eq. (16) to experimental data4). The obvious procedure for estimating such average values is to take the average over spatially resolved values of , that is The SD of is then expected to be a factor reduced compared to the SD for the spatially resolved estimates .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 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 substantially. Larger values of the smoothing length give larger regions of large bias around the vessel (Fig. 4). Thus larger areas around the vessel, parameterized by the diameter , 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 in Eq. (23) and thus a larger value of SD of . 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 becomes for all averaging areas considered, that is, all choices of [cf., in Fig. 9(c)]. With smoothing applied, the SD of becomes even smaller, much [Fig. 9(c)]. We also note that the SD is largest for the largest value of , reflecting that here the averaging area [and thus 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 . The lowest RMSEs are obtained for 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 ( a factor 100 larger, that is, ) is shown in Figs. 9(e)–(g). The spatially resolved RMSE using a smoothing factor of 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 [Fig. 9(f)] is seen to be on the order of 50% for the case without smoothing (), and a smaller RMSE can thus be obtained with smoothing applied [Fig. 9(g)]. The smallest RMSE, less than , is obtained for and . This high-noise example illustrates how accurate estimates of can be obtained even when the spatially resolved estimates for have a large uncertainty. With the parameter values used here, that is, and , a of corresponds to a physical noise level of . [Here we have used that , cf., Eq. (6).] For comparison, the corresponding at the vessel wall in this example would be . 4.DiscussionIn this paper, we have introduced a new method, the diffusion-operator method, to provide spatially resolved maps of estimates based on two-photon measurements of .3,4 The method has two key steps: (i) spatial smoothing of measured 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 is obtained around arterioles assuming circular symmetry.4 4.1.Choice of Inverse-Modeling MethodThe present diffusion-operator method is an approach to the inverse diffusion problem in the context of estimation from high-resolution 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 , or more precisely the variable in principle can be found by applying the Laplace operator on measured maps and (ii) the use of a smoothing routine on 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 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 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 , or more precisely the variable in Eq. (4), is proportional to double spatial derivatives, the smoothing method inherent in csaps effectively penalizes large magnitudes of and thus introduces an unwanted bias. An alternative approach could be to penalize instead changes in the spatial derivatives of , that is, third spatial derivatives of the . Finally, while csaps allows for different weighting of different locations within the map, the weighting functions are restricted to be spatially separable in the and 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 MethodThe noise level and sampling distance in the experimental data reported in Ref. 4 were too large to allow for reliable estimation of spatially resolved maps of (results not shown). Further advancements in engineering of brighter and more sensitive optical probes and further development of optical instrumentation will improve the measurement accuracy4 and facilitate estimation of such maps. Additionally, other inverse-modeling methods may allow for more accurate spatially resolved 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 MethodHere the diffusion-operator method has been applied to estimation of for the case with 2D measurements of (assumed) steady state data. The diffusion-operator method straightforwardly generalizes to the 3D situation and also the nonstationary case where the varies with time. With time-resolved measurements of across a 3D volume of brain tissue, spatiotemporally resolved estimates of 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. DisclosuresThe authors declare that no competing interests exist. Early versions of some of the present work have been published as a Master’s thesis16 and in abstract form.17 AcknowledgmentsWe thank Payam Saisan for useful input in the initial phase of the project and gratefully acknowledge support from the Research Council of Norway via the BIOTEK2021 Digital Life project “DigiBrain,” Grant No. 248828, and the BRAIN Initiative, Grant No. R01MH111359. ReferencesR. Buxton,
“Interpreting oxygenation-based neuroimaging signals: the importance and the challenge of understanding brain oxygen metabolism,”
Front. Neuroenerg., 2 8
(2010). https://doi.org/10.3389/fnene.2010.00008 Google Scholar
A. Devor et al.,
“Neuronal basis of non-invasive functional imaging: from microscopic neurovascular dynamics to bold FMRI,”
Neural Metabolism in Vivo, 433
–500 Springer, Boston, Massachusetts
(2012). Google Scholar
S. Sakadžić et al.,
“Two-photon high-resolution measurement of partial pressure of oxygen in cerebral vasculature and tissue,”
Nat. Methods, 7
(9), 755
(2010). https://doi.org/10.1038/nmeth.1490 1548-7091 Google Scholar
S. Sakadžić et al.,
“Two-photon microscopy measurement of cerebral metabolic rate of oxygen using periarteriolar oxygen concentration gradients,”
Neurophotonics, 3
(4), 045005
(2016). https://doi.org/10.1117/1.NPh.3.4.045005 Google Scholar
A. Krogh,
“The number and distribution of capillaries in muscles with calculations of the oxygen pressure head necessary for supplying the tissue,”
J. Physiol., 52
(6), 409
–415
(1919). https://doi.org/10.1113/jphysiol.1919.sp001839 JPHYA7 0022-3751 Google Scholar
D. Goldman,
“Theoretical models of microvascular oxygen transport to tissue,”
Microcirculation, 15 795
–811
(2008). https://doi.org/10.1080/10739680801938289 Google Scholar
V. Isakov, Inverse Problems for Partial Differential Equations, Springer, New York
(2017). Google Scholar
Y. Kagawa, Y. Sun and O. Matsumoto,
“Inverse solution for Poisson equations using DRM boundary element models—identification of space charge distribution,”
Inverse Prob. Eng., 1 247
–265
(1995). https://doi.org/10.1080/174159795088027583 Google Scholar
A. Farcas et al.,
“A dual reciprocity boundary element method for the regularized numerical solution of the inverse source problem associated to the Poisson equation,”
Inverse Prob. Eng., 11
(2), 123
–139
(2003). https://doi.org/10.1080/1068276031000074267 Google Scholar
J. A. Kołodziej, M. Mierzwiczak and M. Ciałkowski,
“Application of the method of fundamental solutions and radial basis functions for inverse heat source problem in case of steady-state,”
Int. Commun. Heat Mass Transfer, 37
(2), 121
–124
(2010). https://doi.org/10.1016/j.icheatmasstransfer.2009.09.015 Google Scholar
A. Frąckowiak, J. Wolfersdorf and M. Ciałkowski,
“Solution of the inverse heat conduction problem described by the Poisson equation for a cooled gas-turbine blade,”
Int. J. Heat Mass Transf., 54 1236
–1243
(2011). https://doi.org/10.1016/j.ijheatmasstransfer.2010.11.001 Google Scholar
C. Nicholson,
“Diffusion and related transport mechanisms in brain tissue,”
Rep. Prog. Phys., 64
(7), 815
(2001). https://doi.org/10.1088/0034-4885/64/7/202 RPPHAG 0034-4885 Google Scholar
A. Devor et al.,
“Overshoot of is required to maintain baseline tissue oxygenation at locations distal to blood vessels,”
J. Neurosci., 31
(38), 13676
–13681
(2011). https://doi.org/10.1523/JNEUROSCI.1968-11.2011 JNRSDS 0270-6474 Google Scholar
A. Logg, K.-A. Mardal and G. Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Springer Science & Business Media, Berlin
(2012). Google Scholar
L. Wasserman, All of Statistics: A Concise Course in Statistical Inference, Springer Science & Business Media, New York
(2013). Google Scholar
M. J. Sætra,
“Estimation of metabolic oxygen consumption from optical measurements in cortex,”
University of Oslo,
(2016). Google Scholar
L. L. Rubchinsky et al.,
“26th annual computational neuroscience meeting (CNS* 2017): part 2,”
BMC Neurosci., 18
(1), 59
(2017). https://doi.org/10.1186/s12868-017-0371-2 1471-2202 Google Scholar
|