## Abstract

Most oxygen required to support the energy needs of vertebrate tissues is delivered by diffusion from microvessels. The presence of red blood cells (RBCs) makes blood flow in the microcirculation highly heterogeneous. Additionally, flow regulation mechanisms dynamically respond to changes in tissue energy demand. These spatiotemporal variations directly affect the supply of oxygen to parenchymal cells. Due to various limiting assumptions, current models of oxygen transport cannot fully capture the consequences of complex hemodynamic effects on tissue oxygenation and are often not suitable for studying unsteady phenomena. With our new approach based on moving RBCs, the impact of blood flow heterogeneity on oxygen partial pressure (Po_{2}) in the tissue can be quantified. Oxygen transport was simulated using parachute-shaped solid RBCs flowing through a capillary. With the use of a conical tissue domain with radii 19 and 13 μm, respectively, our computations indicate that Po_{2} at the RBC membrane exceeds Po_{2} between RBCs by 30 mmHg on average and that the mean plasma Po_{2} decreases by 9 mmHg over 50 μm. These results reproduce well recent intravascular Po_{2} measurements in the rodent brain. We also demonstrate that instantaneous variations of capillary hematocrit cause associated fluctuations of tissue Po_{2}. Furthermore, our results suggest that homogeneous tissue oxygenation requires capillary networks to be denser on venular side than on arteriolar side. Our new model for oxygen transport will make it possible to quantify in detail the effects of blood flow heterogeneity on tissue oxygenation in realistic capillary networks.

- oxygen transport
- microcirculation
- red blood cells
- hematocrit
- blood flow heterogeneity

the supply of oxygen to tissues is an essential function of the vertebrate circulatory system. Oxygen bound to hemoglobin is carried from the lungs by the blood circulation to the target regions and finally reaches individual cells by diffusive transport from microvessels. Red blood cells (RBCs) make up ∼45% of the blood volume and contain hemoglobin, which is the main oxygen carrier. Gas exchange mostly occurs in the microcirculation, where erythrocytes and vessel diameters are similar in size. In particular, RBCs need to deform to enter capillaries. The particulate nature of blood has profound effects on hemodynamics and hence on oxygen transport. Blood rheological properties and the complex geometry of microvascular networks cause large variations of hematocrit, which are specific to the microcirculation. Additionally, the microcirculation is a dynamic system that adapts to changes in energy metabolism. In the brain, blood flow is controlled by arterioles as well as capillaries (13); in muscles, capillary recruitment increases the surface area for diffusion in response to contractile activity (31). The temporal and spatial variations in the microcirculation render investigations by both experiments and theoretical models challenging. However, new experimental techniques such as two-photon phosphorescence lifetime microscopy were applied to measure in vivo oxygen tensions at depths up to 300 μm (19). In spite of these advances, control of physiological parameters and simultaneous measurements at multiple locations remain difficult to achieve. Theoretical models for oxygen transport ideally complement experiments by providing precise control on all variables and making it possible to isolate their individual influence.

Oxygen modeling started with the seminal work of Krogh (18). For a tissue cylinder with a capillary at its center, the Krogh-Erlang equation yields an estimate of the oxygen gradient that is required to sustain a given rate of oxygen consumption. In the 1970s, Hellums (14) modeled for the first time oxygen transport with individual RBCs and coined the term “erythrocyte-associated transients” (EATs). The presence of EATs in the blood was observed experimentally ∼30 yr later by Golub and Pittman (11) and confirmed with micrometric resolution by Parpaleix et al. (24). Further modeling studies have extended the original Krogh model and considered microvascular networks.

Models for oxygen transport in the microcirculation were reviewed by Goldman (8). Current models for oxygen transport from capillaries to tissue generally employ two distinct approaches. The first class of models focuses on the tissue and does not represent individual RBCs. Instead, they employ a boundary condition at the capillary wall that accounts for oxygen transport from the capillary. While the original Krogh model assumed a constant oxygen tension at the capillary wall, more recent models often use a mass transfer coefficient (MTC) that relates the Po_{2} drop from the RBC to the oxygen flux across the capillary wall (*j* = *k*ΔP). Since these MTCs depend on hematocrit (15, 5), this approach captures the influence of RBC flow on tissue oxygenation. In addition, these models have the advantage that they do not resolve the complex intravascular Po_{2} field with individual RBCs, which makes them applicable to capillary networks (9, 29, 10, 35). However, this first class of models is dependent on other models that compute the oxygen flux out of capillaries.

The second approach models intravascular oxygen transport in more detail and can be used to compute MTCs. Accurate MTC estimates require discrete RBCs to be modeled (14, 6, 15) (as opposed to a continuous hemoglobin solution) and extracapillary oxygen transport to be included (5). Most models with individual RBCs carry out computations in the frame of reference of the erythrocyte, which simplifies the numerical treatment of the reaction between oxygen and hemoglobin in RBCs. In this moving frame, the tissue has an apparent velocity opposite to the RBC velocity and appears to move backwards. This idea was first used by Hellums (14), who used an analytical model with a cylindrical RBC and the adjacent tissue to compute MTCs. Eggleton et al. (5) built on this approach and used a model with concentric layers around the capillary for wall, interstitial fluid, and the tissue. They investigated the dependence of MTCs on hematocrit, RBC velocity, and capillary radius. The resulting MTCs can then be used in simulations of oxygen transport in complex capillary networks (9, 29, 10, 35).

