Effects of Fiber Type and Size on the Heterogeneity of Oxygen Distribution in Exercising Skeletal Muscle

The process of oxygen delivery from capillary to muscle fiber is essential for a tissue with variable oxygen demand, such as skeletal muscle. Oxygen distribution in exercising skeletal muscle is regulated by convective oxygen transport in the blood vessels, oxygen diffusion and consumption in the tissue. Spatial heterogeneities in oxygen supply, such as microvascular architecture and hemodynamic variables, had been observed experimentally and their marked effects on oxygen exchange had been confirmed using mathematical models. In this study, we investigate the effects of heterogeneities in oxygen demand on tissue oxygenation distribution using a multiscale oxygen transport model. Muscles are composed of different ratios of the various fiber types. Each fiber type has characteristic values of several parameters, including fiber size, oxygen consumption, myoglobin concentration, and oxygen diffusivity. Using experimentally measured parameters for different fiber types and applying them to the rat extensor digitorum longus muscle, we evaluated the effects of heterogeneous fiber size and fiber type properties on the oxygen distribution profile. Our simulation results suggest a marked increase in spatial heterogeneity of oxygen due to fiber size distribution in a mixed muscle. Our simulations also suggest that the combined effects of fiber type properties, except size, do not contribute significantly to the tissue oxygen spatial heterogeneity. However, the incorporation of the difference in oxygen consumption rates of different fiber types alone causes higher oxygen heterogeneity compared to control cases with uniform fiber properties. In contrast, incorporating variation in other fiber type-specific properties, such as myoglobin concentration, causes little change in spatial tissue oxygenation profiles.


