Nonstationary Model of Oxygen Transport in Brain Tissue.

The paper addresses the mathematical study of a nonstationary continuum model describing oxygen propagation in cerebral substance. The model allows to estimate the rate of oxygen saturation and stabilization of oxygen concentration in relatively large parts of cerebral tissue. A theoretical and numerical analysis of the model is performed. The unique solvability of the underlying initial-boundary value problem for a system of coupled nonlinear parabolic equations is proved. In the numerical experiment, the tissue oxygen saturation after hypoxia is analyzed for the case when a sufficient amount of oxygen begins to flow into the capillary network. A fast stabilization of the tissue oxygen concentration is demonstrated. The reliability of the results of the numerical simulation is discussed.

1. Introduction

The main requirement of proper functioning of the brain is its sufficient supply with oxygen. Scarcity of oxygen caused, for example, by the decrease of blood circulation can provoke injuring and death of brain cells. Therefore, it is necessary to better recognize all factors affecting oxygen transport in the brain. In this connection, mathematical modeling can be very helpful to understand the cause of impaired oxygen delivery.

There are a great number of mathematical models of oxygen transport in the brain. Among them, the method, where the cerebral substance is regarded as a bifractional (blood and tissue) homogenized material, becomes more and more popular. Mathematically speaking, such models are represented by coupled partial differential equations describing convection, diffusion, and consumption of oxygen in blood and tissue fractions (see, e.g., [1-9]).

As a rule, the concentrations of oxygen in the blood and tissue fractions are distributed state variables of the homogenized models. However, such models should reflect the processes occurring in the prehomogenization phase such as oxygen penetration from blood plasma to tissue through capillary walls, which requires taking into account the plasma oxygen concentration. The nonlinear Hill equation describes a dependence between the oxygen concentrations in plasma and blood, so that the new state variable does not appear in the coupling between the blood and tissue compartments of the resulting homogenized model. It should also be mentioned the presence of a sink term in the tissue-fraction model equation, which describes oxygen consumption in tissue. This term is derived using the Michaelis-Menten equation. Therefore, the above outlined two-compartment models consisting of two coupled quasilinear parabolic equations are difficult for mathematical analysis and numerical implementation. This is the reason why many researchers have tried to simplify the models. Below, an outline of some investigations addressing various models of oxygen transport in the brain is presented.

Mathematical steady-state models of oxygen propagation in a tissue comprising a capillary network are presented in papers [1-4]. The method used there is based on Green's function techniques. In [1], some assumptions such as the neglect of blood plasma oxygen and the independence of metabolic rate on tissue oxygen concentration are imposed. This allows for reducing the degree of nonlinearity in the model. In [2-4], using modifications of the method proposed in [1], the authors consider the case of nonuniform consumption. Moreover, in [4], free oxygen dissolved in blood plasma is taken into account.

Nonstationary models of oxygen transport are the subject of research in many works (see, e.g., [5-11]). Using time-dependent models, it is possible to simulate transition regimes of oxygen transfer such as the saturation of tissues with oxygen.

Paper [5] presents a one-dimensional oxygen transport model containing both nonlinear differential and algebraic equations describing the oxygen concentrations in blood, plasma, and tissue. It should be stressed that the model spatially averages the concentration of oxygen in the tissue fraction.

In [6], a hybrid oxygen transport model is considered. It comprises a one-dimensional equation describing the oxygen transfer in the vessel, a one-dimensional conservation equation for oxygen flux through blood-tissue interface, and a three-dimensional diffusion equation, with a consumption term, governing oxygen propagation in tissue. Thus, there are three computational domains, each for its corresponding equation. The method proposed yields consistent results.

Paper [7] proposes a model of oxygen transportation to living cells. The consideration is done at the scale that allows for taking individual red blood cell into account. The model predicts partial oxygen pressure in capillaries and neighboring tissue areas.

In [8], oxygen transfer from blood to tissue is modeled using a two-compartment model operating with blood and tissue fractions. Numerical tests were conducted for relatively large vessel networks. The tissue oxygenation accounting for hematocrit distribution is computed.