Although the models for intravascular oxygen transport described above are convenient for numerical computations and useful for estimating MTCs, they suffer from limitations that restrict their scope. In the RBC frame of reference, the boundary condition at the distal end of the tissue cylinder has a considerable effect on tissue Po_{2} since the Po_{2} value at that boundary is advected backwards by the apparent tissue motion. Therefore, models that use the RBC frame cannot fully capture the influence of RBC flow on tissue Po_{2}, which is essential in applications such as hypoxia. These models are also inflexible in terms of geometry, since the backward motion of the tissue forces the computational domain to have the same radial cross section along the flow direction. For instance, local capillary dilations, as observed in vivo (13), cannot be simulated with this class of models. Furthermore, the simulation duration is limited to the time that RBCs spend in capillaries [100–300 ms in the cerebral cortex (16)]. For applications that require a larger simulation time (e.g., functional hyperemia), it is also necessary to use the frame of reference of the tissue, as done by models based on MTCs. Unlike other studies, Groebe and Thews (12) modeled individual RBCs in a fixed tissue region. However, their approach is limited to steady-state situations and relies on multiple simplifying assumptions that allow an analytic treatment of the intra-erythrocyte Po_{2} field.

Finally, Goldman (8) pointed out that thorough model validations have yet to be done. For intravascular Po_{2}, this task puts constraints on both the simulation method and the required experimental data. Since Po_{2} is generally measured at one or more fixed locations, a convenient model validation should be performed in the fixed frame of reference of the tissue. In addition, a detailed comparison with measured intravascular Po_{2} requires high spatial and temporal resolution. Pioneering work by Vanzetta and Grinvald (37) has revealed Po_{2} transients related to neuronal activation and oxygen metabolism with the use of phosphorescence lifetime microscopy. Using one photon excitation with a lower excitation volume, Golub and Pittman (11) measured EATs in the rat mesentery. However, until now, only two-photon phosphorescence lifetime microscopy achieved sufficiently high resolution to enable in vivo measurements of the Po_{2} between RBCs in depth. This technique was applied by Parpaleix et al. (24) in the olfactory glomerulus of the rodent brain. Sakadžić et al. (27) used it in the rat cerebral cortex, without reporting details of the intravascular Po_{2} field. Due to the absence of other detailed experimental studies, we compared our simulation results with the data from Ref. 24.

We propose a new model of oxygen transport in the microcirculation that is adapted for validation against experimental data. The main improvement over previous models is the use of overlapping meshes, which simultaneously allows the frame of reference of the tissue to be fixed and individual RBCs to be modeled. Hence, the coupling between intravascular oxygen transport and tissue Po_{2} can be captured together with the details of the Po_{2} field inside and around capillaries. Individual RBCs are followed by moving meshes that are used to compute hemoglobin diffusion and reaction with oxygen. These moving meshes are mapped onto a fixed mesh, where oxygen advection, diffusion, and consumption in the tissue are computed. This approach can capture the influence of heterogeneous RBC flow on tissue oxygenation in a time-dependent manner. Situations with unsteady blood flow such as functional hyperemia can be modeled by adapting blood velocity and hematocrit. A thorough comparison with the experimental data from Parpaleix et al. (24) showed that both intra- and extravascular oxygen transport are accurately simulated. For this comparison, an axisymmetric geometry based on Eggleton et al. (5) with concentric layers for the plasma, the capillary wall, and tissue was used. However, we found that a cone-shaped geometry as used by Hudetz (17) yields a better agreement with the measurements than a cylinder with constant radius. MTCs were also compared with results from previous models.

Although we apply this new model to an axially symmetric geometry, our algorithm is formulated in a general way and can be applied to arbitrary geometries. Therefore, with the use of a model for RBC transport (e.g., Ref. 23) to compute RBC trajectories, oxygen transport can be simulated in arbitrary capillary networks with realistic RBC dynamics. Our efficient time-stepping scheme allows taking large time steps and makes our model tractable in complex geometries. This will enable the investigation of the effects of blood flow heterogeneity during physiologically relevant phenomena such as microstrokes or capillary dilations (13).

## METHODS

#### Mathematical model.

Oxygen transport and consumption was modeled in a domain that consists of four regions: tissue, capillary wall, plasma, and RBCs. Oxygen is consumed only in the tissue; the capillary wall does not consume oxygen and has a lower diffusion coefficient; in both plasma and RBCs, oxygen is convected by the blood flow. Finally, RBCs contain hemoglobin, which carries oxygen in bound form. In fact, due to the low solubility of oxygen in plasma, most oxygen in capillaries is bound to hemoglobin.

Dissolved oxygen can be quantified by its concentration C [mlO_{2}/cm^{3}] and partial pressure P = Po_{2} [mmHg], which are related by Henry's law as
(1)

where α is the solubility coefficient in mlO_{2}·cm^{−3}·mmHg^{−1}. The formulation of the conservation equation for oxygen in terms of C = αP is most convenient for our purposes. Hemoglobin is expressed using the saturation S, which is the concentration ratio of oxyhemoglobin to total hemoglobin.

The reaction between oxygen and hemoglobin in RBCs is most completely described by the Adair equation (3). However, as in many previous studies, here we employ the Hill equation (2)

to describe the equilibrium curve between P and S, where P_{50} is the oxygen partial pressure at hemoglobin half-saturation and *n* is the Hill exponent. This results in a one-step reaction for the four heme groups of the hemoglobin molecule. To model the reaction rates when oxygen and hemoglobin are in nonequilibrium, we followed the approach of Clark et al. (3) and used the function
(3)

