The formation and functional consequences of heterogeneous mitochondrial distributions in skeletal muscle

Mitochondrial density and distribution in skeletal muscle is a plastic property subject to various metabolic factors (Tyler and Sidell, 1984; Howald et al., 1985; Kayar et al., 1986; Hoppeler and Vogt, 2001; Chilibeck et al., 2002; Fluck and Hoppeler, 2003; Guderley et al., 2003; Korzeniewski and Zoladz, 2003; Nyack et al., 2007; Holloszy, 2008; Hardy et al., 2009). The non-uniform distribution of mitochondria in muscle fibers presumably reflects reaction-diffusion constraints associated with aerobic metabolism (Kinsey et al., 2011). Oxygen diffuses from the capillaries into and through the cells to the mitochondria, where ATP is formed, and the ATP subsequently diffuses from mitochondria to cellular ATPases. We have shown that increases in fiber size during hypertrophic growth lead to shifts in the distribution of mitochondria toward the fiber periphery, which reduces diffusion distances for oxygen but increases distances for ATP diffusion (Boyle et al., 2003; Nyack et al., 2007; Hardy et al., 2009; Hardy et al., 2010). Mathematical modeling showed that intracellular oxygen gradients likely play an important role in dictating these distributional shifts, and that the reorganization of mitochondria helps the fibers avoid diffusion limitation (Hardy et al., 2009; Pathi et al., 2011). However, the mechanisms by which mitochondrial distribution is controlled are unknown.

Recently, there have been major advances in our understanding of the molecular pathways that lead to mitochondrial biogenesis and death. Mitochondrial biogenesis entails an increase in mitochondrial mass that is often followed by mitochondrial fission (Terman et al., 2010). At the cellular level, one of the primary mechanisms by which this occurs is through AMP kinase (AMPK). AMPK is considered to be a cellular fuel gauge that is activated when there is an energy deficit, which is specifically manifested as an increase in the [AMP]/[ATP] ratio in mammals (Hardie et al., 2003) and in fish (Jibb and Richards, 2008). AMPK leads to the activation of, among other things, peroxisome proliferator-activated receptor gamma coactivator 1-alpha (PGC-1α), which is considered a master regulator of mitochondrial biogenesis (Puigserver et al., 1998). PGC-1α in mammals induces genes necessary for biogenesis that are present in both the nuclear and mitochondrial genome (Terman et al., 2010). There are several pathways that control biogenesis, but a common linkage is the need to regulate mitochondrial volume density to match ATP supply with demand (Hood et al., 2006).

Mitochondrial degradation was long thought to be due to random autophagic events, but growing evidence indicates that damaged mitochondria are specifically removed in a process referred to as targeted mitophagy (Lemasters, 2005; Kim et al., 2007; Hyde et al., 2010; Gottlieb and Carreira, 2010). Several causative agents have been suggested, with the most widely accepted being reactive oxygen species (ROS) (Terman et al., 2010; Tolkovsky, 2009). ROS are byproducts of oxidative phosphorylation that can damage mitochondria and induce inner membrane permeability transitions that, in turn, stimulate targeted mitophagy or mitoptosis (Wallace, 2005; Kim et al., 2007). Mitophagy therefore serves to eliminate the population of mitochondria producing the most ROS, reducing the oxidative burden and increasing cellular longevity (Kim et al., 2007; Yen and Klionsky, 2008; Zhang et al., 2008; He and Klionsky, 2009). Because mitochondrial ROS production is increased during hypoxia exposure, it is not surprising that low O₂ induces targeted mitophagy (Band et al., 2009). In fact, the hypoxia-inducible factor-1 (HIF-1) pathway that largely regulates the cellular response to hypoxia induces genes that are crucial to mitophagy, and this response is thought to protect the cells from further ROS damage (Lemasters, 2007; Zhang et al., 2008).

Although these and other studies provide a molecular basis for changes in mitochondrial volume density, they do not explain how the heterogeneous distribution of mitochondria arises and is altered. The implication is that mitochondrial heterogeneity is controlled by spatial variation in the signaling pathways highlighted above. Thus, these mitochondrial distributions arise as a consequence of the net effects of the processes governing mitochondrial biogenesis and death. Mitochondrial biogenesis is dependent on gene expression both in the mitochondrion and in the nucleus, as well as the subsequent diffusive transport of nuclear gene products to the mitochondria that will undergo fission (Wallace 1999; Chan, 2006; Hood et al., 2006; Gottlieb and Carreira, 2010; Terman et al., 2010). Thus, the spatial information encoded by signaling molecules that promote biogenesis, such as elevated AMP in the vicinity of a mitochondrion, would seem to be lost as the signaling cascade is propagated through the nucleus. In support of this view, shifts in mitochondrial distribution toward the fiber periphery during hypertrophic fiber growth in a variety of comparative models appear to be independent of opposing shifts in nuclear distribution toward the fiber center (Nyack et al., 2007; Hardy et al., 2009; Hardy et al., 2010; Priester et al., 2011), and experiments with isolated cells indicate that mitochondrial DNA replication is independent of proximity to nuclei (Magnusson et al., 2003). Moreover, reaction-diffusion models indicate that spatial gradients in high-energy phosphate molecules like AMP are very slight in muscle (Jimenez et al., 2008; Dasika et al., 2011; Pathi et al., 2011), and we know of no biogenesis signals that vary spatially across the radius of the fiber. In contrast, as mentioned above, mitochondrial degradation is a regionally variable, selective process, whereby damaged mitochondria are specifically targeted and removed. For instance, a mitochondrion present at the muscle fiber core that is frequently exposed to hypoxic conditions may suffer increased damage by ROS and be targeted for removal as a cytoprotective response and, unlike biogenesis signals such as AMP, there may be steep oxygen gradients across the radius of fibers from mammals (Groebe et al., 2010; Schroeder et al., 2010; Pathi et al., 2011), fishes (Egginton et al., 2002) and crustaceans (Jimenez et al., 2008). Thus, although it is easy to envision how the balance between mitochondrial biogenesis and degradation governs mitochondrial volume density, it appears that only pathways of mitochondrial death, coupled with spatial gradients in key molecules such as oxygen, provide a means of generating a heterogeneous mitochondrial distribution.