In [9], on the base of a new model of the cerebral microvasculature, a nonstationary model of oxygen transport is considered. The vascular and tissue responses to changes in flow and metabolism are studied. Concerning the main assumptions of this approach, the spatial structure of the network is not taken into account, the diffusion within the tissue is neglected, and the metabolic rate is supposed to be independent on tissue oxygen concentration. The advantage of the proposed approach is its ability to be applied to a large enough volume of tissue.

A perspective trend in modeling oxygen transport is related to the so-called continuum models obtained through spatial homogenization of state variables. In [12], such a method is used to simulate heat processes in a tissue comprising a network of blood vessels. In the resulting model, the same spatial domain stands for blood and tissue fractions. A similar ansatz is used in [10, 11, 13, 14], where continuum models governed by coupled partial differential equations are studied. In [10, 11], numerical simulations of a nonstationary continuum model were performed in case of simple domains, and the comparison of results with those obtained for the original (nonhomogenized) model is done. In [13], a theoretical analysis of steady-state oxygen transport model is fulfilled. A priori estimates of solutions, implying the unique solvability of the problem under some conditions, are obtained. An iterative numerical procedure for finding solutions is proposed, along with the proof of its convergence. In [14], the investigation of steady-state continuum models is continued, the existence and uniqueness of solutions for the boundary-value problem are established, and numerical examples that illustrate the theoretical analysis are computed.

Despite the intensive numerical simulations of continuum models of oxygen transport, mathematical analysis of underlying partial differential equations is seldom addressed. The exception is a theoretical analysis of steady-state models considered in [13, 14]. As for nonstationary models, theoretical issues related to the existence and uniqueness of solutions for underlying initial-boundary value problems are still open. The aim of the present paper is to perform accurate mathematical analysis of the underlying nonstationary nonlinear initial-boundary value problem, establish its unique solvability, and implement an iterative solution algorithm based on Finite Element Method. It should be noted that the results of numerical experiments demonstrate a fast stabilization of the oxygen concentration due to diffusion and consumption effects.

Notice that homogenization leads to averaging of oxygen concentrations and smoothing the gradients. In particular, in contrast to the approaches considered in [1–4, 6–8], homogenization of vascular networks does not allow for observing gradients around single capillaries. Nevertheless, this approach allows to simulate other important processes such as diffusion, convection, and consumption of oxygen in relatively large parts of cerebral tissue. Moreover, as demonstrated in numerical experiments, it is possible to estimate the rate of tissue oxygen saturation and the corresponding stabilization of oxygen propagation on the base of the continuum model studied. Additionally, the approach proposed makes it possible to carry out a theoretical analysis of underlying initial-boundary value problems using their weak formulations. The use of continuum models enables to take into account most of the important effects of oxygen transport, for example, the dependence of the metabolic rate on the tissue oxygen concentration, amount of free oxygen dissolved in blood plasma, and diffusion rate of oxygen, which, however, was not always considered in the above cited works.

2. Problem Formulation

We consider the vessel-tissue system as a two-phase flow system, including the blood phase with the volume fraction σ and the tissue phase with volume fraction 1 − σ.

Similar to other two-phase flow models, the oxygen concentrations in the blood and tissue phases are governed by the following coupled parabolic equations (cf. [5, 10, 11]): Here, φ and θ are the blood and tissue oxygen concentrations, respectively, μ describes the tissue oxygen consumption, G is the intensity of oxygen exchange between the blood and tissue fractions, κ = σ(1 − σ)−1, σ is the volumetric fraction of vessels, v is a prescribed continuous velocity vector field in the entire domain G (the averaging of the velocity field of the capillary network), and α and β are diffusivity parameters of the corresponding phases. Notice that, as a result of homogenization, the blood and tissue oxygen concentrations are defined in the same continuum domain Ω.

The Michaelis-Menten equation describes the tissue oxygen consumption rate, μ, as function of θ, the oxygen concentration in tissue, as follows: where μ0 is the maximum value of μ and θ50 is the value of θ at μ = 0.5μ0.