where *k*_ is the dissociation rate. This function satisfies *f* = 0 when oxygen and hemoglobin are in equilibrium (*Eq. 2*). Since no hemoglobin is present in healthy blood plasma, the reaction term *f*(P, S) was only used within RBCs.

Oxygen consumption was modeled using first-order Michaelis-Menten kinetics (8) and assumed to occur only in the tissue, which results in (4)

where M_{0} is the maximal metabolic rate of oxygen consumption in mlO_{2}·cm^{−3}·s^{−1} and P_{crit} is the oxygen level at which consumption is half of M_{0}. Since we compared our results with measurements performed in the rodent brain where no muscles are present, we did not consider myoglobin-facilitated diffusion of oxygen inside the tissue.

Our model is based on a single equation for oxygen for all regions, that is, (5)

where *D* is the diffusion coefficient and *v* the advection velocity. The factor *c* is given by *c* = N_{Hb}V_{mol,O2}, where N_{Hb} is the molar density of heme groups and V_{mol,O2} is the molar volume of oxygen. Hemoglobin saturation is governed by the equation
(6)

where *D*_{Hb} is the diffusivity of hemoglobin in RBCs.

At interfaces between regions with different solubility or diffusion coefficients, continuity of Po_{2} and oxygen flux across the interface have to be satisfied (39). For example, at the wall-tissue interface, the latter condition is
(7)

where the subscripts refer to the wall and the tissue, respectively.

The choice of boundary conditions depends on the computational domain. In this study, we considered representative domains with = 0 at the tissue boundary. At the capillary entrance, a Po_{2} value is required since oxygen is convected into the domain by the blood flow. When a RBC overlaps with the domain boundary, the oxygen tension is interpolated from this RBC to the capillary entrance. When plasma is flowing in, a constant Po_{2} value P_{p,in} was used. At the capillary outlet, the boundary condition = 0 was applied.

Since RBC membranes are impermeable for hemoglobin, the boundary condition for hemoglobin saturation is = 0. Unlike hemoglobin, oxygen is soluble in lipids and can diffuse through cell membranes. The different solubility and diffusion coefficients of oxygen in lipid bilayers were not taken into account since RBC membranes are generally <10-nm thick (33), which is negligible compared with the cell size.

The entry of RBCs into the capillary plays a crucial role, since it determines the amount of oxygen in bound form that enters the domain. The oxygen tension in entering erythrocytes was set to a constant value P_{rbc,in}. The simplest model for capillary spacing is a constant distance between each RBC pair. However, Chaigneau et al. (2) observed large instantaneous fluctuations of the RBC linear density. Moreover, they showed that variations of RBC flow were primarily caused by fluctuations of linear density, whereas instantaneous RBC velocity fluctuations were 2.5 times lower. Therefore, we treated RBC spacings as a random variable and modeled it using a log-normal random variable with independent values for each RBC pair. The parameters were chosen to match experimentally measured mean μ_{LD} and standard deviation σ_{LD} of linear density.

The initial Po_{2} field in RBCs was set to P_{rbc,in}, and hemoglobin saturation was set to equilibrium with oxygen. Outside RBCs, the initial Po_{2} was set to P_{p,in} in the plasma and to 22 mmHg in the tissue.

#### Discretization.

The main objective of this study is to thoroughly compare simulation results with experimental data. To allow an easy comparison with measurements, the numerical model should reflect how experiments are carried out. Our reference data (24) were acquired using two-photon phosphorescence lifetime microscopy. Thus measurements were obtained from the focal plane of the microscope, which may contain both capillaries and tissue. An easy comparison with these data requires a model that focuses on a fixed region. This approach also enables capturing transient phenomena such as local changes in RBC flow or metabolism.

The fixed frame of reference motivated above is problematic when solving *Eq. 6*. Hemoglobin is a large protein that cannot cross erythrocyte membranes. However, the discretization of the advection term would create numerical diffusion, which would in turn cause an unphysical leak of hemoglobin out of RBCs. These problems can be circumvented by solving *Eq. 6* in a Lagrangian frame of reference that follows the moving RBC. This approach enables the no-flux boundary condition for hemoglobin at the RBC membrane to be exactly satisfied.

We therefore used a fixed computational domain for the capillaries and the tissue, denoted by Ω, as well as a moving domain for each RBC, denoted by Ω_{rbc} (Fig. 1). Each domain is covered by its own computational mesh. This overlapping mesh approach was adapted from the overset grid method (26), which has been applied to aerodynamic problems with moving objects. We will also refer to as Ω Eulerian domain and to Ω_{rbc} as Lagrangian domain. To simplify the notation, we omit RBC indexes. Since RBCs are entering and leaving Ω, the Lagrangian domain Ω_{rbc} may be completely or partly inside Ω.

Erythrocytes were assumed to have a fixed shape. While they actually deform, this assumption avoided the expensive treatment of fluid-structure interaction. Therefore, our modeled RBCs behaved similar to solid bodies that follow the plasma flow. As a further simplification, we considered plasma flow to be uniform along radial cross sections of capillaries. Note that the detailed flow field around RBCs is not of importance here, since transport of oxygen is diffusion dominated (see Ref. 36 for a corresponding study about nitric oxide). Consequently, the blood velocity was given by *v* = *Q/A*, where *Q* is the blood volume flow and *A* is the capillary cross section.

