Accounting for Tube Hematocrit in Modeling of Blood Flow in Cerebral Capillary Networks

The aim of this paper consists in the derivation of an analytic formula for the hydraulic resistance of capillaries, taking into account the tube hematocrit level. The consistency of the derived formula is verified using Finite Element simulations. Such an effective formula allows for assigning resistances, depending on the hematocrit level, to the edges of networks modeling biological capillary systems, which extends our earlier models of blood flow through large capillary networks. Numerical simulations conducted for large capillary networks with random topologies demonstrate the importance of accounting for the hematocrit level for obtaining consistent results.


Introduction
Simulation of blood circulation in large capillary networks is a challenging task. Realistic modeling of microvessel structures should take into account not only sophisticated topologies of blood vessel networks but also correct hydraulic resistance of microvessels. e latter is characterized by the apparent blood viscosity which depends on the vessel diameter as well as the discharge and tube hematocrit. e discharge hematocrit is the volume fraction of the red blood cells (RBCs) in the blood delivered by the flow in the vessel. e tube hematocrit is the volume fraction of RBCs that are inside the vessel at a given time instant. e discharge hematocrit is larger than the tube one because the velocity profile in the radial direction is nonuniform; namely, the RBCs velocity is higher than the mean bulk flow speed, which is called the Fåhraeus effect [1,2]. e velocity profile in the radial direction is affected by the presence of the endothelial surface layer (ESL) [1]. e importance of accounting for the hematocrit level in blood flow simulations attracts attention of many researches.
Animal models are utilized, for example, to measure and analyze the distributions of cell velocity and cell flux in the capillary network for different values of systemic hematocrit [3]. By using fluorescent microscopic analysis of rat cerebral capillary networks, the influence of hematocrit on mean RBC capillary velocity and mean arterial pressure can be assessed [4]. Moreover, the effect of hematocrit can be investigated in artificial microvascular branching networks [5]. A combination of an animal model with an effective iterative algorithm allows for finding the distribution of discharge hematocrit and blood flow velocity in a cerebrocortical microvascular network [6]. is approach takes into account the heterogeneity of blood flow and partitioning of red cells at bifurcations. e current paper is related to modeling of computergenerated blood microvessel networks with vessel diameters less than 10 μm. is is motivated by our previous work on simulating cerebral blood flow of preterm infants [7]. In such thin blood vessels, RBCs move in single file. Due to their ability to deform, RBCs can pass through vessels down to 2.7 μm in size, which is less than their diameter, without damaging their membrane [8]. e hematocrit level and the shape of single erythrocytes during their motion in a capillary with diameter less than 8 μm depend on RBC velocity [9].
Two approaches to mathematical modeling of RBC transport through microvessels can be mentioned. e first one is based on continuum models [7,10], where erythrocytes are considered as a homogeneous substance, and the RBC motion is described as a two-phase blood flow; namely, the erythrocyte homogeneous substance moves in the central core of the vessel, whereas the plasma fraction moves in a cell-free layer formation near the wall. Reasonable assumptions allow for deriving an explicit formula for the hydraulic resistance of a single capillary [7], which is very important for the simulation of blood flow in large capillary networks. A model derived in [10] studies the relations between the tube hematocrit level, vessel diameter, and apparent viscosity. e second approach relies on discrete models [11,12], where the effect of each erythrocyte is taken into account.
e capillary blood flow is considered as single-file flow of red cells in blood plasma playing the role of lubrication and filling gaps between the erythrocytes. Such an ansatz [11] is used to develop a model resulting in an efficient algorithm for computing the pressure and flow field as well as the hematocrit distribution in simplified capillary networks.
is approach shows a strong influence of single-file arrangement of RBCs on flow behavior. A coupled model describing the delivery of oxygen to tissue cells is considered [12] at the scale allowing to take into account the size and shape of individual RBCs as well as their deformation. e proposed approach [11,12] takes into account the level of hematocrit in simulations of blood flow.
In the present paper, the approaches described in [7,11] are combined to obtain an analytical formula for the computation of blood flow resistance in microvessels. is formula accounts for gaps between RBCs and, therefore, reflects the dependence on tube hematocrit. e tuning and validation of this formula are performed using hydrodynamical computations based on the representation of the cell-plasma mixture as a fluid with two different viscosities (much larger viscosity for blood cells). A very good consistency of the analytically computed values with the numerical results is obtained. An example of finding the pressure distribution in a relatively large capillary network (the germinal matrix or the whole brain), accounting for the level of tube hematocrit, is presented. e ability to account for the hematocrit level significantly enhances the algorithm proposed in [7] for finding the pressure distribution in the germinal matrix. Numerical experiments show a significant influence of the hematocrit level on the pressure distribution.

