#### Abstract

A computational fluid dynamics (CFD) model is developed to simulate the flow and delivery of oxygen and other substances in a capillary network. A three-dimensional capillary network has been constructed to replicate the one studied by Secomb et al. (2000), and the computational framework features a non-Newtonian viscosity model of blood and the oxygen transport model including in-stream oxygen-hemoglobin dissociation and wall flux due to tissue absorption, as well as an ability to study delivery of drugs and other materials in the capillary streams. The model is first run to compute the volumetric flow rates from the velocity profiles in the segments and compared with Secomb’s work with good agreement. Effects of abnormal pressure and stenosis conditions, as well as those arising from different capillary configurations, on the flow and oxygen delivery are investigated, along with a brief look at the unsteady effects and drug dispersion in the capillary network. The current approach allows for inclusion of oxygen and other material transports, including drugs, nutrients, or contaminants based on the flow simulations. Also, three-dimensional models of complex circulatory systems ranging in scale from macro- to microvascular vessels, in principle, can be constructed and analyzed in detail using the current method.

#### 1. Introduction

The capillary vessels, among other functions, transport oxygen, carbon dioxide, and other materials (nutrients and drugs) to and from cells. These smaller blood vessels in the body (the arterioles, venules, and capillaries) make up the so-called microcirculatory system. Larger blood vessels, that is, arteries and veins with inner diameters greater than approximately 100 *μ*m, are termed “macrovascular.” The two systems have distinguishing characteristics. The macrovascular vessels serve as reservoirs of pressure on the high and low side of the cardiovascular system, as well as a passage for relatively large flow rates of blood between the heart and peripheral organs. The microcirculatory vessels are responsible for regulating flow and also for direct material transport at the cell levels. The transport may involve oxygen, carbon dioxide, nutrients, and other biochemical species. About 80% of the total pressure drop in the circulatory system occurs in these microcirculation networks, as the vessel diameters are small and lengths quite large in total [1]. The microcirculation system inherently includes complex flow patterns such as bifurcations in forward and backward directions, constrictions, and flow turns. It is of both academic and practical biomedical interest to be able to determine the fluid dynamics of the circulatory systems, as well as the transport characteristics of oxygen, carbon dioxide, and other species, and there have been many works devoted to this topic [2–9]. While some of the adverse macrovascular flow conditions and their effects are well known, problems in the microvascular networks can lead to equally malignant conditions. For example, hypoxia in the brain capillaries can lead to irreversible nerve cell damage in a matter of minutes. The critical transports occur in the microvascular network, and understanding of this phenomenon is quite important and may suggest methods to improve and optimize circulation at this level under diverse biomedical conditions of interest.

There have been numerous studies of the flow in circulatory systems, with most early works focusing on the fluid dynamic aspects, for example, flow velocity and pressure distributions. Chen and Lu [10], for example, developed a three-dimensional model that uses a finite-element method to determine flow properties in a single, lager vessel (carotid or aorta) bifurcation. This study revealed the effects of the non-Newtonian fluid property and bifurcation geometry on the wall shear stress. In a paper by Berry et al. [11], a two-dimensional CFD model was used to characterize the blood flow in a coronary artery in the presence of a metallic stent. Dzwinel et al. [12] modeled red blood cells as discreet particles rather than using a continuum approach. They then studied the blood dynamics in capillary vessels with diameters of 10 and 25 *μ*m in a simple, straight segment. Pries et al. [13] devised a method of determining the flow rate and hematocrit in each segment of a capillary network based on mass conservation at node points. First, a linear analysis was used in the computation of the flow in each segment in the network and the pressure at each node, assuming that the apparent viscosity and other rheological parameters are known. A conservation principle was then applied, where the inflows and outflows must sum to zero at a node point where three segments meet. Unsteady and three-dimensional flow simulations in realistic geometry are also well within computational fluid dynamic capabilities, as shown, for example, by Shahcheraghi et al. [14] who investigated the pulsatile flows in the human aortic arch. Transient coronary microcirculation has also been simulated by Lee et al. [15]. These are a small subset of the studies that have been done to serve as illustrations that flow simulation is the first step in studying transport in circulation networks at various scales.