The transfer rate of oxygen from blood to tissue through vessel walls is given by the formula where ψ is the oxygen concentration in plasma (concentration at vessel walls). Note that ψ is expressed through φ so that ψ is not a state variable of the resulting model. Additionally, positive constants a, b, c, and r have the following sense: a defines the oxygen permeability, b characterizes the concentration of tetramer hemoglobin, and c is given by the formula c = ψ, where ψ is the concentration of oxygen in plasma at the hemoglobin level of 50%, and r is the Hill exponent (or coefficient).

We assume that the oxygen concentrations φ and θ satisfy the following conditions on the boundary Γ = ∂Ω: and the following initial conditions: Here, ∂ is the outward normal derivative at points of the domain boundary. Nonnegative functions φ = φ(x), γ = γ(x), δ = δ(x), x ∈ Γ, and the initial functions φ0 = φ0(x) and θ0 = θ0(x), x ∈ Ω, are given. The function g is defined as the inverse of f.

3. Problem Formalization and a Priori Estimates

Let Ω be a bounded Lipschitz domain with the boundary Γ = ∂Ω consisting of finite number of smooth parts. Set Q = Ω × (0, T) and Σ = Γ × (0, T). Denote by L, 1 ≤ p ≤ ∞, the space of p-integrable (essentially bounded if p = ∞) functions. Let H be the Sobolev space W2. The space L(0, T; X) (respectively, C([0, T]; X)) consists of s-integrable on (0, T) (respectively, continuous on [0, T]) functions assuming values in a Banach space X.

Suppose that the following conditions hold for the model data: where γ0 and δ0 are constants.

Denote H = L2(Ω), V = H1(Ω), and V′ the dual of V. Identify H with its dual H′ to get the Gelfand triple V ⊂ H = H′ ⊂ V′. Let ‖·‖, ‖·‖, and ‖·‖∗ denote the norms in H, V, and V′, respectively. Notice that (f, v) is the value of functional f ∈ V′ on an element v ∈ V. If f ∈ H, then (f, v) is the inner product in H.

Introduce the inner product in V by the relation

Define the following space: where y′ = dy/dt. It is well known that W ⊂ C([0, T]; H) is the continuous embedding.

Remembering the problem formulation, introduce strictly increasing odd functions μ : ℝ⟶ℝ and f : ℝ⟶ℝ defined by the formulas

Let g : ℝ⟶ℝ denote the inverse of f. Note that

Define operators A1,2 : V⟶V′ and functionals f1,2 ∈ L2(0, T; V′) using the following relations: where u, w, v ∈ V are arbitrary functions. Notice that the bilinear forms (A1y, z) and (A2y, z) define inner products in V, and the following inequalities hold: where positive constants k1, k2 do not depend on y ∈ V.

Therefore, the problem (1)–(5) can be rewritten as a Cauchy problem in an operator form.

Definition 1 .

A pair {φ, θ} ∈ W × W is a weak solution of the problem (1)–(5) if

Remark 1 .

The above definition can be rewritten in an operator form as follows. A pair u = {φ, θ} is a weak solution iff it solves the equation u′ + Lu = f, where f = {f1, f2}, and L is the operator defined by the left-hand sides of (13) and (14). Unfortunately, as it is shown in [14], the operator L is not monotone for typical values of problem parameters, which prevents from applying standard methods of analysis (cf. [15]).

4. The Existence and Uniqueness of Solution

Theorem 1 .

Let conditions (6) hold. Then, the problem (1)–(5) is unique solvable on any finite time interval [0, T], 0 < T < ∞.

Proof

Define Galerkin approximations φ and θ of solutions of the problem (1)–(5) and derive a priori estimates which are necessary to prove the solvability. To do that, introduce a basis w1, w2, ⋯ of V such that these functions are orthonormal in H. Let satisfy the relations where φ0 and θ0 are H-orthonormal projections of the functions φ0 and θ0 on the subspace V.

