## Abstract

A detailed model of red blood cell (RBC) transport in a capillary network is an indispensable element of a comprehensive model for the supply of the human organism with oxygen and nutrients. In this paper, we introduce a two-phase model for the perfusion of a capillary network. This model accounts for the special role of RBCs, which have a strong influence on network dynamics. Analytical results and numerical simulations with a discrete model and a generic network topology indicate that there exists a local self-regulation mechanism for the flow rates and a global de-mixing process that leads to an inhomogeneous haematocrit distribution. Based on the results from the discrete model, we formulate an efficient algorithm suitable for computing the pressure and flow field as well as a continuous haematocrit distribution in large capillary networks at steady state.

## 1. Introduction

Recent numerical studies of cerebral blood flow (Boas *et al.* 2008; Reichold *et al.* 2009) modelled the perfusion of large vascular trees as complex three-dimensional pipe networks with a viscous flow. The flow and pressure field in such networks can be simulated on the basis of Kirchhoff’s laws for electrical networks. This involves very large systems of linear equations. In order to reduce the total size of these systems, the lowest levels of the vascular tree, i.e. the capillaries, can be treated as a spatial continuum with a given permeability that is based on the local morphology of the capillary network. We will show that this approach falls short of describing the perfusion of a capillary network because it neglects the role of the red blood cells (RBCs). In capillary networks, blood must not be modelled as a single viscous phase. Rather, the transport of RBCs and the plasma phase are to be treated separately. This leads to perfusion characteristics that differ strongly from single-phase models.

The role of RBCs in the dynamics of the microcirculation has been recognized, studied and discussed in a series of papers in the 1980s and 1990s (e.g. Schmid-Schönbein *et al.* 1980; Papenfuss & Gross 1981; Klitzman & Johnson 1982; Dawant *et al.* 1986*a*,*b*; Pries *et al.* 1986, 1989, 1990, 1992; Secomb *et al.* 1986). It was shown that the dynamics in a capillary network is much more than just the ‘sum’ of viscous flows through narrow pipes. The local concentration of RBCs has a strong impact on the local flow resistance, which influences the flow and pressure fields in the whole network. Schmid-Schönbein *et al.* (1980) illustrated local self-regulation mechanisms that tend to balance the flow rates in capillaries diverging from the same parent vessel. Unfortunately, it appears that many of today’s researchers are not aware of these results. Therefore, this paper revisits the basic discrete model for RBC transport in capillary networks (as outlined by Pries *et al.* (1990)). Based on this approach, we will formulate a new continuous algorithm to solve for the pressure, flow and haematocrit distribution in a capillary network. This algorithm is efficient and well suited for today’s applications that often use large morphological datasets of vascular networks (Reichold *et al.* 2009; Tsai *et al.* 2009). We will show for some idealized configurations that the equations for the pressure field can be decoupled from the equations for the flow field and haematocrit distribution. Furthermore, we will present some principal results for capillary networks and discuss their consequences for the perfusion of these networks. In particular, we will stress the differences between the transport of substances in capillary networks with and without RBCs (two phase versus single phase).

We will introduce the theory in a general framework since the presented results apply to the wider field of particulate flows in pipe networks (or even porous media) with the only conditions that the Reynolds number is sufficiently small and that the particle size is comparable to the pipe diameter. In order to reveal the principal phenomena in a clear way, we will idealize and simplify the configuration as far as possible. A more detailed description of the local physical phenomena would not change the global dynamics in a fundamental way, but only take off some of the rough edges of the results discussed here.

## 2. Governing equations for red blood cell transport in a capillary network

We model capillary networks as sets of nodes that are connected by rigid pipes with circular cross sections. A pipe connecting node *i* with node *j* has a length *L*_{ij} and a constant diameter *D*_{ij}. We require that the diameter of the pipes is small, i.e. *L*_{ij}≫*D*_{ij}.

The local Reynolds number in a capillary with diameter *D* and bulk velocity is defined as , where *ρ* is the fluid density, *μ* the dynamic viscosity of the fluid (i.e. the plasma) and *q* the mass flow rate. In capillary networks, this Reynolds number is smaller than unity.

The network transports particles (RBCs) with the same density as the carrier fluid. The RBCs are deformable and incompressible such that the volume is preserved even when their shape changes. Their size is of the order of the capillary diameter such that there is only a narrow lubrication gap between an RBC and the vessel wall. They form a single file as they are advected through a capillary. Further, we assume that the RBCs are flexible enough to pass through each capillary in the network and there is no blockage of capillaries.

Because of the small Reynolds number, the plasma velocity field ** u** and the pressure