In addition to the above purely fluid dynamical aspects, there has been increasing attention on material transport within the circulatory networks involving oxygen, lipoprotein distributions, and drug delivery. An early study by Eggleton et al. [16] considers the transport of oxygen in a linear capillary channel, by calculating the in-fluid advection of oxygen along with absorption of oxygen at the vessel wall. The circulation and uptake of oxygen are an important issue at the capillary scales, since they determine the delivery of oxygen to tissues in complex capillary networks. An overall model of oxygen uptake has also been set up by Wang and Hicks [17] to predict the total oxygen uptake in ectothermic vertebrates. A study by Sinek et al. [18] again focuses on the tissue side transport, with a two-dimensional model of the nanoparticle drug delivery to tumor regions. We can also find mathematical and computational studies of oxygen or drug transport in vascular network [19, 20]. However, thus far these models consider quasi-one-dimensional transport in the network in the sense that the flow in each vessel is considered to be one-dimensional. These studies then focus on the vessel wall and tissue transport, to determine the delivery effectiveness. It should also be noted that most of the studies have used a relatively simple flow balance (see (14)) with flow resistance due to the wall friction and little consideration of the detailed flow patterns through bends, branches, or inlet/exit configurations. It has been shown in a recent study [21] that the flow rates can differ by up to 10% in angled microvascular branches when a full three-dimensional model results are compared with those from a simple one-dimensional one.

For oxygen transport, Popel, Goldman, and coworkers [22–24] have used a quasi-one-dimensional capillary transport model to investigate various issues on the oxygen delivery in hexagonally packed muscle fibers. Secomb and coworkers [25–27] have done extensive work on oxygen transport from capillary networks to tissues. The work can be described in a two-component model, where one component consists of a mathematical prescription of the volumetric flow rate and oxygen partial pressure in the capillary network vessel segments. The oxygen partial pressure at each segment is then used as a source term for distributing oxygen to the tissue in the second component. When all of the source terms are added from all of the segments in the network using Green’s function method, then the tissue oxygen distribution is determined. Thus, the oxygen transport in the vessel is first modeled and then used as a “source term” for perfusion into the tissue region. However, what if the capillary vessels are more complex involving flow turns, divergence, and bifurcations and multiple branches? These will result in static pressure change and therefore have a potential to change the volumetric flow rate distribution in the capillary segments. If the volumetric flow rate and the corresponding available oxygen content are to be accurately modeled, then these effects need to be taken into account so that the flow rates can accurately be calculated. If we further expect to simulate a large-scale oxygen transport from the arteries to the network, then such fluid dynamic and transport effects may deviate even further from simple prescriptions based on a one-dimensional model. Thus, the aim of this work is to provide and demonstrate an approach based on computational fluid dynamics (CFD), involving the calculations of the volumetric flow rates, oxygen content, and potentially other material transports, by solving the fundamental governing equations of convection and diffusion through finite difference schemes (embedded in a commercial software, FLUENT). Such results for the fluid dynamic and transport effects can then be used as the basis for studying delivery, corrective measures, and other applications.

#### 2. Computational Methods

The fluid dynamics is in general described by the Navier-Stokes equations (with appropriate treatment of the non-Newtonian viscosity) in which the fluid momentum (the convection term) is balanced by the pressure and viscous forces. In treatments of biological flows, as noted in the introduction, there have been numerous methods to numerically solve the Navier-Stokes equations. The Navier-Stokes equations can in principle incorporate multiphase, particle-laden flows. Red blood cells, for example, have been modeled as nondeformable particles in continuum “background” flow by Dzwinel et al. [12]. Also, both unsteady- and steady-state processes can be calculated by discretizing the Navier-Stokes equations in time and space for the given geometry and boundary conditions. Many studies have focused on the fluid dynamic aspect, that is, flow patterns arising from pressure and velocity profiles. However, most interesting biological processes involve transport of materials (e.g., oxygen, nutrients, foreign substances, or drugs), and this transport is coupled to the flow motion. That is, the material transport is also governed by a balance equation for convection, diffusion, and boundary conditions at the wall for uptake or tissue absorption. This is precisely our approach in computationally simulating the transport of oxygen (and other materials) in a capillary network. The fluid mechanics is solved using a commercial CFD package, FLUENT, while taking into account the nonlinear viscosity, and in addition transport of oxygen in the stream is also computed along with hemoglobin reactions and wall flux terms. This approach is described in further detail below.

