Abstract
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 .
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, is the intensity of oxygen exchange between the blood and tissue fractions, , is the volumetric fraction of vessels, is a prescribed continuous velocity vector field in the entire domain (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 is the maximum value of and is the value of at .
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 , , , and have the following sense: defines the oxygen permeability, characterizes the concentration of tetramer hemoglobin, and is given by the formula , where is the concentration of oxygen in plasma at the hemoglobin level of 50%, and 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 , , , , and the initial functions and , , are given. The function is defined as the inverse of .
3. Problem Formalization and a Priori Estimates
Let be a bounded Lipschitz domain with the boundary consisting of finite number of smooth parts. Set and . Denote by , the space of -integrable (essentially bounded if ) functions. Let be the Sobolev space . The space (respectively, ) consists of -integrable on (respectively, continuous on ) functions assuming values in a Banach space .
Suppose that the following conditions hold for the model data: where and are constants.
Denote , , and the dual of . Identify with its dual to get the Gelfand triple . Let , , and denote the norms in , , and , respectively. Notice that is the value of functional on an element . If , then is the inner product in .
Introduce the inner product in by the relation
Define the following space: where . It is well known that is the continuous embedding.
Remembering the problem formulation, introduce strictly increasing odd functions and defined by the formulas
Let denote the inverse of . Note that
Define operators and functionals using the following relations: where are arbitrary functions. Notice that the bilinear forms and define inner products in , and the following inequalities hold: where positive constants do not depend on .
Therefore, the problem (1)–(5) can be rewritten as a Cauchy problem in an operator form.
Definition 1. A pair 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 is a weak solution iff it solves the equation , where , and is the operator defined by the left-hand sides of (13) and (14). Unfortunately, as it is shown in [14], the operator 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
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 of such that these functions are orthonormal in . Let satisfy the relations where and are -orthonormal projections of the functions and on the subspace .
Assuming and in (18) and (19), adding these equations, taking into account properties (10) and (12) as well as the nonnegativity of the products and , yield the following inequality: The terms of the last inequality are being estimated using the relation holding for all . Taking , , and yields Along with the integration over , this yields the following estimate: Here,
Applying the Gronwall inequality implies the claim
Obtain now an estimate guaranteeing the compactness of the sequences in . For this end, set in (18), integrate the result over , and integrate that over where is sufficiently small. Finally, where
Notice that, accounting for the continuity of embedding , the following estimate holds: Here and below, is a positive, independent on constant. The terms in (27), that depend on , can be estimated by changing the integration order in (25). Along with the boundedness of and in , 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 , to prove that the limiting functions and satisfy equations (13) and (14) in the sense of distributions on () and the initial conditions hold. Note that the passage to the limit in terms containing and can be easily done due to the following inequalities: which follows from the estimates of the derivatives of the functions and . Note also that estimate (24) implies the inclusions which means that the time derivatives and belong to the space and satisfy equations (13) and (14) almost everywhere on . Thus, the pair is a solution of (1)–(5) in the weak sense.
Show now the uniqueness of weak solutions. Let and be two solutions of the problem (1)–(5), and . Then, the following equalities hold:
Setting and , omitting the nonnegative terms, and taking into account the properties (10) and (12) yield where ;
Adding the last inequalities yields the estimate
This estimate, along with the Gronwall inequality, implies that 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: (see [10]), , (see [7]), (see [5, 17]), (see [5]), (see [5]), (see [18]), (see [18]), and (see [5]).
To specify the boundary conditions, we set at the inlets and at the outlets and at the edges of square. Moreover, we set , . The initial distributions of the concentrations are the following: and , which simulate the effect of tissue hypoxia.
The velocity field 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: . 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]).
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 (the inverse of ) 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.
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 , and the weight of the brain equals (see [16]), and setting the density of the brain equals , 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).
Data Availability
The data used can be found in the references cited in this paper.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The work was supported by the Klaus Tschira Foundation, Buhl-Strohmaier Foundation, and Würth Foundation.