*p*are governed by the Stokes equation2.1and the continuity equation2.2These equations assume implicitly that the flow field is stationary or at least changing so slowly that the term

*ρ*∂

**/∂**

*u**t*of the unsteady Stokes equation is negligible with respect to the other terms in equation (2.1). A dimensional analysis of the situation in the microcirculation confirms that this is the case.

### (a) Flow without red blood cells

For a one-dimensional flow through a pipe, it follows directly from equation (2.1) that the mass flow rate *q*_{ij} is proportional to the pressure drop between the two adjacent nodes, i.e.2.3where *p*_{i} and *p*_{j} are the pressure values at the nodes *i* and *j*, respectively. The nominal resistance *R*_{ij} is given by2.4where *r*_{ij} is the specific resistance of the capillary. According to Poiseuille’s law (Sutera & Skalak 1993), it is given by2.5The linear relationship between pressure drop and mass flux drastically simplifies the flow calculations. In analogy to Kirchhoff’s first rule for electrical networks, we formulate an equation for the mass balance at each node *i*,2.6where *q*_{i} is a source term describing a flux into or out of the network at node *i* (except at the in- and outflow of the network, we set *q*_{i} typically to zero). This results in a system of linear equations. To close the system, we have to fix one node pressure by a Dirichlet boundary condition (e.g. by setting the pressure at the outflow to zero). The linear system of equations (2.6) is the discrete equivalent to Darcy’s law coupled with mass balance. Its solution yields all node pressure values *p*_{i}, and the flow rates *q*_{ij} can then be computed from equation (2.3).

### (b) Flow with red blood cells

The situation becomes more complex with the presence of RBCs. Because of the small Reynolds numbers and the quasi-stationary flow field, we may neglect the RBC inertia. This means also that the RBC velocities adjust ‘immediately’ to balance the pressure and the viscous forces.

In general, the flow resistance of a capillary will increase if RBCs are present. This increase is related to the altered flow field in the vicinity of an RBC (figure 1). Instead of a Poiseuille flow profile, we find a narrow shear flow (similar to Couette flow) in the narrow gap between an RBC and the vessel wall. This shear flow leads to an increased wall shear stress and, thus, to a higher resistance. The additional pressure drop owing to an RBC can be written as2.7where is the effective length of the section in which the velocity profile is influenced by the RBC, and is the specific resistance within that section.

The additional pressure drop is illustrated in figure 1 for a capillary with two RBCs. The total pressure drop in the pipe *i*–*j* with *N*_{ij} RBCs is then2.8where we have introduced the effective resistance2.9with the tube haematocrit2.10and2.11The variable *β*_{ij} is typically of order 1 and is known as the apparent intrinsic viscosity (Secomb *et al.* 1986). Several numerical results illustrate the dependence of *β*_{ij} on various flow parameters (Pozrikidis 2003, 2005*a*,*b*) and on the endothelial surface layer (ESL) (Secomb *et al.* 2003). In the present work, we simply assume that *β*_{ij} is a constant property of a capillary and of order 1. This idealization is legitimate because we will see in §3*b* that one of the central mechanisms in capillary networks does not depend on the exact form of equation (2.9), but only requires that the effective resistance increases *monotonically* with the haematocrit. Similar to equation (2.6), we can use the effective resistance to derive a linear system of equations2.12which yields the node pressure values in a network with RBCs.

The velocity of an RBC is given by2.13where the coefficient corrects the RBC velocity with respect to the mean velocity . Very tightly fitting RBCs will be advected with the mean fluid velocity (), whereas RBCs at the centreline of very large blood vessels (e.g. arteries) will be advected nearly with the peak fluid velocity (). The velocity correction coefficient is directly related to the Fåhræus effect in capillaries (Albrecht *et al.* 1979).

### (c) Bifurcation rule

If the RBC locations in the network are given, the pressure and flow field can be computed from equations (2.9) and (2.12). Equation (2.13) can then be used to compute the RBC velocities with which the RBCs are advected. As long as the RBCs remain within their capillary, this propagation is well defined and the pressure and flow field do not change. Once an RBC arrives at a bifurcation, it will enter a different capillary such that the effective resistance changes and the pressure and flow fields have to be recomputed.

