1 Introduction

Due to its highly specialized function, the brain is one of the organs with the highest basal energy demand. With essentially no substantial energy reserves, it is thus extremely vulnerable to sudden interruptions in oxygen and nutrients delivery by the blood, which can induce neuronal death within minutes with devastating consequences, e.g., for stroke victims [1]. It is also highly sensitive to chronic cerebral hypoperfusion, which can lead to progressive neurodegeneration and cognitive decline, not only in hypoperfusion dementia [2] but also, as increasingly accepted, in Alzheimer’s disease [3–6]. However, despite its critical role in the transition between health and disease, many aspects of oxygen transport and metabolism in the brain remain poorly understood.

This motivated the development of high-resolution brain imaging techniques [7]. Together with the increased sophistication of experimental protocols, which enabled the brain of living rodents to be studied in various conditions including sleep, resting and awake states, these provide an unprecedented window on microvascular dynamics (e.g. diameters, red blood cell velocities, blood and tissue oxygenation, neural activity) [7–10]. However, due to the intrinsically heterogeneous and non-local nature of network flows [11–13], the results obtained in different conditions have been difficult to interpret. As we shall see next, this contributed to casting doubt on previously accepted ideas, including the fundamental idea that both structure and function of the brain microcirculation are subservient to cerebral metabolic demand.

With regard to brain function, the physiological role of neurovascular coupling, i.e. local surges in blood flow driven by increased neuronal activity (also referred to as functional hyperamia), has been questioned. On the one hand, even the baseline level of blood flow is indeed globally sufficient to supply oxygen to neurons with elevated levels of activity [12]. On the other hand, in the words of Drew [12], “low-flow regions are an inescapable consequence of the architecture of the cerebral vasculature” and “cannot be removed by functional hyperemia”. In fact, “increases in blood flow—whether local or global—will serve only to move the location of the low-blood-flow regions, not eliminate them [13]”.

With regard to structure, the local variations of vascular density have been believed for decades to reveal underlying differences in neuronal organization and metabolic load [14–16], as a result of cerebral angiogenesis being driven by their oxygen requirements [17–19]. Recent breakthroughs in brain-wide vascular network imaging and reconstruction in rodents, associated to scaling analyses, support this vision at the scale of the whole brain [20]. However, detailed measurements of periarteriolar oxygen profiles across cortical layers in awake mice, associated with estimates of the corresponding cerebral metabolic rate of oxygen, recently suggested that baseline oxygen consumption may decrease with cortical depth, from Layer I to Layer IV [10], in contrast to the known increase of capillary density [20, 21].

Solving these apparent contradictions requires the development of models integrating the non-local nature of microvascular blood flow [11, 13, 22], which account for the complex architecture of brain microvascular networks but simplify or neglect transport and metabolism within the tissue, with models of oxygen dynamics going beyond the geometrical oversimplifications associated to Krogh-type analytical descriptions [10, 23–28].

However, the computational cost of simulating oxygen transport and consumption in the brain parenchyma by standard numerical methods, such as finite volume or finite element methods, is prohibitive. In fact, they imply to finely mesh the extravascular tissue so as to resolve the strong oxygen concentration gradients building up in the vicinity of each vessel (e.g. [29]), not to mention the technical challenge of automatically meshing its complex three-dimensional volume. A popular alternative, specifically designed to solve oxygen transport in the microcirculation, formulates the problem using Green’s functions [30–34]. The non-local nature of this formulation allows the description of concentration gradients around microvessels while circumventing the need for meshing the intricate geometry of the extravascular space. However, it generally relies on the infinite domain form of the Green’s function, making difficult the application of boundary conditions at the limits of the tissue domain (e.g. periodic boundary conditions). Additionally, oxygen metabolism exhibits non-linear behavior [23, 35], which is challenging to describe using Green’s function and generally requires additional meshing [30, 36]. This, coupled with the non-local formulation at the core of the approach, requires the creation of large and dense matrices that are computationally costly to invert. Therefore, solving oxygen transport whether using standard methods or the Green’s function approach limits the size of the regions that can be considered and hinders the potential of such methods to be used in inverse problems, where measured spatial oxygen dynamics are used to deduce local metabolic rate constants or permeability coefficients, which requires to run the direct problem many times. In the latter case, the spatial resolution of the solver is much higher than that of the measurements, which requires averaging of the numerical results, deviating from an optimal allocation of computational resources.

Such challenges have been bypassed by introducing dual mesh techniques, where the extravascular domain is coarsely meshed independently of microvessel locations [37], or by simplifying the mesh structure, e.g. based on cartesian grids, to approximate the extravascular domain [38]. These approaches decrease the computational cost, but do not leverage recent progresses in other fields, where analytical solutions to similar problems (analogous form of equations with same underlying mathematical structure) could be used to capture the smallest features of the extravascular oxygen field (perivascular gradients). This would circumvent the need of mesh refinement around the sources. In geosciences (well or fractured reservoir modelling), for example, coupling models are often used where analytical functions help provide a relationship between the highly conductive slender structures (commonly modeled as 1D sources) and the 3D simulation domain [39–45]. In particular, in operator-splitting approaches [46–48], the scalar field (concentration, pressure, heat, etc.) is decomposed into a slowly varying contribution and a rapidly varying contribution. The former can be solved numerically over a coarse cartesian mesh, while the later can be approximated analytically, thus enabling a precise estimation of exchanges at the vessel-tissue interface as well as an a posteriori highly-resolved reconstruction of the concentration field in each mesh cell.