*Equation 5* for oxygen was solved in the Eulerian domain Ω, whereas the hemoglobin *Eq. 6* was solved in the Lagrangian domain Ω_{rbc}. Since Ω_{rbc} moves with the velocity *v*_{rbc}, the coordinate transformation *x*′ = *x* + *v*_{rbc}*t* cancels the advection term and yields
(8)

Since this equation is discretized in Ω_{rbc}, the oxygen partial pressure is also needed in that domain. This field, denoted by P_{rbc}, is obtained by interpolation from Ω to Ω_{rbc}. Likewise, since *Eq. 5* is solved in Ω, values of S in the Eulerian domain, denoted by S_{Euler}, have to be interpolated from Ω_{rbc} (Fig. 1). The interpolation method may considerably affect simulation results, since most oxygen in the blood is bound to hemoglobin. Thus interpolation errors that cause inaccurate values of S_{Euler} may have a large effect on the resulting Po_{2}. A conservative interpolation scheme is therefore crucial.

To obtain P_{rbc} and S_{Euler}, we used a volume-based interpolation scheme that is discretely conservative in the sense that the integral of the interpolated field on any subset of the target mesh is conserved. For grid cells V_{I} and V_{rbc,J} in Ω and Ω_{rbc}, respectively, interpolation weights were defined by
(9)

and (10)

The interpolation formulas for P_{rbc} and S_{Euler} are then given by
(11)

and (12)

The discrete conservation property for the interpolated field S_{Euler} is shown as follows. Consider a subdomain Ω′ = V_{Ik} that consists of *m* grid cells V_{Ik}. The integral of S_{Euler} on Ω′ is given by
(13)
(14)
(15)
(16)
(17)

The same argument can be used for the integral of P_{rbc} on a subset of Ω_{rbc}, which shows that the interpolation scheme given by *Eqs. 11* and *12* is discretely conservative.

Grid cells in Ω that overlap with the RBC border require special care. If the intersection of a grid cell V_{I} with Ω_{rbc} occupies a small volume, S_{Euler,I} will be also small. This fact has to be accounted for in the discretization of the reaction term *f*(P,S). We introduce the RBC volume fraction
(18)

In V_{I}, we consider that the chemical reaction between hemoglobin and oxygen only occurs in a fraction of V_{I} with volume γ_{I}|V_{I}| where all the hemoglobin is contained. Since this volume fraction has hemoglobin saturation S_{Euler,I}/γ*I*, the discretized reaction term in is given by
(19)
(20)

Continuity of the oxygen flux at interfaces between regions with different solubility or diffusion coefficient (*Eq. 7*) is enforced by adequately interpolating the Krogh diffusion coefficient *D*α. At cell faces, mass conservation is enforced by using the harmonic average of *D*α in both neighboring grid cells (25). The boundary condition at the capillary inlets of Ω also requires interpolation. If a RBC overlaps a cell face at the capillary inlet, the Po_{2} value at that face is obtained by bilinear interpolation of the RBC Po_{2} at the corresponding location. Otherwise, the boundary Po_{2} is set to the constant value P_{p,in}.

The governing equations were discretized using a finite-volume method with the central scheme for the divergence operator. For the Laplace operator, Gauss integration, centered differences for the surface normal gradient and harmonic interpolation for the diffusion coefficient were used. Time stepping and coupling between *Eqs. 5* and *6* are addressed in appendix a. The algorithm was implemented using the open source software package OpenFOAM v.2.1.1.

#### Model parameters.

Our main goal is the validation of the method explained above against the experimental data from Parpaleix et al. (24). These data were acquired in the rodent olfactory glomerulus, which is an area with a high capillary density.

We used an axially symmetric geometry with a capillary at its center, similar to the classical Krogh model (18). Instead of a cylinder, we employed a cone-shaped domain with different radii at the proximal (arteriolar) and distal (venular) ends. Due to symmetry, Ω can be represented by a two-dimensional domain. As shown in Fig. 2, Ω consists of three regions, that is, the plasma, the capillary wall, and the tissue region.

In the olfactory glomerulus, the average distance from any point to the nearest capillary is 10.8 μm (20). In a hexagonal array of Krogh cylinders with a capillary diameter of 4 μm, this corresponds to a radius of 16 μm. Therefore, unless stated otherwise, the radii on the arteriolar and venular sides were set to *r*_{t,a} = 19 μm and *r*_{t,v} = 13 μm, respectively. The length of the capillary was set to 100 μm.

The RBC shape was taken from Secomb et al. (30) for a RBC velocity of 1 mm/s. This shape (computed for human RBCs) was scaled down to the size of mouse erythrocytes with volume V_{rbc} = 59.0 fl (32). We used the mean RBC velocity *v*_{rbc} = 0.57 mm/s measured in the olfactory glomerulus by Chaigneau et al. (2).

The cerebral metabolic rate of oxygen consumption (CMRO_{2}) is an essential model parameter. To our knowledge, no measurement of CMRO_{2} in the olfactory glomerulus has been performed. Therefore, we chose the value CMRO_{2} = 197 μM/s to obtain Po_{2} values in the tissue between 15 and 20 mmHg approximately (using the perfect gas law at 36.9°C, this corresponds to M_{0} = 5 × 10^{−3} mlO_{2}·cm^{−3}·s^{−1}). The resulting values of Po_{2} in the plasma agree well with the results of Parpaleix et al. (24).

## RESULTS