The present study tested the hypotheses that, in a diversity of fibers: (1) mitochondrial death alone, with dependence on the prevailing oxygen concentration, can adequately explain observed skeletal muscle mitochondrial distribution and plasticity, and (2) these non-homogeneous distributions alleviate diffusion constraints, raise the cellular energy state and increase permissible rates of aerobic metabolism. We tested these hypotheses in red and white muscle from dolphinfish (Coryphaena hippurus Linnaeus 1758), white muscle that undergoes extreme hypertrophic growth from black sea bass [Centropristis striata (Linnaeus 1758)] and four fiber types from two muscles in mouse (Mus musculus Linnaeus 1758). Light and electron microscopy, coupled with stereological analyses, were used to characterize fiber size, capillarity, and mitochondrial density and distribution. These data were then compared with mitochondrial distributions that were predicted from a coupled reaction-diffusion/cellular automata (CA) model of mitochondrial function. This allowed us to link a reaction-diffusion model of aerobic metabolism for a population of mitochondria to CA rules that govern the mitochondrial life cycle based on general principles of the signaling pathways highlighted above. This simple approach was effective at predicting mitochondrial distribution, and the model output indicated that these non-homogeneous distributions of mitochondria facilitate higher rates of aerobic metabolism than fibers with a homogeneous mitochondrial distribution.

Fiber size, capillaries around a fiber (CAF), and mitochondrial volume density and distribution were examined in a diversity of skeletal muscle fiber types from three species. Red (aerobic) and white (anaerobic) muscle, which differ dramatically in mitochondrial volume density, was analyzed in dolphinfish (C. hippurus), a highly active pelagic species that attains a large adult body mass. White muscle that undergoes extreme hypertrophic fiber growth leading to a large increase in fiber size during animal growth was examined in juvenile (small fibers) and adult (large fibers) black sea bass (C. striata), a relatively inactive benthic species. Finally, four distinct fiber types (1a, 2a, 2x and 2b) that differ in mitochondrial volume density, capillarity and fiber size were examined in skeletal muscle of mice (M. musculus).

All data from epaxial white muscle from juvenile and adult black sea bass were taken from Nyack et al. (Nyack et al., 2007), and these data were collected using the same procedures described herein. A 7.7-kg adult dolphinfish was caught using a hook and line 60 miles off the coast of Wilmington, NC, immediately placed on ice, and weighed and dissected several hours later upon return to land. Epaxial white and lateral red muscle was harvested from beneath the anterior portion of the dorsal fin. Inbred strains of BALB/cByJ (C) and C57BL/6J (B6) mice obtained from The Jackson Laboratory (Bar Harbor, ME, USA) were maintained in a breeding colony and crossed to produce F1 hybrids, which were used in the present study. Individual breeding pairs and their offspring were housed in polycarbonate cages in a laminar flow hood at 23±1°C with a 12:12 h light:dark cycle. Animals were provided Agway 3000 Mouse Chow (Agway Inc., Syracuse, NY, USA) and water ad libitum. Breeding pairs were observed weekly and once the female was judged to be pregnant the male was removed. Litters were weaned at 17 days and separated into same-sex sibling groups with no more than four mice per cage. Mice were euthanized with an overdose of CO₂ and the soleus (SOL) and extensor digitorum longus (EDL) hindlimb muscles were removed and tied to a small wooden rod at resting length. All procedures were approved by the University of North Carolina Wilmington Institutional Animal Care and Use Committee.

All tissues to be used for transmission electron microscopy (TEM) to measure mitochondrial distribution were immediately fixed in Karnovsky’s fixative (3.0% glutaraldehyde, 3.0% paraformaldehyde in 0.2 mol l^–1 sodium cacodylate buffer, pH 7.6) for a minimum of 24 h at room temperature. The tissues were then rinsed in 0.2 mol l^–1 sodium cacodylate, pH 7.6, followed by secondary fixation in 0.5% osmium tetroxide in 0.2 mol l^–1 sodium cacodylate buffer, pH 7.6. Samples were rinsed and dehydrated after which they were embedded in Spurr epoxy resin (Electron Microscopy Sciences, Hatfield, PA, USA). Embedded muscle samples were cut into a trapezoid with beveled edges. Muscle samples were sectioned at 90 nm with a diamond knife on a Reichert Ultracut E ultramicrotome (Leica Microsystems Inc., Bannockburn, IL, USA). Series of sections were collected using a systematic random sampling method (Howard and Reed, 1998) to ensure complete representation of the mitochondrial distribution throughout the muscle. Five sections were collected from a random starting point and then a distance of 500 nm was skipped before collecting another five sections. Sections were collected and mounted on Formvar-coated (0.25 g Formvar in 100 ml ethylene dichloride) or Collodion-coated (1.0% Collodion in amyl acetate), 200 μm mesh, high-transmission copper grids. In this way the tissue samples were sectioned so that five grids, with five sections per grid, were obtained. Samples were stained with 2% uranyl acetate in 50% ethyl alcohol, and then stained with Reynolds’ lead citrate (Reynolds, 1963) for 15 min each.