First, a geometrical model for the capillary network needs to be constructed. To assemble this network with some level of realism and also to compare with the existing work, we use the geometry that has been studied in Secomb et al. [25], as shown in Figure 1. The model is a three-dimensional section of a rat brain (more detailed description is referred to in [25]); however, in our model all of the capillary vessels are made of a uniform diameter of 5 microns and the flow rates are matched with Secomb’s model rather than velocity. The connections between the branches of different diameters were difficult to implement without creating sharp edges in the computational meshes, particularly when there are triple junctions as shown in Figure 1. In the Secomb et al.’s model [25], the mean vessel diameter is 5.77 *μ*m with a standard deviation of 0.77 *μ*m; therefore the variations in the vessel diameters are about 13%. The effect of varying vessel diameters is currently being investigated in this laboratory. Aside from that aspect, the current geometrical model is a reconstruction of Secomb’s model and consists of seven inlets and two exits as in the original work [25].

In Figure 1, each of the segments is numbered for later identification purposes, in an identical manner as in the work of Secomb et al. [25]. Using the meshing scheme available in FLUENT, cylindrical meshes were constructed for straight sections while tetrahedral meshes were created for junctions, with a total number of nodes being 673,479. Thus, all the flow and transport parameters are computed at 673,479 nodes embedded in the capillary network. Using a PC with Intel Core 2 Quad 2.5 GHz processor with 2 GB SDRAM at 1066 MHz, the nominal CPU time was a moderate 48 hours with the oxygen-hemoglobin reactions included.

Blood consists of red blood cells (RBCs), white blood cells, and platelets suspended in plasma. RBCs can be described as being biconcave discs of about 8 *μ*m in diameter [28]. Capillary vessels often have smaller diameters than this; hence RBC membranes are viscoelastic, allowing the RBC to deform quite drastically. Plasma itself is considered (and measured to be) a Newtonian fluid; however, blood as a whole displays non-Newtonian effects because of the RBC deformations. Due to their relatively low concentrations, white blood cells and platelets do not typically make a significant contribution to the blood viscosity [29]. The effective blood viscosity that includes the effect of RBC content is a function of several different parameters. Blood viscosity increases with increased hematocrit (volume percentage of RBCs). This is due to the deformability of the RBCs. Blood viscosity also increases with decreasing temperature, is very high at low shear rates (less than 100 s^{−1}) with a drop to an asymptotic value as shear rate increases, and decreases with a decrease in capillary diameter.

In this work, we use the Carreau-type viscosity model that has been shown to work well for blood flows [30]: where is the non-Newtonian viscosity, is the power law index which depends on the hematocrit (as well as other blood constituents), is the relaxation time constant, and and are the zero and infinite shear viscosities, respectively. The local shear rate, , is defined as The rate-of-deformation tensor, , is defined as where is the local fluid velocity vector, and the superscript “” means the transpose of the vector.

Shibeshi and Collins [30] compared several different viscosity models using FLUENT software with parameter values gathered from the literature. The resulting velocity profiles (radial and axial), shear stress, and vortex length were investigated and compared. The results of that study indicated that there is a significant variation in viscosity (as a function of shear rate) depending on which model is used and that the appropriateness of the viscosity model depends on the flow conditions. On the other hand, Hsu et al. [31] have considered flow in a T-shaped junction and conclude that Carreau-type model better reproduces the real blood flow characteristics out of several others in use (power-law and Casson models). The Carreau model has been used quite frequently in computational simulations of blood flows and also can readily be incorporated in the FLUENT computational framework (unlike another often used model, Carreau-Yasuda). The values for the Carreau model parameters that were used in the present work have been obtained from Shibeshi and Collins [30] and are presented in Table 1.