We now show simulated oxygen tensions inside the sample capillary and the surrounding tissue region shown on Fig. 2. Whenever possible, we compare our results with the data measured by Parpaleix et al. (24) using two-photon phosphorescence lifetime microscopy in the rodent olfactory glomerulus. They characterized intracapillary oxygen tensions by the three following quantities: RBC Po_{2}, mean Po_{2}, and inter-RBC Po_{2}. RBC Po_{2} is the maximal oxygen tension in the plasma, which is attained at the erythrocyte membrane. Mean Po_{2} is the average Po_{2} value between two erythrocytes, and inter-RBC Po_{2} is the minimal Po_{2} between two RBCs. The EAT amplitude is the difference between RBC Po_{2} and inter-RBC Po_{2}. Throughout this section, the coordinate *x* denotes the axial direction.

Using the parameters listed in Table 1, we obtained an averaged EAT amplitude of 29.7 mmHg (RBC Po_{2} = 50.8 mmHg, inter-RBC Po_{2} = 21.1 mmHg, and mean Po_{2} = 27.4 mmHg). These values were obtained by sampling Po_{2} on the capillary centerline at nine evenly spaced longitudinal locations (between *x* = 10 and 90 μm). The maximal Po_{2} in the plasma was attained on the rear side of the RBC membrane. Parpaleix et al. (24) also observed significant differences between these quantities [RBC Po_{2} = 57.1 ± 1.3 mmHg (means ± SE), inter-RBC Po_{2} = 23.6 ± 0.7 mmHg, and mean Po_{2} = 30.8 ± 0.9 mmHg]. Since they performed 241 measurements, the results for our sample capillary differ from these average values by less than one-third of a standard deviation.

Figure 3 shows instantaneous longitudinal profiles on the capillary centerline and at various radial distances from the capillary wall. In RBCs close to the arteriolar end of the domain, the intracellular Po_{2} variation exceeds 30 mmHg and decreases to 15 mmHg at the venular end. These strong intravascular oxygen variations extend to the nearby tissue. At 1 μm from the outer side of the wall, the amplitude of these fluctuations ranges from 12.7 to 4.2 mmHg. Away from the capillary entrance, these values agree well with the mean pulse amplitude of 5.0 mmHg reported by Parpaleix et al. (24) outside the vessel (<2 μm). At 5 μm from the endothelium, these pulses are almost entirely smeared out. The influence of instantaneous linear density fluctuations on inter-RBC Po_{2} is clearly illustrated by the second and the third RBC spacings. Since short RBC spacings cause higher inter-RBC Po_{2} values, the EAT amplitude drops when the instantaneous linear density increases.

We then investigated longitudinal variations of Po_{2} along our sample capillary. Figure 4 shows time-averaged oxygen partial pressures for the cone-shaped geometry (Fig. 2) and for a cylinder with equal tissue volume. Since RBC Po_{2} declines faster than inter-RBC Po_{2}, the EAT amplitudes also decrease along the capillary. Parpaleix et al. (24) reported longitudinal variations of Po_{2} in single capillaries over a mean distance of 49.7 μm. Table 2 contains these values as well as our simulated Po_{2} variations in the conical and cylindrical geometries. The maximal gradients in the cone-shaped geometry are a consequence of the high RBC Po_{2} at the capillary entrance. However, the gradients away from the arteriolar end of the domain correspond very well to the experimental data, while in the cylinder geometry the gradients of mean Po_{2} and inter-RBC Po_{2} are significantly higher than in the reference data. A better match could not be obtained in a cylindrical geometry by changing CMRO_{2}, since this would considerably decrease the agreement of RBC Po_{2} and inter-RBC Po_{2} with experimental data. The chosen geometry with *r*_{t,a} = 19 μm and *r*_{t,v} = 13 μm had the smallest taper that yielded a good match with the measured longitudinal Po_{2} variations. These results suggest that a cylindrical geometry is not a suitable model for capillaries, at least in the brain region considered in this study.

Our model includes instantaneous variations of linear density similar to those observed by Chaigneau et al. (2). Figure 5 shows values of RBC Po_{2} and inter-RBC Po2 that were collected during 3 s at 30 μm from the capillary entrance. The linear density on the horizontal axis was quantified by the length occupied by RBCs over a given capillary segment divided by the segment length. As previously observed in Fig. 3, inter-RBC Po_{2} is correlated with the linear density. The dependency of inter-RBC Po_{2} on linear density agrees very well with the experimental data, but the simulated RBC Po_{2} is almost constant, while the reference data exhibit a positive correlation between linear density and RBC Po_{2}. Our simulations did not reproduce this trend, since a single capillary with constant RBC Po_{2} at its arteriolar end was used. However, Parpaleix et al. (24) measured EAT properties in 42 capillaries, which limits the scope of this comparison. This difference between the pooled experimental data and our computations in a single capillary indicates that capillaries with high average linear density also have a higher Po_{2}. Besides, Parpaleix et al. (24) have observed that inter-RBC Po_{2} attains similar values as Po_{2} in the neuropil. Figure 5 also shows the difference between inter-RBC Po_{2} and tissue Po_{2} at 10 μm from the capillary wall as a function of linear density. For linear densities lower than 0.25, this difference stays below 2.0 mmHg. For high hematocrit values, this gap exceeds 10 mmHg. Thus our results indicate that inter-RBC Po_{2} may significantly exceed tissue Po_{2} for high linear densities.