Assuming v = φ and w = θ in (18) and (19), adding these equations, taking into account properties (10) and (12) as well as the nonnegativity of the products g(φ)φ and μ(θ)θ, yield the following inequality: The terms of the last inequality are being estimated using the relation uv ≤ εu2/2 + v2/2ε holding for all ε > 0. Taking ε = k1/2, ε = k2, and ε = 1 yields Along with the integration over t, this yields the following estimate: Here,

Applying the Gronwall inequality implies the claim

Obtain now an estimate guaranteeing the compactness of the sequences φ, θ in L2(0, T; H). For this end, set v = φ(t) − φ(s) in (18), integrate the result over t ∈ (s, s + h), and integrate that over s ∈ (0, T − h), where h > 0 is sufficiently small. Finally, where

Notice that, accounting for the continuity of embedding V ⊂ H, the following estimate holds: Here and below, C is a positive, independent on m constant. The terms in (27), that depend on t, can be estimated by changing the integration order in (25). Along with the boundedness of φ and θ in L2(0, T; V), this yields the following equicontinuity estimate:

Similarly, the equicontinuity of the sequence θ can be shown:

The estimates (24)–(29) allow us to pass, up to subsequences, to the limit. There are functions φ and θ such that

Convergences listed in (30) allow us to pass to the limit in (18) and (19), as m⟶∞, to prove that the limiting functions φ and θ satisfy equations (13) and (14) in the sense of distributions on (0, T) and the initial conditions hold. Note that the passage to the limit in terms containing µ(θ) and g(φ) can be easily done due to the following inequalities: which follows from the estimates of the derivatives of the functions g and μ. Note also that estimate (24) implies the inclusions which means that the time derivatives φ′ and θ′ belong to the space L2(0, T; V′) and satisfy equations (13) and (14) almost everywhere on (0, T). Thus, the pair {φ, θ} ∈ W × W is a solution of (1)–(5) in the weak sense.

Show now the uniqueness of weak solutions. Let {φ1, θ1} and {φ2, θ2} be two solutions of the problem (1)–(5), and φ = φ1 − φ2, θ = θ1 − θ2. Then, the following equalities hold:

Setting v = φ and w = θ, omitting the nonnegative terms, and taking into account the properties (10) and (12) yield where C3 = ((1/4k1)‖v‖2 + a/2);

Adding the last inequalities yields the estimate

This estimate, along with the Gronwall inequality, implies that φ = θ = 0, which means the uniqueness of solutions.

5. Numerical Simulation

Numerical example involves a 2D square domain with the area of 3.24 mm2. It contains 64 holes corresponding to 32 inlets and 32 outlets that are interpreted as arteriolar and venular ends of the capillary network. The density of inlets and outlets is chosen in accordance with cerebral physiological characteristics reported in [16].

The following parameter values were used: σ = 0.03 (see [10]), α = 2.2 · 10−3mm2/s, β = 2.4 · 10−3mm2/s (see [7]), a = 39 s−1 (see [5, 17]), b = 9.2 mM (see [5]),  ψ = 3.6 · 10−2 mM (see [5]), µ0 = 0.08 mM/s (see [18]), θ50 = 5 · 10−5 mM (see [18]), and r = 2.73 (see [5]).

To specify the boundary conditions, we set φ = 9.2 mM, θ = 0.16 mM at the inlets and φ = 8.2 mM, θ = 0.076 mM at the outlets and at the edges of square. Moreover, we set γ = 1000α, δ = 1000β. The initial distributions of the concentrations are the following: φ0 = 4 mM and θ0 = 0.01 mM, which simulate the effect of tissue hypoxia.

The velocity field v is computed in advance using the Stokes equation. To solve the Stokes equation, we set the following boundary conditions at the edges of the square: v = (0,0.6) mm/s. Moreover, following [19], velocities of 3.4 mm/s and 1.7 mm/s are set at the ends of arterioles and venules (boundaries of holes shown in Figure 1), respectively. Note that we use the Stokes equation to obtain an example of velocity field satisfying the specified boundary conditions in the considered perforated domain. Nevertheless (see Figure 1), in most of the domain, the velocity norm computed lies in the range of acceptable values (from 0.3 to 1.7 mm/s), which is necessary for normal functioning of brain cells (see [20]).