To decide which capillary the RBC will enter, we need a *bifurcation rule*. At first glance, the RBC dynamics at a bifurcation appears to be a very complex matter that does not only depend on the local flow field but also on the exact bifurcation geometry and RBC properties. Intuitively, we would assume that an RBC prefers the capillary that requires the least change in direction. This reasoning is motivated by the inertia of the RBC. In our case, however, we have a low-Reynolds-number configuration governed by the Stokes equation (2.1) that does not contain any inertial terms. Therefore, we postulate that the RBCs will follow the path of the steepest pressure gradient (regardless of their inertia). Further, we assume that the work required to deform an RBC (e.g. to squeeze it into a capillary with a smaller diameter) is negligible (Secomb & Hsu 1996). We formalize the bifurcation rule as follows: an RBC at node *i* will enter the pipe *i*– with the highest negative pressure gradient2.14where the local pressure gradient *G*_{ij} of the pipe *i*–*j* is given by2.15Note that *G*_{ij} does not depend on the effective resistance , but rather on the specific resistance *r*_{ij}, because the RBC at the bifurcation only ‘feels’ the local pressure gradient and does not ‘see’ what is behind the bifurcation (cf. figure 1). Therefore, the bifurcation dynamics is not directly influenced by any obstacle, constriction or dilation behind the bifurcation.

Note that this formulation of the bifurcation rule is different from the bifurcation rule used in older studies (e.g. Klitzman & Johnson 1982; Pries *et al.* 1990), where the models for the bifurcation dynamics are based on empirical results (e.g. Fung 1973; Pries *et al.* 1989). Our bifurcation rule is motivated by physical arguments and follows from a strict analysis of an idealized configuration. It may neglect some of the more subtle effects present at a bifurcation (e.g. RBC deformation), but has the advantage that it can be used to derive analytical results for capillary networks (as will be shown, for instance, in §4*a*).

In idealized configurations where the local specific resistances *r*_{ij} are equal for all diverging vessels, an RBC will prefer the capillary with the highest flow rate. In this case, our bifurcation rule is consistent with the bifurcation rule used by Schmid-Schönbein *et al.* (1980). Therefore, it can be considered an idealization of the experimental findings by Fung (1973) and Pries *et al.* (1989), and it is consistent with the more recent numerical results by Barber *et al.* (2008) if only narrow vessels are considered.

## 3. Discrete model: red blood cell transport in capillary networks

The mass balance at the bifurcations, the effective resistance and the bifurcation rule constitute a complete model for the RBC transport in a capillary network. The basic simulation algorithm reads as follows:

Choose initial RBC positions.

Compute effective resistance using equation (2.9).

Set up a linear system of equations according to equation (2.12).

Solve the linear system of equations for node pressure values

*p*_{i}.Compute flow rates

*q*_{ij}from equation (2.8).Propagate the RBCs with the RBC velocities (2.13) until the first RBC reaches a bifurcation.

Apply the bifurcation rule (2.14) to assign the RBC to a new capillary.

Continue with step 2.

Schmid-Schönbein *et al.* (1980) was the first to use such an algorithm to simulate blood flow in the microcirculation. Although they studied only a small network with 12 vessels, their results revealed a self-regulation mechanism for flow rates in the vessels diverging from a bifurcation (see discussion below).

We use this algorithm in the following to study the dynamics in a generic square network. It consists of nodes with 1≤*i*_{1},*i*_{2}≤*n*, which are arranged in a square. Each node is connected to its neighbours, i.e. to the nodes , , and . The inflow is at the node with a constant flow rate *q*_{0}>0 and the outflow with −*q*_{0} at the node . The pressure at the outflow node is set to zero.

Such a network is the simplest possible configuration for studying basic mechanisms of RBC transport. Although it cannot be found in nature, it is well suited for illustrating the dramatic differences between a single- and a two-phase model.

For the following numerical experiments, a network with *n*=11 is seeded randomly with *n*_{p}=550 RBCs (uniformly distributed), which corresponds to an average haematocrit of 0.4. The network is homogeneous in the sense that all capillaries have the same physical properties (*D*=5 μm, *L*=50 μm, m, *β*=2.7, *θ*=1). The number of RBCs is held constant by reinjecting all RBCs that exit the network through the outflow. The mass flow rate at the inflow is set to *q*_{0}=*ρπD*^{2}/4×10^{−3} ms^{−1}. The fluid (i.e. the plasma) is assumed to have density *ρ*=1050 kg m^{−3} and dynamic viscosity *μ*=1.5×10^{−3} kg ms^{−1}.

### (a) Transition towards a steady-state configuration