Since linear density affects tissue Po_{2}, we investigated the influence of the standard deviation σ_{LD} of linear density on tissue Po_{2}. Figure 6 shows tissue Po_{2} at 10 μm from the capillary wall and *x* = 50 μm for two different values of σ_{LD}. The same random numbers were used and the parameters of the log-normal distribution for RBC spacings were adjusted to obtain an average linear density of 0.28 over 4 s and the desired standard deviation. Only the last second of the simulation is shown. Random fluctuations of linear density led to large Po_{2} oscillations. For σ_{LD} = 0.08, the difference between minimal and maximal Po_{2} was 5.7 mmHg, and for higher fluctuations (σ_{LD} = 0.16), it increased to 10.9 mmHg. This is a consequence of RBC groups that are close to or far away from each other. Occasionally, a large RBC spacing resulted in a sudden drop of tissue Po_{2} by several millimeters of mercury. Therefore, if linear density fluctuations as reported by Chaigneau et al. (2) are present, Po_{2} in the tissue cannot be considered to be constant.

Finally, we compare our results with previous works by examining the intracapillary resistance to oxygen transport. MTCs were computed using a constant linear density and compared with previously published values. The MTC may be defined as *k* = *j*/(P* − P_{w}), where *j* is the oxygen flux (mlO_{2}·cm^{−2}·s^{−1}), P* is the oxygen tension in equilibrium with the mean hemoglobin saturation in the RBC, and P_{w} is the average oxygen tension at the capillary wall around a RBC. For a tube hematocrit of 0.25, we obtained *k* = 1.67 × 10^{−6} mlO_{2}·cm^{−2}·s^{−1}, which exactly matches the results of Eggleton et al. (5) for the same hematocrit and capillary radius (*r*_{p} = 2.0 μm). This consistency was expected, since the same equations as in Ref. 5 were solved (except myoblogin-facilitated diffusion in the tissue) and similar diffusion and solubility coefficients were chosen.

Comparison with earlier works can also be performed using the Nusselt number, which is defined by (21)

where d_{p} is the capillary diameter. For tube hematocrit values between 0.15 and 0.36, we obtained Nusselt numbers from 0.48 to 1.7. Hellums et al. (15) summarized Nusselt numbers from various studies. For a diameter of 3.6 μm and a tube hematocrit of 0.28, Secomb and Hsu (28) calculated Nu = 1.22 using a solid cylinder model. Our computed value for this tube hematocrit is 1.17. Therefore, our model reproduces oxygen fluxes from previous studies in steady-state situations.

## DISCUSSION

Oxygen transport from a capillary with moving RBCs to the surrounding tissue has been simulated in an axisymmetric cone-shaped geometry. Oxygen partial pressure in the capillary and the tissue was compared with experimental data (24). Longitudinal oxygen variations and the influence of linear density were investigated. As an application of our model, we studied the impact of instantaneous hematocrit fluctuations on tissue oxygenation.

Our simulations reproduced a number of results from Parpaleix et al. (24). Their average measured EAT amplitude was 33.5 mmHg, and similar amplitudes were obtained in the first section of our sample capillary (Fig. 4). At 30 μm from the capillary entrance, the simulated EAT amplitude was 33.6 mmHg. Close to the venular end, RBC Po_{2} was lower due to oxygen consumption in the tissue, which gave rise to smaller EATs (<25 mmHg). Therefore, our average EAT amplitude of 29.7 mmHg over the nine sampled positions is slightly lower than that from Parpaleix et al. (24). Since the experimental data were collected independently of the measurement position in the vascular bed, it is difficult to further interpret these differences. However, the dependency of EAT values on the distance from the arteriolar side could for example be studied experimentally in the brain cortex.

The relationship between intracapillary oxygen tensions and tissue Po_{2} was also examined. For linear densities <0.25, simulated inter-RBC Po_{2} exceeds tissue Po_{2} at 10 μm from the capillary wall by <2.0 mmHg (Fig. 5), while this difference is >10 mmHg for higher hematocrit values. These findings are only in partial agreement with the observation by Parpaleix et al. (24) that inter-RBC Po_{2} attains similar values as in the neuropil. However, measurements in capillaries and tissue were not performed simultaneously and results were averaged over several seconds, which filtered out Po_{2} fluctuations, whereas we report instantaneous snapshots. Moreover, the influence of hematocrit fluctuations was not examined in this part of the experiment. Therefore, our simulations indicate that inter-RBC Po_{2} is similar to tissue Po_{2} only close to the capillary or at low linear densities. Since concentration gradients drive molecular diffusion, we suggest that inter-RBC Po_{2} is on average higher than tissue Po_{2} far away from capillaries, provided they are not close to an arteriole. This hypothesis can be tested in vivo by measuring the dependency of tissue Po_{2} on the distance to the nearest capillary.

Our simulation setup with RBCs moving through a fixed capillary allows the computation of longitudinal oxygen gradients. Motivated by the fact that capillary segments with high oxygen tensions can supply a correspondingly large tissue volume, we used a cone-shaped geometry (Fig. 2) similar to Hudetz (17). We compared results obtained with this geometry and with a simple cylindrical domain to the data (24), where longitudinal Po_{2} variations were measured in individual capillaries. While gradients of mean Po_{2} and inter-RBC Po_{2} in the classical Krogh cylinder geometry are much higher than in the reference data (Table 2), the cone-shaped domain leads to a very good agreement. Although a conical geometry is idealized, it appears to be a suitable model to reproduce in vivo intracapillary oxygen gradients in the brain. This finding may imply that capillary density increases along RBC paths through capillary networks. In other words, we suggest that an evenly distributed tissue Po_{2} requires denser capillary networks on venular side. However, one should examine whether these simulation results hold in realistic networks, where capillary interactions and tortuosity are present.