As blood flows through the capillary network, oxygen is delivered through the capillary wall to the surrounding tissue. As this occurs, the concentration of free O_{2} (O_{2} dissolved in the stream) decreases [32]. However, there is another mechanism at work which seeks to replenish the O_{2} lost to the tissue. Hemoglobin is an oxygen-carrying protein found in the red blood cells of mammals and many other animals. The reaction between oxygen and the hemoglobin is reversible and can be described as follwos:
Clark Jr. et al. [33] described a way of modeling the oxygen-hemoglobin binding kinetics by considering three distinct species: free O_{2}, oxygenated hemoglobin (), and nonoxygenated hemoglobin (). As the molar concentration of the free O_{2} decreases, so will the molar concentration of as the concentration of Hb increases (due to hemoglobin’s loss of oxygen). The species conservation equation for two of the species can be expressed as follows:
where is the net rate at which the oxygen and hemoglobin bind to form . An equilibrium exists, where there is sufficient free O_{2} in the blood and no reaction (dissociation or association) takes place. In the limit, as ,
where is an experimentally determined equilibrium constant. According to Clark Jr. et al. [33], can be expressed as
where is the dissociation rate, and is the association rate. The dissociation rate depends upon the concentration of oxygenated hemoglobin and can be written as
where is an experimentally determined dissociation rate constant. If we define the total molar concentration of hemoglobin to be
then the association rate is
where and are parameters found from experimental data [33]. Implementation of this “reaction” set is relatively straightforward in the FLUENT computational framework, as is the wall loss rate of oxygen to mimic the oxygen absorption into the surrounding tissue.

To impose oxygen flux at the vessel wall to simulate the delivery of oxygen to the surrounding tissue, a constant O_{2} flux condition at the capillary wall surface is desired. However, within FLUENT computational framework, the mass flux boundary condition cannot easily be set up. Therefore, a wall surface reaction was used in place of a mass flux condition. This reaction would consume the O_{2} at the capillary wall surface at a determined rate that corresponds to the wall loss rate. In the species conservation equation of oxygen, this would appear as a negative source term which occurs only at the wall surface. To avoid errors involving an apparent mass imbalance in the stoichiometric equation, a dummy species was created and called simply “zero-oxygen.” The basic stoichiometric equation describing this surface reaction was
where represents the rate constant of this reaction. Thus, at the wall the oxygen will be depleted to a dummy species called “zero-oxygen” at the rate specified by . For example, an oxygen depletion rate through the vessel wall of 0.01448 (cm^{3}O_{2}(100 g)^{−1 }min^{−1}) is converted, after many unit and geometrical considerations, to the surface reaction rate of
The oxygen depletion rate and therefore the surface reaction rate depend on the physiological conditions, and it is possible to examine how the variations in the oxygen transport rate simply by varying the above parameter, . However, the surface oxygen depletion rate was set to be constant as the baseline case. Depending on physiological conditions, the oxygen uptake rate can certainly be different, and if there were such conditions with available estimates on the oxygen uptake rate data, the current simulations would allow for input of such altered distribution of oxygen in the capillary network.

The aforementioned geometry, governing equations, and reaction kinetic equations are solved using the upwind scheme in FLUENT code with double precision. Various pressure, inlet/exit conditions were run along with introduction of blockages in the capillary network, as will be discussed in the next section. The resistance to blood flow decreases with decreasing vessel diameter in small vessels (<0.33 mm diameter), which is known as the Fahraeus-Lindqvist effect [34]. This effect could be due to a number of factors, such as, the Fahraeus effect itself, since the lowering of the viscosity of a red cell suspension would result from a lowering of the hematocrit in the vessel [30]. In the present study, the capillary network consists of segments of equal diameter, so the Fahraeus-Lindqvist effect is neglected.

