Research Article  Open Access
Computational Simulations of Flow and Oxygen/Drug Delivery in a ThreeDimensional Capillary Network
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 threedimensional capillary network has been constructed to replicate the one studied by Secomb et al. (2000), and the computational framework features a nonNewtonian viscosity model of blood and the oxygen transport model including instream oxygenhemoglobin 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, threedimensional 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 socalled 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 threedimensional model that uses a finiteelement method to determine flow properties in a single, lager vessel (carotid or aorta) bifurcation. This study revealed the effects of the nonNewtonian fluid property and bifurcation geometry on the wall shear stress. In a paper by Berry et al. [11], a twodimensional 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 threedimensional 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 influid 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 twodimensional 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 quasionedimensional transport in the network in the sense that the flow in each vessel is considered to be onedimensional. 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 threedimensional model results are compared with those from a simple onedimensional one.
For oxygen transport, Popel, Goldman, and coworkers [22–24] have used a quasionedimensional 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 twocomponent 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 largescale 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 onedimensional 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 NavierStokes equations (with appropriate treatment of the nonNewtonian 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 NavierStokes equations. The NavierStokes equations can in principle incorporate multiphase, particleladen 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 steadystate processes can be calculated by discretizing the NavierStokes 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 threedimensional 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 oxygenhemoglobin 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 nonNewtonian 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 Carreautype viscosity model that has been shown to work well for blood flows [30]: where is the nonNewtonian 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 rateofdeformation 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 Tshaped junction and conclude that Carreautype model better reproduces the real blood flow characteristics out of several others in use (powerlaw 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, CarreauYasuda). 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 oxygencarrying 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 oxygenhemoglobin 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 “zerooxygen.” 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 “zerooxygen” 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 FahraeusLindqvist 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 FahraeusLindqvist 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 threedimensional model results are compared with those from a simple onedimensional 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, fluiddynamically 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 threedimensional model and onedimensional 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 1920 involve a triple junction and a bend and are far downstream inlets 1 and 3, as do segments 2122, 3940, and 3940. 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 threedimensional capillary networks. One of the advantages of the current fully threedimensional 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 onedimensional 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 nonNewtonian viscosity model such as the one shown in (1), whereas in onedimensional 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 threedimensional 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, threedimensional 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 steadystate 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 steadystate results.
The steadystate 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 lowactivity or lowoxygendemand 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 oxygencontaining 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 lowoxygen 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 seveninlet/threeexit 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 3inlet/7exit 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 5inlet/5exit 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 instream 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 threedimensional 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 hemoglobinoxygen kinetics along with oxygen instream 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, threedimensional 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.
References
 A. S. Popel and P. C. Johnson, “Microcirculation and hemorheology,” Annual Review of Fluid Mechanics, vol. 37, pp. 43–69, 2005. View at: Publisher Site  Google Scholar
 J. Aroesty and J. F. Gross, “Convection and diffusion in the microcirculation,” Microvascular Research, vol. 2, no. 3, pp. 247–267, 1970. View at: Publisher Site  Google Scholar
 A. G. Hudetz, A. S. Greene, G. Feher, D. E. Knuese, and A. W. Cowley Jr., “Imaging system for threedimensional mapping of cerebrocortical capillary networks in vivo,” Microvascular Research, vol. 46, no. 3, pp. 293–309, 1993. View at: Publisher Site  Google Scholar
 A. G. Hudetz, G. Fehér, and J. P. Kampine, “Heterogeneous autoregulation of cerebrocortical capillary flow: evidence for functional thoroughfare channels?” Microvascular Research, vol. 51, no. 1, pp. 131–136, 1996. View at: Publisher Site  Google Scholar
 A. G. Hudetz, “Blood flow in the cerebral capillary network: a review emphasizing observations with intravital microscopy,” Microcirculation, vol. 4, no. 2, pp. 233–252, 1997. View at: Google Scholar
 W. Malkusch, M. A. Konerding, B. Klapthor, and J. Bruch, “A simple and accurate method for 3D measurements in microcorrosion casts illustrated with tumour vascularization,” Analytical Cellular Pathology, vol. 9, no. 1, pp. 69–81, 1995. View at: Google Scholar
 E. D. F. Motti, H.G. Imhof, and M. G. Yasargil, “The terminal vascular bed in the superficial cortex of the rat. An SEM study of corrosion casts,” Journal of Neurosurgery, vol. 65, no. 6, pp. 834–846, 1986. View at: Google Scholar
 A. Krogh, The Anatomy and Physiology of Capillaries, Hafner, New York, NY, USA, 1959.
 J. Prothero and A. C. Burton, “The physics of blood flow in capillaries. I. The nature of the motion,” Biophysical Journal, vol. 1, pp. 565–579, 1961. View at: Google Scholar
 J. Chen and X. Lu, “Numerical investigation of the nonNewtonian blood flow in a bifurcation model with a nonplanar branch,” Journal of Biomechanics, vol. 37, no. 12, pp. 1899–1911, 2004. View at: Publisher Site  Google Scholar
 J. L. Berry, A. Santamarina, J. E. Moore Jr., S. Roychowdhury, and W. D. Routh, “Experimental and computational flow evaluation of coronary stents,” Annals of Biomedical Engineering, vol. 28, no. 4, pp. 386–398, 2000. View at: Publisher Site  Google Scholar
 W. Dzwinel, K. Boryczko, and D. A. Yuen, “A discreteparticle model of blood dynamics in capillary vessels,” Journal of Colloid and Interface Science, vol. 258, no. 1, pp. 163–173, 2003. View at: Publisher Site  Google Scholar
 A. R. Pries, T. W. Secomb, P. Gaehtgens, and J. F. Gross, “Blood flow in microvascular networks. Experiments and simulation,” Circulation Research, vol. 67, no. 4, pp. 826–834, 1990. View at: Google Scholar
 N. Shahcheranhi, H. A. Dwyer, A. Y. Cheer, A. I. Barakat, and T. Rutaganira, “Unsteady and threedimensional simulation of blood flow in the human aortic arch,” Journal of Biomechanical Engineering, vol. 124, no. 4, pp. 378–387, 2002. View at: Publisher Site  Google Scholar
 J. C. J. Lee, A. J. Pullan, and N. P. Smith, “A computational model of microcirculatory network structure and transient coronary microcirculation,” in Proceedings of the 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC '04), pp. 3808–3811, September 2004. View at: Google Scholar
 C. D. Eggleton, T. K. Roy, and A. S. Popel, “Predictions of capillary oxygen transport in the presence of fluorocarbon additives,” The American Journal of Physiology—Heart and Circulatory Physiology, vol. 275, no. 6, pp. H2250–H2257, 1998. View at: Google Scholar
 T. Wang and J. W. Hicks, “An integrative model to predict maximum O_{2} uptake in animais with central vascular shunts,” Zoology, vol. 105, no. 1, pp. 45–53, 2002. View at: Publisher Site  Google Scholar
 J. Sinek, H. Frieboes, X. Zheng, and V. Cristini, “Twodimensional chemotherapy simulations demonstrate fundamental transport and tumor response limitations involving nanoparticles,” Biomedical Microdevices, vol. 6, no. 4, pp. 297–309, 2004. View at: Publisher Site  Google Scholar
 D. Goldman, “Computational modeling of drug delivery by microvascular networks,” in Proceedings of the IEEE 29th Annual Bioengineering Conference, pp. 321–322, March 2003. View at: Google Scholar
 A. Stéphanou, S. R. McDougall, A. R. A. Anderson, and M. A. J. Chaplain, “Mathematical modelling of flow in 2D and 3D vascular networks: applications to antiangiogenic and chemotherapeutic drug strategies,” Mathematical and Computer Modelling, vol. 41, no. 10, pp. 1137–1156, 2005. View at: Publisher Site  Google Scholar
 M. Z. Pindera, H. Ding, M. M. Athavale, and Z. Chen, “Accuracy of 1D microvascular flow models in the limit of low Reynolds numbers,” Microvascular Research, vol. 77, no. 3, pp. 273–280, 2009. View at: Publisher Site  Google Scholar
 D. Goldman and A. S. Popel, “A computational study of the effect of capillary network anastomoses and tortuosity on oxygen transport,” Journal of Theoretical Biology, vol. 206, no. 2, pp. 181–194, 2000. View at: Publisher Site  Google Scholar
 W. J. Federspiel and A. S. Popel, “A theoretical analysis of the effect of the particulate nature of blood on oxygen release in capillaries,” Microvascular Research, vol. 32, no. 2, pp. 164–189, 1986. View at: Google Scholar
 A. Vadapalli, D. Goldman, and A. S. Popel, “Calculations of oxygen transport by red blood cells and hemoglobin solutions in capillaries,” Artificial Cells, Blood Substitutes, and Immobilization Biotechnology, vol. 30, no. 3, pp. 157–188, 2002. View at: Publisher Site  Google Scholar
 T. W. Secomb, R. Hsu, N. B. Beamer, and B. M. Coull, “Theoretical simulation of oxygen transport to brain by networks of microvessels: effects of oxygen supply and demand on tissue hypoxia,” Microcirculation, vol. 7, no. 4, pp. 237–247, 2000. View at: Google Scholar
 A. R. Pries, B. Reglin, and T. W. Secomb, “Structural adaptation of vascular networks: role of the pressure response,” Hypertension, vol. 38, no. 6, pp. 1476–1479, 2001. View at: Google Scholar
 T. W. Secomb, R. Hsu, E. Y. H. Park, and M. W. Dewhirst, “Green's function methods for analysis of oxygen delivery to tissue by microvascular networks,” Annals of Biomedical Engineering, vol. 32, no. 11, pp. 1519–1529, 2004. View at: Publisher Site  Google Scholar
 F. J. H. Gijsen, F. N. van de Vosse, and J. D. Janssen, “The influence of the nonNewtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model,” Journal of Biomechanics, vol. 32, no. 6, pp. 601–608, 1999. View at: Publisher Site  Google Scholar
 Q. Chen, M. Schlichtherle, and M. Wahlgren, “Molecular aspects of severe malaria,” Clinical Microbiology Reviews, vol. 13, no. 3, pp. 439–450, 2000. View at: Publisher Site  Google Scholar
 S. S. Shibeshi and W. E. Collins, “The rheology of blood flow in a branced arterial system,” Applied Rheology, vol. 15, no. 6, pp. 398–405, 2005. View at: Google Scholar
 C.H. Hsu, H.H. Vu, and Y.H. Kang, “The rheology of blood flow in a branched arterial system with threedimensional model: a numerical study,” Journal of Mechanics, vol. 25, no. 4, pp. N21–N24, 2009. View at: Publisher Site  Google Scholar
 Q. H. Gibson, F. Kreuzer, E. Meda, and F. J. W. Roughton, “The kinetics of human haemoglobin in solution and in the red cell at 37°C,” The Journal of Physiology, vol. 129, no. 1, pp. 65–89, 1955. View at: Google Scholar
 A. Clark Jr., W. J. Federspiel, P. A. A. Clark, and G. R. Cokelet, “Oxygen delivery from red cells,” Biophysical Journal, vol. 47, no. 2 I, pp. 171–181, 1985. View at: Google Scholar
 H. L. Goldsmith, G. R. Cokelet, and P. Gaehtgens, “Robin Fåhraeus: evolution of his concepts in cardiovascular physiology,” The American Journal of Physiology, vol. 257, no. 3, pp. H1005–1015, 1989. View at: Google Scholar
 M. L. Ellsworth, C. G. Ellis, D. Goldman, A. H. Stephenson, H. H. Dietrich, and R. S. Sprague, “Erythrocytes: oxygen sensors and modulators of vascular tone,” Physiology, vol. 24, no. 2, pp. 107–116, 2009. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2014 T.W. Lee et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.