Figure 2 shows the network with the RBCs at different times. Starting with the uniformly distributed initial condition (figure 2*a*), the network goes through a transient phase (figure 2*b*,*c*), during which the uniform haematocrit distribution becomes increasingly inhomogeneous. Eventually, this de-mixing of the RBC and plasma phases reaches a statistically steady state (figure 2*d*) in the sense that the ensemble averages of the flow field, the pressure field and the haematocrit distribution are constant in time. The duration of the transient phase is comparable to the total network turnover time *τ*=11 *s*, which is defined as the ratio of the total network volume to the volumetric flow rate *q*_{0}/*ρ* at the inflow. The particular steady-state configuration is a direct consequence of the bifurcation rule and of the effective resistance that depends on the local haematocrit. For the given network topology, the steady-state configuration exhibits a low haematocrit and a relatively high flow rate in the centre of the network. An ‘RBC highway’ emerges along the main diagonal, while the RBCs in the outer corners accumulate and stagnate.

In order to illustrate the relevance of the bifurcation rule and the effective resistance, we show in figure 3 four different steady-state solutions. Only figure 3*a* has been obtained with the full model. In the three other plots, the effective resistance and/or the bifurcation rule have been replaced by the nominal resistance *R*_{ij} and/or by a random bifurcation (biased by the local flow rate). Only the full model exhibits a low particle concentration and a relatively high flow rate in the centre of the network. If the bifurcation rule is switched off (figure 3*b*,*d*), the particles maintain their uniform distribution across the whole network. The flow rates in figure 3*c*,*d* are identical because the effective resistance is replaced by the nominal resistance such that the flow field is independent of the particles (single-phase flow). In figure 3*c*, the RBCs follow the path of the steepest pressure gradient. In figure 3*d*, the RBCs act as passive tracers that are carried along with the fluid without influencing the flow.

### (b) Self-regulation mechanism

The steady-state configuration in a network is the global result of a local phenomenon, i.e. the self-regulation mechanism. Schmid-Schönbein *et al.* (1980) described this mechanism for a bifurcation with two identical diverging vessels. Here, we discuss the mechanism in the context of our more general bifurcation rule. According to this rule, an RBC arriving at a bifurcation moves into the capillaries with the highest pressure gradient. This increases the effective resistance and thus reduces the flow rate and the pressure gradient. This process continues until another capillary has the highest pressure gradient. On average, an equilibrium is attained where all capillaries (with *q*_{ij}>0) have the same pressure gradient. For capillaries with equal specific resistances *r*_{ij}, this process is tantamount to a self-regulation mechanism that leads to equal flow rates in all divergent capillaries, as reported by Schmid-Schönbein *et al.* (1980). Their results show also that the self-regulation mechanism is able to adjust the flow rates very rapidly, i.e. the time necessary to attain a statistical equilibrium is of the order of the network turnover time.

There are cases, however, where the self-regulation mechanism breaks down. A simple example for such a situation is the case of two parallel capillaries connecting the same two nodes. If we assume that one of these capillaries has a nominal resistance *R*_{1}, which is substantially higher than the other resistance *R*_{2}, then there exists a minimal number *N* of RBCs necessary to sustain the self-regulation mechanism, i.e.3.1where *L*_{2}, and *β*_{2} characterize the capillary with the lower nominal resistance. If there are fewer particles, the resistance *R*_{2} cannot be increased sufficiently, i.e. such that reaches the level of *R*_{1}. Nevertheless, the system will reach a steady state in which all RBCs are in the capillary with the lower nominal resistance and the other capillary contains no RBCs at all. For the remainder of this paper, we neglect this situation and assume that there are always sufficiently many RBCs.

### (c) Unsteady red blood cell transport and transit times

Apart from the self-regulation mechanism, we can observe several phenomena that have also been observed *in vivo*. Figure 4 illustrates the RBC transport in a space–time diagram of RBC positions in five different capillaries of our generic network. Note that these diagrams show a striking similarity to fig. 3*a*–*c* and fig. 4*a*–*d* of Kleinfeld *et al.* (1998), who show space–time diagrams obtained with two-photon laser scanning microscopy in a rat neocortex.

Although the network is statistically in steady state, figure 4 shows that the motion of the RBCs is neither homogeneous nor uniform. Apparently, most RBCs prefer to travel in groups, and the velocities of the RBCs are not constant. It has been suspected (e.g. Lawrence & Clapper 1972) that the flow ‘may pulsate with the heart beat’. Here, we see that unsteady behaviour (e.g. pulsating velocity) can be caused by the nonlinear dynamics of a capillary network. The same authors and others (e.g. Kleinfeld *et al.* 1998) also report cases of flow reversal, i.e. situations where the RBCs change direction intermittently. This phenomenon can also be observed with our model (see diagram 4 in figure 4 at *t*≈15 s), and is directly related to the regulation of the flow rates by the self-regulation mechanism and the fluctuation of the flow rates around the regulated average flow rate. In capillaries with a very small average flow rate, this statistical fluctuation may lead to intermittent changes of the sign of the flow velocity. Finally, we would like to point out that there are also capillaries in our capillary network that (at times) do not contain any RBCs at all. This inhomogeneous haematocrit distribution (cf. fig. 6 in Secomb *et al.* 1995) leads to the *network Fåhræus effect* (Pries *et al.* 1992).

The results shown in figure 4 illustrate the local dynamics of RBC transport. To study the transport on a global level, we measure the transit time *T* of the RBCs, i.e. the time taken by an RBC to pass through the whole network. We observe that some RBCs pass rapidly along the main ‘highway’ towards the outlet, while many RBCs take a slower route via the periphery of the network. This is reflected in the histogram of the transit times (figure 5). For comparison, we also show the transit times of passive tracer particles that are equivalent to a single-phase model. While both histograms peak at around *T*=*τ*, we observe that the transit time varies less for passive tracer particles than for RBCs (long tail in the histogram). At the same time, the fastest RBCs traverse the network slightly faster than the fastest passive tracers. It is also worthwhile to note that the two peaks around *T*/*τ*=5 and 6.5 appear to be genuine in the sense that the statistics have converged. The preference for these two transit times is probably related to the specific network topology, but is not fully understood at the moment.