Sections were examined on a Philips CM 12 TEM (Philips Research, Briarcliff Manor, NY, USA), operated at 80 kV, and photographed using a 3.25×4.25 inch plate camera. Mitochondria in fish white muscle are very sparse, so in fishes mitochondria were characterized as either subsarcolemmal (SS), mitochondria that were either in direct contact with the sarcolemma or part of a cluster that was in contact with the sarcolemma, or as intermyofibrillar (IM), mitochondria that were not in contact with the sarcolemma. In mouse muscle, the much higher mitochondrial volume density and capillary number permitted a more detailed analysis that included measurement of SS mitochondrial volume density as well as the radial profile of IM mitochondria from a capillary to the core of the fiber.

For fish muscle, micrographs were collected by a systematic random sampling method where sections were first observed on a TEM grid under low magnification of either 710×, 530× or 400×, after which the section was zoomed in to a magnification of 3000× and the first micrograph was taken and utilized as a random starting point for two subsequent micrographs. In mouse muscle, images were collected at a magnification of 1850×. The field of view was then moved to the right by a distance equal to three times the field of view and then a second micrograph was taken. If any position was obstructed by a TEM grid bar or was no longer on a section, the previous field of view was reestablished and a new field of view was established three fields down from this position. If this position was also not appropriate, the previous field of view was reestablished and a new field of view was established two fields down and two to the right of the previous position. If this also did not produce a suitable field of view, the grid was repositioned under low magnification and a new arbitrary starting point was utilized. This process was repeated until three micrographs were taken from each grid.

Negatives were developed and images were then digitized using a Microtek Scanmaker 4 scanner (Microtek Inc., Carson, CA, USA). Digitized micrographs were then processed using Adobe Photoshop version 7.0 (Adobe Systems, San Jose, CA, USA). To calculate mitochondrial volume density (V_vMT), digital images were analyzed with ImagePro Plus (version 4.1.09; Media Cybernetics, San Diego, CA, USA) using a stereological point-counting method (Howard and Reed, 1998). A point grid was applied to the images and the points falling on extracellular space were subtracted from the total number of points. Points that fell on SS or IM mitochondria were divided by the number of points that fell on the sarcoplasm to yield volume density of each population, and the IM and SS mitochondrial points were added to obtain the total volume density (Nyack et al., 2007).

In mouse muscle, V_vMT and distributions were stereologically analyzed using a modification of the methods described by Kayar et al. (Kayar et al., 1986). To determine the radial distribution of IM mitochondria, a set of three grids each with an area of 7 μm² was applied to all micrographs. One of the capillaries surrounding the fiber was randomly selected and a transect from the fiber centroid to the center of the capillary was created. One grid was placed at the fiber centroid, which was calculated in ImagePro Plus after the fiber periphery was manually traced, one was placed midway between the fiber centroid and the capillary, and one was placed at the fiber periphery adjacent to the capillary. Point counting was then applied as above to calculate volume density at each radial position (excluding SS mitochondria). The volume density of SS mitochondria was determined by selecting the SS mitochondrial clusters using the threshold function in ImagePro Plus (as well as manual selection when necessary) and dividing this area by the total fiber area.

Fiber size, fiber type and CAF in C. hippurus were measured using light microscopy because the relatively large size and low capillarity of the white fibers made measurement of these properties difficult using TEM. Tissue was fixed for 24 h in 4% paraformaldehyde in phosphate buffered saline (PBS). Fixed tissue was rinsed in a 5% sucrose solution followed by an overnight infiltration in 30% sucrose. Tissue was trimmed, embedded in Tissue-Tek OCT compound (Sakura Finetek Inc., Torrance, CA, USA) and frozen for 20–40 s in isopentane (2-methylbutane; Fisher Scientific, Pittsburgh, PA, USA) cooled with liquid nitrogen. Frozen tissue was cut in 30 μm sections at –19°C with a Reichert-Jung Cryocut 1800 cryostat (Leica Microsystems). Sections were placed on Superfrost PLUS slides (Fisher Scientific) and allowed to air dry for approximately 15 min. Sections to be examined from each tissue preparation were selected using a systematic random sampling method (Howard and Reed, 1998) as described previously (Boyle et al., 2003; Hardy et al., 2009; Priester et al., 2011).