Instantaneous variations of hematocrit as observed by Chaigneau et al. (2) can be accounted for by our model, which overcomes a limitation of the models based on MTCs. We treated linear density as a random process governed by a log-normal RBC spacing distribution. The resulting dependency of inter-RBC Po_{2} on linear density agrees very well with the data from Parpaleix et al. (24) (Fig. 5). On the other hand, RBC Po_{2} stayed constant, which means that the drop in hemoglobin saturation along RBC paths was not influenced by instantaneous hematocrit fluctuations. Since our results were produced in one sample capillary and the data from Parpaleix et al. (24) were pooled from 42 capillaries, we propose the following interpretation of this discrepancy: while fast fluctuations of linear density do not influence RBC Po_{2}, capillaries with high average hematocrit have a higher RBC Po_{2}. This explanation should be investigated by measuring RBC Po_{2} in capillaries that have different average linear densities. Additionally, these hematocrit fluctuations also affect tissue Po_{2} (Fig. 6). With a RBC length of 7.27 μm, the standard deviation of linear density reported by Chaigneau et al. (2) is 0.12. Our results show that for this value, oscillations of oxygen tension in the tissue approach 10 mmHg. During transient periods of low RBC density and/or velocity, it therefore seems possible that tissue oxygenation drops at times below the critical level for oxidative phosphorylation, although the average tissue Po_{2} remains above this level. Since the geometry of complex capillary networks affects tissue Po_{2}, it will be essential to further study the influence of linear density fluctuations.

Although multiple experimental results could be reproduced, the simulation setup presented here has several limitations, in particular the axisymmetric geometry. While such a geometry is most relevant for parallel capillary arrays in muscles, Krogh cylinder models fail to capture the minimal tissue Po_{2} in the capillary beds of the brain cortex (29). Accordingly, our conclusions on the relationship between inter-RBC Po_{2} and tissue Po_{2} will certainly need to be refined for realistic networks. The hypothesis that capillary networks are denser on venous side should also be verified in such networks. Nevertheless, the simulated oxygen tensions in the plasma mainly depend on hemoglobin saturation in nearby erythrocytes and should not be directly affected by diffusive interactions between capillaries. This is confirmed by the good agreement between the simulated inter-RBC Po_{2} and experimental data (Fig. 5).

Other limitations include constant blood velocity, the absence of shifts of the oxygen-hemoglobin dissociation curve, and the uncertainty in the choice of parameters. While RBC velocity undergoes fluctuations, their amplitude is lower than that of linear density (2); hence, we chose to keep it constant. However, RBC velocity is an important factor for tissue oxygenation and should be realistically modeled. In addition, variations of carbon dioxide concentration and pH are known to shift the equilibrium curve modeled by *Eq. 2*. This may be significant in regions with low Po_{2} and high CO_{2} concentration (4). The inclusion of these shifts would require further modeling efforts. Finally, tissue oxygenation highly depends on CMRO_{2}, which is difficult to measure experimentally. Our chosen value (197 μM/s) is almost three times as high as the CMRO_{2} in the cortex of awake rats (73.5 μM/s), which was obtained using the value 420 μmol (100 g)^{−1}·min^{−1} (7) and a brain density of 1.05 g/cm^{3} (22). Based on estimates by Nawroth et al. (21), the neuron density in the olfactory glomerulus of the rat is 6.9 × 10^{5} cells/mm^{3}, whereas this value is 1.17 × 10^{5} in the mouse neocortex (34). The high density of neural elements (possibly in combination with a high steady-state firing rate) in the olfactory glomerulus may explain why a high CMRO_{2} value was needed to reproduce the tissue Po_{2} observed by Parpaleix et al. (24). However, using a theoretical energy budget, Nawroth et al. (21) obtained a CMRO_{2} value of 75 μM/s for the olfactory glomerulus, which is lower than our chosen value. Further interpretation of this discrepancy would require actual measurements of CMRO_{2} in the olfactory bulb.

In addition to model limitations, the comparison with experimental data is also limited. To the author's knowledge, only the data from Ref. 24 allowed a detailed comparison of simulated intracapillary Po_{2}. A good agreement was obtained by adapting CMRO_{2} and the initial Po_{2} in RBCs on the arteriolar side, and by choosing a tapered cylinder. Further data on intracapillary Po_{2} and its relationship to tissue Po_{2} should be obtained and compared with (24). The parameters mentioned above will most likely need to be modified to reproduce further experiments. The computational model presented in this study will be a useful tool to interpret possible differences between future experimental data.

Although our model for oxygen transport was applied to a simple axisymmetric geometry, the numerical algorithm is independent of the domain topology and can be extended to realistic capillary networks provided velocities of single RBCs are known. This can be achieved by coupling our method with a detailed model of RBC transport such as that of Obrist et al. (23). This combined approach will remove the need for separately computed mass transfer coefficients and is suitable for investigating unsteady scenarios. For example, Hall et al. (13) recently observed that capillary pericytes participate in the regulation of cerebral blood flow. Our model will enable quantifying the influence of capillary dilations on tissue oxygenation. Therefore, our present study is a first step toward an oxygen transport model that can capture a wide range of dynamic physiological phenomena while taking into account the complex properties of RBC flow.