## 4. Continuous model: dynamics in steady state

We have seen in §3*b* that the self-regulation mechanism regulates the pressure gradients in the diverging blood vessels of a bifurcation. This local mechanism leads to steady-state configurations of the whole network.

To make further progress in investigating the dynamics of steady-state configurations, we assume that and . It is then a valid idealization to consider the haematocrit *H*_{ij}≥0 to be a continuous real number. In contrast to the previous section, where we tracked discrete RBCs in a transient process, we will now study the dynamics of a stationary system with continuous variables. We use the haematocrit to formulate an equation for the conservation of RBCs within the network, i.e. we require that the sum of all RBC mass flow rates at a bifurcation is zero. This is similar to the application of Kirchhoff’s first rule in equation (2.6), except that we are now directly involving the haematocrit in our equations. Further, we formulate a reordering approach for determining the flow rates. This will allow us to determine the steady-state solution in an iterative process. A similar scheme was used by Papenfuss & Gross (1981) and Dawant *et al.* (1986*a*,*b*) to study the haematocrit distribution in microvascular networks. Here, we will formulate an algorithm that makes use of our specific formulation of the bifurcation rule. This algorithm is able to solve the governing equations in an efficient manner such that it is well suited for studying very large networks. We will use this algorithm to illustrate some basic phenomena. Furthermore, we will derive some analytical results for idealized network configurations.

### (a) Red blood cell mass balance

The RBC mass balance for node *i* reads4.1where is the rate at which RBCs flow into or out of the network at node *i* (typically zero, except at the in- and outflow). The RBC flow rate in the pipe *i*–*j* is defined as4.2We can use equations (2.8) and (2.9) to express the haematocrit by4.3such that the RBC flow rate is given by4.4With this expression, the RBC mass balance becomes4.5Note that equation (4.5) relates the flow rates *q*_{ij} to the node pressure values *p*_{i} and does not contain any terms that depend explicitly on the haematocrit *H*_{ij}. To determine the statistically stationary haematocrit distribution, an additional equation is required.

### (b) Reordering approach

The self-regulation mechanism ensures that all pressure gradients *G*_{ij} (for the capillaries with *q*_{ij}>0) attain the same value *G*_{i} in steady state. If all node pressure values were known, we could calculate the RBC distribution by a reordering approach. To this end, we start at the node *i* with the highest pressure value (where the inflow rate *q*_{i} is known). In steady state, the mass flow rate for all diverging capillaries is4.6with4.7These flow rates provide the inflow rates for a new network that we obtain by eliminating the node *i*. Exactly the same procedure can be applied to this new network. We repeat this scheme until the flow rates in all branches are known.