Introduction
Oxygen transport from capillaries to muscle fibers plays an essential role in the maintenance of physiological functions of skeletal muscle across a wide range of conditions and various forms of exercise. Oxygen transport is regulated by convection in blood vessels and by diffusion across the vessel walls and into the interstitial space and parenchymal cells. Hemoglobin transports oxygen in the blood, while myoglobin serves as an oxygen storage and facilitates its diffusion inside myocytes. At resting conditions, oxygen consumption is less than the available oxygen delivered by the microvasculature, with nearly half of the oxygen returned to the venous circulation; thus, muscle at rest is typically oversupplied. During exercise, the oxygen consumption rate in skeletal muscle can increase as much as 50 fold compared to resting conditions [1]. The increased oxygen demand is partially compensated by the increased blood flow in tissue microcirculation (10 to 25 fold) [2], and the increase of oxygen extraction by the tissue. The mismatch between oxygen demand and oxygen supply may lead to regional tissue hypoxia, and prolonged hypoxia can result in angiogenesis (capillary growth from pre-existing vasculature), an adaptive response that leads to a decrease of oxygen diffusion distances [3]. Angiogenesis may also result from the elevated shear stress often associated with functional hyperemia [4]. Insufficient oxygen supply is a major downstream pathological effector of chronic ischemic diseases such as coronary artery disease and peripheral artery disease [5]. In those ischemic diseases, the obstruction of upstream blood vessels limits blood flow and convective oxygen transport.
Theoretical aspects of oxygen transport processes have been extensively studied. The goal has been to describe the non-uniform oxygen distribution in skeletal muscle tissue. In the pioneering work by August Krogh [6], an oxygen transport model (also called the Krogh-Erlang model) was developed to describe a single capillary supplying oxygen to a cylindrical volume of surrounding tissue. It was the first theoretical model used to understand oxygen transport to tissue and was built on a number of idealized assumptions including: constant oxygen consumption rate; uniform oxygen diffusivity; and homogeneous capillary distribution. Since then, many alternative oxygen transport models have been developed (reviewed in [7][8][9]) with a number of additional physiological features neglected in the Krogh-Erlang model, including: variable oxygen consumption rate [10]; myoglobinfacilitated diffusion [11]; flow redistribution [12]; axial oxygen diffusion in the tissue [13]; intravascular oxygen transport resistance [14], and pre-and post-capillary transport [15,16]. In addition, oxygen transport models have been extended to include complex microvascular network geometry with capillary tortuosity and anastomoses [17][18][19]. The Krogh-Erlang model is informative; however because of the complexity and heterogeneity of tissue microarchitecture, it would be difficult to establish a single Krogh-type parameter to describe the fiber types or the vessels adjacent to them. Instead, we construct a spatially detailed heterogeneous model based on the experimentally reported values of capillary:fiber ratio for specific fibers and report the distribution of oxygen across the tissue being simulated.
Skeletal muscle is composed of an array of fibers with capillaries running primarily parallel to the fiber direction. Muscles in the body have different functions and different fiber compositions; fiber type composition within a muscle is closely correlated to the function of that muscle. Fibers are commonly categorized as type I (slow oxidative twitch) and type II (fast twitch) based on contractile properties and oxidative capacity [20][21][22]. Type I fibers are mainly involved in aerobic activities and endurance exercise; they contain ample amounts of mitochondria and utilize mostly oxidative phosphorylation, making them fatigue-resistant. In addition, they are myoglobin-rich with a red appearance. Type II fibers, appearing whiter than the 'red' Type I fibers due to lower myoglobin expression, can be classified into three subtypes, Type IIa, IIb and IId/x. Type IIa (or oxidative fast twitch) fibers generate ATP through the glycolytic cycle but also have a high mitochondrial count, allowing them to obtain ATP through oxidative metabolism. They have higher contraction velocities compared to Type I fibers but are not fatigue-resistant. Type IIb (or glycolytic fast twitch) fibers have far fewer mitochondria and thus depend on glycolytic metabolism to generate ATP. They are activated when short and powerful bursts of contraction are required. Type IId/x fibers have properties intermediate between Type IIa and Type IIb. Type I and Type IIa fibers have higher oxidative capacities and therefore consume more oxygen than Type IId/x and Type IIb. These muscle fibers are distributed throughout the skeletal muscle tissue in varying proportions, depending on the muscle group involved. In addition, there is a correlation between the fiber type and the number of adjacent blood vessels: Type IIb fibers generally have fewer adjacent blood capillaries compared to Type I and Type IIa. In addition to histological and metabolic differences, muscle fibers also differ in their size and their components, with variations in the proportions of aqueous cytosol and lipid-rich membranes and droplets, allowing varying oxygen diffusivity in different fiber types [9].
The objective of our study is to assess the impact of fiber type composition on the heterogeneity of oxygen distribution in tissue using a computational oxygen transport model. Our working hypothesis was that fiber-type dependent characteristics may balance with the number of nearby capillaries to locally 'match' oxygen supply and demand. We further hypothesized that this matching may be more relevant in exercise conditions. One goal of the study was to use computational modeling techniques as a hypothesis testing tool to understand the role of each parameter. The computational model is an extension of our previously developed oxygen model (see Methods). We began by building several tissue geometries that incorporate fiber type-specific properties: at the tissue level this includes the proportion of each type of fiber present; at the fiber level, each fiber type has a different characteristic fiber size, oxygen consumption rate (M c ), oxygen diffusivity (D O2 ), myoglobin concentration (C Mb ), and number of surrounding capillaries. We constructed a set of tissue geometries; each geometry has a different combination of parameters set to be uniform and parameters set to be fiber type-specific. This allows us to capture the different effects of each parameter on spatial heterogeneity in the muscle. For each geometry, the oxygen transport model is used to simulate the spatial profiles of oxygen under exercise conditions of varying intensities. Our simulation results suggest a marked increase of spatial heterogeneity of oxygen due to fiber size distribution in a mixed muscle. The difference in oxygen consumption rates of different fiber types also causes higher oxygen heterogeneity compared to control cases with uniform fiber properties. In contrast, incorporating variation in other fiber type-specific properties, such as myoglobin concentration, causes little change in spatial tissue oxygenation profiles. To the best of our knowledge, this is the first computational model used to examine the importance of fiber type and size in tissue oxygen transport.

Methods
We adapted our multiscale computational model of oxygen transport [18,23,24] to study the effects of fiber type and fiber size on the oxygen distribution and oxygen gradients in working muscle. Two new features were introduced into the model. First, we incorporated six variables representing experimentally-measured physiological properties of specific fiber types in skeletal muscle. These variables are fiber type composition, fiber size, D O2 , M c , C Mb and capillary distribution. Second, the modified model describes the oxygen transport process using coupled partial differential equations in three regions: inside the muscle fibers, in the interstitial space, and inside the vascular space. Previously, fibers and interstitial space had not been treated separately.
The model consists of three components: the muscle geometry module, the microvascular blood flow module, and the oxygen transport module. We will describe each module, with the focus on the new model features in this study.