The goal of the present paper is to revisit the problem of oxygen transport and metabolism in the brain parenchyma to introduce a simple, scalable and accurate numerical method for its direct resolution. By simplicity, we mean the ability to use cartesian mesh cells independent of vessel locations, thus avoiding meshing the extravascular space, as well as the ability to impose various boundary conditions at the outer limits of the computational domain. By scalability, we refer to a mathematical formulation of the problem at the core of which is a low-bandwidth linear system of equations, so that the numerical resolution can be fully and efficiently parallelized. By accuracy, we mean the ability to control the numerical errors even in the case of a coarse mesh. Here, we present the associated concepts in two dimensions (Section 2), so as to increase the readability of the mathematical developments. This also permits to exploit current commercial finite element solvers, which enable to obtain reference solutions of the initial boundary value problem. This enables to carefully study how the underlying simplifications translate into numerical errors in idealized test cases that sequentially challenge these assumptions (Section 3). We then show how this model helps understanding the recent counter-intuitive experimental results on cortical oxygenation and metabolism [10, 24, 26] (Section 4). Finally, we discuss how this novel approach compares to previous work and how it will provide the groundwork for computationally affordable oxygen transport and metabolic simulations, fully coupled with intravascular transport in large microvascular networks.

2 Model and methods

We first focus on the diffusive transport of oxygen in the brain parenchyma, i.e. the brain tissue except for blood vessels, denoted Ω_σ in Fig 1A, for which we present the general three-dimensional formulation in Section 2.1. We then restrict ourselves to a 2D configuration, where vessels are reduced to a collection of circular sources, as schematized in Fig 1B. This enables to maintain the readability of the mathematical developments, introduced from Section 2.2 onwards, without significant loss of generality. We finally consider oxygen consumption in Section 2.4.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Terminology and notations for parenchyma and vessel spaces.

Panel A represents a 3D region Ω of the brain tissue, which includes the parenchyma Ω_σ and the vessel space Ω_β. The external boundary is denoted by ∂Ω, the vessel walls by ∂Ω_β, the vessels center-lines by Λ, the curvilinear coordinate system for the vessels by (s, r, θ) and the outer normal to the vessel walls by n. Panel B illustrates a 2D geometry, as used in the present paper to establish the modeling framework in Section 2, with two sources. The source walls are denoted by ∂Ω_β,1 and ∂Ω_β,2. Panel C displays one example of a tesselation of space Ω_σ into 4 sub-spaces V_k for k = {1, 2, 3, 4}. Here, only two sub-spaces contain sources, i.e., E(V₁) = E(V₄) = ⌀, E(V₂) = {2}, and E(V₃) = {1}.

https://doi.org/10.1371/journal.pcbi.1011973.g001

2.3 Assembly of a system of discrete algebraic equations

For that purpose, we set tessellation to match a cartesian grid of cell side-length h = |∂V_k,m|, where m is a direct neighbour of k (i.e ) with: (14) as defined in Fig 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Terminology and notations for the FV discretization (A) and sub-grid interpolation (B).

Panel A displays the current cell k of the cartesian mesh in green and its direct neighbours in blue. is the value of the slow term at the center of cell k. The dummy variables and represent the values of the slow term on both sides of the interface ∂V_k,e. Generally, due to the jump introduced by Eq 12. Panel B displays the dual mesh used for sub-scale interpolation in red. This dual mesh is constructed by joining the centers of the FV grid. Its cells are denoted by numbers (0, 1, 2 and 3).

https://doi.org/10.1371/journal.pcbi.1011973.g002

From this point forward, we use symbol ^~ to represent the discrete average values of a field on each FV cell k. Noteworthy, for harmonic functions such as , from Gauss’s harmonic function theorem [59], if we neglect the small volume occupied by the vasculature (≈3% [60]), this average value can be approximated at second-order by the value at the cell’s center.

Therefore, the unknowns of the system are the values of the slow term at the center of each FV cell , and the vessel-tissue flux for each source q_j. These are represented by two vectors of discrete variables and q = {q₁, q₂, q₃, …, q_S}, respectively.

2.3.1 FV discretization for the slow term.

The gradient of the slow term is approximated by the Two Point Flux Approximation (TPFA): (15) where h is the size of the FV cell face h = |∂V_k,m| and are dummy variables, to be eliminated by substitution from the final system, which represent the values of the slow term on interfaces ∂V_k,m.

Additionally, we use the classic FV formulation by integrating Eq 13a over each FV cell k (see Section A in S1 Methods for more detail). This yields: (16)

The discrete versions of boundary conditions 13e and 13f are:

From the above equations, we can express the dummy variables as follows: (18) with: (19) where J_k,m is a function of the sources q as we shall see in Section 2.3.2, and it accounts for the discontinuities of the rapid term across the interfaces of the FV. We can express Eq 16 as a function of the unknowns of the system: (20) where and are found under the vector .