Continuous Model of Red Cell Transport
In [7], the transport of red cells in capillaries was modeled as a continuum flow with spatially variable viscosity. A high viscosity was assigned to the central part of the capillary, RBC substance, whereas the layer between the RBC substance and capillary wall was assigned with a small viscosity typical for blood plasma (Figure 1.) Such an approach is motivated by the results of [13] claiming that a rigid body moving in a fluid can be replaced with another fluid whose viscosity tends to infinity. Ignoring gaps between the red cells was explained by small length of gaps compared to a high velocity of RBCs. us, the resistance of a capillary was assumed to be independent on the hematocrit in a first approximation. e viscosity as function of the vessel's radius was chosen as follows: where r c is the radius of the capillary, r 0 is the radius of the RBC substance, and μ 1 and μ 2 are the viscosities of blood plasma and the RBC substance, respectively. Typically, μ 1 � 0.001 Pa·s and μ 2 � 0.1 Pa·s. e last high value is used here to make the RBC column effectively rigid. It should be noted that further increasing μ 2 has practically no effect. e velocity profile corresponding to (1) was derived in [7,14], under some assumptions, as follows: where A � Δp/L and Δp is the pressure drop at the capillary of the length L. e total flux is given by the following formula: and the resistance reads as follows: In [7] (see the end of Section 2 and references there), the following relation between r 0 and r c was proposed: and arguments for the reliability of this formula were presented. Notice that formula (4) with r 0 � 0 (i.e., the whole capillary is filled with blood plasma) transforms into the Poiseuille formula: Plasma RBC substance Plasma Figure 1: Schematic representation of continuum stream of erythrocytes in a capillary (the axial section is shown). e gaps between erythrocytes are ignored. e motion occurs due to the lubrication plasma layer between erythrocytes and vessel wall.
It should be noted that the model described in this section is a rough approximation of RBC motion in a capillary. Obviously, plasma recirculation between cells should affect the flow properties. Section 3 considers an extended model, assuming the presence of gaps between RBCs, which means accounting for tube hematocrit.

Discrete Model of Red Cell Transport
In [11], RBCs in capillaries are considered as separate rigid bodies immersed in blood plasma, as it is sketched in Figure 2. e RBCs are schematically shown as cylindershaped objects (their projections onto the axial section are depicted). e following formula for the effective resistance R e of a capillary, assuming that RBCs move in single-file, is proposed in [11]: where H � L/L is the value of tube hematocrit, and β � (9/9) − 1, where L � i l i is the part occupied by RBCs, 9 is the specific resistance for the parts occupied by erythrocytes, and 9 is the specific resistance for the parts filled with blood plasma. Moreover, R is the Poiseuille resistance of the vessel filled only with blood plasma, i.e., R � L9. Obviously, Note that the Obrist et al. [11] neglected the plasma layer between RBCs and the capillary wall, which yields some error. For example, the values 0.1, 0.4, and 1 of H correspond to the values 0.082, 0.33, and 0.82 of exact tube hematocrit. Nevertheless, the definition of [11] will be kept for consistency. Note that single-file flow of RBCs is typical for small vessel diameters (down to 10 μm [15]), which holds for the capillary network of the brain. Formula (7) can be transformed into the following one: which means that R e is the resistance of serially connected parts corresponding to the intervals l i and gaps between them. Note that paper [11] does not propose a precise definition of β in (7). In contrast, we are able to compute β directly using the expressions (4), (5), and (6); that is, 9 and 9 will be set as follows:

Finite Element Model of Red Cell Transport
Similar to the modeling method described in Section 2, the RBCs and blood plasma are considered as one flow with two different viscosities. As in Section 2, the viscosity of blood plasma is assumed to be μ 1 � 0.001 Pa·s, whereas the viscosity of RBCs is set to be μ 2 � 0.1 Pa·s to make RBCs effectively rigid. Moreover, it is assumed that the flow is steady state, without transition effects. erefore, the model is described by the steady state Stokes equation with spatially variable viscosity. As usually, Euler's reference system is used; that is, the flow velocity is computed at each spatial point of the unmovable simulation domain. Assuming that the flow is axisymmetric and all variables depend only on the radial and longitudinal coordinates r and z, the problem is reduced to a two-dimensional one ( Figure 3). Using the linear size and volume of erythrocytes reported in [16,17], the average length of erythrocytes was estimated to be 3 μm.
e FE method is implemented with triangle linear finite elements. e simulation domain is partitioned into 50 × 1000 rectangles in the radial and longitudinal directions, respectively, and each rectangle was divided into two triangles.
Functions v � (v 1 , v 2 ) T and q are the test ones. e viscosity distribution μ(x 1 , x 2 ) is the following discontinuous function: Figure 2: Schematic representation of a discrete stream of erythrocytes in a capillary (the axial section is shown). e RBCs are separated by blood plasma. ere is also a plasma lubrication layer between the RBCs and vessel wall.
us, the RBCs are modeled as fluid parts with a high viscosity to make them effectively rigid. e model (10)-(11) is equivalent to the following one: where ∇p � (zp/zx 1 , zp/zx 2 ). e system (13)-(14) is implemented with finite element method using the FreeFEM++ package [19]. Figure 4 shows the z-component of the velocity. e simulation is done for r c � 2.8 μm, L � 60 μm, p 0 � 0.63 mmHg, and H � 0.4. Only a part (20 μm) is presented in Figure 4 for better illustration. Figure 5 shows the dependence of specific resistance on the value of tube hematocrit for r c � 2.8 μm. An almost linear behavior holds in the range H ≤ 0.6. e dashed line, 1.28H + 0.41, is a good linear approximation for this range, which includes practically all realistic values of tube hematocrit.

Adjustment of the Formula for Capillary Flow Resistance
Fitting formula (7), or equally (8), to the numerical results is performed by increasing the blood plasma viscosity in (9) to the value μ 1 � 0.0011 Pa·s for H � 0.1 and to the value μ 1 � 0.0012 Pa·s for H � 0.4. In this case, a very good agreement of this formula with numerical results is shown in Figures 6 and 7. us, for fitting the formula (7) to the results of the numerical simulation, we need to increase the viscosity of the plasma when the hematocrit is increased. Notice that the interval from 0.1 to 0.4 covers a large part of the normal tube hematocrit range [16]. us, we can interpolate the above values of μ 1 to obtain μ 1 (H) � 10 − 3 (16 + 5H)/15. Obviously, μ 1 (0.1) � 0.0011 and μ 1 (0.4) � 0.0012. e influence of hematocrit on the specific resistance is shown in Figure 8. Here, the specific resistances computed by the Finite Element method for different values of  Finally, the effective specific resistance is given by the following formula: where ρ and ρ are defined by (9) with μ 1 � μ 1 (H) ≔ 10 − 3 (16 + 5H)/15. Simulation results presented in Figures 6 and 7 correspond to cylinder-shaped erythrocytes with the radius r 0 and length l i � 3 μm. To check the robustness of the model with respect to the change of RBC shape, simulations with a modified RBC shape were performed. e new shape, with parabolic front and back parts, is shown in Figure 9. In the simulations, l � 3 μm and h � l/3 were used. Figures 10 and 11 show the specific resistance versus capillary radius for the hematocrit values H � 0.1 and H � 0.4, respectively. e results were computed both by the FE method and analytic formula (15) with μ 1 � μ 1 (H) ≔ 10 − 3 (16 + 5H)/15. A good agreement between the results delivered by the two methods is observed.     Computational and Mathematical Methods in Medicine