Skeletal Muscle Geometry Module
We constructed a set of skeletal muscle tissue geometries. These constructions were done computationally, using a similar algorithm to that applied in previous studies [23][24][25][26][27][28], with modifications described below. The muscle in our model is the rat extensor digitorum longus (EDL) that has been used extensively in experimental studies of activity-induced angiogenesis. To assess the effects of fiber type, fiber size and also fiber-type dependent microvascular structures on oxygen distribution profiles, we specifically introduced four types of skeletal muscle geometry. Each geometry (designated G1-G4) incorporates an experimentally-observed number (or percentage) of each fiber type. However they differ as follows: G1, uniform fiber size and uniform capillary distribution; G2, uniform fiber size and fiber-type-specific capillary distribution; G3, non-uniform fiber size and uniform capillary distribution; and G4, non-uniform fiber size and fiber-type-specific capillary distribution.
Computational generation of these 3D muscle geometries (Figures 1 and 2) is implemented as a two-step method, with fiber geometry constructed first followed by microvascular network insertion. The outline of the design algorithm is as follows: (1) In the rectangular cuboid tissue volume, muscle fibers are represented as right circular cylinders; (2) Microvessels, specifically capillaries, are placed between the muscle fibers; the vessels are divided into short segments and each segment is represented as a right circular cylinder; tortuosity and anastomoses are introduced to mimic physiological vascular networks by varying the orientation of adjoining segments; (3) Pre-capillary arterioles and postcapillary venules are placed in a staggered pattern to connect with the microvessel network; the capillaries originate from an arteriole and terminate at a venule approximately 350 mm away.
For fiber geometry with uniform-size fibers (G1 and G2), we place the fibers on a regular grid [23,[25][26][27][28]. To avoid edge effects, a number of fibers are bisected by the edges, but appear on the other side in a periodic fashion. The simulations apply periodic boundary conditions for oxygen diffusion. We then use a randomizing algorithm, combined with experimentally-observed rules, to assign fiber types to fibers. The fiber type composition is dependent on the muscle type, animal species, and individual differences. The proportion of fiber type I in the EDL ranges from 2 to 9%, fiber type IIa from 11 to 27%, and fiber type IIb from 61 to 82% [29]. For this model, we only consider fiber type I, fiber type IIa and fiber type IIb due to incomplete information about fiber type IId/x (for example, values for the myoglobin concentration and oxygen diffusion coefficient in this fiber type are not available). In our constructed geometries, we consider 48 fibers in the tissue volume (35862336800 mm), of which 4 fibers are type I (8.3%), 11 are type IIa (22.9%) and 33 are type IIb (68.8%). Second, based on histochemical staining reported in the literature [29], there are two rules for fiber type distribution in rat EDL muscle: (1) type I fibers are not typically found next to each other; (2) type IIa fibers are typically found close to type I fibers.
For the microvascular network in the tissue, capillary placement is also governed by a stochastic algorithm, augmented by several rules including: (1) the average capillary-to-fiber ratio (C:F) is 1.1 (i.e., the total number of capillaries is about 52). The capillary density is 623/mm 2 , which is close to 650 (recalibrated value,  An example of the 3D structure of skeletal muscle including nonuniform-size fibers and microvasculature. B) 3D structure of the microvasculature. C) 2D cross section of skeletal muscle with nonuniform size fibers and uniform capillary distribution (G3). D) 2D cross section of skeletal muscle with non-uniform size fibers and fiber-typedependent capillary distribution (G4). Fiber type I is shown in blue, type IIa in green and type IIb in yellow. Microvessels are shown in red. doi:10.1371/journal.pone.0044375.g002 number per fiber area including interstitial space, reported in [35]); (2) the radii of capillary, venule and artery segments are set as uniform in the vascular network (3, 6, 6 mm); and (3) the number of capillaries around a fiber (NCAF) is either specific (G2, G4) or non-specific (G1, G3) to fiber type, depending on the constructed geometry (i.e. uniform or heterogeneous fiber-typespecific microvascular distribution. For the uniform vascular network geometries (G1, G3), the average NCAF values for three fiber type groups are the same (approximately 3.24) and we also seek to minimize the variability of NCAF between individual fibers in the geometry ( Fig. 1 and 2). For fiber-type-specific vascular geometries (G2, G4), we follow experimental observations in rat skeletal muscle that muscle fibers with a high oxidative potential are generally associated with a denser capillary network [36,37], although not all studies agree with this finding [38]. In this type of fiber-type-specific geometry ( Fig. 1 and 2), fiber type I has the largest number of surrounding capillaries, followed by type IIa and then type IIb. NCAF values for each fiber varies from 1 to 6. In addition, we also designed the geometry with a significant degree of heterogeneity in each fiber type group.

Flow Module
We apply an in vivo hemorheological model [39] to calculate the blood flow rate (Q) and discharge hematocrit (H D ) in each of the blood vessel segments in the muscle, which incude terminal arterioles, capillaries and collecting venules. The vascular network is digitized as a set of nodes (vascular bifurcations) and a set of vessel segments linking those nodes together. The governing equations for volumetric blood flow rate and red blood cell flow rate at the j th node are: where Q ij and H Dij are the volumetric flow rate and discharge hematocrit in vessel segment ij, a cylindrical segment whose ends are the ith and jth nodes of the network. Using the governing equations combined with empirical equations describing red blood cell-plasma separation at vascular bifurcations [39], we obtained a set of nonlinear algebraic equations for all N segments and the equations were solved for pressure and hematocrit, from which flow in each segment is calculated.

Oxygen Module
We use a modified version of our previously published oxygen transport model to compute oxygen distribution in the tissue. In previous studies we assumed that muscle fibers and the interstitial space were a single tissue phase [18,23]. Here the model consists of three partial differential equations, governing: intravascular oxygen transport; oxygen diffusion in interstitial space; and oxygen diffusion and consumption inside fibers.
Local oxygen tension in the fibers, P f , is governed by free oxygen diffusion, myoglobin-facilitated diffusion, and oxygen consumption by myofibers: where D O2,f and D Mb are the diffusivities of oxygen (O 2 ) and myoglobin (Mb) inside the fiber; S Mb is the oxygen-myoglobin saturation; a f is the oxygen solubility in the fiber; C Mb bind is the binding capacity of myoglobin for oxygen; M c is the Michaelis-Menten-type maximal oxygen consumption rate; P crit is the critical P f at which oxygen consumption rate equals 50% of M c . S Mb is defined as P f /(P f +P 50,Mb ), assuming local binding equilibrium between oxygen and myoglobin, where P 50,Mb is the P f corresponding to 50% myoglobin saturation with oxygen. It should be noted that the myoglobin-facilitated diffusion formulation is simplified in that it does not consider diffusion of myoglobin as a separate molecular species and assumes that myoglobin and oxygen are locally in chemical equilibrium; however, this assumption is justified a posteriori by the small effect of myoglobin under steady-state conditions.
Local oxygen tension in the interstitial space, P i , is governed only by free oxygen diffusion (we assume negligible consumption in this space): Oxygen transport within the blood vessels, P b , is governed by hemoglobin binding and blood convection: Here S RBC Hb is the oxygen-hemoglobin saturation in blood vessel; a RBC and a pl are oxygen solubility in red blood cell and plasma, respectively; P b is the oxygen tension in blood plasma; n b is the mean blood velocity (n b = Q/(pR 2 )); H T is the vessel (tube) hematocrit calculated from blood flow model; C RBC bind is binding capacity of hemoglobin with oxygen; j is the distance along a vessel's longitudinal axis; J wall is the capillary wall flux. S RBC Hb is defined as P 50,Hb ) assuming binding equilibrium between oxygen and hemoglobin, where P 50,Hb is the P b at 50% hemoglobin oxygen saturation.
Continuity of oxygen flux at the interface between blood vessels and interstitial space, and between muscle fiber and interstitial space, yields: where n is the unit normal vector, and k 0 is the mass transfer coefficient. This system of nonlinear partial differential equations was solved using the finite difference method, with a grid size of 1 micron as described in [23].

Model parameters and simulation platform
Parameters used in the model are listed in Table 1 (fiber-type specific parameters) and Table 2 (non-fiber-type-specific parameters). Most of the model parameters were taken from experimental data for rat EDL; some were based on measurements in other muscle types; theoretical estimates were used when parameters are unavailable from literature. We use the multiscale modeling platform developed previously [24] to combine the three modules described above into an integrated model. This platform allows integration of different modules written in different programming languages and using different modeling methodologies. While the flow module remained unchanged from the previous version, the geometry and oxygen transport modules were updated based on the changes described above. The simulation experiments were run on an IBM cluster and each simulation was run on one node with 8 cores and 32 Gbyte memory. The oxygen transport module was implemented with OPENMP support and can be parallelized. JDK (Java Development Kit) 1.6.16 (Oracle, Redwood Shores, CA) is used as Java compiler, and Intel Fortran/C++ compiler suite 11.1 (Intel, Santa Clara, CA) is used as Fortran/C++ compiler.

Results
We computationally evaluated the contribution of fiber type composition and fiber type-specific parameters to spatial heterogeneity of tissue oxygenation. These parameters include: fiber type-dependent oxygen consumption (M c ), oxygen diffusivity in fiber (D O2 ), myoglobin concentration (C Mb ), fiber size, and number of capillaries around a fiber (NCAF). To understand the effects of each factor separately as well as their combined effects, we ran a series of computational experiments to examine the heterogeneity of oxygen distribution in skeletal muscle during exercise in different muscle geometries (G1-G4) and with scenarios varying a number of the other factors under consideration (S1-S5). The annotations for each simulation are summarized in Table 3. For example, G3S1 refers to G3 geometry (non-uniform fiber size and uniform capillary distribution) with scenario S1 (uniform M c , D O2 , and C Mb ). The corresponding microvascular and fiber structures are shown in Figs. 1C (G1), 1D (G2), 2C (G3) and 2D (G4), representing the same location in networks (z = 400 mm). The three-dimensional geometry of G1 is shown in Fig. 1A,B and that of G3 in Fig. 2A,B. The number of capillaries around a fiber (NCAF) for each of the simulation geometries is shown in Table S1. In the following sections, we report the mean and the coefficient of variation (CV) of the tissue oxygen profile as a measure of the heterogeneity in the tissue for each simulation.
It is worth noting that microvascular hemodynamics is a wellknown important determinant of tissue oxygenation and it is not the focus of this study. Our blood flow simulations suggest that these vascular networks in two geometries (G1 and G2) share similar distribution patterns of blood flow velocity and hematocrit, and have similar total blood volume flow rates (see Table S2 and Fig. S1). In addition, vascular networks G3 and G4 share similar microvascular blood flow and red blood cell distribution (Table S2 and Fig. S2).

Capillary distribution is an important determinant of oxygen distribution
For the G1 and G2 geometries (uniform fiber size; Fig. 1) we first computed the steady-state tissue oxygen distribution without consideration for fiber type composition, not varying O 2 consumption rates, O 2 diffusivities, and myoglobin concentrations in different fiber types. In other words, uniform fibers with the same parameters: volume-averaged M c (3.34, 6.68, 10.02610 24 mlO 2 ml 21 s 21 for low, moderate, high intensity exercise), D O2 (1.72610 25 mlO 2 ml 21 s 21 ), and C Mb (5.7610 23 mlO 2 ml 21 s 21 ) were used for every fiber in the simulation. In this case, with light intensity exercise, an increase in the heterogeneity of microvascular structure (G2 compared to G1), and therefore in the heterogeneity of oxygen supply, is predicted to lead to a slight increase in the heterogeneity of O 2 (CV G1S1 : 0.09, CV G2S1 : 0.10) and a slight lowering of mean tissue PO 2 level (P G1S1 : 28.9 mmHg, P G2S1 : 27.6 mmHg) (Fig. 3 and Table 4). In exercise of moderate or high intensity, the increase in oxygen consumption levels causes even greater heterogeneity. When capillaries are non-uniformly distributed around the fibers (Fig. 3, solid lines), the variation further increases; mean PO 2 levels decrease from 17.4 to 15.1 mmHg (moderate intensity exercise), and from 8.84 to 6.82 mmHg (high intensity exercise). These results indicate that capillary distribution in muscle tissue affects non-uniform oxygen distribution and heterogeneity.

Fiber type composition does not contribute significantly to oxygen distribution
We then computed steady-state tissue oxygen distribution in geometries G1 and G2 when fiber type composition and varying M c , D O2 , and C Mb in different fiber types (fiber type-specific parameters shown in Table 1) are considered; in other words, all fiber-type-specific parameters except for size. The simulation results with fiber type composition considered (scenario S2) were compared to uniform properties (scenario S1) ( Table 4). Fig. 3 shows the histogram of muscle fiber PO 2 in two muscle geometries (G1 and G2). Table 4 summarizes the key characteristics from all the simulation results of PO 2 distribution under different intensities of exercise in G1 and G2.
At all levels of exercise, the difference in oxygen distribution between the heterogeneous fiber-type specific parameters (S2) and the uniform fiber parameters (S1, control cases) is minimal (Fig. 3, red vs. black lines). This includes the mean, range and variance of oxygen levels and the portion of tissue that is hypoxic, i.e. with PO 2 ,1 or 2% (Table 4). For example, compared to the control cases with uniform fiber properties (i.e. G1S1, G2S1), the effect of fiber type composition on oxygen spatial heterogeneity and mean PO 2 level is small when muscles are stimulated with light or moderate exercise intensity (e.g., for light exercise, CV for the geometry with uniform capillary distribution, CV G1S2/S1 : ,0.09, CV for the geometry with fiber-type-specific capillary distribution, CV G2S2/S1 : ,0.10, P G1S2/S1 : ,28.9 mmHg, P G2S2/S1 : ,27.6 mmHg). Under high intensity exercise conditions, oxygen spatial heterogeneity changes slightly compared to control cases (i.e. G1S1, G2S1) (CV G1S2/S1 : 0.56 vs. 0.59; CV G2S2/S1 : 0.65 vs. 0.68) and mean PO 2 remains at the same levels (PO 2 G1S2/S1 : 8.8 vs. 8.5 mmHg, high;PO 2 G2S2/S1 : 6.8 vs. 6.6 mmHg, high). The oxygen distribution is much more dependent on capillary distribution (Fig. 3, solid vs. dotted lines) than on heterogeneity of fiber properties (Fig. 3, red vs. black lines).

Non-uniform fiber size significantly enhances the heterogeneity in oxygen distribution
To investigate the effect of fiber size on oxygen distribution in the EDL, we constructed geometries G3 and G4 (Fig. 2C,D), with non-uniform size fibers within the same dimension cuboid (35862336800 mm 3 ) and the same total fiber volume (79%) as G1 and G2 (Fig. 1C,1D). We first compared the tissue PO 2 profile in geometry G3S1 with heterogeneous fiber sizes to the profile in geometry G1S1 with homogeneous fiber sizes. Both of them use uniform fiber type properties (scenario S1, i.e., M c , D O2 , C Mb ) and have uniform capillary distribution. Our oxygen transport simulation results suggest that fiber size distribution in muscle geometry plays an important role in determining tissue oxygen profile and spatial heterogeneity. Figure 4 shows that PO 2 in heterogeneous fiber size geometry is much more broadly distributed compared to the control case (i.e., G1S1) with uniform size fiber under all exercise conditions (G3S1 vs. G1S1, black dashed lines in Fig. 4 vs. Fig. 3; 3D graphical representation shown in Fig. 5A We further examined the effects of fiber-type-specific properties on oxygen distribution in the non-uniform size fiber geometries (geometries G3 and G4, results in Figure 4 and Table 5). Conclusions similar to those described above can be drawn from these data. First, oxygen distribution is sensitive to local capillary distribution around the fibers. However, their effects are different from previous case with uniform size fiber distribution (Fig. 4, solid vs. dashed lines). With fiber-type-dependent capillary distribution (smaller fibers have more adjacent capillaries), the coefficients of variation of the tissue PO 2 profiles under light exercise conditions are slightly decreased compared to the tissue with uniform capillary distribution (CV: 0.14 vs. 0.12, light; 0.42 vs. 0.42, moderate), and their mean PO 2 levels shift to slightly lower values   Fig. 4 and three-dimensional graphical representation in Fig. 5B,C. There are only slight differences between black and red lines (including dashed and solid cases) at all exercise intensities, with one exception: under high intensity exercise and with fibertype-specific capillary distribution, heterogeneous fiber type composition increases spatial heterogeneity (CV: 0.81 vs. 1.00) and reduces its tissue PO 2 level (PO 2 : 5.93 vs. 3.22 mmHg). We further studied the distribution of average PO 2 for each fiber from three fiber type groups under various cases (Fig. 6). This shows that under light intensity exercise, the difference in average fiber PO 2 level among the three fiber type groups is low. With the increase of exercise intensity (to moderate and high intensity), fibers of type I have higher local oxygenation levels than types IIa and IIb. With fiber-type specific capillary distribution (G4 vs. G3), the oxygenation levels increase for oxidative fibers type I and type IIa, but decrease for type IIb. Therefore, at high intensity exercise conditions, the difference of PO 2 levels between fiber type IIb and fiber type I increased significantly, suggesting that local oxygen supply matches local demand well for oxidative fiber but not glycolytic fiber in exercise. Greater heterogeneity of average PO 2 level in each fiber type was also observed for higher exercise intensities and with capillary heterogeneity.
Fiber-type-specific M c and D O2 , but not C Mb , cause greater heterogeneity in oxygen distribution We next investigated the contribution of the heterogeneity of each fiber-type-dependent variable by itself (M c , D O2 , and C Mb ) to tissue PO 2 spatial heterogeneity. To do this, we performed simulations for the non-uniform fiber size geometry with alterations of the factor to be considered (scenarios S3-S5; see Table 3 for a detailed description). Control cases were simulated with uniform fiber properties for comparison. Figures 7, 8, 9 and Table 6 show the simulation results for the effects of M c , D O2 , C Mb heterogeneity in different fiber types. The results indicate that the heterogeneity of M c in different fiber types (Fig. 7,   Black lines are the simulation scenarios S1 with uniform fiber-type properties (M c , D O2 , C Mb ) and red lines are the simulation scenarios (S2 using fiber-type-specific parameters (except size) for different fiber types. doi:10.1371/journal.pone.0044375.g003 Table 4. Effects of fiber type properties on oxygen (PO 2 ) distribution.
Tissue Geometry* Exercise** P mean P min P max CV(P) %P,1*** %P,2  Table 5). Fig. 7 shows that the tissue oxygen levels are distributed in a much broader range. Incorporation of heterogeneity of oxygen diffusivity in different fiber types (Fig. 8, black vs. red lines) also affects tissue oxygen profiles. The PO 2 spatial heterogeneity increased significantly for both G3 (dashed lines) and G4 (solid lines) geometries. However, incorporation of heterogeneity of oxygen diffusivity in different fiber types causes a small increase in the mean PO 2 level (e.g., G3S1vs G3S4, 26.5 vs. 27.5 mmHg, light; 13.6 vs. 15.04 mmHg, moderate; 6.9 vs. 7.6 mmHg, high; Table 5). Lastly, comparison of simulation results for G3 and G4 geometries with incorporation of myoglobin concentration variation in different fiber types (Fig. 9, black vs. red lines) to control cases indicates that the contribution of C Mb heterogeneity to PO 2 spatial heterogeneity is small ( Fig. 9 and Table 6). Further  investigation of the sensitivity of PO 2 profile to myoglobinfacilitated diffusion by using a much smaller D Mb (2610 214 cm 2 /s, i.e. 7 orders of magnitude lower) shows that the tissue PO 2 profile with no myoglobin-facilitated diffusion does not deviate significantly from control cases at all exercise conditions (Fig. S3).

Discussion and Conclusions
Insufficient oxygen supply from microvascular networks under conditions of increased metabolic demand may lead to new capillary growth from pre-existing microvasculature during physiological conditions (e.g., exercise). Spatial heterogeneity is present in the angiogenic response within exercising muscles. The study of inherent geometric and cellular heterogeneities at the tissue level is critical to a fundamental understanding of organ function and may provide answers for drug resistance and suggestions for efficient therapeutic strategy. The sources of tissue heterogeneity can be attributable to complex microvascular structure and geometry (such as irregular vessel structure, vessel length and diameter distributions), and correspondingly heterogeneous hemodynamic variables such as blood flow rate, shear stress, and hematocrit [40,41]. These factors have been shown to play an important role in determining heterogeneous oxygen perfusion in tissue using both experimental and theoretical approaches [18,40,41]. The effect of the heterogeneity can be difficult to predict without a computational model; for example, inclusion of interstitial space as a separate non-oxygen-consuming volume decreased oxygen heterogeneity from 0.12 to 0.09 in light exercise. However, the heterogeneity in local capillary supply is often correlated with the difference of oxygen consumption capacity in different types of muscle fibers to maintain tissue homeostasis. In this study, we applied a multiscale oxygen transport model to quantitatively evaluate the impact of incorporation of fiber typespecific properties including fiber size on the oxygen distribution. Using rat EDL skeletal muscle as an example, we simulated the distribution of oxygen partial pressure under three exercise conditions (light, moderate, and high intensities) and studied the contributions of fiber type-dependent properties and fiber size to tissue oxygenation profiles.
Our results show that combined effects of fiber type properties other than size do not contribute significantly to tissue oxygen spatial heterogeneity and do not significantly change oxygen profiles compared to the control cases with uniform fiber distribution. By itself, variation of oxygen consumption ratio across the three fiber types has a significant effect (Fig. 7). But that Table 5. Effects of fiber type properties and fiber size on oxygen (PO 2 ) distribution.
Tissue Geometry* Exercise** P mean P min P max CV(P) %P,1*** %P,2  heterogeneity can be compensated by fiber type-specific oxygen diffusivity based on experimentally-measured physiological parameters (Fig. 8); the heterogeneity of myoglobin concentration in different fiber types does not significantly affect PO 2 profiles (Fig. 9). These three effects are thus collectively small. However, fiber size heterogeneity increases oxygen spatial heterogeneity significantly under all exercise conditions. These data are consistent with experimental observations that oxygen supply is closely correlated with fiber size [42]. In addition, it was shown in rat skeletal muscle that the degree of hypoxia-induced angiogenesis can be due to a fiber size effect [43]. Our results also suggest a small contribution of myoglobin-facilitated diffusion, and consequently differences of myoglobin concentration in different fiber types cause small deviations in oxygen distribution. Its minor contribution at steady state is consistent with other studies [23,44,45], although it might play an important role in regulating oxygen PO 2 at non-equilibrium states and high-intensity exercise conditions [46,47]. Our model has several limitations. First, current heterogeneous capillary distribution used in the geometry construction is based on the findings that capillary density is primary correlated with fiber oxidative capacity. Other studies have shown that capillary distribution is also a function of fiber size, muscle position Figure 7. Effects of variation of oxygen consumption rate in different fiber types on PO 2 histograms. Fiber PO 2 probability distribution profiles in muscle tissues under A) light intensity exercise; B) moderate intensity exercise; and C) high intensity exercise, for simulation cases G3S1, G3S3, G4S1 and G4S3. Dotted lines are simulation results from G3 geometry, solid lines are results from G4 geometry. Black lines are the simulation scenarios S1 and red lines are the simulation scenarios (S3) when M c is fiber-type-dependent while others use uniform values for different fiber types. doi:10.1371/journal.pone.0044375.g007 Figure 8. Effects of variation of oxygen diffusivity in different fiber types on PO 2 histograms. Fiber PO 2 probability distribution profiles in muscle tissues under A) light intensity exercise; B) moderate intensity exercise; and C) high intensity exercise, for simulation cases G3S1, G3S4, G4S1 and G4S4. Dotted lines are simulation results for G3 geometry, solid lines are results for G4 geometry. Black lines are the simulation scenarios S1 and red lines are simulations scenarios S4 when D O2 is fiber-type-dependent while others use uniform values for different fiber types. doi:10.1371/journal.pone.0044375.g008 (superficial or deep region), muscle type and animal species [48]. These factors can be considered in the future. Second, oxygen distribution in tissue is dependent on the balance of oxygen consumption, oxygen supply and oxygen transport. Thus the parameters related to these variables are critical for tissue PO 2 profiles. Most of these parameters used in our model are from rat EDL skeletal muscle and thus the conclusions we draw here may only hold true for the EDL. These hypotheses can be tested in other muscle types when muscle-specific parameters are available. Third, in these simulations, we hold blood flow constant (mean capillary velocity at approximately 1 mm/sec, [49]) in order to investigate and dissect the effect of the heterogeneity of fiber type. Fourth, depending on the level of exercise, different muscle fibers are likely recruited. While we do not consider these physiological effects in the current study, the model makes it possible to incorporate such regulatory changes in the future. Fifth, our oxygen transport model is limited to simulation of endurance exercise (primarily aerobic metabolism), but not resistance exercise. Sixth, in our model, exercise is simulated as an increase of oxygen consumption, without taking account the physical contractions and consequent fluid displacement and compartment deformation. Seventh, we do not include remodeling in this study. However, clearly oxygen perfusion is a dynamic process that is regulated by microvascular network geometry, which itself is regulated by oxygen-dependent vascular remodeling. Models that include such a feedback system have been proposed [50], and our current model can be extended to include remodeling. Finally, specific vascular networks may have specific heterogeneous structure that causes blood flow to be imbalanced. Our tissue geometry models are somewhat idealized and actual distribution of hypoxia may be sensitive to the exact placement of vessels -for example, a fiber with six rather than five capillaries nearby may skew the local distribution. However, we hypothesize that significant outliers would be removed by the physiological processes of angiogenesis and vessel regression. We generated two additional sets of tissue geometries and ran all the simulations for these to determine whether the results could be generalized. The conclusions presented here hold also for those additional geometries (see Additional Supplemental Material S2).
In summary, we conclude that differences in fiber size and capillary arrangement (but not other fiber type properties) contribute significantly to oxygen tension profiles in EDL skeletal Figure 9. Effects of variation of myoglobin concentration in different fiber types on PO 2 histograms. Fiber PO 2 probability distribution profiles in muscle tissues under A) light intensity exercise; B) moderate intensity exercise; and C) high intensity exercise, for simulation cases of G3S1, G3S5, G4S1 and G4S5. Dotted lines are simulation results for G3 geometry with uniform capillary distribution, solid lines are results for G4 geometry with fiber type-dependent distribution. Black lines are simulation scenarios S1 and red lines are simulations scenarios S5 when C Mb is fiber-type-dependent while others use uniform values for different fiber types. doi:10.1371/journal.pone.0044375.g009 Table 6. Effects of heterogeneities of M c , D O2 , and C Mb on oxygen (PO 2 ) distribution.
Tissue Geometry* Exercise** P mean P min P max CV(P) %P,1*** %P,2 muscle. Our extended heterogeneous-fiber oxygen transport model can be used to understand local oxygen distribution in other muscle types when fiber-type-specific parameters become available.