Moreover, if the current mesh cell k belongs to a boundary, the boundary condition 10c is used instead of 17, yielding: (21)

Therefore, the discretized version of BVP 13 can be assembled from Eqs 20 and 21 into an algebraic system with as many equations as grid-cells: (22) where matrix contains the classic diffusion stencil, and the vector J contains the values of J given by Eq 19. Therefore, for each row k, contains one diagonal value and 4 off-diagonal values associated to its neighbours, while the vector b_∂Ω contains the entries relevant to enforce the BCs 21.

We have constructed a system of algebraic equations that enforces mass balance of the concentration field in each FV cell through Eq 20. To go further, we must specify the choice of the rapid term that will allow the the entries of J to be deduced from Eq 19. We note could be obtained numerically as in [56, 57], or approximated analytically based on the Green’s function formulation, as detailed in the next Section.

2.3.2 Potential-based localized formulation for the rapid term .

We first recall that, as written in Section 2.2, must be harmonic functions that satisfy Eq 9 for all k, ensuring to consider, a minima, the rapid contribution of all sources within V_k.

Straightforward analytical approximations for in 2D are therefore: (23) where represents any extension of V_k, i.e. any region of space containing V_k, and P_j is the single-source potential associated to source j.

According to potential theory [61–63], P_j can be written as (see Section B in S1 Methods): (24)

With this explicit definition of the potentials, the expression of the rapid term (Eq 23) only depends on the vessel-tissue exchanges (q) and on the position ||x−x_j||, so that . From this expression of the rapid term, now we have an explicit definition of J from Eq 22 as a function of the vessel-tissue exchanges q. We can thus assemble the discrete system of equations with and q as follows: (25)

Noteworthy, strictly fulfills the constraints corresponding to Eq 9 if and only if there is a single source i in , for which x_i ∈ V_k. Any additional source j in induces a perturbation : (26) of the normal flux around source i (see Section C in S1 Methods), which is not accounted for in the model. However, the integral contribution of these errors is null so the model remains conservative (Section C in S1 Methods). The impact of this perturbation will be examined in the Results Section 3.3.

Inspired by our previous work in [64], we set to correspond to a finite number n² of cells in the finite-volume mesh, with n ≥ 3 to avoid a special treatment for sources lying on the interface between two mesh cells.

By this way, n sets up the characteristic size of the region in which we account for the contribution of nearby sources to , while the contribution of sources outside is only implicitly treated through , as illustrated in Fig 3. Thus, increasing n leads to a better approximation of the concentration field (see Section 2.3.4), but at the same time increases the density of matrix in Eq 25. In the limit case where , we would obtain an element-wise non-zero , leading to a non-sparse system similar to [30, 51, 65] where the boundary integrals of the classic Green’s function formulation are estimated by . Since the goal here is to obtain a sparse linear system, () is chosen to be small in comparison to the domain of computation, but large enough to include the near source gradients.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Localization strategy illustrated for two sources with neighbourhood of size 3x3 grid-cells (i.e., n = 3).

Panel A: Cells where the rapid term accounts for source 1 and source 2 are displayed in blue and green respectively, with superposition in cells h and k, so that ∀w ∈ {a, b, c, d, e, f, m}, ∀w ∈ {g, p, q, r, t, u, v} and ∀w ∈ {h, k}; cells lying further from sources 1 and 2 are displayed in white. In these cells, ; Panel B: Flux balance for all cells highlighted in dark blue in panel A. Green arrows represent the contributions of slow terms while grey arrows those of rapid terms. The latter may exhibit jumps, e.g. at interfaces V_km and V_nm due to the localization-induced discontinuities in the rapid term. Panel C: Concentration field decomposition (Eq 7) along the x-axis crossing the center of source 2 (dashed axis in Panel A). We show in red the fine-grid reconstructed solution through Eq 30, in green the coarse-grid slow term and in grey the rapid term.

https://doi.org/10.1371/journal.pcbi.1011973.g003

The estimation of the single source potential based on the Green’s integral formulation has a natural extension to 3D. The circular sources that appear in the 2D model provide a simple formulation to the potential since the double layer potential is null (see Section B in S1 Methods). In contrast, an open cylinder provides a non-null value for the double layer integral resulting in a second potential in Eq 24 [51] which accounts for the axial variations. The rest of the developments presented in Section 2.3, including FV discretization and localization of the slow term, can be simply extrapolated to 3D.

2.3.3 Sub-grid reconstruction to estimate vessel-tissue exchanges (q).

The vessel-tissue exchanges are governed by Eq 3, which in 2D translates into: (27) for each source j ∈ E(Ω). In 3D, this equation should be coupled to an intravascular transport problem that introduces a discrete 1D description of average intravascular concentrations along vessel centerlines [22] as additional unknowns (see e.g. [45, 66]). In our 2D case, however, sources are disconnected, so that the values of 〈C_v〉_j are provided as boundary conditions. The average wall concentration is given by: (28)

However, the numerical model only provides an approximation of the concentration field at the grid-cell centers x_k: (29) and, from Eq 21, at the boundary nodes.