Sections were rinsed for 5 min in PBS, stained with 5 μg ml^–1 fluorescein-labeled Griffonia simplicifolia lectin (GSL-I) in PBS at room temperature for 1–2 h, and rinsed with PBS twice for 15 min. Sections were then incubated with 25 μg ml^–1 wheat germ agglutinin lectin (WGA) labeled with Alexa Fluor 594 in PBS for 30 min at room temperature, rinsed with PBS twice for 15 min, and coverslipped with 1:9 tris-glycerol. WGA binds carbohydrates on the basement membrane of the sarcolemma, whereas GSL-I binds endothelial cells in capillaries. Sections were examined with an Olympus FV1000 laser scanning confocal microscope (Olympus, Center Valley, PA, USA). One-micrometer-thick optical slices were simultaneously collected on three channels to observe the signal from labels on GSL-I and WGA, as well as the differential interference contrast signal. Stacks of 20–25 images of fiber cross-sections were collected and viewed using Olympus Fluoview v1.6a software. Images of muscle fibers were outlined with an Intuos 3 tablet (Wacom Co., Ltd., Vancouver, WA, USA) in Adobe Photoshop (v. 7.0) and resulting polygons were analyzed with ImagePro Plus (v. 6.1.0.346). The mean diameter was determined by taking the average of the distances across the cell through the centroid in 2 deg increments around the perimeter of the cell, and capillaries were counted to yield CAF. In mouse muscle, the same analytical approach was used to measure fiber diameter and CAF, but was applied to the TEM images.

In the F1 mice used in the present study, >99% of the EDL cross-sectional area is comprised of two fiber types: fast-twitch oxidative 2x fibers (≈25%) and fast-twitch glycolytic 2b fibers (≈75%), where the mean fiber diameter of 2x fibers is 32 μm and that of 2b fibers is 53 μm (Luedeke et al., 2004). Therefore, in the TEM images, fibers that were <40 μm in diameter and had a V_vMT >0.2 were considered to be 2x fibers, and those that were >40 μm in diameter and had a V_vMT <0.2 were considered to be 2b fibers. In the SOL from the F1 mice, >95% of the cross-sectional area is comprised of slow-twitch oxidative 1a and fast-twitch oxidative 2a fibers, each of which comprises an equal volume of the fiber (Luedeke et al., 2004). These fibers do not differ in mean diameter, but the 1a fibers have a higher V_vMT. Therefore, fibers in the SOL with a V_vMT >0.2 were considered to be type 1a, whereas those with a lower value were considered to be type 2a. Although this approach for characterizing fiber types is not definitive, it is likely to be highly accurate because there are only two dominant fiber types in each muscle and they are very different in either fiber size and V_vMT (or both). However, even if there are some mischaracterized fibers, our primary interest is in evaluating fibers that differ in size, capillarity and aerobic capacity, and this approach effectively does this.

Both our previous (Pathi et al., 2011) and the present work use a two-dimensional hybrid CA reaction-diffusion model to analyze the behavior of individual mitochondria in response to both the local environment and the metabolic state of the cell. The reactions and model equations used are described in detail in our previous work (Pathi et al., 2011). In brief, each fiber is represented as a circle with capillaries placed at the periphery (the number of capillaries varies with type of fiber), and the fiber is divided into a grid where each square in the grid represents either a mitochondrion or contains only cytoplasm. Capillaries were placed on the periphery of the cell such that 5 μm of each individual capillary surface was adjacent to the muscle sarcolemma (Mathieu-Costello et al., 1991). The diameter of each mitochondrion was assumed to be 1 μm and equivalent to a single grid space (i.e. 1×1 μm). For small fibers (e.g. juvenile black sea bass white fibers and EDL 2x fibers), the size of the grid spaces (and therefore the size of the mitochondria) was reduced to 0.7 μm. Use of smaller-sized mitochondria allowed a large enough number of mitochondria to realistically represent the population in these smaller fibers. Because mitochondria in skeletal muscle are typically 0.5–1 μm in diameter, both sizes used were reasonable, but the larger grid size was used for the large fibers to reduce computation time. For the very large fibers in adult black sea bass, only a quadrant of the cell was modeled in order to reduce computational time, whereas for all other cases the whole cell was considered. An initial distribution of mitochondria was generated by randomly placing mitochondria in the grids using a random number generator and a counter for the total mitochondrial number until the desired mitochondrial volume density was achieved for the specified fiber type. Initial metabolite concentrations inside the fiber and the oxygen concentration in the boundary condition (capillary) (35 μmol l^–1) were set. The parameters for the dolphinfish and black sea bass were taken from Nyack et al. (Nyack et al., 2007) and the corresponding parameters for mouse muscle fibers from Overton et al. (Overton et al., 2009).

Oxygen consumption and ATP production at each mitochondrion were based on a rate law derived from the model of oxidative phosphorylation developed by Beard (Beard, 2005). The full model of Beard (Beard, 2005) is too computationally intensive to apply to our population of mitochondria. Therefore, Dasika et al. (Dasika et al., 2011) developed a rate law that closely approximates the full model of Beard (Beard, 2005) that expresses ATP production in the mitochondria as a function of the concentration of oxygen, ADP and inorganic phosphate (P_i). This rate law was used to describe reactions in the mitochondria. All the mitochondria in the present study from the different species i.e. mammalian and fish (both red and white) muscle are treated the same and it is assumed that mitochondrial respiration per mitochondrial volume is the same for all fiber types. Therefore, mitochondria from these species used the same rate law for ATP production and the same maximum rate of ATP production (V˙_max) per mitochondrion. This is based on prior evaluation of maximal (state 3) respiration in isolated mitochondria per volume of mitochondria (Schwerzmann et al., 1989; Burpee et al., 2010). Schwerzmann et al. (Schwerzmann et al., 1989) evaluated respiration per mitochondrial volume in mammalian muscle, and Burpee et al. (Burpee et al., 2010) examined both the body mass and species dependence of state 3 respiration per volume of mitochondria in white muscle of several species of fishes. When rates were adjusted to physiological temperature, the mitochondrial O₂ consumption in mammals was the same as in fishes, and in fishes this was independent of body mass (Burpee et al., 2010). The ATP consumption by the ATPases in the cytoplasmic spaces between mitochondria was modeled as a single Michaelis–Menten expression that is dependent on ATP concentration. Creatine kinase (CK) also was present in the cytoplasmic spaces and was described by the expression of Morrison and James (Morrison and James, 1965).