Supporting Information
Supplemental Material S1 Supporting information text and table S1-2. (DOC) Additional Supplemental Material S2 Additional supporting information for two additional sets of geometries.
(DOC) Figure S1 Distribution of blood flow rate and hematocrit distribution of the two computed geometries G1 and G2. Histograms of velocity (A, C) and hematocrit (B,D) distributions for the vascular network G1 and G2. G1 refers to the geometry with uniform fiber size and uniform capillary distribution; G2, uniform fiber size and fiber type-dependent capillary distribution.
(PNG) Figure S2 Distribution of blood flow rate and hematocrit distribution of the two computed geometries G3 and G4. Histograms of velocity (A, C) and hematocrit (B,D) distributions for the vascular networks G3 and G4. G3 refers to the geometry with non-uniform fiber size and uniform capillary distribution; G4, non-uniform fiber size and fiber type-dependent capillary distribution. (PNG) Figure S3 Effects of fiber-type specific myoglobin diffusivity on oxygen distribution in geometries of nonuniform-size fibers. Fiber PO2 under A) light intensity exercise, volume-averaged m c = 3.64610 24 mlO2 ml 21 s 21 ; B) moderate intensity exercise, volume-averaged m c = 6.68610 24 mlO2 ml 21 s 21 ; and C) high intensity exercise, volume-averaged m c = 1.02610 23 mlO2 ml 21 s 21 in muscle tissue of G3S1,G3S6,G4S1,G3S6. Dashed lines are simulations for G3 geometry with uniform capillary considered; solid lines are results from G4 geometry with fiber type-dependent capillary distribution. Black lines are the simulation cases (S1) with uniform fibertype properties (m c , D O2 ,C Mb ) and red lines are simulations when D Mb uses low value 3610 214 instead of 3610 27 (S6). (PNG)