To estimate from Eq 28, we must reconstruct the concentration field everywhere in Ω_σ. For that purpose, we interpolate the slow term from its values at x_k and x_k,∂Ω using a classical set of linear shape functions γ_i associated to these points, as defined in Section D in S1 Methods. We also introduce a extended rapid term bridging the discontinuities across the interfaces of the FV cells , as detailed in Section D in S1 Methods. The resulting interpolation function reads: (30) where represents the set of FV grid-cell centers x_k and of boundary nodes x_k,∂Ω.

Since both and are harmonic functions in V_k, the average needed to estimate in Eq 28 can be deduced from Gauss’s harmonic function theorem, yielding so that Eq 27 becomes: (31)

We can now assemble Eq 31 into a discrete linear system of S equations with and q as vectors of unknowns: (32)

Noteworthy, the interpolation function uses the nearby sources as well as the four nearest FV unknowns of the mesh grid (see Section D in S1 Methods). Therefore, matrix is sparse, with 4 non-zero entries per line, while the density of matrix depends on the size of . Additionally, vector contains the values of intravascular concentrations (〈C_v〉_j) treated here as boundary conditions.

2.3.4 Full discrete system and error induced by localization.

The full discrete system is therefore: (33) with a total of F + S unknowns, being the total number of FV grid-cells and S the number of sources. This general form is independent of the specific choice made for the size of extensions used to define as linear combinations of potentials based on the Green’s formulation. This size, however, strongly influences the densities of matrices and . Nevertheless, matrices and always remain sparse since is the classic FV diffusion matrix with only 5 non-zero terms per line and only depends on the interpolation function , resulting in 4 non-zero elements per line.

Of course, the global error resulting from approximating BVP 13 by the above system depends on the size of . This global error induced by the localization strategy can be estimated by considering the neglected contribution of sources outside of to the concentration field ().

The error associated to FV methods is commonly given by [44]: (34) where C₀ is bounded by the norm of the second derivative of the estimated field. Using Eq 24 and considering that the minimal distance between a source in V_k and one outside is of order (n−1)h/2, we get an upper-bound of : (35)

Simplifying, we obtain: (36)

Therefore the localization error is expected to decrease with n², i.e., .

3 Results: Error estimation

In this Section, we first consider test cases involving a single source (Section 3.1) and a single dipole, i.e., the combination of a single source and a single sink (Section 3.2). Then, in Section 3.3, we turn to multiple sources and sinks. Noteworthy, we generically designate by “source” any vessel j whose concentration is greater than the local tissue concentration, i.e., for which the resulting flux q_j will be positive. In the same way, we use “sink” for any vessel j whose concentration is lower than the local tissue concentration, i.e., for which the resulting flux q_j will be negative. This enables “diffusional shunts” in the parenchyma, which have been evidenced experimentally between arterioles and venules [69], to be considered.