Verification of the Method on Large Capillary Networks
To check the consistency of formula (15), the total resistance of the brain capillary network is computed. Similar to [7], it is assumed that capillaries of the brain are connected in a network having a variable random topology, i.e., a random number of incident edges are generated for each node. e topology is characterized by the average number of incident edges for each node so that 2-edge, 3-edge, 4-edge, etc. random topologies can be considered. Let us explain the network generation in the case of 4-edge topology. First, the set of all nodes, (i, j, k), is generated. For a node (i, j, k), the successive closest nodes, e.g., (i + 1, j, k), (i, j + 1, k), and (i, j, k + 1), are considered, and two of them are randomly chosen. ese two nodes are to be connected to (i, j, k) with the corresponding edges (capillaries). e probability that the current node (i, j, k) already has two incident edges, constructed when (i, j, k) was considered as a successive node with respect to the nodes (i − 1, j, k), (i, j − 1, k), and (i, j, k − 1), is close to 1. erefore, the network constructed in such a way has in average 4-edge topology. Additionally, it is supposed that the length and radius of capillaries are random values distributed according to data reported in [20]. Moreover, a tube hematocrit level is assigned to all capillaries, and formula (15) is applied to compute the resistance of each capillary. Similar to [7], it is supposed that the network contains blood sources and sinks (inlets and outlets) distributed over the network. ey are associated with the arteriolar and venular endpoints, respectively. e calculation of the total resistance is performed by the direct computation of the total blood flux through the capillary network by analogy with electric circuits, i.e., using Kirchhoff's law and solving a large sparse system of linear algebraic equations. e results yielded by the above sketched procedure are compared with data obtained from a model by Piechnik et al. ( [21]) fitted to experimental data of papers [22][23][24], where intravascular pressures have been measured on animals. Since the model from [21] assumes the parallel connection of all capillaries, which reduces the global hydraulic resistance, the length of capillaries has apparently been increased by factor 10 to be equal 600 μm there. is allowed the authors of [21] to fit the modeled pressure drop in the capillary compartment to the data reported in [22][23][24]. However, experimentally retrieved parameters of capillaries reported in [20] are the following: the mean length is equal to 57.4 μm and mean diameter equals 5.9 μm (M1-mosaic). We use these realistic values and, additionally, different levels of tube hematocrit in the test runs. e computation, using different random topologies, of the hydraulic resistance of the entire adult cerebral capillary network containing 756 million capillaries is presented in Figure 12. Here, the horizontal solid black line corresponds to the value of 9.9 · 10 7 (Pa·s/m 3 ) computed with the Piechnik et al. model ([21]) based on experimental data. e red stars show values of the total brain resistance calculated by the authors' method assuming 3-edge random topology and different values of tube hematocrit. e dashed line connecting the red stars demonstrates a linear dependence of the result on the hematocrit level. e blue hollow circle dots stand for 4-edge random topology. It is seen that the case of 3-edge random topology (along with the hematocrit level of 0.35) yields the best approximation of the value delivered by Piechnik's model that, generally speaking, interpolates experimental data. erefore, it is established that 3-edge random topology should be more consistent with the structure of the brain capillary network, which is in agreement with physiological data. Note that 4-edge (not quite realistic) random topology was found in [7] to be the  Figure 12: e scaled total resistance of the brain capillary network for different values of tube hematocrit: solid line-Piechnik's model immediately based on experimental data; red asterisks-the case of 3-edge random topology; blue hollow circles-the case of 4-edge random topology. e dashed lines demonstrate linear dependency of the resistance on the hematocrit level. e scale factor for the vertical axis is equal to 10 − 7 . Note that this simulation shows the sensitivity of the method against the change of the hematocrit level and random topology. 6 Computational and Mathematical Methods in Medicine best one because of ignoring the tube hematocrit level (it was always equal to 1 by default).