Because the blood flow within a physiological capillary network is driven by the pressure gradient, the flow model can be set up where we specify the pressure at the inlets and exits. The pressure drop () across a capillary network in a rat varies from 5 to 20 mmHg depending on location, geometry, and other physiological factors [1], and the pressure drop was set at the inlet with the exit pressure being set to 0 mmHg (gauge pressure). Three different values of were also tested: 5, 10, and 20 mmHg. For example, in the case where = 5 mmHg, each of the seven inlets was given a pressure condition of 5 mmHg, while each of the three outlets was set at pressure of 0 mmHg.

#### 3. Results and Discussion

First, we compare our volumetric flow rate data with those of Secomb et al. [25], as shown in Figure 2(a). As noted above, the capillary network geometry and inlet conditions in the current work were adopted from Secomb et al. [25], and we use the same segment numbers as in that work in Figure 2(a). Table 2 gives the inlet boundary conditions for the volumetric flow rates, which again are the same between Secomb et al.’s [25] and this work. In Secomb et al.’s work [25], the measured diameters vary from one inlet to the other, which explains the various inlet flow velocities in Table 2. At the exit, a fixed pressure boundary condition of 25 mmHg is applied in the current computations. The oxygen mass fraction is nominally set at , with small variations due to static pressure differences corresponding to different flow rates at the inlets. We use the above inlet and exit boundary conditions for all the subsequent simulations in this work unless otherwise noted.

**(a)**

**(b)**

Figure 2(a) shows that the flow rates in the capillary streams in the most part are quite similar to Secomb’s data, except in a small number of segments (e.g., 19–22, 30, 41–45, and 50). The agreement between the full computational simulations contained in this work and those in Secomb et al. [27] is noteworthy in the sense that Secomb et al. [25] use a relatively simple flow balance method to estimate the flow rates. For example, at a node point where three segments meet, the inflows and outflow(s) must sum to zero [22, 25]: where are the segment flow rates (negative for outflows), are the segment conductances (given by Poiseuille’s law), is the pressure at the node, and are the pressures in the adjacent nodes. When this equation was applied at every node in the network, a linear system is found that can be solved for the unknown pressures. The segment conductance according to the Poiseuille’s law takes into account the pressure loss due to the wall viscous stresses in a circular duct. As noted earlier, it has been shown in a recent study [21] that the flow rates can differ by up to 10% in angled microvascular branches when a full three-dimensional model results are compared with those from a simple one-dimensional one. The reason for this difference is attributed to the fact that there can be nonuniform pressure distributions around bends and bifurcations even in microvascular channels. Although the magnitude of the pressure nonuniformity may be small, fluid-dynamically they can alter the flow streamlines and therefore how the flow is distributed in a bifurcation or a triple junction (see Figure 2(b)). These discrepancies will further be amplified in downstream junctions, further enlarging the difference between full three-dimensional model and one-dimensional model results. Indeed, in Figure 2 the largest differences, although rare, occur in segments far downstream the inlets after having gone several flow bends and bifurcations. For example, the segments 19-20 involve a triple junction and a bend and are far downstream inlets 1 and 3, as do segments 21-22, 39-40, and 39-40. Segment 30 is straight, but it is downstream a triple bend and far removed from inlets 6 and 7 as well. These are locations where the differences in the data are more pronounced in Figure 2(a). Thus, it appears that flow bends, constrictions, and bifurcations can cause nonuniform pressure distributions and flow patterns that can alter downstream flow rates in complex three-dimensional capillary networks. One of the advantages of the current fully three-dimensional computational framework is that most of these complex effects can be captured and accounted for. As shown in Figure 1, there are indeed many flow turns and branches in the current model, and some of the discrepancies in Figure 2(a) are attributed to this effect as noted above.