Thus, for all simulations we assign 〈C_v〉_j = ϕ_max to all sources and 〈C_v〉_j = 0 to all sinks, where ϕ_max represents the oxygen concentration in penetrating arterioles at the inlet of the brain cortex. We also use a diffusion coefficient D = 2 × 10⁻⁵cm² ⋅ s⁻¹ [25, 30, 70], an effective permeability for the capillaries of K_eff = 2 × 10⁻⁵cm² ⋅ s⁻¹ [29, 30] and a maximum metabolic consumption of M = 2.4 μmol ⋅ cm⁻³ ⋅ min⁻¹ which falls within physiological range (see Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Parameter values.

Radii R: see Fig 7; ϕ_max: oxygen concentration in penetrating arterioles at the inlet of the brain cortex; D: diffusion coefficient in the parenchyma; K_eff: effective diffusive permeability of the capillary walls; α: oxygen solubility in water at atmospheric pressure; M: maximum metabolic rate of oxygen; K: concentration where consumption is half of its maximum.

https://doi.org/10.1371/journal.pcbi.1011973.t001

Moreover, for all test cases considered in this Section, we purposely put ourselves in Case 1.1 above by considering relatively small domains of side L = 240μm, i.e., only 50 times larger than the source/sink radii (R = 4.8μm). In doing so, we aim at providing reasonable estimates for the upper bounds of the numerical errors.

Errors are estimated by comparison with a fine mesh finite element (FE) solution of the original BVP (Eq 1 for the linear problem or Eq 37 for the non-linear problem) without any additional modeling assumptions, in the same spirit as [29]. This reference FE solution, ϕ_ref, was obtained with COMSOL Multiphysics using a triangular mesh fine enough to accommodate the contours of the circular sources, to handle the azimuthal variations of the concentration field around the sources and to ensure convergence in the estimation of q_ref, obtained by integrating the normal derivative of ϕ_ref along the vessel wall.

We define the following metrics to compare our multiscale model with this reference solution. The local errors on the vessel-tissue exchanges for each source (q_j) are given by: (41) and the local errors on the concentration field at the center of each grid-cell are given by: (42) where is given by Eq 29. We then define the global errors as the average of the local errors: (43) and: (44) where where F is the number of discrete grid-cells in the cartesian mesh and S is the total number of sources, i.e., S = Card(E(Ω)).

We consider the error on vessel-tissue exchanges () as the main metric to assess the model’s accuracy, since proper estimation of q_j’s relies on an accurate evaluation of the microscale dynamics and provides crucial information on oxygen exchanged between blood and tissue. The error on the concentration field serves as a secondary metric, offering valuable insights into the interactions between sources.

We compare these errors with the errors resulting from a coarse-grid FV approach without multiscale coupling, in the same spirit as [37]. Such an approach solves the simplified BVP (Eq 6 for the linear problem or Eq 37 for the non-linear problem) by approximating the average concentration on the vessel wall () by the value of the concentration field in the nearest FV cell k ( for Ω_β,j ∈ V_k). As a result, the exchange term is given by (45) where k is the grid-cell containing source j. This coupling condition is not a multiscale coupling condition, as it doesn’t integrate any description of the near source concentration gradients that could compensate for the scale gap with the coarse-grid for the estimation of q_j. At its core, it assumes a well-mixed concentration within each mesh cell, i.e., it neglects the effect of concentration gradients near sources when using a coarse grid, generating significant errors in the estimation of q_j (see Figs 4–6). On the one hand, increasing mesh discretizations can solve this issue and allow to (asymptotically) recover the influence of such gradients [38], with the significant trade-off of increased computational cost. On the other hand, including a multiscale component by reconstructing analytically the local concentration near sources as done in Section 2, allows to capture the influence of the gradients whilst allowing to use a coarse-grid discretization of the tissue space. The FV solution with resolution matching that of the coarse-grid is thus useful to illustrate the interest of the multiscale coupling at the core of the developments presented in Section 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Error estimation for a single source: Impact of discretization and boundary effects.

A: schematics of the configuration under study, highlighting the detail of the boundary conditions. The domain size is L = 240μm, the source radius is R = 4.8μm, and the neighbourhood size is (30R)², i.e., n = 3 for a 5x5 grid (h/L=0.2); B: evolution of global errors as a function of grid size for the linear and non-linear problems, and for both the multiscale and the coarse-grid FV model (see legend); C: schematics of the boundary test, for which we use a mesh size h/L = 0.2, i.e., a 5x5 grid; D: evolution of global errors as a function of d, from d = 0 where the source is in contact with the no-flux boundary, to d = 1.2h where the source lies in the contiguous grid-cell. The dashed vertical line illustrates the limit of the boundary cell.

https://doi.org/10.1371/journal.pcbi.1011973.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Error estimation for a single dipole: Interplay between source/sink separation distance and neighborhood size.

A: evolution of the local errors on vessel-tissue exchanges for the source (filled symbols) and for the sink (empty symbols) as a function of distance d between source and sink; B: schematics of the smaller-neighborhood configuration (n = 3); C: schematics of the large-neighborhood configuration (n = 5); D: reconstruction of the sub-grid concentration field for the case n = 3 and d = 60 ⋅ R. The value of the concentration (ϕ) is non-dimensionalized by the value of the intravascular concentration in the source. The dashed vertical line in Panel A illustrates the transition between a situation where, for n = 3, the intersection of the source and sink neighborhoods contains both of them to a situation where the source and sink lie outside each other’s neighborhood.

https://doi.org/10.1371/journal.pcbi.1011973.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Error estimations for a realistic source distribution.

A: synthetic capillary bed generated by the method of [60] with cutting plane highlighted in blue; B: intersections of each capillary vessel in A with the cutting plane (red dots) and coarse cartesian grid (dashed lines); C: coarse-grid solution for the concentration field with metabolic consumption; D: sub-grid reconstruction of ϕ using from Section 2.3.3. E: global errors on vessel-tissue exchanges estimations for a grid size h = L/16, therefore 256 FV cells in total. The curve with κ = 0.1 is represented by the dashed black line; F: evolution of the errors in the estimation of the vessel-tissue exchanges and the concentration field with decreasing mesh cell size and with n chosen so that the size of is approximately 3L/5 = 30R.

https://doi.org/10.1371/journal.pcbi.1011973.g006

For the sake of comparison, the following conventions are used in all figure legends in this Section:

• Blue lines are used for the present multiscale method while red ones are used for the coarse-grid FV model.
• Continuous lines are used for the linear, non-reactive model (Eq 33) while discontinuous ones are used when metabolism is considered, i.e. reactive model (Eq 39).
• Square markers are used to display the global errors on the vessel-tissue exchanges while triangular ones are used to display errors on the concentration field.

4 Results: Periarteriolar oxygen concentration gradients

Now that we have shown the ability of our model to efficiently solve for the oxygen concentration field, including around vessels where gradients are the strongest, we turn to its exploitation in the context of brain metabolism. We specifically ask if variations of the radial peri-arteriolar concentration profiles that were recently measured across cortical layers in awake mice [10] could result from the layer-specific (laminar) increase of capillary density with cortical depth rather than from variations of baseline oxygen consumption.

For that purpose, we consider the typical case of a single penetrating arteriole (PA) and its surrounding tissue, as illustrated in Fig 7A. To account for the capillary-free space that encircles the PA, we include a cylindrical tissue region devoid of capillaries, with typical radius of 100 μm [10, 24]. Further away, we generate a random spatial distribution of sources with densities approximately matching the capillary density in cortical layer II (Table 2). We deduce the equivalent 2D source density (E2DSD) by using synthetic capillary networks similar to Section 3.3. We then create a randomized but statistically homogeneous distribution of sources following [60, 73]. We assign concentrations at the outer walls of the PA and capillaries, by using an asymptotically large value of K_eff and following experimental measurements in layer II [8]. For the capillaries, we assign a random distribution of normalized concentrations at capillary walls ∂Ω_β,j, drawn from a Gaussian distribution with mean ϕ_cap = 0.45ϕ_PA and standard deviation σ = 0.1ϕ_PA, approximately corresponding to experimental measurements in layer II (Table 2). We also impose periodic boundary conditions on the limit of the domain to mimic the larger cortical space. Finally, the maximum metabolic rate of oxygen is chosen to be M = 2.4μmol ⋅ cm⁻³ ⋅ min⁻¹ (Table 1), an intermediate value within the physiological range. Other transport parameters used to solve the non-linear problem (Eq 39 and 40) are deduced from reference values in the literature (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Effect of capillary density and intravascular concentration on radial periarteriolar concentration profiles.

A: sketch of simulated configuration, with a random and homogeneous capillary bed for r > R_cyl and a capillary-free region around the central PA. The cartesian grid of size 20x20 matches that of experimental sampling used in [10]. B: corresponding coarse-grid partial pressure deduced by linear transformation from the concentration field using the solubility of oxygen in brain tissue [10, 30, 37]. Capillary density and concentration correspond to layer II (Table 2) and M = 2.4μmol ⋅ cm⁻³ ⋅ min⁻¹. Note that all simulations (panels B, D, E, F) use the same value of M and n = 10; C: example of an experimentally sampled oxygen partial pressure field around a PA at 100μm under cortical surface, i.e. at the interface between layer I and II; D: estimated metabolic consumption deduced from Panel B; E: radial concentration profiles predicted in layers I to IV, each obtained by averaging the results of 30 simulations; Inset: result obtained when only variations of the capillary density are considered; F: resulting spatial average of the tissue concentration for r > R_cyl, as a function of capillary density for four values of the average capillary intravascular concentration (ϕ_cap/ϕ_PA from 0.4 to 0.55) corresponding to the four layers in Table 2; the black dots represent the resulting spatial average of the tissue concentration for r > R_cyl for the four different layers, i.e., varying both capillary density and concentration. The black dashed line represents the associated linear fit ().

https://doi.org/10.1371/journal.pcbi.1011973.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Layer-specific (laminar) variations of capillary density and average intravascular capillary PO₂.

The capillary length density (CLD) and depth of the corresponding layers are approximated from data in [20]. The equivalent two-dimensional source density (E2DSD) is deduced using synthetic capillary networks from [60] (Fig 6A). The ratio between arteriole and capillary concentration is approximated from data in [9].

https://doi.org/10.1371/journal.pcbi.1011973.t002

A typical realization of the resulting coarse-grid concentration field, converted to partial pressure (PO₂) using the solubility of oxygen in brain tissue [10, 30, 37], is displayed in Fig 7B for a 20x20 grid (h ∼ 20μm). This matches the experimental sampling used in [10] and results in a very good qualitative agreement with the field measured in a 100μm-deep plane perpendicular to a PA (see Fig 7C, data acquired following [10]). This field can be used to deduce the sub-grid concentration dynamics by interpolation (Eq 30, see e.g. Figs C and D in S1 Figures), as well as the local cerebral metabolic rate of oxygen (CMRO₂, see Eq 40 and Fig 7D). Interestingly, the cerebral metabolic rate of oxygen exhibits ∼±10% variations in the outer region (around capillaries) and up to ∼±20% in the periarteriolar region, in contrast to the common assumption of a spatial homogeneity [10, 23, 24, 74].

Moreover, the above results can be post-processed to deduce the azimuthal average of the normalized concentration around the PA, as displayed in Fig 7E as a function of the radial distance r to the PA center. In this figure, the plain orange line corresponds to the mean over 30 realizations of source distributions for layer II, while the faint orange areas shows the associated standard deviation. This radial concentration dynamics exhibits three regimes (Fig 7E): 1/ a constant value for r ≤ R_PA, i.e. within the PA, consistent with the source potential (Eq 24); 2/ a fast decrease for R_PA ≤ r ≤ ∼0.8R_cyl corresponding to the inner part of the region devoid of capillaries around the PA (see Fig 7A) and 3/ a re-increase followed by a slowly-varying region for larger values of r. The presence of a local minima, which can also be observed in the measurements (Fig 7C and Fig E in S1 Figures, dashed lines) suggests that the outer region of the capillary-free cylinder is both fed by the PA and the capillary bed.

Next, as the capillary density approximately increases linearly with depth in the cortex from layer I to layer IV [20], we increased the source density from 250 to 475 mm⁻² (Table 2). This results in an increase in size of regions with high oxygen concentration around capillaries (see panel B vs. A in Fig C in S1 Figures) and therefore 1/ in a slight decrease of the steepness of radial periarteriolar PO₂ gradients averaged over 30 realizations and 2/ in a slight increase of the partial pressure in the plateau region (r ≥ R_cyl), see inset in Fig 7E. This increase can be quantified by plotting the spatial average (see isocolor variations in Fig 7F).

Increases of the average intravascular PO₂ within the capillary bed (Table 2), which can be speculated based on depth-resolved experimental measurements of vascular oxygen within the cortex [9], result in a higher increase in size of regions with high oxygen concentration around capillaries (see panel C vs. A in Fig C in S1 Figures) and thus to higher increase of (i.e. black dashed line vs. colored dashed lines in Fig 7F).

Combined together, the increase in capillary density and intravascular PO₂ that has been reported experimentally from layer I to layer IV in the cortex of living rodents leads to an even faster decrease of the perivascular concentration gradient (Fig 7E). This results in a faster increase of from layer I to layer IV, see black dots in Fig 7F. For a constant value of M, this yields an increasing metabolic rate of oxygen from layer I to layer IV (panel D vs. A in Fig C in S1 Figures), consistent with the increased density of mitochondrial cytochrome oxidase through these layers [10, 75, 76].

Noteworthy, similar results have been obtained for different values of the maximal metabolic rate of oxygen M within the physiological range (Table 1), as illustrated in Fig D in S1 Figures. The only notable difference is that the amplitude of the dip in the radial concentration profile decreases with decreasing values of M (Fig E in S1 Figures). For the set of parameters representative of layers I and II, this leads to a monotonous decrease of the average radial concentration followed by a region with nearly constant oxygen concentration, similar to experimental measurements reported in [10, 24], as soon as M ≤ 1.2μmol ⋅ cm⁻³ ⋅ min⁻¹ (Fig E in S1 Figures). As a result, the oxygen dynamics in the close vicinity of the arteriole (r < R_cyl/2) may be highly similar in different cortical layers for different values of M (e.g. M = 2.4μmol ⋅ cm⁻³ ⋅ min⁻¹ in layer IV, see red line in Fig 7E vs. M = 1.6μmol ⋅ cm⁻³ ⋅ min⁻¹ in layer I, see blue dashed-dotted line in Fig E in S1 Figures).

Altogether, the present model suggests that laminar variations of the capillary density may be sufficient to explain the differences in periarteriolar radial oxygen profiles measured at different depths within the cortex [10], without any variation of the maximal cerebral metabolic rate of oxygen M. If laminar variations of intravascular capillary PO₂ are also considered, the predicted differences are even larger than the experimental ones. This demonstrates the interplay between metabolic consumption, capillary density and intravascular availability of oxygen in the capillary bed to determine the radial oxygen gradient in the vicinity of PAs. This makes it difficult to consider the steepness of the radial periarteriolar oxygen profile as a surrogate for the baseline oxygen consumption, with the potential to reconcile recent experimental measurements with the idea that laminar variations of capillary density could reveal underlying differences in metabolic load.

5 Discussion

In this paper, we revisited the problem of oxygen transport and metabolism in the brain parenchyma, with the goal to introduce a simple, scalable and accurate numerical method for its direct resolution. Getting inspiration from previous work on blood flow and oxygen transport in the brain [22, 30, 45, 64] and on mixed-dimensional problems in applied mathematics for geosciences [39, 56, 57] and biology [48], we applied the notion of operator-splitting, which allowed us to describe the oxygen concentration field in the parenchyma as the sum of a slow and a fast varying contributions. The slow contribution was treated using a classic finite volume approach on a coarse grid, while the fast contribution was described using Green’s functions that allowed to analytically capture the sub-grid perivascular concentration gradients. This made it possible to locally bridge the scale gap between sources and the coarse-grid with higher flexibility than [39, 41] regarding the position of the sources within the coarse grid-cells, including multiple sources within a single grid-cell, the proximity of the boundaries, and the control of the matrix sparsity thanks to the the size of the neighbourhood (). Similarly to singularity removal approaches [45, 48, 57], this also made it possible to mix the dimensionality reduction of the Green’s functions approaches [30, 32, 51] with the versatility of FV methods [56, 57, 72]. Moreover, solving for a slowly varying background concentration field, thanks to a change of variable, offers a huge advantage with respect to the Green’s functions resolution since it allows for a localization of the source potentials, thereby providing a much sparser system. In addition, this enable the use of a much coarser mesh that reduces considerably the size of the system compared to FV or FE methods, but without loss of precision thanks to the sub-grid reconstruction of the concentration field (Figs 5D and 6D). This, in turn allows the addition of non linear volume terms (metabolism) without significant loss of accuracy (Figs 4B, 4D and 6E).

To provide rigorous but still intelligible mathematical developments, we focused on a two-dimensional version of the problem (Section 2). In this way, we were able to introduce the localization scheme enabling us to control the bandwidth of the associated linear system of equations, by manipulating the size of the region over which the interactions with nearby vessels are accounted for analytically (scalability). We demonstrated the existence of an optimal size for this region, above which the errors induced by localization are smaller than those induced by deviations from local azimuthal symmetry of the concentration field around each vessel (see Fig 6E). This emphasizes the importance of comparing the results of any simplified model for oxygen transport in the brain parenchyma with a reference solution that is able to fully resolve these deviations. This is neither the case if only single vessels with Krogh-type configurations are considered for validation, as in [10, 23, 27, 37], or if the discrete version of the problem is compared to the corresponding continuous version (i.e., comparing the solution of Eq 33 to the solution of Eq 6 instead of Eq 1), as in [58, 66, 72, 77]. To our knowledge, such comparisons had never been performed before in this context. Crucially, they enabled to provide careful estimates of the numerical errors associated to the use of coarse meshes, demonstrating the unprecedented balance between reduction of problem size and minimization of errors associated to our method, compared to previous strategies in the literature (accuracy). With this regard, it is worth insisting that we designed test cases that enabled the origins of errors to be understood by purposefully choosing configurations with deviations from the model underlying assumptions (Section 3). Thus, all errors provide upper-bounds of the errors expected when considering larger, physiological-like problems. Moreover, the mathematical groundwork provided by the Green’s function framework (Section B in S1 Methods) permitted to trace back the source of errors to specific modeling assumptions, which in turn offers a rationale for choosing the model parameters, including discretization and neighborhood size. Finally, the method makes use of a cartesian mesh independent of vessel locations, thereby belonging to the class of mesh-less approaches [38]. In contrast with the widespread semi-analytical methods [30, 32, 34, 51, 65] that require computationally intensive Fourier transforms to enforce conventional boundary conditions such as Neumann and Dirichlet [30], it also shows remarkable versatility with regard to the boundary conditions that can be handled (simplicity).

Of course, the two-dimensional version of the problem we considered doesn’t enable coupling oxygen transport in the parenchyma with oxygen transport in blood vessels, because intersecting the vascular network with a plane yields disconnected vascular sources, as schematized in Fig 1. Thus, in the present work, intravascular concentrations have been treated as inputs, while in a three-dimensional version they should be treated as unknowns, with additional blocks in the final system accounting for intravascular transport, as highlighted in [45, 52, 66]. Together with previous work by our group that focused on revisiting intravascular transport [22], these will provide the foundations for an extension to three dimensions. The fact that the integrated potential arising from a circular source can be written analytically (Eq 24 and Section B in S1 Methods) is a peculiar characteristic of the 2D model. In 3D, additional errors may arise due to the approximations needed to estimate the potential of a small cylindrical element used to discretize the vascular networks, that will require careful evaluation in the spirit of [78, 79].

However, considering two-dimensional situations already offers the opportunity to investigate important physiological questions, as illustrated in Section 4. Thanks to the computational efficiency of our method, we were able to easily generate the large number of synthetic data (600 simulations per layer and 30 measurements per simulation) needed to incorporate statistical information about the capillary bed available in the literature [9, 20]. By this way, we shed new light on the interpretation of experimentally measured variations of periarteriolar oxygen profiles [10] in the context of laminar variations of capillary density and their relationships with baseline cerebral metabolic load. Due to the non-local nature of blood flow, understanding the measured variations of average intravascular concentration with depth will require a fully coupled three-dimensional analysis of oxygen exchange in the brain. In this regard, it is worth noting that the matrix assembly process only depends on the location of vessels in the considered domain and of the size of grid-cells and neighbourhood , and requires to be performed only once when the structure of the vascular network is known. Besides parametric analyses, this will pave the way for the inverse modelling of brain metabolism from three-dimensional oxygen measurements. Inverse modelling indeed requires fast and precise forward model resolutions, overcoming the geometric simplifications of the Krogh cylinder type, which are at the basis of all current work that aims at measuring the cerebral metabolic rate of oxygen [10, 24, 26].

6 Conclusion

We developed a multi-scale model describing the spatial dynamics of oxygen transport in the brain that is simple, conservative, accurate and scalable. Our strategy was to consider that the oxygen concentration in the parenchyma is the result of a balance between contributions at local (cell metabolism, delivery by neighbouring arteriole) and larger (capillary bed) spatial scales. This allowed us to split the oxygen concentration field into a slow and a fast varying terms, underlying the separation of scales between local and distant contributions. Doing so allowed us to combine a coarse-grid approach for the slow-varying term to a Green’s function approach for the fast-varying terms. This resulted in a computationally efficient model that was able to capture precisely gradients of concentration around microvessels and to describe boundary condition with flexibility (Dirichlet, Neumann, and periodic) along with the non-linear metabolic activity by the cells.

We then compared our model with reference solutions of the oxygen transport problem in scenarios of increasing complexity. We showed that our model was able to maintain small errors, even in scenario where the separation of scales was challenged or when azimuthal variations of flux were no longer negligible, demonstrating the robustness of the model.

While the present multi-scale model focused on two-dimensional problems, it has been designed to be easily extended to three-dimensional problems by adapting the expression for the source potential and by including an intravascular description of oxygen transport, both of which being available in the literature (see e.g. [22, 30, 45]).

Despite this limitation, we showed that the model was already capable of generating synthetic data reproducing the heterogeneous distribution of oxygen in the brain parenchyma. Doing so, we showed that periarteriolar gradients were the result of the balance between local cellular oxygen consumptionand supply by not only the neighbouring arteriole but also distant capillaries, thus reconciling recent measurements of periarteriolar oxygen gradients across cortical layers with the fundamental idea that variations of vascular density within the depth of the cortex may reveal underlying differences in neuronal organization.

Supporting information

Acknowledgments

We gratefully acknowledge Maxime Pigou for technical help with the development and maintenance of the code and Franca Schmid for carefully reading our manuscript.