Computing the Pressure Distribution in the Germinal Matrix Accounting for Tube Hematocrit
As it is seen from Figure 8, the hydraulic resistance of a capillary essentially depends on the level of tube hematocrit. us, the pressure distribution in a cerebral capillary network should be significantly affected by the hematocrit level. is is especially important in the case of the subependymal germinal matrix, which is a specific part of the immature brain with high vascularity and a fragile capillary network [25], because most hemorrhages originate from this structure [26].
To calculate the pressure distribution in the germinal matrix, accounting for the hematocrit level, the approach proposed in [7] is used in the current paper with some modifications. Recall first that the total vascular system of the infant brain was described in [7] using a model proposed in [21]. is model comprises 19 levels of brain vessels: 9 levels of arterioles, 9 levels of venules, and 1 level of capillaries. All vessels at each level i, i � 1, . . . , 19, are connected in parallel. In [7], to adapt the model to infants, the numbers of vessels, their lengths, and their diameters were reduced by dividing the original values over 12 − 1.22|i − 10|, 1 + 0.14|i − 10|, and 1 + 0.11|i − 10|, respectively, which fits the CBF to a typical value of the infants brain corresponding to the age of 25 weeks [27,28]. Resistances of noncapillary vessels are given by Poiseuille's formula (6) with the apparent viscosity set to be equal to 0.003 Pa·s ( [21]). e capillary level is assumed to consist of two networks corresponding to the germinal matrix and the rest part of the brain, respectively.
In the current paper, it is supposed that these networks have 3-edge random topology, and the capillary mean length and diameter are equal to 57.4 μm and 5.9 μm, respectively ( [20]). e total resistances R GM (of the germinal matrix) and R B (of the rest part of the brain) are computed as it is outlined in Section 6, assuming that formula (15) is used for computing capillary resistances. To compute the pressure drop in the germinal matrix, the capillary level is replaced with two parallel connected lumped objects having the resistances R GM and R B , which yields the total resistance of the whole 10th level ( Figure 13): For the other levels, the resistances are calculated as where R i is the resistance of a single vessel of i th level and N i the number of vessels of i th level ( [21]). e total resistance is given as and the total flow is given by the following formula: where p A and p V are the arterial and venous (intracranial) pressures for infants. Finally, the pressure drops in the 10th level of the model, and therefore, in the germinal matrix, it is given as Remember that Δp is the difference of pressures exerted on the blood sources and sinks (inlets and outlets) distributed in the germinal matrix, and, therefore, Δp is the driving force of flow.
Parameters of a brain corresponding to infants of 25 weeks' gestational age (Table 1) are used.
Two tube hematocrit values, 0.1 and 0.4, that are close to extreme levels [16] are taken. In both cases, the resistances of the germinal matrix and the rest part of the brain, the pressure drop in the germinal matrix, and the cerebral blood  Figure 13: A vascular brain model used in [7] for the calculation of the pressure drop in the germinal matrix. e germinal matrix and the rest part of brain are considered as lumped objects. e model is a modification of that developed in [21].
Computational and Mathematical Methods in Medicine 7 flow are presented in Table 2. It is seen that the pressure drop in the germinal matrix varies by factor 1.6 if the hematocrit level increases from 0.1 to 0.4. e corresponding pressure drop distributions in a longitudinal cross-section of the germinal matrix are shown in Figure 14. e elliptic shape is chosen for better visualization. e gradient is caused by a heterogeneous distribution of inlets and outlets in the germinal matrix (see the middle of Section 4 of [7]). An essential difference in the pressure drop distributions is clearly seen there.

Conclusion
When creating a model of cerebral capillary network, it is necessary to randomly assign diameter and length to almost billion of capillaries and to compute the RBC flow resistance for each capillary. erefore, it is important to have a simple analytical formula for the resistance accounting for the hematocrit level. In the current paper, such a formula is proposed and numerically verified through finite element simulations. Accounting for hematocrit will justify and enhance models of cerebral capillary networks, allowing us to study the danger of vessel rupture in the germinal matrix, dependence on the hematocrit level.
Bearing in mind that the main goal of our study is simulation of large capillary networks, we neglect some fine effects. For example, we assume that the hematocrit level is constant over the capillary network, the plasma lubrication layer between the erythrocyte flow and the capillary wall depends on the capillary diameter only, the effect of bifurcation at network nodes is dropped because of approximate homogeneity of capillaries, and the effect of ESL (endothelial surface layer) is not accounted for. We believe that the effect of such simplifications is not damaging in context of modeling large capillary networks. Our future intentions include enhancements of the current model, in Table 1: Average data corresponding to the age of 25 gestational weeks.

Parameter
Name Value Source N A Number of capillaries of the adult brain 756 · 10 6 [21] w A Weight of the adult brain 1.2 kg [21] w Weight of the infant brain of 25 gestational weeks 0.1 kg [29] p GM Weight part of the germinal matrix 0.023 [30] w GM Weight of the germinal matrix for infants of 25 gestational weeks, w GM ≔ p GM w 2.3 g -   particular through including the bifurcation and ESL effects. We are planning the investigation of fluid interaction with ESL using homogenization theory for partial differential equations.
Data Availability e data used can be found in the references cited in this paper.

Conflicts of Interest
e authors declare that they have no conflicts of interest.