An example of the flow patterns that can be visualized using the current computational model based on the numerical solutions of the fundamental governing equations is illustrated in Figure 2(b), where we plot the detailed velocity vectors at a junction between segments 1, 2, 3, and 24. The junction walls are “hidden” in Figure 2(b) to show the internal velocity vectors. Segment 1 is an exit, so that the streams from the three other segments join flow out through segment 1. Right at the junction, the fluid streams merge and the velocity field can involve some complex flow patterns. Also, it can be observed that the velocity magnitudes are much smaller close to the walls, due to the viscous effects. All of the flow velocities and therefore the volumetric flow rates in the segments, which are integrated from the velocity data, are computed in this manner, instead of relying on one-dimensional models. Again, the flow rates and also the flow patterns will affect the transport of oxygen and other materials in the capillary network.

A full computational model is able to take into account all of the pressure nonuniformities and their effects on the local volumetric flow rates and flow patterns across a large range of scales and geometries. In addition, the current approach allows us to use a more realistic non-Newtonian viscosity model such as the one shown in (1), whereas in one-dimensional network models as in Secomb et al. [25] a constant viscosity is used in a Poiseuille’s law regardless of the shear rates. These aspects may be of importance since the transport characteristics depend on the volumetric flow rates, which in turn depend on the exact viscosity and the fully three-dimensional geometry. Using the fundamental governing equations to solve the fluid mechanics of the capillary flow, the full flow effects are accounted for and quite accurate flow distributions in complex, three-dimensional geometry at arbitrary scales can be computed. The primary function of the capillary network is to transport and distribute oxygen and other essential substances. Since the capillary diameter is quite small, this transport often is determined by the flow rate in the longitudinal direction. The net amount available for oxygen transport across the capillary wall, for example, is determined by the wall flux and the flow rate. For this reason, it is important to obtain correct flow rate in complex capillary network geometry, and based on this information we can assess the effects of changes in the conditions within the capillary network and examine essential biotransport processes such as oxygen distribution and drug delivery.

For the above inlet and exit conditions, we have run the flow and oxygen transport calculations for both unsteady and steady states. The unsteady calculations were an attempt to observe the time evolution of the flow and oxygen distributions, with the inlet and exit boundary conditions suddenly applied at = 0 in a stepwise manner. Therefore, it does not yet fully simulate the pulsating pressure field representative of physiological processes. Using this framework, however, more realistic pressure fluctuations as a function of time can in principle be applied, as in pulsating flows. Figure 3 illustrates the development of the oxygen partial pressure at various time steps after the application of the boundary conditions. It can be observed that steady state is achieved in less than 10 ms. Due to the short distances involved in a capillary network, both the flow and oxygen transport are established in a short amount of time in comparison to, say, the time scale of pulsatile pressure fluctuations. For this reason, it appears that steady-state computations suffice to characterize the main features of flow and transport in capillary networks of this scale. In the remainder of this work, we confine our interests to steady-state results.

The steady-state data from Figure 3 can be visualized as a contour plot of oxygen mass fraction in the capillary network as shown in Figure 4, which shows the spatial distribution more directly. As shown in Table 2, the volumetric flow rate conditions are different from one inlet to the other, and Figure 4 shows that the higher volumetric flow rates are associated with higher oxygen mass fractions by a roughly proportional amount. This is due to the higher available supply of oxygen at larger flow rates, at constant wall loss conditions that exist throughout the capillary network. This also affects the oxygen distribution in the rest of the capillary network. At low-activity or low-oxygen-demand conditions, the oxygen is well supplied [35]. It is only during hyperactive states or stenosis conditions that the oxygen delivery becomes critical and regions of hypoxia can develop. Hypoxia can lead to irreversible brain cell damage in minutes. The ability to compute and visualize the exact distributions of oxygen in the capillaries is, therefore, an important tool to investigate the above and other specific conditions that may arise.