Figure 1

The absolute velocity (mm/s).

To approximate the initial-boundary value problem (1)–(5), the difference approximation of the time derivative and the first-order Finite Element spatial approximation were used. To resolve the nonlinearities at each time instant, the simple iterative procedure was applied. The function g (the inverse of f) was interpolated using cubic splines. It was observed that 20 steps of a simple iterative procedure provide a good accuracy in each time step. The Finite Element package FreeFEM++ (see [21]) was used to implement the solution method. To build a computational mesh, the following partition of the domain boundaries was specified: 60 segments for each edge of the square and 8 segments for each inlet and outlet. The time step length was taken to be equal to 0.25 s because of a good stability of the numerical scheme. Note that the Finite Element Method is suitable for solving the PDE models in complex geometries, for example, in the perforated domain considered here. Also, this approach is popular and developed enough for solving the diffusion-convection models. For example, for natural convection models in [22, 23], the stability and convergence analysis of the Finite Element Method is performed and its efficiency is demonstrated.

The oxygen concentration in tissue at different time instants is shown in Figures 2–5. The dynamics of oxygen distribution is basically defined by the convection of oxygen in the blood fraction and its diffusion and consumption in tissue. Notice that a rapid stabilization (within 6-7 seconds) of oxygen distribution in tissue is observed. Nevertheless, more fast stabilization (within 3-4 seconds) occurs in the blood fraction.

Figure 2

Brain tissue oxygen saturation after hypoxia: 1 second after oxygen begins to flow into the capillary network (mM).

Figure 3

Brain tissue oxygen saturation after hypoxia: 2 seconds after oxygen begins to flow into the capillary network (mM).

Figure 4

Brain tissue oxygen saturation after hypoxia: 3 seconds after oxygen begins to flow into the capillary network (mM).

Figure 5

Brain tissue oxygen saturation after hypoxia: 7 seconds after oxygen begins to flow into the capillary network (mM).

6. Discussion

In the numerical example considered, the initial oxygen concentration in tissue corresponds to the effect of hypoxia, which is damaging for brain cells. Nevertheless, a quick supply of a sufficient amount of oxygen from the ends of arterioles into the homogenized capillary domain rapidly improves the situation so that the oxygen concentration is being stabilized at a safe level. This is in agreement with the estimation of time needed for oxygen to propagate in the capillary network and surrounding tissue. Indeed, assuming the number of inlets (ends of a1-arterioles) and outlets (ends of v1-venules) in an adult brain equals N = 76.8 · 106, and the weight of the brain equals m = 1200 g (see [16]), and setting the density of the brain equals ρ = 1.04 g/cm3, we obtain a mean distance between inlets and outlets estimated by the value of . Moreover, accounting for the speed of blood flow in the capillary network (see [20]), the time for which the blood flow passes from inlets to outlets is estimated by 2 seconds. Also, note that the maximal distance between brain cells and capillaries is in average of 0.02 mm (see [24]). With the diffusion coefficient of 0.0024 mm2/s (see [7]), oxygen molecules travel this distance in 0.1 second. Thus, in about 2.1 seconds after entering the capillary network, oxygen molecules begin to arrive the furthest parts of the brain. This estimation is in agreement with the results of our numerical experiment, which points out to the consistency of the continuum oxygen transport model considered.

A faster stabilization of oxygen concentration in the blood fraction compared with that in tissue is explained by the Hill equation (see [25]) determining the level of oxygen in plasma. The plasma oxygen level specifies the amount of oxygen penetrating from capillaries into tissue. According to the Hill equation, a small change in plasma oxygen level may lead to a relatively large change in tissue oxygen concentration if the oxygen concentration in blood is sufficiently large.

Thus, the considered continuum model of oxygen transport can be used to study the rate of oxygenation of brain tissue in dependence on the initial level of hypoxia as well as on blood oxygen concentrations at the inlets of capillary network (the ends of arterioles).