In conclusion, we have developed a new model of oxygen transport from capillaries with moving RBCs based on overlapping grids. We successfully validated it against experimental data acquired in the rodent brain. EATs and longitudinal gradients of Po_{2} could be reproduced using a cone-shaped geometry. Instantaneous variations of hematocrit were shown to cause considerable fluctuations of oxygen tension in the tissue. Further work includes the extension of the model to realistic capillary networks. The coupling of RBC dynamics with oxygen transport will eventually allow simulations of blood flow regulation mechanisms in health and disease with unprecedented detail.

## GRANTS

This research was funded by the Swiss National Science Foundation Grant No. 140660.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: A.L., B.W., and P.J. conception and design of research; A.L. performed the numerical simulations, analyzed the data, interpreted the simulation results, prepared figures, and drafted manuscript; A.L., B.W., and P.J. edited and revised manuscript; B.W. and P.J. approved final version of manuscript.

## ACKNOWLEDGMENTS

We thank S. Charpak for providing raw data presented in Parpaleix et al. (24). The authors are also grateful for valuable discussions with F. Schmid and for M. Barrett's useful comments on the manuscript.

## Appendix A: TIME INTEGRATION

Generation of Po_{2} maps in realistic capillary network may require simulations with at least hundreds of RBCs during several seconds. The ability to use large time steps is therefore crucial to keep the computational time sufficiently low. Special care is required to achieve this within our framework based on overlapping meshes. The nonlinear reaction term *f*(P, S) (*Eq. 3*) combined with RBC displacements prevents from using an explicit scheme. As observed by Clark et al. (3), the boundary layer inside erythrocytes is a region of chemical nonequilibrium, such that large explicit time steps inevitably cause overshooting. Another requirement is that the coupling between hemoglobin and oxygen equations conserves the total of free and bound oxygen.

To achieve this, we use Godunov splitting for *Eq. 5* and linearization of the reaction and consumption terms using Picard's method. While the equation for oxygen can be integrated without Godunov splitting, this unsplit approach would severely limit the maximal stable time step, since the linearization of the reaction term requires Po_{2} values in Ω to vary moderately. If RBCs undergo large displacements during one time step, the resulting large Po_{2} variations would lead to instabilities.

Let the superscript *k* indicate the current time *t*^{k}. To integrate *Eqs. 5* and *6* from *t*^{k} to *t*^{k} + Δ*t*, an intermediate solution P^{*} is obtained by integrating only the advection term:
(22)

Here, the solubility α* corresponds to RBC positions after their displacement. The reaction term *f*(P, S) and the consumption term M(P) were both linearized and their linear part is treated implicitly as
(23)

and (24)

where ν is the iteration number and P^{(0)} = P*. The coupling between both equations conserves the total oxygen amount, if the integral of both terms in square brackets are equal. Although the volume-based interpolation method (*Eqs. 11* and *12*) conserves P and S, it does not exactly conserve the integral of *f*(P, S) since the reaction term is nonlinear in P. However, this only causes a minimal amount of oxygen loss in the domain (<0.2% for total RBC discharge).

The moving meshes Ω_{rbc} are displaced during each time step by the increment *v*_{rbc}Δ*t*. When a RBC leaves the domain and no longer overlaps it, the corresponding mesh is moved to the front of the RBC queue and placed at a distance to the next RBC, which is randomly generated based on a log-normal distribution. In the plasma, the coefficients α and *D* have to be updated to reflect RBC motion. In a grid cell V_{I}, the discretized coefficients are given by
(25)
(26)

where the subscripts “rbc” and “p” refer to values in the RBCs and in the plasma. The algorithm is summarized in Table A1.

The domain Ω was discretized using a Cartesian grid with constant grid spacing Δ*x* = 0.1 μm in the axial direction. In the radial direction, the grid cell spacing in the capillary was constant (Δ*y* = 0.1 μm) and decreasing in the tissue region, since oxygen gradients decrease away from capillaries. The ratio between the height of the top-most grid cell to the bottom-most in the tissue was set to four. This results in a grid with 333 × 29 grid cells.

The RBC domain Ω_{rbc} consists of those Cartesian grid cells that lie entirely inside the RBC shape, which results in a “staircase” geometry (Fig. 1). A curvilinear shape-conforming mesh is not necessary for such an advection diffusion problem. In addition, the computation of the interpolation coefficients defined in *Eqs. 9* and *10* is easier for Cartesian grids.

The tolerance (tol) in the algorithm shown in Table A1 was set to 10^{−4}. A smaller tolerance affected results by <0.1 mmHg. Unless stated otherwise, the time step Δ*t* was set to 0.5 ms. All our simulations were run for 4 s. After 1 s, the influence of the initial condition disappeared. The results were collected during the following 3 s.

The accuracy of the algorithm with a coarser Eulerian grid and larger time steps was also examined. Table A2 shows absolute and relative errors on the capillary centerline and in the tissue against a baseline case with Δ*t* = 0.1 ms and Δ*x* = Δ*y* = 0.1 μm in the capillary. The relative error was normalized by the maximum Po_{2} value in the considered longitudinal profile. When multiplying the grid spacing and the time step by three, the relative error stays <2.5%. With a 50 times larger time step (Δ*t* = 5 ms), the absolute error in the tissue is still <1 mmHg, while the computational time is divided by 10. This is an indication that our numerical algorithm is very robust in terms of time step size and grid spacing. This property will allow for simulations of oxygen transport in larger capillary networks.

- Copyright © 2015 the American Physiological Society