Abnormal pressure conditions are frequently observed in circulatory systems and are referred to as hyper- or hypotension. They are known to be the causes of a variety of malignant physiological conditions. In this work, we consider the consequences of altered pressure conditions on flow and transport phenomena, although based on the pressure changes local mechanical wall stresses or possibly deformations, for example, can also be calculated. Figure 5 shows the changes in the oxygen mass fraction partial pressure in the capillary network at three difference pressure differences from the inlet to exit: 5, 10, and 20 mmHg. Again, the exit pressure is fixed and the inlet pressures are adjusted to result in the above pressure differences from the inlets to the exits. Higher pressure difference between the inlet and exit, obviously, causes an increase in the flow speed, and the consequence is that the residence time for the oxygen-containing stream in the capillary network is reduced. Substantial differences in the oxygen mass fraction can be seen at different pressure conditions, with highest oxygen levels at 20 mmHg pressure drop. For example, about equal levels of oxygen are found close to the inlet, with the levels diverging further into the network at different pressure conditions. Again, the higher flow speed and lower residence time conspire to reduce the transport through the vessel walls leaving more oxygen in the stream at higher pressure conditions. From further calculations, it is found that for each 25% increase in the pressure drop from the normal condition there is a 15% increase in the oxygen that is leaving the network on the average. That essentially amounts to 15% decrease in oxygen being delivered through the vessel walls. Higher residence time associated with low volumetric flow rates at low pressure conditions would indicate a more successful delivery of oxygen through the vessel walls. On the other extreme, however, hypotension or excessively low pressure difference may cause insufficient fluid momentum to reach all of the capillary network branches leading to localized spots of low-oxygen delivery. In general, higher pressure at the inlet results in increased static pressure through the capillary vessels, leading to increased vessel stress in a proportionate manner.

Stenosis refers to an abnormal narrowing of arterial blood vessels, caused by a number of factors including inflammation and accumulation of lipoprotein. Stenoses in capillary networks can be particularly problematic, as they may not be subject to as obvious early symptoms or diagnosis as in coronary vessels. Although flow constrictions associated with stenosis can be circumvented in the current capillary network due to redundant branches, stenosis can in general cause serious medical conditions ranging from reduced oxygen delivery to heart attacks even from those in capillary vessels. Moreover, hypertension is believed to cause narrowing of the capillary vessels, which will yet increase the local pressure initiating a vicious cycle of increasing flow restrictions. In this work, we simulate the stenosis points as local pressure increases at the segments 25 and 32, that is, downstream the inlets 5 and 7 (see Figure 1). Local pressure increase is fluid dynamic manifestation of a flow constriction such as stenosis. For example, an orifice placed in a pipe acts to increase the downstream pressure device in fluid flows, according to where the subscripts “1” and “2” refer to the upstream and downstream conditions, respectively. Then, the severity of the stenosis can be equated with the magnitude of the pressure increase at that point. In principle, any arbitrary area ratio, , can be simulated by setting the pressure increase appropriately. Figure 6 shows the static pressure distribution involving these two stenosis points (here, set up as complete blockage of the flow) in the capillary network. The most pronounced reductions in the static pressures are observed at segments 21 through 24. These segments are, as expected, points directly downstream the stenosis points. Overall, however, the pressure reduction in the other parts of the network is relatively minor, and the reason for this is the redundancy in the flow channels in this particular network, with multiple inlets and exits along with interconnections. The large number of inlets, in comparison to the number of stenosis points, allows for nearly normal circulation in the network. This obviously will not be the case in isolated networks with small number of inlets. The oxygen mass fraction shows a similar effect in Figure 7 where the lowered oxygen levels are found upstream the stenosis points as the flow is not able to proceed in these segments and the oxygen content is stagnantly depleted through the vessel walls.