The CA rules that were used were based on the hypothesis that mitochondrial biogenesis signals are uniform across the fiber because of the lack of evidence for gradients across the radius of the fiber in AMP or other known biogenesis signals, and because if spatial variation of an appropriate scale does exist, at least some of that information is lost when that signal is propagated to the nucleus and then back to the mitochondria via diffusion. This is supported by shifts in mitochondrial distribution toward the fiber periphery during hypertrophic fiber growth that appear to be independent of nuclear distribution (Nyack et al., 2007; Hardy et al., 2009; Hardy et al., 2010; Priester et al., 2011), and experiments indicating that mitochondrial DNA replication is independent of proximity to nuclei (Magnusson et al., 2003). In contrast, mitochondrial death is considered to be spatially variable and more likely in regions with lower oxygen (Lemasters, 2005; Kim et al., 2007; Hyde et al., 2010; Gottlieb and Carreira, 2010; Pathi et al., 2011).

The model was solved by first generating an initial random mitochondrial distribution with known V_vMT in a fiber with a specific size and number of capillaries. Reaction rate constants were set based on the presence/absence of mitochondria. Grid spaces containing mitochondria consume oxygen, ADP and P_i, produce ATP, and had their ATPase reaction rate constants set to zero. Grid spaces lacking mitochondria consume ATP, produce ADP and P_i (cellular ATPases), and had their mitochondrial reaction rate constants set to zero. In each CA iteration, all of the grid spaces were evaluated to determine whether a mitochondrion lives or dies (if the space had a mitochondrion) or whether a new mitochondrion was formed (if it was a cytoplasmic grid space adjacent to an existing mitochondrion). For each CA iteration, the simulation was run until the concentrations of all metabolites reached a temporal steady state, because the time scale for metabolite diffusion (on the order of seconds) is much shorter than the time scale for the mitochondrial life cycle (on the order of days). A series of CA iterations were run in sequence until a steady-state mitochondrial pattern was reached.

The system of second-order parabolic partial differential equations describing the metabolite species concentrations (Eqn 1) were solved using an alternating direction implicit finite difference scheme (Peaceman and Rachford, 1955; Madzvamuse and Maini, 2007; Yang et al., 2005). The model was implemented using MATLAB release 7.9 (The Math Works, Inc., Natick, MA, USA).

The mitochondrial distributions in the model simulations were compared with experimental data. Unlike the experimental data, where the V_vMT was derived from sampling a three-dimensional muscle preparation, the simulations are two-dimensional. Therefore, the mitochondrial densities in the model are actually fractional areas rather than volume densities. However, we will refer to the simulation data as V_vMT for consistency in terminology with the experimental data. The total IM V_vMT was calculated from the area of IM mitochondria divided by the total fiber area. To calculate the radial distribution of IM mitochondria (for comparison with the mouse experimental data), the V_vMT was calculated by determining the area of mitochondria present in 4-μm-wide bands of the fiber starting from the center of the fiber and moving to the fiber periphery normalized by the area of the bands containing those mitochondria. The SS V_vMT for the different fiber types was calculated from the model following the same procedure used for the experimental analysis, where mitochondria touching the sarcolemma or touching mitochondria in contact with the sarcolemma were counted as SS mitochondria. The sum of the area of all of these mitochondria normalized by the area of the cell gives the V_vMT of the SS mitochondria in the cell.

Table 1 summarizes the measurements of fiber size, CAF and V_vMT for the dolphinfish red and white muscle, black sea bass white muscle from juveniles and adults, and mouse 2x and 2b fibers (EDL) and 1a and 2a fibers (SOL).

Fig. 1 shows an example of the general modeling approach and some of the model outputs, where the initial random distribution of mitochondria (Fig. 1A) changes during the CA iterations as mitochondria shift toward the capillary boundaries until a steady-state pattern is established (Fig. 1B). A decrease in the oxygen gradient is observed after redistribution (Fig. 1D) in comparison to before redistribution (Fig. 1C). A noticeable increase was observed in the PCr concentration profile after the CA iterations (Fig. 1F) compared with the profile prior to CA iterations (Fig. 1E), indicative of a higher energy state for the fiber when mitochondria are distributed closer to capillaries. The ADP concentrations did not change noticeably after redistribution (Fig. 1G,H). The absence of gradients shown for PCr and ADP is consistent with results for all of the high-energy phosphate molecules, including ATP, which is in contrast to the steep gradients for oxygen. Thus, the peripheral distribution of mitochondria helps lessen oxygen gradients without incurring gradients for high-energy phosphate molecules that must diffuse greater distances when mitochondria are non-uniformly distributed.