Equations (4.6) and (4.7) constitute implicitly the bifurcation rule in the continuous model because they enforce the self-regulation mechanism. This becomes intuitively clear for the special configuration where the resistances of all diverging vessels are the same, *r*_{ij}=*r*. In that case, equations (4.6) and (4.7) simplify to *q*_{ij}=*q*_{i}/*N*_{diverging} (where *N*_{diverging} is the number of diverging vessels) such that the inflow *q*_{i} is distributed equally to all diverging vessels.

### (c) Iterative solution of the pressure and flow fields

In order to solve for all node pressure values, mass flow rates and tube haematocrits, we combine the reordering approach with the system of equations (4.5). This yields the following iterative procedure:

Set

*ν*=0 and .Solve the linear system of equations4.8for the node pressure values .

The new node pressure values are used for reordering and the new mass flow rates are obtained from equation (4.6).

Increment the iteration count

*ν*and continue with step 2 until the algorithm has converged, i.e. when .Finally, the converged node pressure values and mass flow rates are used to calculate the mean tube haematocrits

*H*_{ij}according to equation (4.3).

In the following, we use this algorithm to study a number of network configurations similar to the networks studied in §3.

#### (i) Homogeneous networks

First, we study the steady-state configuration of a homogeneous network where all capillaries have the same diameter *D*_{ij}=*D* and where and *β*_{ij}=*β*. In that idealized case, it is not necessary to apply the full iterative process introduced above. Instead, we find that equation (4.5) reduces to the linear system4.9where we expressed the RBC inflow by *q*_{i} and the inflow haematocrit *H*_{i}. If the injection flow rates *q*_{i} and the corresponding injection haematocrits *H*_{i} are specified, the pressure distribution can be computed directly from equation (4.9), regardless of how the RBCs get distributed in the network.

For the special case where the injected haematocrit is constant, i.e. *H*_{i}=*H*, we can derive an important and surprising result. We find that a network without RBCs yields the same system of equations as described in equation (4.9), except that the source terms on the right-hand side of equation (4.9) are scaled with the factor 1+*βH*. Therefore, we find that the steady-state pressure field in a capillary network is the same as the pressure field in a network without RBCs up to the factor 1+*βH*.

Once we have computed the node pressure values, we can compute the mass flow rates through the reordering approach. For the generic square network introduced in §3, the flow rate into node can be determined recursively by4.10where *δ*_{ij} is the Kronecker delta. The result of this algorithm for a network with *n*=4 is illustrated in figure 6. It has been confirmed that the steady-state results for the discrete RBC dynamics from §3 converge towards the continuous result as the number of discrete RBCs increases.

Figure 7*a* shows the pressure field, mass flow rate, haematocrit and passive tracer distribution in a homogeneous 30×30 network with and without RBCs. As expected, we find that the pressure fields are the same up to a multiplicative constant. At the same time, we observe strong differences in the transport characteristics. The distribution of a passive tracer after the dimensionless time *t*=0.4*τ* illustrates these differences most clearly. Without RBCs, the tracer front assumes roughly a circular shape. With RBCs, most of the tracer moves along a preferential highway directly from the source to the sink. The haematocrit distribution is qualitatively the same as the distribution in the steady-state configuration in figure 2*d*. The highest haematocrits are found in areas where the flow rates are small, i.e. a large fraction of the RBCs is used to ‘plug’ the network in these places.

To illustrate this ‘plugging’, figure 7*b* shows results for a truncated network. In such a network, we eliminate all connections for which the effective resistance is larger than a certain threshold. In other words, we assume that the capillaries with high haematocrits are ‘plugged’ such that these capillaries might as well be omitted from the network. Surprisingly, we find that the truncated network yields almost the same network dynamics as the original configuration in figure 7*a*, whereas the haematocrit distribution in the truncated network is more homogeneous.

We suspect that the truncated network is closer to a natural network topology than the original square network because the truncated network does not contain any plugged vessels that are almost useless for transporting oxygen or nutrients. Therefore, the remaining numerical experiments will use the truncated network.

#### (ii) Heterogeneous networks

Next, we study heterogeneous networks with different capillary diameters *D*_{ij}. In this more general configuration, it is not possible to use the idealized equations (4.9) and (4.10) for homogeneous networks such that we have to apply the iterative solution algorithm.

Figure 7*c* shows results for the truncated network from the previous experiment but with randomized diameters, i.e. *D*_{ij}=*D*[0.25+1.5*ξ*_{ij}]^{1/4}, where *ξ*_{ij}∈[0,1] is a uniform random number. Aside from some localized differences, the results look very similar to the results shown in figure 7*b*. Apparently, the overall network dynamics are very robust with respect to statistical fluctuations in the vessel morphology.

#### (iii) Networks with local constrictions or dilations

As a final generalization, we consider networks with constant diameters *D*_{ij}=*D* at the bifurcations (and with *β*_{ij}=*β*, ), but with constrictions or dilations between the nodes *i* and *j*. Because our mathematical framework has been developed for capillaries with a constant diameter, we have to modify our network configuration. To this end, we split a capillary into three subsections by introducing two auxiliary nodes. The two outer sections have the diameter *D*, whereas the interior section of length *L*′ has the diameter *D*′ and an apparent intrinsic viscosity *β*′. Each subsection can now be treated as a separate capillary. The application of the bifurcation rule at the auxiliary nodes is trivial because there is only a single capillary to choose from. It can even be shown that equation (2.9), for computing the effective resistance, remains valid, albeit with modified parameters and . The effective resistance *R*^{e} of the whole capillary is given by4.11where *H*_{ij} is the tube haematocrit in the unaltered vessel sections,4.12aand4.12bWith these modifications, we can use again the iterative solution algorithm to find the steady-state solution of such networks. For the corresponding numerical experiments, we use the truncated network and dilate or constrict capillaries in the centre of the network.

Figure 7*d* shows the relative changes in the haematocrit (top row) and flow rate (bottom row) with respect to the truncated network shown in figure 7*b*. We study four different cases.

*— Horizontal dilation*. The diameter of the horizontal capillary in the centre of the network is increased such that the resulting specific resistance is reduced by 50 per cent.*— Vertical dilation*. The diameter of the vertical capillary in the centre of the network is increased such that the resulting specific resistance is reduced by 50 per cent.*— Horizontal and vertical dilation*. The diameters of the horizontal and the vertical capillaries in the centre of the network are increased such that both specific resistances are reduced by 50 per cent.*— Horizontal and vertical constriction*. The diameters of the horizontal and the vertical capillaries in the centre of the network are decreased such that both specific resistances are doubled.

For the horizontal dilation (leftmost column in figure 7*d*), we find that the haematocrit increases locally by more than 5 per cent to compensate for the lower specific resistance, i.e. there are more RBCs in the dilated capillary to ‘readjust’ the effective resistance to its undilated value. At the same time, we find a region just above the dilated capillary with reduced haematocrit, which is probably due to a local shortage of RBCs because of the higher RBC demand in the dilated vessel. The mass flow rates change as well: there is a region of lower flow rate in the wake of the dilated vessel and a region of increased flow rate on the opposite side. Because of the symmetry of the network, the results for the vertical dilation are, of course, identical to those for the horizontal dilation, except that they are mirrored with respect to the main diagonal.

If we dilate the horizontal and the vertical vessels simultaneously (third column in figure 7*d*), the wakes behind the dilations disappear. There is only a very localized increase in the haematocrit (and a small decrease in the neighbouring vessels). The flow rate is not noticeably affected by this dilation. For a constriction of the horizontal and vertical vessels, we obtain the same result, but with opposite sign. Again, we find that local modifications lead to local changes in the haematocrit and barely affect the flow rates.

## 5. Discussion and concluding remarks

We have presented a model for RBC transport in capillary networks, with a focus on the formulation of efficient algorithms and on the fundamental mechanisms of the network dynamics. We have shown that single-phase models for capillary networks lead to wrong results. The unphysiological build-up in the corners of our generic network, for instance, would not be detected by a single-phase Kirchhoff approach. This finding is further supported by results as in figure 7*a*, which make it clear that the presence of RBCs in a capillary network has a dramatic impact on the network dynamics. Figure 3 emphasizes the importance of using a full model with a haematocrit-dependent flow resistance as well as a bifurcation rule.

Of course, the proposed bifurcation rule may be questioned. Results by Desjardins & Duling (1987) suggest that the distribution of the discharge haematocrit throughout a capillary network is quite homogeneous, whereas the variations in the tube haematocrit are mostly due to an ESL rather than due to phase separation at bifurcations. If this were the case, we would need to replace our bifurcation rule by a random bifurcation biased by the flow rate. Apart from the missing dynamical effects of the ESL (see the discussion in subsequent sections), we would then obtain results as in figure 3*b*. For the following reasons, however, we are confident that our bifurcation rule is appropriate: First, it is derived on the basis of strict physical arguments for a low Reynolds number regime. Second, experimental observations at bifurcations (e.g. Fung 1973; Pries *et al.* 1989) are consistent with our numerical results. It can be regarded as a sign of strength of the present formulation that it does not just duplicate the experimentally observed behaviour by design, but that this behaviour follows as a result of the bifurcation rule.

An efficient model for the simulation of the fluid dynamics in the microcirculation is particularly critical for the investigation of organs in which the level of perfusion is high. This can either be because the organ has a high metabolic rate (e.g. heart or brain) or is involved in the clearance of metabolites and toxins (e.g. liver and kidney). The impact of the present model to study cerebral blood flow is evident. However, there are several other fields where our model is relevant. In cancer therapy, for instance, a detailed understanding of the perfusion of densely vascularized tumours (Pries *et al.* 2009) is very important for general advances in tumour biology and an efficient administration of pharmaceuticals. The dynamics at bifurcations can be used to design sophisticated filters (e.g. blood plasma separation; Yang *et al.* 2006).

### (a) Limitations of the model

We have limited our studies to idealized square networks with a high level of symmetry. Certainly, it will be interesting to investigate the dynamics in more irregular and realistic networks (e.g. Wieringa *et al.* 1982; Kassab *et al.* 1999). The local mechanisms (e.g. self-regulation) will also be found in these networks, but it is not clear *per se* what effect the local mechanisms will have on the global dynamics of the whole network.

Further, the present formulation of the model neglects the distensibility of the vessels (Kassab *et al.* 1999) and the role of the ESL. The ESL interacts with the plasma and the blood particles, and has an influence on the local haematocrit and flow resistance (Desjardins & Duling 1990; Vink & Duling 1996; Constantinescu *et al.* 2001). The effect of an ESL could be introduced to our model by a velocity-dependent resistance of the RBC, i.e. the parameter in equation (2.7) is then also a function of the RBC velocity . This will have a considerable impact on the transient dynamics of RBC transport. Nevertheless, the self-regulation mechanism will remain intact as long as the relations between the intrinsic viscosity and the RBC velocity and the haematocrit are monotonic. Fig. 4 in Secomb *et al.* (2003) and fig. 4 in Pries & Secomb (2005), indicate that these relations are indeed monotonic. Therefore, we expect that the distribution of the flow rates in the steady-state configuration would remain unchanged (up to a constant multiplicative factor), whereas the haematocrit distribution would adapt to the changed conditions.

The present model neglects the presence of other blood particles apart from RBCs. Possible implications of leucocytes on the RBC transport have been discussed by Schmid-Schönbein *et al.* (1980) with respect to the self-regulation mechanisms and by Han *et al.* (2006) with respect to the endothelial glycocalyx.

Although we believe that our formulation represents an appropriate idealization of the actual dynamics at a bifurcation, we would like to point out that our bifurcation rule does not prevent an RBC from entering a capillary that is already filled with RBCs. Such a situation can be found in diagram 1 of figure 4, which shows more than 10 RBCs in a single capillary section of 50 μm length. We can prevent this from happening by adding another condition to the bifurcation rule that allows RBCs to enter a capillary only if is satisfied. We have performed simulations with this modified bifurcation rule and found no fundamental differences in the presented results. We only noticed that the RBC distribution is somewhat more homogeneous and that the ‘highway’ along the main diagonal is less pronounced.

### (b) Relevance for the study of cerebral blood flow

A better knowledge of the cerebral microvascular system and the regulation of cerebral blood flow is fundamental for understanding the dynamics of many physiological and pathological processes in the brain. Furthermore, most of the functional neuroimaging methods, including the widely applied functional magnetic resonance imaging, rely on the haemodynamic response following neural activation. As a consequence, in order to correctly interpret the results of neuro-imaging experiments, mechanisms governing cerebral blood flow regulation must be well understood. Despite its small relative size, the brain uses approximately one-quarter of the body’s total glucose turnover and about one-fifth of the oxygen. The continuous delivery of oxygen is particularly critical and already temporary shortages can lead to transient failures or irreversible damage of nervous tissue. The present study provides the basis for studying oxygen transfer from blood to tissue, particularly when applied to realistic vascular networks. More specifically, the haematocrit and RBC residence time fundamentally determine the tissue oxygen extraction. Recent studies have demonstrated that pericyte constrictions can regulate vascular resistance on the capillary level (Peppiatt *et al.* 2006; Yemisci *et al.* 2009). Our simulations provide first evidence that the effect of a local increase in capillary diameter is a concomitant increase in the local haematocrit. Recent advances in high-resolution imaging (e.g. two-photon laser scanning microscopy) should be exploited to confirm this result *in vivo*. Another interesting finding is the fact that RBC flow stagnates in specific areas of our artificial network. When looking at this phenomenon from a developmental viewpoint, it becomes evident that the vascular topology undergoes a kind of selection process. Flow in ‘inefficient’ capillaries would stall and would eventually lead to the degeneration of these vessels. In fact, there is evidence from histological studies that capillaries do indeed degenerate, which results in basal membrane remnants known as ‘string vessels’ (Challa *et al.* 2002).

The presented model for RBC transport in capillary networks is a building block for a comprehensive model of human physiology. In such a framework, the model can serve as an algorithmic tool for simulating blood flow in a spatially resolved capillary network. It can be linked through in- and outflow boundary conditions to a model for flow in larger blood vessels. The link to the perfused organ is established through the haematocrit distribution (oxygen supply) and/or through a passive tracer distribution (e.g. glucose concentration). Alternatively, the model may be used in a multi-scale approach where basic results for the network dynamics are upscaled to a next higher level of spatial abstraction. In that sense, numerical results for isolated prototypical networks (e.g. figure 7) could be used to estimate the permeability tensor for Darcy’s law (Reichold *et al.* 2009), such that capillaries can be treated as a spatial continuum between larger blood vessels.

## Footnotes

One contribution of 13 to a Theme Issue ‘The virtual physiological human: computer simulation for integrative biomedicine II’.

- © 2010 The Royal Society