The baseline geometry for the capillary network again follows that provided by Secomb et al. [25]. While retaining the basic geometry, however, we can alter number and designations of the inlets and exits, to see the effects on the pressure and oxygen distributions. As can be seen in the previous plots thus far, both the static pressure and oxygen levels are highly nonuniform in the current network configuration. To be sure, the highest peaks of both the pressure and oxygen levels are found at the seven inlets. Outside these peaks and low levels at the exits, there is still a fairly high level of nonuniformity. The optimum capillary network configuration would be one where the oxygen undergoes a linear, monotonic decrease with minimum pressure loss. This will be an indication that the oxygen is being delivered through the vessel walls in a nearly uniform manner and that the flow resistance is the minimum. Biological systems exhibit naturally evolved configurations, and it is of interest to see how optimum configurations may be reached in capillary networks. A full permutation of the inlet and exit configurations would require much computational time, and we defer such complete investigation to a subsequent work using a systematic optimization algorithm. Here, we test three alterations to the inlet and exit configurations to check their effects on the pressure and oxygen level distributions. The variable parameters are the number and location of the inlets and exits, altered from the seven-inlet/three-exit configuration shown in Figure 1. First, we consider a configuration with an equal number of inlets and exits, five each. Next, we rearrange the inlets and exits so that there are three inlets and seven exits. The volumetric flow rates and the oxygen mass fraction for the two modified configurations are shown in Figures 8 and 9, respectively, along with those from the original configuration.

It can be observed in Figure 8 that a substantial change in flow rates occurs, with two configurations showing the most uniform. It should be noted that in order to make a meaningful comparison the total volumetric flow rate is fixed in all of the configurations. Thus, the larger number of exits in the original configuration is associated with generally lower flow rates in the inlets and large flow rates at the three exits (e.g., segment 1 in Figure 8). The converse is true for 3-inlet/7-exit configuration. Overall, an equal number of inlets and exits are expected to yield the most uniform flow rates, as verified in Figure 8. It is of interest to see how the uniformity of the flow rates translates to the oxygen distribution, as in Figure 9. It can be seen that the extreme levels, both high and low, are found for the uneven number of inlets and exits, and again the 5-inlet/5-exit configuration exhibits the most moderate variations. This is indicative of the gradual depletion in the stream or more uniform diffusion through the vessel walls.

Transport of other materials, such as drugs, nutrients, or other foreign substances, can also be tracked using the current computational framework, so long as their in-stream and wall diffusion properties are known. The wall diffusion in most instances is a strong function of the local physiological conditions. We take a liquid drug, digoxin (C_{41}H_{64}O_{14}), as an example case, by “injecting” it at the inlet(s) and show its distribution in the capillary network under two different injection conditions. Digoxin is a widely used drug for treatment of heart conditions, including arterial fibrillation and flutter. First, digoxin is injected uniformly at all of the seven inlets, so it simulates a case where the injection occurs far upstream the capillary. This is compared with a capillary injection where the same total amount is injected at only one of the inlets (inlet 5 in Figure 1), in Figure 10. It can be observed that the localized capillary injection results in digoxin diffusing only to segments connected to inlet 5, which is to be expected. The multiple injection scheme results in a drug distribution in the capillary which is quite similar to that of oxygen shown in previous plots. The reason is that the stream flow patterns are identical and the capillary distribution is primarily due to the flow patterns minus the wall diffusional losses. Almost any drug and its delivery can be simulated in the capillary and other networks using this framework, as long as its diffusive and biochemical reactive properties, if applicable, are known. Dissolution and other phase changes can also be handled using standard thermodynamic models, as needed.

#### 4. Conclusions

A realistic three-dimensional model of a capillary network is used to simulate the flow processes and the transport of oxygen. The governing equations are from the fundamental conservation principles, and the solutions are obtained numerically via a use of a commercial software called FLUENT. Additional hemoglobin-oxygen kinetics along with oxygen in-stream diffusion and delivery through the vessel walls are incorporated in the above computational framework. The results for the volumetric flow rate have been compared with the data of Secomb et al. [25], and good agreement is found. In addition, current approach allows for inclusion of oxygen and other materials transport, including drugs, nutrients, or contaminants. Moreover, three-dimensional models of complex circulatory systems ranging in scale from macro- to microvascular vessels in principle can be constructed and analyzed in detail using the current method. Potential applications are delivery methods for drug or nutrients tailored for individual circulatory network geometry, design of surgical alterations for optimized microcirculations, and so on. A number of flow and oxygen transport conditions were examined using the current method for this particular example of microvascular network, including pressure variations, unsteady effects, stenosis conditions, and changes in the inlet and exit configurations.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.