Abstract
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.
1. 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. The latter is characterized by the apparent blood viscosity which depends on the vessel diameter as well as the discharge and tube hematocrit. The discharge hematocrit is the volume fraction of the red blood cells (RBCs) in the blood delivered by the flow in the vessel. The tube hematocrit is the volume fraction of RBCs that are inside the vessel at a given time instant. The 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]. The velocity profile in the radial direction is affected by the presence of the endothelial surface layer (ESL) [1].
The 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]. This approach takes into account the heterogeneity of blood flow and partitioning of red cells at bifurcations.
The current paper is related to modeling of computer-generated blood microvessel networks with vessel diameters less than 10 μm. This 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]. The 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.
The 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.
The second approach relies on discrete models [11, 12], where the effect of each erythrocyte is taken into account. The 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. This 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. The 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. This formula accounts for gaps between RBCs and, therefore, reflects the dependence on tube hematocrit. The 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. The 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.
2. 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. Thus, the resistance of a capillary was assumed to be independent on the hematocrit in a first approximation.
The viscosity as function of the vessel’s radius was chosen as follows:where is the radius of the capillary, is the radius of the RBC substance, and and are the viscosities of blood plasma and the RBC substance, respectively. Typically, and . The last high value is used here to make the RBC column effectively rigid. It should be noted that further increasing has practically no effect.
The velocity profile corresponding to (1) was derived in [7, 14], under some assumptions, as follows:where and is the pressure drop at the capillary of the length L. The 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 and was proposed:and arguments for the reliability of this formula were presented.
Notice that formula (4) with (i.e., the whole capillary is filled with blood plasma) transforms into the Poiseuille formula:
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.
3. 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. The RBCs are schematically shown as cylinder-shaped objects (their projections onto the axial section are depicted).
The following formula for the effective resistance of a capillary, assuming that RBCs move in single-file, is proposed in [11]:where is the value of tube hematocrit, and , where is the part occupied by RBCs, is the specific resistance for the parts occupied by erythrocytes, and 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., . Obviously, if , and if . 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 is the resistance of serially connected parts corresponding to the intervals 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, and will be set as follows:
4. 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 , whereas the viscosity of RBCs is set to be to make RBCs effectively rigid. Moreover, it is assumed that the flow is steady state, without transition effects. Therefore, 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 . The FE method is implemented with triangle linear finite elements. The simulation domain is partitioned into 50 1000 rectangles in the radial and longitudinal directions, respectively, and each rectangle was divided into two triangles.
Let and be radial and longitudinal flow velocities, respectively, p the pressure, and .
The model is mathematically formulated in [18] in a weak form, which allows for using spatially discontinuous viscosity functions. With , , , , , , and , the weak formulation reads in cylindrical coordinates as follows:where and .
Functions and q are the test ones. The viscosity distribution is the following discontinuous function:
Thus, the RBCs are modeled as fluid parts with a high viscosity to make them effectively rigid.
The model (10)-(11) is equivalent to the following one:where .
The system (13)-(14) is implemented with finite element method using the FreeFEM++ package [19]. Figure 4 shows the z-component of the velocity. The simulation is done for , , , and . Only a part () is presented in Figure 4 for better illustration.
Figure 5 shows the dependence of specific resistance on the value of tube hematocrit for . An almost linear behavior holds in the range . The dashed line, , is a good linear approximation for this range, which includes practically all realistic values of tube hematocrit.
5. 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 for and to the value for . In this case, a very good agreement of this formula with numerical results is shown in Figures 6 and 7. Thus, 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]. Thus, we can interpolate the above values of to obtain . Obviously, and .
The 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 hematocrit (, , and ) are presented. The value corresponds to the model used in [7], where the gapless flow of RBCs was assumed. For the tube hematocrit values of 0.1 and 0.4, the corresponding specific resistances approximately differ by factor 2, which points out to the importance of incorporating tube hematocrit in models of capillary blood flow.
Finally, the effective specific resistance is given by the following formula:where and are defined by (9) with .
Simulation results presented in Figures 6 and 7 correspond to cylinder-shaped erythrocytes with the radius and length . To check the robustness of the model with respect to the change of RBC shape, simulations with a modified RBC shape were performed. The new shape, with parabolic front and back parts, is shown in Figure 9. In the simulations, and were used. Figures 10 and 11 show the specific resistance versus capillary radius for the hematocrit values and , respectively. The results were computed both by the FE method and analytic formula (15) with . A good agreement between the results delivered by the two methods is observed.
6. 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. The 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, , is generated. For a node , the successive closest nodes, e.g., , , and , are considered, and two of them are randomly chosen. These two nodes are to be connected to with the corresponding edges (capillaries). The probability that the current node already has two incident edges, constructed when was considered as a successive node with respect to the nodes , , and , is close to 1. Therefore, 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. They are associated with the arteriolar and venular endpoints, respectively. The 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.
The 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–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. This allowed the authors of [21] to fit the modeled pressure drop in the capillary compartment to the data reported in [22–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.
The 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 computed with the Piechnik et al. model ([21]) based on experimental data. The 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. The dashed line connecting the red stars demonstrates a linear dependence of the result on the hematocrit level. The 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. Therefore, 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 best one because of ignoring the tube hematocrit level (it was always equal to 1 by default).
7. 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. Thus, the pressure distribution in a cerebral capillary network should be significantly affected by the hematocrit level. This 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]. This 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, , 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 , , and , 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 ([21]). The 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]). The total resistances (of the germinal matrix) and (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 and , which yields the total resistance of the whole 10th level (Figure 13):
For the other levels, the resistances are calculated aswhere is the resistance of a single vessel of i th level and the number of vessels of i th level ([21]). The total resistance is given asand the total flow is given by the following formula:where and 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 is the difference of pressures exerted on the blood sources and sinks (inlets and outlets) distributed in the germinal matrix, and, therefore, 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 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. The corresponding pressure drop distributions in a longitudinal cross-section of the germinal matrix are shown in Figure 14. The elliptic shape is chosen for better visualization. The gradient is caused by a heterogeneous distribution of inlets and outlets in the germinal matrix (see the middle of Section 4 of [7]).
(a)
(b)
An essential difference in the pressure drop distributions is clearly seen there.
8. 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. Therefore, 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 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
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
This work was supported by the Klaus Tschira Foundation (Grant 00.302.2016), Buhl-Strohmaier Foundation, and Würth Foundation.