Fig. 1I shows the changes in concentrations over time of several metabolites during the first CA iteration, where mitochondria are randomly distributed (Fig. 1A). The concentrations change from the initial conditions (based on resting muscle values) until a new steady state is reached at the elevated rates of aerobic ATP demand, similar to a rest-to-work transition. Each CA iteration is of sufficient duration to achieve steady state, as in Fig. 1I. After a number of CA iterations, the mitochondrial distribution shifts and reaches a new steady state, which also leads to a new steady state in the metabolite concentrations. Fig. 1J shows an example of the changes in concentrations during 10 CA iterations, where the increase in PCr and ATP and the decrease in ADP indicate an elevated cellular energy state.

As in other biological models using a hybrid CA and reaction-diffusion approach (Ferreira et al., 2002; Alarcon et al., 2003), the metabolites are assumed to be in a temporal steady-state distribution throughout the domain for each CA iteration because the time scale for metabolite diffusion is known to be much shorter than the time scale for the mitochondrial life cycle. Based on a half-life of mitochondria in muscle of approximately 7 days, each CA iteration corresponds to approximately 1 week (Booth and Holloszy, 1977; Hepple, 2010). Based on this, examination of Fig. 1J indicates that a new steady state in metabolite concentrations (and therefore in mitochondrial distribution) would be reached in approximately 8 weeks.

In the example in Fig. 1 and in other simulations below, the V˙_max of the ATPase reaction is adjusted so that rates of aerobic ATP turnover are maintained within the physiological range. This is particularly apparent in the relative PCr concentration, where it is clear that PCr is only partially depleted (to approximately 40% of the resting value during steady-state ATP demand after the CA iterations are complete), and in the relative ATP concentration, which is not depleted at all (Fig. 1J). This approach for all simulations is consistent with experimental results of relative PCr and ATP depletion determined using nuclear magnetic resonance (NMR) in mammalian muscle under conditions of elevated steady-state ATP demand (Kushmerick and Meyer, 1985; Kushmerick et al., 1992).

Model outputs for the dolphinfish red and white fibers are presented in Fig. 2, which illustrates steady-state mitochondrial distributions, corresponding TEM images and model predictions of SS and IM mitochondrial volume densities. Dolphinfish red fibers operate at relatively high aerobic ATP turnover rates during steady-state contraction, leading to a substantial mitochondrial spatial shift toward the capillaries (Fig. 2A). Comparison of the simulation to the TEM image reveals that the simulations capture the clustering behavior of mitochondria near the capillary, as well as other characteristics such as the occurrence of linear arrays of IM mitochondria. The simulation also predicts an IM and SS mitochondrial volume density that is similar to that observed experimentally (Fig. 2B).

Dolphinfish white fibers have very low V_vMT and generally operate at very low aerobic ATPase rates. Because of the low ATPase rate and relatively modest fiber size in white muscle of dolphinfish, the redistribution of mitochondria toward the capillary is relatively slight (Fig. 2C). Nevertheless, the simulations again capture the basic pattern observed in the TEM images, including SS mitochondria near capillaries and small clusters of a few IM mitochondria (Fig. 2C). As in red muscle, the predicted SS and IM V_vMT are in good agreement with the observations (Fig. 2D).

Fig. 3 represents the corresponding results for hypertrophic growth in the black sea bass white fibers (juvenile and adult fibers). Our mathematical model results were compared with previously observed shifts in mitochondrial distribution that occur during hypertrophic fiber growth in the white muscle of black sea bass (Nyack et al., 2007). Juvenile black sea bass white muscle operates at low R_ATPase of approximately 1 mmol l^–1 min^–1 (Nyack et al., 2007). Because of the low aerobic ATPase rate, small size (35 μm in diameter) and low V_vMT, oxygen in the cell is plentiful across the fiber. Hence, the steady-state distribution of mitochondria after applying the CA rules is only modestly shifted toward the capillaries (Fig. 3A). However, adult black sea bass white muscle operates at a similarly low aerobic ATPase rate of 0.9 mmol l^–1 min^–1, but the fibers grow to an average diameter of 240 μm (Nyack et al., 2007). Because of the larger fiber size, even at such low ATPase rates depletion of oxygen in the cell is much higher than that of the juvenile black sea bass fibers. Hence, a dramatic redistribution of mitochondria is observed in this case in comparison to the juvenile muscle (Fig. 3C). In the transition from juvenile to adult, the simulations compare favorably to the TEM images (Fig. 3A,C), and the predicted SS and IM V_vMT values are similar to the experimentally measured values (Fig. 3B,D).

The spatial mitochondrial distribution patterns were also analyzed in EDL and SOL muscles from mouse. Fig. 4 illustrates the predicted steady-state mitochondrial distributions, TEM images and model predictions of SS and IM V_vMT of EDL (i.e. 2b and 2x fibers) and Fig. 5 represents the corresponding results for SOL (i.e. 1a and 2a fibers).

The EDL 2b fibers have the lowest V_vMT and were the largest in size amongst all the mouse fibers considered in the present study. In contrast, the EDL 2x fibers are the smallest of the fibers containing the largest volume density of mitochondria. The SOL 2a and 1a fibers were both nearly of the same size, but the 1a fibers contained a larger volume density of mitochondria. All these mouse fibers demonstrated varying amounts of shifts in mitochondrial populations towards the fiber capillaries following CA iterations that were consistent with the TEM images (Fig. 4A,C and Fig. 5A,C).

