1. Introduction {#sec1}
===============

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 \[[@B1], [@B2]\]. The velocity profile in the radial direction is affected by the presence of the endothelial surface layer (ESL) \[[@B1]\].

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 \[[@B3]\]. 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 \[[@B4]\]. Moreover, the effect of hematocrit can be investigated in artificial microvascular branching networks \[[@B5]\]. 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 \[[@B6]\]. 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 \[[@B7]\]. 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 \[[@B8]\]. 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 \[[@B9]\].

Two approaches to mathematical modeling of RBC transport through microvessels can be mentioned.

The first one is based on continuum models \[[@B7], [@B10]\], 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 \[[@B7]\], which is very important for the simulation of blood flow in large capillary networks. A model derived in \[[@B10]\] studies the relations between the tube hematocrit level, vessel diameter, and apparent viscosity.

The second approach relies on discrete models \[[@B11], [@B12]\], 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 \[[@B11]\] 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 \[[@B12]\] at the scale allowing to take into account the size and shape of individual RBCs as well as their deformation. The proposed approach \[[@B11], [@B12]\] takes into account the level of hematocrit in simulations of blood flow.

In the present paper, the approaches described in \[[@B7], [@B11]\] 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 \[[@B7]\] 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 {#sec2}
=========================================

In \[[@B7]\], 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](#fig1){ref-type="fig"}.)

Such an approach is motivated by the results of \[[@B13]\] 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:$$\begin{matrix}
{\mu\left( r \right) = \begin{cases}
{\mu_{1},} & {{if}\, r_{0} \leq r \leq r_{\text{c}},} \\
{\mu_{2},} & {{if}\, r < r_{0},} \\
\end{cases}} \\
\end{matrix}$$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. The 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.

The velocity profile corresponding to ([1](#EEq1){ref-type="disp-formula"}) was derived in \[[@B7], [@B14]\], under some assumptions, as follows:$$\begin{matrix}
{v\left( r \right) = \begin{cases}
{\frac{A\left( {r^{2} - r_{\text{c}}^{2}} \right)}{4\mu_{1}},} & {{if}\, r \geq r_{0},} \\
{\frac{A\left( {r_{0}^{2} - r_{\text{c}}^{2}} \right)}{4\mu_{1}} + \frac{A\left( {r^{2} - r_{0}^{2}} \right)}{4\mu_{2}},} & {{if}\, r < r_{0},} \\
\end{cases}} \\
\end{matrix}$$where *A*=Δ*p*/*L* and Δ*p* is the pressure drop at the capillary of the length *L*. The total flux is given by the following formula:$$\begin{matrix}
{q = - \frac{\pi\Delta p}{8L}\left( {\frac{r_{\text{c}}^{4} - r_{0}^{4}}{\mu_{1}} + \frac{r_{0}^{4}}{\mu_{2}}} \right),} \\
\end{matrix}$$and the resistance reads as follows:$$\begin{matrix}
{\hat{R} = \frac{8L}{\pi}\left( {\frac{r_{\text{c}}^{4} - r_{0}^{4}}{\mu_{1}} + \frac{r_{0}^{4}}{\mu_{2}}} \right)^{- 1}.} \\
\end{matrix}$$

In \[[@B7]\] (see the end of [Section 2](#sec2){ref-type="sec"} and references there), the following relation between *r*~0~ and *r*~c~ was proposed:$$\begin{matrix}
{r_{0} = 0.3\,\mu\text{m} + 0.8r_{\text{c}},} \\
\end{matrix}$$and arguments for the reliability of this formula were presented.

Notice that formula ([4](#EEq4){ref-type="disp-formula"}) with *r*~0~=0 (i.e., the whole capillary is filled with blood plasma) transforms into the Poiseuille formula:$$\begin{matrix}
{R = \frac{8L}{\pi}\left( \frac{r_{\text{c}}^{4}}{\mu_{1}} \right)^{- 1}.} \\
\end{matrix}$$

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](#sec3){ref-type="sec"} considers an extended model, assuming the presence of gaps between RBCs, which means accounting for tube hematocrit.

3. Discrete Model of Red Cell Transport {#sec3}
=======================================

In \[[@B11]\], RBCs in capillaries are considered as separate rigid bodies immersed in blood plasma, as it is sketched in [Figure 2](#fig2){ref-type="fig"}. The RBCs are schematically shown as cylinder-shaped objects (their projections onto the axial section are depicted).

The following formula for the effective resistance *R*^e^ of a capillary, assuming that RBCs move in single-file, is proposed in \[[@B11]\]:$$\begin{matrix}
{R^{\text{e}} = R\left( {1 + H\beta} \right).} \\
\end{matrix}$$where $H = {\hat{L}/L}$ is the value of tube hematocrit, and $\beta = \left( {\hat{\varrho}/\varrho} \right) - 1$, where $\hat{L} = {\sum_{i}l_{i}}$ is the part occupied by RBCs, $\hat{\varrho}$ 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., *R*=*Lϱ*. Obviously, *R*^e^=*R* if *H*=0, and $R^{\text{e}} = \hat{R}\, ≔ \,\hat{\varrho}L$ if *H*=1. Note that the Obrist et al. \[[@B11]\] 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 \[[@B11]\] will be kept for consistency. Note that single-file flow of RBCs is typical for small vessel diameters (down to 10 *μ*m \[[@B15]\]), which holds for the capillary network of the brain.

Formula ([7](#EEq7){ref-type="disp-formula"}) can be transformed into the following one:$$\begin{matrix}
{R^{\text{e}} = \hat{L}\hat{\varrho} + \left( {L - \hat{L}} \right)\varrho,} \\
\end{matrix}$$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 \[[@B11]\] does not propose a precise definition of *β* in ([7](#EEq7){ref-type="disp-formula"}). In contrast, we are able to compute *β* directly using the expressions ([4](#EEq4){ref-type="disp-formula"}), ([5](#EEq5){ref-type="disp-formula"}), and ([6](#EEq6){ref-type="disp-formula"}); that is, $\hat{\varrho}$ and *ϱ* will be set as follows:$$\begin{matrix}
{\hat{\varrho} = \frac{8}{\pi}\left( {\frac{r_{\text{c}}^{4} - r_{0}^{4}}{\mu_{1}} + \frac{r_{0}^{4}}{\mu_{2}}} \right)^{- 1},} \\
{\varrho = \frac{8}{\pi}\left( \frac{r_{\text{c}}^{4}}{\mu_{1}} \right)^{- 1},} \\
{r_{0} = 0.3\,\mu\text{m} + 0.8r_{\text{c}}.} \\
\end{matrix}$$

4. Finite Element Model of Red Cell Transport {#sec4}
=============================================

Similar to the modeling method described in [Section 2](#sec2){ref-type="sec"}, the RBCs and blood plasma are considered as one flow with two different viscosities. As in [Section 2](#sec2){ref-type="sec"}, 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. 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](#fig3){ref-type="fig"}). Using the linear size and volume of erythrocytes reported in \[[@B16], [@B17]\], the average length of erythrocytes was estimated to be 3 *μ*m. 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 *u*~*r*~ and *u*~*z*~ be radial and longitudinal flow velocities, respectively, *p* the pressure, and *Ω*=(0, *r*~c~) × (0, *L*).

The model is mathematically formulated in \[[@B18]\] in a weak form, which allows for using spatially discontinuous viscosity functions. With *x*~1~=*r*, *x*~2~=*z*, *u*~1~=*u*~*r*~, *u*~2~=*u*~*z*~, *u*=(*u*~1~, *u*~2~)^*T*^, *p*(*r*, 0)=*p*~0~, and *p*(*r*, *L*)=0, the weak formulation reads in cylindrical coordinates as follows:$$\begin{matrix}
{{\int_{\Omega}{x_{1}\left( {2\mu\left( {x_{1},x_{2}} \right){\sum\limits_{i,j = 1}^{2}{D_{ij}\left( u \right)D_{ij}\left( v \right)}} + \frac{u_{1}v_{1}}{x_{1}^{2}}} \right)dx}} - {\int_{\Omega}{x_{1}p\text{ div}\left( v \right)dx}} = {\int_{\Gamma_{0}}{x_{1}p_{0}v_{2}\, dx}},} \\
\end{matrix}$$$$\begin{matrix}
{\varepsilon{\int_{\Omega}{x_{1}pq\, dx}} - {\int_{\Omega}{x_{1}\text{ div}\left( u \right)q\, dx = 0}},} \\
 \\
{\left. u \right|_{\Gamma_{2}} = 0,} \\
{\left. v \right|_{\Gamma_{2}} = 0,} \\
\end{matrix}$$where *D*~*ij*~(*u*)=1/2((∂*u*~*i*~/∂*x*~*j*~)+(∂*u*~*j*~/∂*x*~*i*~)) and div(*u*)=(*u*~1~/*x*~1~)+(∂*u*~1~/∂*x*~1~)+(∂*u*~2~/∂*x*~2~).

Functions *v*=(*v*~1~, *v*~2~)^*T*^ and *q* are the test ones. The viscosity distribution *μ*(*x*~1~, *x*~2~) is the following discontinuous function:$$\begin{matrix}
{\mu\left( {x_{1},x_{2}} \right) = \begin{cases}
{\mu_{1} = 0.001\text{ Pa} \cdot \text{s},} & {{if}\left( {x_{1},x_{2}} \right) \in \text{plasma part},} \\
{\mu_{2} = 0.1\text{ Pa} \cdot \text{s},} & {{if}\left( {x_{1},x_{2}} \right) \in \text{RBCs  part}.} \\
\end{cases}} \\
\end{matrix}$$

Thus, the RBCs are modeled as fluid parts with a high viscosity to make them effectively rigid.

The model ([10](#EEq10){ref-type="disp-formula"})-([11](#EEq11){ref-type="disp-formula"}) is equivalent to the following one:$$$$$$$$where ∇*p*=(∂*p*/∂*x*~1~, ∂*p*/∂*x*~2~).

The system ([13](#EEq13){ref-type="disp-formula"})-([14](#EEq14){ref-type="disp-formula"}) is implemented with finite element method using the FreeFEM++ package \[[@B19]\]. [Figure 4](#fig4){ref-type="fig"} shows the *z*-component of the velocity. The 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](#fig4){ref-type="fig"} for better illustration.

[Figure 5](#fig5){ref-type="fig"} 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. The dashed line, 1.28*H*+0.41, 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 {#sec5}
==========================================================

Fitting formula ([7](#EEq7){ref-type="disp-formula"}), or equally ([8](#EEq8){ref-type="disp-formula"}), to the numerical results is performed by increasing the blood plasma viscosity in ([9](#EEq9){ref-type="disp-formula"}) 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](#fig6){ref-type="fig"} and [7](#fig7){ref-type="fig"}. Thus, for fitting the formula ([7](#EEq7){ref-type="disp-formula"}) 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 \[[@B16]\]. Thus, we can interpolate the above values of *μ*~1~ to obtain *μ*~1~(*H*)=10^−3^(16+5*H*)/15. Obviously, *μ*~1~(0.1)=0.0011 and *μ*~1~(0.4)=0.0012.

The influence of hematocrit on the specific resistance is shown in [Figure 8](#fig8){ref-type="fig"}. Here, the specific resistances computed by the Finite Element method for different values of hematocrit (*H*=1, *H*=0.4, and *H*=0.1) are presented. The value *H*=1 corresponds to the model used in \[[@B7]\], 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:$$\begin{matrix}
{r^{\text{e}} = \rho + \left( {\hat{\rho} - \rho} \right)H,} \\
\end{matrix}$$where *ρ* and $\hat{\rho}$ are defined by ([9](#EEq9){ref-type="disp-formula"}) with *μ*~1~=*μ*~1~(*H*)≔10^−3^(16+5*H*)/15.

Simulation results presented in Figures [6](#fig6){ref-type="fig"} and [7](#fig7){ref-type="fig"} 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. The new shape, with parabolic front and back parts, is shown in [Figure 9](#fig9){ref-type="fig"}. In the simulations, *l*=3 *μ*m and *h*=*l*/3 were used. Figures [10](#fig10){ref-type="fig"} and [11](#fig11){ref-type="fig"} show the specific resistance versus capillary radius for the hematocrit values *H*=0.1 and *H*=0.4, respectively. The results were computed both by the FE method and analytic formula ([15](#EEq15){ref-type="disp-formula"}) with *μ*~1~=*μ*~1~(*H*)≔10^−3^(16+5*H*)/15. A good agreement between the results delivered by the two methods is observed.

6. Verification of the Method on Large Capillary Networks {#sec6}
=========================================================

To check the consistency of formula ([15](#EEq15){ref-type="disp-formula"}), the total resistance of the brain capillary network is computed. Similar to \[[@B7]\], 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, (*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. These two nodes are to be connected to (*i*, *j*, *k*) with the corresponding edges (capillaries). The 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. 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 \[[@B20]\]. Moreover, a tube hematocrit level is assigned to all capillaries, and formula ([15](#EEq15){ref-type="disp-formula"}) is applied to compute the resistance of each capillary. Similar to \[[@B7]\], 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. (\[[@B21]\]) fitted to experimental data of papers \[[@B22]--[@B24]\], where intravascular pressures have been measured on animals. Since the model from \[[@B21]\] 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 \[[@B21]\] to fit the modeled pressure drop in the capillary compartment to the data reported in \[[@B22]--[@B24]\]. However, experimentally retrieved parameters of capillaries reported in \[[@B20]\] 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](#fig12){ref-type="fig"}. 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 (\[[@B21]\]) 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 \[[@B7]\] 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 {#sec7}
============================================================================================

As it is seen from [Figure 8](#fig8){ref-type="fig"}, 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 \[[@B25]\], because most hemorrhages originate from this structure \[[@B26]\].

To calculate the pressure distribution in the germinal matrix, accounting for the hematocrit level, the approach proposed in \[[@B7]\] is used in the current paper with some modifications. Recall first that the total vascular system of the infant brain was described in \[[@B7]\] using a model proposed in \[[@B21]\]. 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*, *i*=1,..., 19, are connected in parallel. In \[[@B7]\], 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 \[[@B27], [@B28]\]. Resistances of noncapillary vessels are given by Poiseuille\'s formula ([6](#EEq6){ref-type="disp-formula"}) with the apparent viscosity set to be equal to 0.003 Pa · s (\[[@B21]\]). 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 (\[[@B20]\]). The 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](#sec6){ref-type="sec"}, assuming that formula ([15](#EEq15){ref-type="disp-formula"}) 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](#fig13){ref-type="fig"}):$$\begin{matrix}
{R_{10} = \left( {R_{\text{GM}}^{- 1} + R_{\text{B}}^{- 1}} \right)^{- 1}.} \\
\end{matrix}$$

For the other levels, the resistances are calculated as$$\begin{matrix}
{R_{i} = \frac{{\overline{R}}_{i}}{N_{i}},\quad i = \overline{1,19},\, i \neq 10,} \\
\end{matrix}$$where ${\overline{R}}_{i}$ is the resistance of a single vessel of *i* th level and *N*~*i*~ the number of vessels of *i* th level (\[[@B21]\]). The total resistance is given as$$\begin{matrix}
{R_{\text{T}} = {\sum\limits_{i = 1}^{19}R_{i}},} \\
\end{matrix}$$and the total flow is given by the following formula:$$\begin{matrix}
{Q = \frac{\left( {p_{\text{A}} - p_{\text{V}}} \right)}{R_{\text{T}}},} \\
\end{matrix}$$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$$\begin{matrix}
{\Delta p = QR_{10}.} \\
\end{matrix}$$

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](#tab1){ref-type="table"}) are used.

Two tube hematocrit values, 0.1 and 0.4, that are close to extreme levels \[[@B16]\] 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](#tab2){ref-type="table"}. 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](#fig14){ref-type="fig"}. 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](#sec4){ref-type="sec"} of \[[@B7]\]).

An essential difference in the pressure drop distributions is clearly seen there.

8. Conclusion {#sec8}
=============

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.

This work was supported by the Klaus Tschira Foundation (Grant 00.302.2016), Buhl-Strohmaier Foundation, and Würth Foundation.

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.

![Schematic representation of continuum stream of erythrocytes in a capillary (the axial section is shown). The gaps between erythrocytes are ignored. The motion occurs due to the lubrication plasma layer between erythrocytes and vessel wall.](CMMM2019-4235937.001){#fig1}

![Schematic representation of a discrete stream of erythrocytes in a capillary (the axial section is shown). The RBCs are separated by blood plasma. There is also a plasma lubrication layer between the RBCs and vessel wall.](CMMM2019-4235937.002){#fig2}

![Schematic drawing of the FE domain in cylindrical coordinates; only three RBCs are shown. The RBC subdomain is shown in pink, and the plasma subdomain is colorless.](CMMM2019-4235937.003){#fig3}

![The *z*-component of flow velocity (mm/s). The velocity is smaller in the parts occupied by red cells.](CMMM2019-4235937.004){#fig4}

![The dependence of specific resistance (red dots) on the value of tube hematocrit for *r*~c~=2.8 *μ*m, computed with the FM model. The dashed line is a good linear approximation of this dependence for *H* ≤ 0.6. The scale factor for the vertical axis is equal to 10^−20^.](CMMM2019-4235937.005){#fig5}

![Comparison between analytically (solid line) and FEM-computed (dots) specific resistances in the case of *H*=0.1. The scale factor for the vertical axis is equal to 10^−20^.](CMMM2019-4235937.006){#fig6}

![Comparison between analytically (solid line) and FEM-computed (dots) specific resistances in the case of *H*=0.4. The scale factor for the vertical axis is equal to 10^−20^.](CMMM2019-4235937.007){#fig7}

![Specific resistances computed with FEM for different values of the hematocrit: *H*=1 (solid red), *H*=0.4 (dashed blue), and *H*=0.1 (dashed-dot green). The scale factor for the vertical axis is equal to 10^−20^.](CMMM2019-4235937.008){#fig8}

![Noncylindrical (deformed) RBCs.](CMMM2019-4235937.009){#fig9}

![Comparison between analytically (solid line) and FEM-computed (dots) specific resistances for noncylindrical (deformed) RBCs in the case *H*=0.1. The scale factor for the vertical axis is equal to 10^−20^.](CMMM2019-4235937.010){#fig10}

![Comparison between analytically (solid line) and FEM-computed (dots) specific resistances for noncylindrical (deformed) RBCs in the case *H*=0.4. The scale factor for the vertical axis is equal to 10^−20^.](CMMM2019-4235937.011){#fig11}

![The 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. The dashed lines demonstrate linear dependency of the resistance on the hematocrit level. The 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.](CMMM2019-4235937.012){#fig12}

![A vascular brain model used in \[[@B7]\] for the calculation of the pressure drop in the germinal matrix. The germinal matrix and the rest part of brain are considered as lumped objects. The model is a modification of that developed in \[[@B21]\].](CMMM2019-4235937.013){#fig13}

![Pressure drop distributions in the longitudinal cross-section of the germinal matrix: (a)  *H*=0.1; (b)  *H*=0.4.](CMMM2019-4235937.014){#fig14}

###### 

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^     \[[@B21]\]
  *w* ~A~          Weight of the adult brain                                                                                          1.2 kg          \[[@B21]\]
  *w*              Weight of the infant brain of 25 gestational weeks                                                                 0.1 kg          \[[@B29]\]
  *p* ~GM~         Weight part of the germinal matrix                                                                                 0.023           \[[@B30]\]
  *w* ~GM~         Weight of the germinal matrix for infants of 25 gestational weeks, *w*~GM~≔*p*~GM~*w*                              2.3 g           ---
  *κ*              Capillary density correction factor for the germinal matrix                                                        1.5             \[[@B31]\]
  *N* ~GM~         Number of capillaries of the germinal matrix for infants of 25 gestational weeks, *N*~GM~≔*κw*~GM~*N*~A~/*w*~A~    2.17 · 10^6^    ---
  *N* ~B~          Number of capillaries of the rest part of the brain of 25 gestational weeks, *N*~B~≔(*w* − *w*~GM~)*N*~A~/*w*~A~   61.55 · 10^6^   ---
  *p* ~A~/*p*~V~   Arterial/venous pressure                                                                                           34/5 mmHg       \[[@B28], [@B32]\]

###### 

Dependence of computed network characteristics on the tube hematocrit level.

  Characteristic                                                  *H*=0.1       *H*=0.4
  --------------------------------------------------------------- ------------- --------------
  Resistance of the germinal matrix, *R*~GM~(Pa · s/m^3^)         150 · 10^8^   229 · 10^8^
  Resistance of the rest part of the brain, *R*~B~(Pa · s/m^3^)   8.2 · 10^8^   13.2 · 10^8^
  Pressure drop in the germinal matrix, Δ*p* (Pa)                 129.2         203.3
  CBF (ml/min)                                                    9.97          9.77

[^1]: Academic Editor: Michele Migliore