Comparisons of the radial distribution of IM V_vMT and the SS V_vMT for the four mouse muscle fiber types are shown in Fig. 4B,D and Fig. 5B,D. The spatial gradients in IM V_vMT seen experimentally were also apparent in the simulations, and the model predicted radial profiles well for EDL 2b and 2x and SOL 1a and 2a fibers (Fig. 4B,D, Fig. 5B,D). The model outputs corresponding to the SS V_vMT after redistribution also agree better with experimental observations than the initial distributions prior to the CA iterations (EDL 2b and 2x, SOL 1a and 2a (Fig. 4B,D, Fig. 5B,D).

In all of the fibers analyzed, the heterogeneous distribution of mitochondria led to higher average O₂ and PCr concentrations compared with a random distribution similar to that shown in Fig. 1. The effects of the non-uniform mitochondrial distributions observed experimentally and predicted by the model on R_ATPase, η and free energy of ATP hydrolysis (ΔG) are shown in Fig. 6. Increases in R_ATPase after redistribution were observed in all the fiber types (Fig. 6A,B). Similar trends were observed in the η values, indicating a reduction in diffusion limitation of aerobic metabolism (Fig. 6C,D). Moreover, the fibers had a lower (more negative) ΔG following the CA iterations, indicating that the fibers are able to maintain a higher energy state during elevated rates of aerobic ATP demand (Fig. 6E,F).

The major findings of the present study were that the spatial heterogeneity in mitochondrial distribution observed in different skeletal muscle fibers could be predicted from mitochondrial death governed by the local oxygen concentration. Analysis of the mitochondrial redistributions in the different fiber types showed an increase in the average concentrations of essential metabolites such as oxygen, PCr and ATP, and a decrease in the concentration of ADP. This indicates that the non-uniform mitochondrial distribution observed experimentally and predicted with our model benefits the muscle fiber by increasing the cellular energy status (more negative ΔG), increasing sustainable rates of aerobic ATP demand (higher R_ATPase) and decreasing diffusion constraints (higher η).

The hybrid CA reaction-diffusion modeling approach used here has been successfully used to predict dynamics associated with tissue engineering (Chung et al., 2010), tumor growth (Dormann and Deutsch, 2002; Ferreira et al., 2002; Mallett and de Pillis, 2006) and microbial biofilm growth and function (Picioreanu et al., 1998a; Picioreanu et al., 1998b; Picioreanu et al., 2004). However, to our knowledge, this is the first time it has been applied to intracellular organization.

In the model, mitochondrial biogenesis required a nearest neighbor, simulating fission, and biogenesis signals were considered to be spatially uniform across the fiber. However, because of the presence of more mitochondria near the periphery of the fiber after several CA iterations, there was more biogenesis near the periphery than the fiber core. The dependence of mitochondrial biogenesis on gene expression in both the nuclear and mitochondrial genomes, as well as the subsequent transport of nuclear encoded proteins to the growing mitochondria (Wallace 1999; Chan, 2006; Hood et al., 2006; Gottlieb and Carreira, 2010; Terman et al., 2010), make it unlikely that biogenesis can lead to the observed spatial heterogeneity. In addition, mitochondrial biogenesis is mainly dictated by the energy status of the cell, i.e. the AMP/ATP ratio (Jager et al., 2007). As the ratio of AMP/ATP drops, protein kinases responsible for mitochondrial biogenesis are activated (Hood, 2001; Zong et al., 2002; Hardie, 2003; Hardie et al., 2006; Hood et al., 2006; Lopez-Lluch et al., 2006). AMP was not explicitly evaluated in the model, but it is in equilibrium with ADP, which was evaluated, through the adenylate kinase reaction. ADP, and hence AMP, concentration was uniform across all of the fibers (Fig. 1G,H), as seen previously (Jimenez et al., 2008; Locke and Kinsey, 2008; Kinsey et al., 2011), and so the model results are consistent with the idea that spatial variation in biogenesis signals will have little effect on the mitochondrial distribution.

Reaction-diffusion models used in our previous studies also suggest very slight spatial gradients in ADP and hence AMP (Dasika et al., 2011; Pathi et al., 2011). However, it is certainly possible that there are spatial gradients in biogenesis signals, but to our knowledge, none have been identified. Also, Ca²⁺ stimulates biogenesis (also possibly mediated through AMPK) and spatial gradients in Ca²⁺ may play a role in heterogeneous biogenesis (Hood et al., 2006). However, we didn’t find any evidence that suggested spatial gradients in Ca²⁺ that span the radius of the fiber. Recent studies suggest that exercise-induced ROS enhance mitochondrial biogenesis signaling markers (Kang et al., 2009; Powers et al., 2011a; Powers et al., 2011b). However, the exercise-induced ROS that stimulate mitochondrial biogenesis are primarily non-mitochondrial and are produced in other cellular regions (Powers and Jackson, 2008). Similar to Ca²⁺, we know of no studies that suggest substantial spatial gradients in ROS across the fiber that are similar to those seen in mitochondrial distribution. Although there is no evidence of spatial gradients in Ca²⁺ and ROS, it is certainly a possibility that there is a regional variation in these species that leads to targeted biogenesis. Another possibility is that higher oxygen concentration at the fiber periphery leads to higher mitochondrial biogenesis. Mitochondria are the site of heme biosynthesis (Mense and Zhang, 2006) and heme proteins function as oxygen sensors, perhaps providing a mechanism for oxygen-dependent mitochondrial biogenesis.

In contrast to biogenesis, the CA rules dictated that mitochondrial death depends on prevailing local oxygen concentration, and mitochondria deficient in oxygen (i.e. subjected to hypoxia) experienced greater probability of death. Such a fate of mitochondria was observed in fiber cores of very large fibers or in fibers working at a very high ATPase rate. The mitochondrial death rules were consistent with previous studies concerning the selective (targeted) death of mitochondria based on ROS production levels (He and Klionsky, 2009; Kim et al., 2007; Yen and Klionsky, 2008; Zhang et al., 2008) and hypoxia (Band et al., 2009). Also, in contrast to low spatial variation in ADP, substantial spatial variation is observed in the oxygen concentration (Fig. 1C), as seen previously (Jimenez et al., 2008; Kinsey et al., 2011). Thus, significant spatial variation in oxygen concentration is consistent with the view that there is a non-uniform signal for mitochondrial death (i.e. targeted mitophagy) in muscle fibers.

Model results depicting steady-state mitochondrial spatial distributions post CA iterations are consistent with experimental observations in dolphinfish red and white fibers (Fig. 2), juvenile and adult black sea bass white fibers (Fig. 3) (Nyack et al., 2007), as well as four putative fiber types in mouse (Figs 4, 5). These distributions of mitochondria in the wide range of fiber types considered in the present study are crucial to fiber function because the shifts result in increases in concentrations of vital metabolites such as oxygen and PCr compared with a random distribution. In all our model simulations, the ATP turnover rates considered were such that changes in high-energy phosphate molecules were within the physiological range. In all cases, the ATP depletion during simulated elevations in ATP demand was minimal, and the PCr concentration during contraction ranged from approximately 25 to 60% of the resting value. These results are consistent with changes in high-energy phosphate molecules that have been measured in muscle using NMR (Kushmerick and Meyer, 1985; Kushmerick et al., 1992).

Increases in η after mitochondrial redistribution were also observed, which indicate that the non-uniform distribution allows fibers to operate at higher rates of steady-state aerobic metabolism (higher R_ATPase) without being diffusion limited compared with the initial random distribution of mitochondria. It is interesting to observe that the maximum increase in η is seen in mammalian muscle fibers (Fig. 6D) operating at very high values of R_ATPase. This is because high rates of ATP demand lead to greater depletion of oxygen in the fiber core, resulting in more redistribution towards the fiber periphery (i.e. oxygen source). Finally, muscle fibers with a non-uniform mitochondrial distribution following CA iterations showed a significantly more negative ΔG, indicating a higher energy state.

The coupled CA/reaction-diffusion model used in the present analysis could adequately predict the mitochondrial distribution for a variety of fiber types as well as during hypertrophic growth in black sea bass (Figs 2, 3). However, a θ_det value of 0.35 was used for fish white fibers and 0.55 was used for fish red fibers and mouse muscle fibers, as a single value could not adequately describe the experimentally observed distributions in both aerobic and anaerobic fiber types. Initially, a single θ_det value of 0.35 was used for both aerobic and anaerobic fibers (results not shown here). The use of this value in aerobic fibers led to a greater migration of mitochondria towards the fiber periphery, which did not match experimental observations. In contrast, a value of 0.35 as θ_det for anaerobic fibers yielded distributions consistent with experimental findings. Also, a value of 0.55 was used for both aerobic and anaerobic fibers, and spatial distributions of mitochondria in aerobic fibers matched well with experimental observations, but distribution patterns observed in anaerobic fibers were more homogeneous and did not align with experimental observations. It is possible that two different values of θ_det may be needed because the sensitivity of mitochondrial death to low oxygen is different in aerobic and anaerobic fibers. It is also possible that white and red fibers are exposed to different boundary oxygen concentrations, which is not accounted for in the present analysis. Alternatively, the simple probability function we used may not fully capture the oxygen dependence of mitochondrial degradation. Although the need to use different θ_det values likely reflects some shortcomings in the current version of the model, it is not at odds with the basic principles examined in this study.

The mathematical model in the present study treated all mitochondria the same and did not differentiate between SS and IM mitochondria. Experimental findings suggest that there are differences between these two subpopulations (Kuznetsov and Margreiter, 2009; Wikstrom et al., 2009; Lombardi et al., 2000). However, there is not enough information to incorporate these in the birth and death rules of mitochondria used in the mathematical model in the present analysis.

The present study coupled experimental analyses of mitochondrial distribution with a predictive CA/reaction-diffusion model to evaluate a potential mechanism of controlling mitochondrial distribution as well as the functional advantage of a non-uniform arrangement of mitochondria in skeletal muscle. We found that in dolphinfish red and white muscle, black sea bass white muscle fibers during hypertrophic growth, and in four mammalian muscle fiber types, spatially uniform signals for mitochondrial biogenesis and oxygen-dependent mitochondrial death can explain observed mitochondrial distributional patterns. The non-uniform distributions alleviated diffusion constraints in the fibers, leading to an increased energy status and higher permissible rates of aerobic ATP turnover, which are findings consistent with our previous work (Pathi et al., 2011). The model results also compare favorably with independent experimental observations for the different fiber types considered in this study. Thus, the present study provides a plausible mechanism by which muscle fibers attain non-uniform mitochondrial distributions, and illustrates the importance of these distributions on muscle fiber structure and function. Experimental evidence of regionally variable ROS production and mitochondrial death are needed to further evaluate the proposed mechanism.

We thank Mark Gay from the University of North Carolina Wilmington for assistance with image processing.