## Abstract

The study of cardiac energetics commonly involves the use of isolated muscle preparations (papillary muscles or trabeculae carneae). Their contractile performance has been observed to vary inversely with thickness. This inverse dependence has been attributed, almost without exception, to inadequate diffusion of oxygen into the centers of muscles of large diameter. It is thus commonly hypothesized that the radius-dependent diminution of performance reflects the development of an anoxic core. We tested this hypothesis theoretically by solving a modification of the diffusion equation, in which the rate of oxygen consumption is a sigmoidal function of the partial pressure of oxygen. The model demonstrates that sufficiently thick muscles, operating at sufficiently high rates of oxygen demand or sufficiently low ambient partial pressures of oxygen, will indeed show diminished energetic performance, whether indirectly indexed as stress (force per cross-sectional area) development or as the rate of heat production. However, such simulated behavior requires the adoption of extreme parameter values, often differing by an order of magnitude from their experimental equivalents. We thus conclude that the radius-dependent diminution of muscle performance in vitro cannot be attributed entirely to an insufficient supply of oxygen via diffusion.

- heat production
- stress development
- oxygen diffusion limitation
- critical radius
- anoxic core

all animal life depends on oxygen. Hence, when experimentalists study an animal tissue preparation in vitro they strive to supply oxygen at a rate commensurate with its rate of metabolic consumption. When the in vitro preparation is a whole organ, such as the isolated heart, oxygen is usually supplied to the tissues via the coronary circulation, whether directly, in gaseous form (58), in simple saline solution (57), or reversibly bound in whole blood (21, 22, 101), erythrocytes (80), or perfluorocarbon-supplemented media (11, 87). When the in vitro cardiac preparation is isolated tissue, be it an excised “strip” (69), a papillary muscle (43), a trabecula (82), or a single cell (41), it is commonly merely superfused. Under this circumstance, the adequacy of oxygenation depends critically on the prevailing Po_{2}, the metabolic demand for oxygen, and the radius of the preparation.

There is abundant evidence that the performance of isolated cardiac preparations varies inversely with size. Indeed, an early investigation (15) showed that even the basal rate of oxygen consumption of excised feline papillary muscles depends inversely on muscle diameter. Even much smaller guinea pig right ventricular trabeculae are not immune when the oxygen demand is increased by the induction of K^{+} contractures or the addition of agents chosen to uncouple oxidative phosphorylation (16). Comparable examples exist in “mechanics” literature. Thus, Schouten and ter Keurs (86) documented an inverse relationship between peak isometric twitch force production and the diameter of rat ventricular trabeculae and small papillary muscles. More recently, Raman et al. (81) revealed a striking radius-dependent diminution of isometric stress (force per cross-sectional area) development by rat trabeculae in the face of high metabolic demand (8-Hz stimulation at 37°C). It is these latter results, in particular, that have motivated us to reexamine the issue of insufficiency of diffusive supply of oxygen to an isolated trabecula in vitro. Our primary investigative tool is mathematical modeling.

Mathematical modeling of the adequacy of tissue oxygenation requires solution of the diffusion equation. This approach has a lengthy history, commencing with the still much-referenced work of A. V. Hill (42), in which analytic solutions were developed for the steady-state Po_{2} within muscles of various cross-sections consuming oxygen at various Po_{2}-independent rates. Numerous variations on this theme have been published over the years. Some have restricted consideration to simple diffusion (63, 68, 92, 93), whereas others have included myoglobin-facilitated diffusion (14, 26, 33, 55, 59, 71, 73, 74, 83, 99, 100).

To explore the effects of time-varying metabolic demands, Barclay (3) solved the full partial differential diffusion equation allowing both time- and radial location-dependent solutions for papillary muscles and trabeculae carneae. Subsequently, Cairns et al. (10) used a comparable approach to probe the role of oxygen insufficiency in the course of a fatiguing stimulation protocol imposed on isolated fast- and slow-twitch skeletal muscle preparations. More recently, Goo et al. (34) developed a novel method of solution of the diffusion equation for cardiac trabeculae of any arbitrary cross-section in either the presence or absence of capillary (arteriolar) sources or (venous) sinks of oxygen. However, none of these recent reports has examined the radius-dependent decline of cardiac performance (whether indexed directly as the rate of oxygen consumption or indirectly as stress production or as the rate of heat production).

In the present study, we extend the use of Hill's diffusion equation to simulate the radius-dependent diminution of cardiac muscle performance in vitro using well-characterized parameters values. The cross-section of the muscle preparation is assumed to be circular. We further assume that active stress production, active rate of heat production, and basal rate of oxygen uptake are linearly proportional to the cross-sectional area that remains oxygenated. We demonstrate that preparations of sufficiently large diameter must, indeed, develop anoxic cores; however, there remains an irreconcilable quantitative disparity between experimental and simulated results.

## METHODS

### Models

#### The one-dimensional model.

A theoretical analysis of the adequacy of oxygen supply was made by estimating the distribution of steady-state Po_{2} (0 ≤ p ≤ P; in Pa) as a function of muscle radius (0 ≤ *r* ≤ *R*; in μm) throughout the circular cross-section of the muscle (3, 4, 10, 34, 59, 97). Our analysis followed the general approach described by Hill (42) but included a modification introduced by Loiselle (63) to account for a more realistic relationship between Po_{2} and the metabolic rate of oxygen consumption (*m*; in mol·m^{−3}·s^{−1}) and that incorporated myoglobin-facilitated oxygen diffusion (3, 14, 26, 33, 55, 59, 71, 73, 74, 83, 99, 100). The governing equation is as follows:
*D* are the solubility (in mol·m^{−3}·Pa^{−1}) and diffusion constant (in m^{2}/s) of oxygen in muscle tissue, respectively; *t* is time; C_{M} and *D*_{M} are the concentration (in mol/m^{3}) and diffusion constant (m^{2}/s) of myoglobin in muscle tissue, respectively; and *S* is the fraction of myoglobin that is saturated with oxygen.

At steady state, ∂p(*r*,*t*)/∂*t* = 0, so *Eq. 1* reduces to the following:
*K* is Krogh's constant (equivalent to the product of *D* and σ, in mol·m^{−1}·Pa^{−1}·s^{−1}).

*S*, as a function of local Po_{2}, is given by the following:
_{M50} (in Pa) is the Po_{2} that achieves half-saturation of myoglobin by oxygen and *n*_{M} is the Hill coefficient for myoglobin saturation.

*m* (in mol·m^{−3}·s^{−1}) is also a function of local Po_{2}, given by the following:
*m*_{b} and *m*_{a} are muscle basal and active rates of oxygen consumption, respectively; Po_{50} (in Pa) is the Po_{2} that yields the half-maximal value of *m*; and *n*_{o} is the Hill coefficient constraining the sigmoidal dependence of *m* on Po_{2}.

Following Cairns et al. (10), the fraction of the muscle cross-section that remains oxygenated (*A*_{oxy}) is calculated as follows:
*Eq. 4*) accounts for the sigmoidal diminution of *m* as Po_{2} passes through Po_{50} and approaches zero.

The one-dimensional (1-D) model was implemented in Matlab R2008a (The MathWorks) using the built-in “pdepe” function. The model was solved as a time-dependent problem (*Eq. 1*) and was run until steady-state values of p were reached (*Eq. 2*). Initially, the entire cross-section of the muscle was assigned the same value of Po_{2} as prevails at its surface, i.e., p(*r*,0) = P. At all times, Po_{2} at the muscle surface was constant, i.e., p(*R*,*t*) = P, and a flux symmetry condition was applied at the muscle center: ∂p(0,*t*)/∂*t* = 0.

#### Validation of the 1-D model.

Results were validated against the analytic solutions of the steady-state Hill diffusion equation (42), in which *m* is constant, independent of *r*, and myoglobin-facilitated oxygen diffusion is ignored. The critical radius [*R*_{crit(Hill)}], defined as the thickest cylinder which oxygen will fully penetrate (42), is as follows:
*r*′/*R*)^{2} and *r*′ is the radius at which p = 0 (42). Note the difference between θ and *A*_{oxy} (*Eq. 5*), where *A*_{oxy} takes into account the sigmoidal dependence of *m* on local Po_{2}.

#### Model parameters.

The values (expressed in both SI and conventional units) of parameters adopted for the model, together with their sources, are shown in Table 1. Literature values are further documented in Tables A1⇓–A3 (in appendix a), where those adopted for the model have been highlighted.

We adopted the active heat rate values, measured at 37°C, reported by Loiselle et al. (66) despite the fact that those data arise from only a single superfused rat trabecula. We converted the heat rate values from milliWatts per gram dry weight to milliWatts per meters cubed by taking into account the radius, length, and dry weight of the muscle and by converting oxygen consumption to the rate of heat production using an energetic equivalent of oxygen of 20 kJ/l (91). At 2 and 6 Hz, the equivalent rates of oxygen consumption were 0.044 and 0.150 mol·m^{−3}·s^{−1}, respectively. These values are consistent with those reported by Hütter et al. (46) for the rat heart over the same range of frequencies: 0.053 and 0.159 mol·m^{−3}·s^{−1}, respectively. Similarly, values at 3 and 4 Hz (0.069 and 0.094 mol·m^{−3}·s^{−1}, respectively) were in good agreement with the average values at 3 and 4 Hz obtained by Schenkman (85) in guinea pig hearts (0.092 mol·m^{−3}·s^{−1}). Finally, the value measured by Loiselle et al. (66) at 1 Hz was 0.022 mol·m^{−3}·s^{−1}, which again compared well with that of 0.024 mol·m^{−3}·s^{−1} obtained by Holmes et al. (44) in rabbit papillary muscles at optimal muscle length. With such excellent correspondence, we felt justified in adopting the values reported by Loiselle et al. (66) for our simulation of experimental data at 37°C.

The output of the model (Po_{2} as a function of *r*) is clearly a function of its parameter values. To explore this functional dependence, we performed a full parameter sensitivity analysis (13). The results of this analysis appear in Table B1 (appendix b), where it can be seen that the influence of each of the myoglobin-related parameters is negligible, being completely dominated by three parameters: namely, P, *m*, and *K* (*D*σ).

#### The two-dimensional model.

The 1-D diffusion equation is suitable for modeling preparations under a high rate of flow of superfusate through a muscle chamber since it allows the assumption that Po_{2} at the muscle surface is constant along the length of the muscle. However, it is necessary to invoke a two-dimensional (2-D) analog to describe the advective-diffusive supply of oxygen to an isolated cardiac trabecula in our low-flow micromechanocalorimeter, since Po_{2} is expected to decline progressively along the length of the muscle. For this simulation, we adopted the same parameter values as in the 1-D case. The 2-D case requires a description of the geometry of the measurement chamber. It was modeled as a cylinder of 500-μm radius whose length equals that of the trabecula. The upstream and downstream hooks, which attached the preparation, were modelled as cylinders of radius 200 μm and length 150 μm.

The 2-D model was implemented in COMSOL Multiphysics 3.5a (COMSOL) using the “convection and diffusion” module coupled with the “incompressible Navier-Stokes fluid dynamics” module. The model was solved in cylindrical coordinates as a time-dependent problem until steady state was achieved.

In cylindrical coordinates, the governing equation of advection and diffusion of oxygen is as follows:
^{3}), **u** is the velocity field (in m/s), ∇ is the standard del operator, and *m* and *D* have the meanings defined above. Oxygen concentration is the product of its partial pressure and its solubility, i.e., c(*r*,*l*,*t*) = σ × p(*r*,*l*,*t*), where *l* is location along the muscle or chamber length (0 ≤ *l* ≤ *L*; in m).

In the muscle domain, **u** = 0, and *m* is as follows:
_{o50} is the oxygen concentration yielding the half-maximal value of *m*.

In the superfusate domain, *D* was taken to be 2.07 × 10^{−9} m^{2}/s, which was calculated from van Stroe and Janssen (95) at 22°C and 130 mM NaCl. This value agreed well with 2.1 × 10^{−9} m^{2}/s as measured by Evans et al. (23) and is within the range calculated by Ju and Ho (47) (1.88 × 10^{−9} m^{2}/s) and by Hung and Dinius (45) (2.86 × 10^{−9} m^{2}/s). The solubility of the superfasate at 22°C was calculated to be 1.32 × 10^{−5} mol·m^{−3}·Pa^{−1} using the value for physiological saline (22.7 ml·l^{−1}·atm^{−1} at 37°C and Q_{10} of 0.83) as reported by Loiselle (61). The calculated value at 22°C was in excellent accord with a value of 3.93 × 10^{−2} ml·l^{−3}·Torr^{−1} (or, equivalently, 1.31 × 10^{−5} mol·m^{−3}·Pa^{−1}) as measured by Graham (35) and 1.32 × 10^{−5}·mol·m^{−3}·Pa^{−1} as measured by Evans et al. (23) at the same temperature.

The velocity field **u** of the superfusate was calculated as the following ratio of the superfusate flow rate to the annular area of the superfusate domain: **u** = 1/π(*R*_{c}^{2} − *R*^{2}), where *R*_{c} is the radius of the superfusate chamber.

The Navier-Stokes model of incompressible Newtonian flow (constant viscosity) was used to calculate the velocity field *u* that results from advection in the superfusate domain as follows:
_{s} is pressure (0 Pa), and η is the dynamic viscosity of the superfusate, approximated by that of water

(1 × 10^{−3} Pa·s).

In the advection-diffusion model, the following boundary conditions were assumed: the advective flux, *n* × (−*D*∇c) = 0 for the muscle center and the side walls, there is Po_{2} continuity at the muscle-superfusate interface, and a constant c = P × σ at the superfusate inlet. In the incompressible Navier-Stokes model, the boundary condition at the superfusate inlet was assumed to obey parabolic Poiseuille flow. At the superfusate outlet, the outflow condition was defined as pressure with no viscous stress, using a value of 0 Pa. At the chamber wall and at the muscle-superfusate interface, a nonslip (**u** = 0) boundary condition was adopted. Triangular elements were used, and refinement of the mesh produced some 5,000 elements. The model was solved using the UMFPACK direct solver.

The Hagen-Poiseuille equation, which can be derived from the Navier-Stokes equations, governs laminar and incompressible flow through a circular cross-section. The axial velocity (*v*_{z}; in m/s) distrubution along the muscle or chamber length (0 ≤ *l* ≤ *L*; in m) for Poiseuille flow with a nonslip boundary condition at the muscle wall, *v*_{z}(*R*) = 0, and at the chamber wall, *v*_{z}(*R*_{c}) = 0, is given by the following:
*L*) was assumed to be a constant whose value is obtained by equating Q = 1 μl/s (or 10^{−9} m^{3}/s). For example, for our largest trabecula preparation (*R* = 150 μm), dP/d*L* was calculated to be 68.8 Pa/m within the measurement chamber (*R*_{c} = 500 μm).

The oxygenated fraction of muscle volume (V_{oxy}) was computed as follows:

#### Validation of the 2-D model.

Results from the 2-D model were validated against those from the 1-D model *1*) by setting c = P × σ constant at all points along the length of the muscle surface and *2*) by increasing the superfusate flow rate 100-fold. Computed values of c(*r*, *l*, *t*) were converted to p(*r*, *l*, *t*). The former criterion is essentially the same as solving the 1-D model in which *c* (or P) is constant along the muscle length. The latter criterion, likewise, sets constant *c* (or P) along the muscle length as a result of an extremely high superfusate flow rate, implying a high rate of oxygen supply. The results arising from the two models were indistinguishable (results not shown).

We also compared differences between experimentally measured values of c at two different stimulus frequencies with those predicted by the 2-D model. Experimentally, we simultaneously measured the active rate of heat production, active stress production, and downstream oxygen concentration of the superfusate for a trabecula (radius: 125 μm and length: 2.8 mm) contracting at 0.2, 2, and 4 Hz. The superfusate oxygen concentration was measured at the downstream end of the trabecula using a commercially available oxygen probe (Foxy-18G, Ocean Optics). Using the model, we computed the average oxygen concentration in the annular region of superfusate flow at the downstream end of the muscle as follows:
^{3}) agreed well with that predicted by the 2-D model using the measured active rate of heat production (8.7 mmol/m^{3}). Likewise, between 0.2 and 4 Hz, the difference between the measured values of c (14.9 mmol/m^{3}) was in excellent agreement with that predicted by the 2-D model (14.6 mmol/m^{3}), thereby giving us confidence in the numerical output of the model.

### Experimental Measurements

We simultaneously measured the steady-state active rate of heat production and active stress production of geometrically uniform right ventricular trabeculae from Wistar rats, as previously described (39, 40), at two stimulus frequencies: 0.2 and 4 Hz. The composition of the superfusate was (in mM) 130 NaCl, 6 KCl, 1 MgCl_{2}, 0.5 NaH_{2}PO_{4}, 2 CaCl_{2}, 10 HEPES, and 10 glucose. pH was adjusted to 7.4 using Tris. The superfusate was vigorously bubbled with 100% O_{2} at room temperature (22°C) and was maintained at a flow rate of 1 μl/s. In total, 15 trabeculae of radii between 60 and 150 μm were examined. Subsequently, the extracellular Ca^{2+} concentration ([Ca^{2+}]_{o}) was halved, i.e., to 1 mM, and the thermal and mechanical performance of 11 of the trabeculae were reexamined. All experimental protocols were approved by the Animal Ethics Committee of The University of Auckland.

To extend our observations made at 22°C, we explored the radius dependence of stress production in rat trabeculae under more physiological conditions (37°C and 1.5 mM [Ca^{2+}]_{o}) and at the substantially higher superfusate flow rate of 6–8 μl/s. We examined the performance of 16 trabeculae at stimulus frequencies of 5 and 10 Hz; 11 preparations were also examined at 1 Hz.

## RESULTS

### Experimental Results

Figure 1 shows, as a function of *R*, steady-state active twitch stress (force per cross-sectional area) production (*A*) and active rate of heat production (*B*) at 22°C and 2 mM [Ca^{2+}]_{o} of 15 trabeculae in response to both 0.2 Hz (open circles) and 4 Hz (solid circles). As is evident, in neither case was there radius dependence at 0.2 Hz, but diminution of both stress and heat rate with increasing muscle radius was evident at 4 Hz.

Given the negative relationship between active heat rate and *R* at 4 Hz (Fig. 1*B*), which is qualitatively consistent with the conjecture of insufficient diffusive supply of oxygen in large muscles, we extrapolated the linear regression line to zero radius. We assume that this extrapolated value (79.4 kW/m^{3}), which we define as the “extreme heat rate,” represents the heat rate achievable by a single myocyte unfettered by diffusion insufficiency. The extreme heat rate was converted to the “extreme rate of O_{2} consumption” using an energetic equivalent of oxygen of 20 kJ/l and muscle density of 1.06 × 10^{3} kg/m.

Figure 2*A* shows that the decline of stress with *R* at room temperature and 4-Hz stimulus frequency observed in Fig. 1*A* was not a consequence of elevated Ca^{2+} concentration; reduction of [Ca^{2+}]_{o} to 1 mM produced comparable behaviour. Similarly, at 37°C and physiological Ca^{2+} concentration appropriate for the rat (1.5 mM), stress declined with *R* at all stimulus frequencies examined (Fig. 2, *B–D*). We note, once again, that each behavior is qualitatively consistent with a radius-dependent limitation of oxygen diffusion.

### Modeling Results

#### Negligible contribution of myoglobin.

We solved the time-dependent 1-D diffusion equation (*Eq. 1*) until steady state (*Eq. 2*) with and without the myoglobin term and plotted Po_{2} as a function of *r* (note the difference between *r* and *R*: 0 ≤ *r* ≤ *R*). As shown by the comparison in Fig. 3, *A* and *B*, the contribution of myoglobin was found to be negligible, consistent with the experimental observation by Merx et al. (70) that ambient Po_{2} must be extremely low (∼0.01 atm in isolated cardiomyocytes) before such a contribution can be detected. It is also consistent with the theoretical results reported by Loiselle (59) and confirmed by Barclay (3) as well as with the result of our parameter sensitivity analysis (Table B1). Because of its negligible contribution to the diffusive flux of oxygen, the myoglobin term was ignored in subsequent modeling results.

#### The relationship between p and r.

In Fig. 3, we show the calculated Po_{2} profiles for the thinnest (60 μm) and thickest (150 μm) muscles examined (at 22°C) in the presence (*A*) and absence (*B*) of a contribution from myoglobin. For these simulations, we used both the experimentally measured active heat rates (red) as well as the (extrapolated) extreme heat rate (black). In addition, using the extreme heat rate, we explored the profiles of Po_{2} by doubling the literature value of *m* (or, equivalently, halving the value of *K* in the absence of a contribution from myoblobin; see *Eq. 2*; blue), or halving the value of P (green). Under all scenarios, the thinnest muscle was predicted to have a fully oxygenated cross-section. The cross-sectional area of the thickest muscle, in contrast, was predicted to be fully oxygenated except for those cases where *m* was doubled or P was halved.

#### Validation.

We validated our model by comparing its output with analytic solutions arising from the Hill formulation. For the thickest muscle, when *m* was doubled, the radius of the anoxic core (the point at which the Po_{2} profile first intersects the abscissa) was predicted to be 51.0 μm and *A*_{oxy} (*Eq. 5*) was estimated to be 0.877. Similar results were obtained when P was halved, i.e, the radius of the anoxic core was calculated to be 49.8 μm and *A*_{oxy} was 0.875. In both cases, Hill's analytic calculations yielded values of 52.2 μm for *r*′ and 0.879 for θ (*Eq. 7*). The differences between our values and those arising from solution of the Hill equation reflect the Po_{2} independence of *m* assumed by Hill (42) vis-à-vis the Po_{2} dependence that we assume (*Eq. 4*). The differences were very small because Po_{50} is only 0.2% of P (i.e., the ratio of 213 to 1.01 × 10^{5} Pa; see Table 1).

#### Simulated A_{oxy} as a function of R.

As shown in Fig. 4 (solid lines), the calculated *A*_{oxy} is independent of *R* up to a point where *A*_{oxy} decreases monotonically. We define this point as the “critical radius” [*R*_{crit}; note that *R*_{crit} differs from R_{crit}_{(Hill)} (*Eq. 6*), in which *R*_{crit} takes into account the Po_{2} dependence of *m*] and divided the relationship between *A*_{oxy} and *R* into two portions. The first portion started from zero radius and extended to *R*_{crit}. Muscles whose radii fell in this segment were predicted to be fully oxygenated and, hence, were predicted to develop maximum stress. The second portion of the relationship extended beyond *R*_{crit}. Muscles in this region were predicted to develop anoxic cores and to develop submaximal stresses, in strict proportion with *A*_{oxy}.

Using the parameter values for 22°C and the extreme heat rate value, converted to its equivalent value of *m*, *R*_{crit} was computed to be 165.9 μm (Fig. 4*A*), in agreement with that calculated using Hill's analytic solution (167.5 μm; *Eq. 6*). As shown in Fig. 4*A*, *insets*, a muscle of radius 100 μm was predicted to be fully oxygenated and, thus, capable of developing maximal active stress. Muscles of radii 200, 400, or 1,000 μm were predicted to have *A*_{oxy} of 0.912, 0.529 and 0.226, respectively. We assumed that such muscles would develop correspondingly submaximal active stresses. When the extreme heat rate value (*m*) was doubled (or the literature value of *K* was halved), the simulated value of *R*_{crit} was reduced to 117.9 μm (cf. Hill's analytic solution of 118.4 μm). Under this scenario, we would expect some of the muscles examined in this study to be hypoxic. The thickest muscle (radius 150 μm), for example, would be expected to developed an anoxic core such that its *A*_{oxy} would be 0.877 (Fig. 4*B*), thereby producing only 0.877 of its potential maximum active stress. When *m* was arbitrarily increased 10-fold (or *K* was reduced 10-fold), *A*_{oxy} of the thickest muscle was 0.454 and would be expected to produce 0.454 of its maximum active stress (Fig. 4*B*).

#### Simulation of experimental data.

Following the same approach used in Fig. 4, we simulated our experimentally measured steady-state active stress-radius data at 4 Hz (22°C and 2 mM [Ca^{2+}]_{o}) and superimposed the predicted *A*_{oxy}. As shown in Fig. 5*A*, using the literature value for *m*/*K*, *R*_{crit }was computed to be 165.9 μm. All muscles examined (radii 60–150 μm) were smaller than this, such that they were expected to have been fully oxygenated. The predicted profile of *A*_{oxy} completely failed to mimic the experimental data, despite using the extreme heat rate (Fig. 1*B*). To mimic, quantitatively, the decline of experimentally measured active stress with muscle radius, it would be necessary to increase the value of *m*/*K* by 10-fold. Under this implausible 10-fold increase, *R*_{crit} would become 52.4 μm, such that all muscles examined would be predicted to be hypoxic, thereby developing stress quantitatively in line with the predicted *A*_{oxy} (under the assumption that active stress production is linearly proportional to the fractional cross-section oxygenated). As shown in Fig. 5*B*, a comparable 10-fold reduction of P was found to be required. It should be emphasised that a Po_{2} value of 10 kPa would be only one-half of that of a room air-equilibrated superfusate.

We applied the same simulation procedure to our data recorded at 37°C and 10 Hz. As shown in Fig. 6, all the trabeculae examined were predicted to be fully oxygenated. That is, their radii were smaller than the predicted value of *R*_{crit}. To mimic the observed decline of stress production, it was necessary to elevate the literature value of *m*/*K* by at least 10-fold (or to decrease the value of P by a comparable factor).

In addition to our own data, we also simulated those of others. We applied the same simulation procedure to the digitized experimental data of Raman et al. (81), as shown in Fig. 7. These data were of particular interest since they showed a “plateau” of active stress in trabeculae of small radii. They also showed that active stress plummeted in trabeculae of ∼60- to 90-μm radius, allowing the authors to infer that the *R*_{crit} is 75 μm in trabecula contracting at 8 Hz and 37°C. However, the value of *R*_{crit }predicted by our 1-D model was 157.2 μm using literature values of *m*/*K*. Once again, quantitative simulation of these data could be approximated only by the substitution of extreme parameter values (at least 10 times those reported in the literature) for either *m*/*K* or P.

Similar results were obtained for the digitized experimental data of Schouten and ter Keurs (86), as shown in Fig. 8. Using the literature value of *m*/*K*, *R*_{crit} was 191.7 μm, suggesting that those preparations of radii greater than *R*_{crit} (all of which were papillary muscles) were hypoxic. However, they developed stresses often much less than predicted. For those preparations of radii smaller than *R*_{crit}, the pronounced radius-dependent decline of stress production remains a mystery, as it cannot be attributed to inadequate diffusive oxygen supply. By increasing the value of *m*/*K* by 10- or 20-fold, *R*_{crit} was reduced to 60.6 and 44.0 μm, respectively. Only under this improbable scenario could acceptable quantitative agreement between predicted and experimentally measured stress be achieved.

To avoid the assumptions invoked when relating measurements of either heat or stress to oxygen consumption, we simulated the experimentally measured rates of basal oxygen uptake of quiescent papillary muscles of the cat, as reported by Cranefield and Greenspan (15). As shown in Fig. 9, using literature values for *m*/*K*, the predicted value of *R*_{crit} was 389.0 μm. *R*_{crit} calculated by the authors (320 μm) was roughly comparable. About one-half of the preparations were predicted to be hypoxic (i.e., to have had radii greater than *R*_{crit}), but their measured rates of oxygen uptakes were much less than predicted. It was necessary to multiply *m*/*K* by a factor of five or to decrease the value of P to 20 kPa (room air equilibrated) to mimic the observed decline of basal rate of oxygen uptake with radius.

To underscore our discomfort with the hypothesis that radius-dependent diminution of performance (i.e., active stress) is exclusively due to oxygen diffusion limitation, in Fig. 10 we show the digitized data of Fisher and Kavaler (24), which were measured from canine papillary muscles in situ. Instead of being isolated preparations superfused with an external medium, these preparations were blood perfused by their own coronary circulation and thus very unlikely to have developed hypoxic cores. Nevertheless, a similar radius-dependent decline of active stress was evident.

#### The 2-D model result.

As a final test of the lack of correspondence between simulation and experiment, we determined the steady-state results of 2-D modeling of the advection-diffusion equation (*Eq. 8*) applied to a specific case (22°C, 2 mM [Ca^{2+}]_{o}, and 4-Hz stimulation), chosen for its extreme conditions: the largest muscle that we examined (radius: 150 μm and length: 4.0 mm) subjected to the extrapolated extreme heat rate value (Fig. 1*B*) at a superfusate flow rate of 1 μl/s.

Figure 11 shows steady-state Po_{2} contours (blue: 101 kPa and red: 0 kPa) within the superfusate and muscle domains. The geometries of the superfusate chamber and hooks mimic those of the measurement chamber (90) of our micromechanocalorimeter (39). Figure 12*A* shows Po_{2} as a function of *r* (radial distance from the muscle center), whereas Fig. 12*B* shows Po_{2} as a function of *l* (axial distance from the upstream end). At the upstream end, the Po_{2} of the superfusate domain was set to 101 kPa (1.24 mol/m^{3}, boundary condition), whereas within the muscle, it dropped to 9.3 kPa at the center. At the downstream end, in contrast, the superfusate Po_{2} was lower (average of 97.9 kPa or, equivalently, 1.21 mol/m^{3}; *Eq. 14*) such that at the muscle surface and muscle center, Po_{2} was 62.3 and 0.04 kPa, respectively. Note that the latter value is in the vicinity of Po_{50} (0.213 kPa; Table 1). *A*_{oxy} monotonically decreased along the muscle length (Fig. 12*C*), from a value of 1.000 to 0.962, such that V_{oxy} (*Eq. 13*) was estimated to be 0.984. We predict that this muscle could have been producing heat at a rate 0.984 times the extreme heat rate. However, its measured rate of heat production was 45.8 kW/m^{3}, only ∼0.58 times the predicted value.

## DISCUSSION

In isolated cardiac preparations, an inverse relationship between muscle stress production (6, 19, 24, 25, 27, 30, 31, 50, 53, 60, 63, 65, 81, 86), basal rate of heat production (64, 98), or basal rate of oxygen uptake (15) and muscle radius has repeatedly been reported. Many experimentalists have used the observed inverse relationship between muscle contractile performance and muscle radius to demonstrate the inadequacy of oxygenation in preparations that exceed some particular size and ipso facto to declare a “critical radius” above which diffusion fails to supply sufficient oxygen to the core of the preparations to prevent the diminution of performance (6, 15, 27, 50, 81, 86). In the present report, we tested the conjecture that the radius-dependent decline of muscle twitch stress production, observed experimentally, reflects inadequacy of the diffusive oxygen supply. Our “test” consisted of comparing the results of mathematical models of diffusive oxygen transport with experimental data, based on the assumption that muscle stress production is linearly proportional to the cross-sectional area of the muscle that remains oxygenated. We further assume that isolated cardiac preparations (papillary muscles and trabeculae carneae) are circular in cross-section. The implications of these assumptions require examination.

### Assumptions of the Model: Implications

#### Cross-sectional geometry.

The elliptical geometry of trabeculae in cross-section is now well recognized (29, 54, 66, 81, 84, 96). Indeed, Goo et al. (34) recently applied an index of eccentricity (ε) to quantify the extent of departure from circularity in a variety of specimens from both ventricles of the rat, emphasizing the potential error of estimating stress development under the “circularity” assumption. However, no comparable risk applies to oxygen diffusion. The reason is straightforward. For any given cross-sectional area, the diffusion distance from the perimeter to core will necessarily be less for an ellipse than for a circle. This means that an even larger “radius” (semi-major axis) can be tolerated before an anoxic core would develop in a trabecula of elliptical cross-section. This issue can be quantified by comparing two trabeculae of equal cross-sectional areas: *A* = π*R*^{2} (circular) and *A* = π*ab* (elliptical; where *a* and *b* are the semi-major and semi-minor axes, respectively). ε is given by the following (34):
*Eq. 15* yields the following:
*b* < *R*, thereby favoring diffusion into an ellipse in proportion to the extent of its eccentricity, as given by *Eq. 16*.

#### Linear dependence of stress on A_{oxy.}

This is an admittedly gross assumption that denies any significant contribution by anaerobic metabolism to the production of ATP. It does, however, allow for active stress development to decline progressively in the annulus between the fully oxygenated region and the anoxic core (as shown in the shaded regions in Fig. 4). This annulus is of small extent under the parameters (Table 1) adopted to describe the dependence of the consumption of oxygen on its partial pressure. More troublesome is the implicit assumption that the ability of a muscle to develop active stress is unaffected by an anoxic region in which rigor bridges may develop, attended by an increase of stiffness. At the very least, passive force would be expected to increase, thereby probably diminishing the ability of the muscle to produce active stress. We make no attempt to accommodate this scenario.

### Simulated Versus Experimental Results: an Order of Magnitude Discrepancy

Experimentally, at room temperature, we observed that both the stress and rate of heat production of cardiac trabeculae remain reasonably independent of muscle radius at low (0.2 Hz) stimulus frequency (Fig. 1). At the relatively high stimulus frequency of 4 Hz, muscle performance declined with muscle radius, at both 2 mM (Fig. 1) and 1 mM (Fig. 2*A*) [Ca^{2+}]_{o}. Under physiological conditions (37°C, 1.5 mM Ca^{2+}) and at 1-, 5-, or 10-Hz stimulus frequency, a radius-dependent decline of contractile performance was also observed (Fig. 2, *B–D*). We simulated these experimental data using literature values for model parameters and predicted a critical radius that was greater than that of any of the trabeculae examined (Figs. 5 and 6). This suggests that all of our preparations were fully oxygenated and that their radius-dependent decline of stress production cannot be attributed to a limited diffusive supply of oxygen. Indeed, to mimic the observed decline quantitatively, it was necessary to increase the metabolic rate of oxygen consumption, to decrease Krogh's diffusion constant, or to decrease superfusate Po_{2} by at least an order of magnitude.

We created a 2-D model to examine the consequences of a comparatively very low rate of flow (on the order of 100-fold less than that commonly used in organ baths). Despite the consequent greatly reduced rate of provision of oxygen, the model predicted that even our thickest (diameter: 300 μm) and longest (4 mm) preparation, consuming oxygen at the “extreme rate,” would not have experienced an anoxic core at its downstream end. [Note that the predicted minimum value reached (120 Pa) was of the same order of magnitude as the P_{50} for oxygen consumption: 213 Pa (Table 1).] Nevertheless, the actual heat rate observed experimentally was only 0.58 of the predicted extreme heat rate, a discrepancy that cannot be attributed to insufficient oxygen supply under the assumptions of either our 1-D or 2-D mathematical models.

Schouten and ter Keurs (86) and Raman et al. (81) inferred *R*_{crit} values of 100 and 75 μm, respectively, based on their experimental findings that muscle stress production drops dramatically in preparations that exceed these observed radii. Our simulation results suggest that muscle cross-sectional area remains fully oxygenated up to radii of 160–220 μm, whereas experimentally measured stress production declines at considerably lower values (Figs. 7 and 8). We could simulate the *R*_{crit} values reported by Schouten and ter Keur (86) and by Raman et al. (81) by increasing the metabolic rate of oxygen consumption (or by decreasing either Krogh's diffusion constant or the Po_{2} of the superfusate) by at least an order of magnitude. However, in each case, the required parameter values were wholly improbable.

Using parameter values gleaned from the literature, our model predicted that some of the large papillary muscle preparations of Schouten and ter Keur (86) were hypoxic since their radii were greater than the predicted *R*_{crit} (Fig. 8). However, the stresses produced by those putatively hypoxic preparations were still substantially less than the predicted values. A comparable mismatch was obtained when we examined the data of Cranefield and Greenspan (15) for quiescent papillary muscles from the cat, where the preparations predicted to be hypoxic had rates of oxygen uptake substantially lower than our predicted values (Fig. 9). Clearly, these results imply that inadequate diffusive supply of oxygen cannot entirely explain the lower stress production or lower rate of basal oxygen consumption of hypoxic muscles.

For completeness, we also simulated (data not included) the negative radius-dependent contractile performance observed by several other investigators (6, 19, 27, 30, 31, 53, 67). In every case, the smallest preparations were at least 250–300 μm in radius, which precludes the estimation of suitable values for maximum active stress production by small, and presumably fully oxygenated, preparations. As is evident from our data, as well as from those of Raman et al. (81) and Schouten and ter Keurs (86), even small preparations (radii < 200 μm) experience decreased stress production with increasing radius.

Clearly, there is an unbridgeable disparity between measurement and simulation if the hypothesis that the inverse relationship between active stress development and muscle radius is to be upheld. But, that is not to deny that stress development varies inversely with radius. Indeed, Fisher and Kavaler (24) showed that this much-observed phenomenon prevails even in blood-perfused papillary muscles in situ. It thus appears that we must turn elsewhere to seek an explanation for the radius-dependent decline of contractile performance of small preparations and something more than oxygen diffusion limitation for large and hypoxic preparations. What factor(s) might contribute? At least six distinct possibilities, in three different categories, arise.

#### Structural differences.

*Possibility 1* is that trabeculae contain a variable amount of noncontractile (collagenous) extracellular tissue. However, the hyperbolic dependence of collagen content on radius (84) means that specimens of larger cross-sectional area have a greater proportion of their cross-sections occupied by contractile tissues and are thus capable of developing proportionately greater stress production. However, this is contrary to what is observed. *Possibility 2* is that perhaps there exists some radius dependence of intracellular components that could provide an explanation. To investigate that possibility, Delbridge and Loiselle (19) sought ultrastructural differences in rabbit papillary muscles of disparate cross-sections. Contrary to their expectations, they found that the relative proportion of cross-sectional area occupied by mitochondrial and contractile matrix elements was invariant between large and small specimens.

#### Contractile pattern.

Sarcomere inhomogeneity (*possibility 3*) has been proposed to play a role in reducing the active stress developed by thicker muscles. Several authors (19, 24, 53) have invoked this explanation, which would require that the degree of sarcomere inhomogeneity increases with increasing muscle diameter, resulting in a reduced capacity of larger muscles to operate at optimal filament overlap and, accordingly, less active stress development. In support of this contention, macroscopic structural inhomogeneity, in the form of nonuniform segmental length changes, has indeed been observed (24, 56, 79). Nevertheless, we remain sceptical since it is difficult to imagine how sarcomere inhomogeneity could produce a comparable radius-dependent decline of basal metabolism, as revealed by Cranefield and Greenspan (15) in quiescent papillary muscles.

#### Ions and metabolites.

Although reduced intracellular K^{+} concentration (*possibility 4*) has an inhibitory effect on stress production (8, 9), no correlation between muscle size and intracellular K^{+} concentration was found by Page and Solomon (75), even in preparations as large as 650 μm in radius. We therefore consider it unlikely that a radius-dependent depletion of transmembrane K^{+} gradient arises. *Possibility 5* involves glucose. In all of the experiments quoted, the exogenous metabolic substrate was glucose (5–10 mM). Perhaps the radius dependence of performance reflects the limited diffusion of glucose rather than oxygen. However, consideration of stoichiometries and concentration gradients argues to the contrary. The stoichiometry of glucose oxidation is 6 mol oxygen/mol glucose. But, the diffusion constant of glucose is about one-sixth that of oxygen (12, 38). Thus, the difference between the diffusion constants is offset by the difference in stoichiometry. Furthermore, at a bath concentration of 10 mM glucose, the maximal difference in concentration, from the surface of the muscle to its core, is 10 mol/m^{3}, whereas that of oxygen is 1.4 mol/m^{3} (Table 1). Hence, the relative diffusivity (i.e., the ratio of the products of the diffusion constants and their concentration gradients) favors glucose by a factor of seven. Finally, Han et al. (40) demonstrated that varying glucose concentrations between 5 and 30 mM had a negligible effect on either the active stress or active heat production of rat trabeculae. Furthermore, given our modeling results, which showed an adequacy of oxygen supply (Figs. 5⇑⇑–8), it seems unlikely that metabolites such as lactate or P_{i}, which become elevated under anaerobic conditions, could be responsible (*possibility 6*). Nevertheless, we have examined this possibility in more detail, following the modeling approaches adopted by Stuyvers et al. (89) and Stienen et al. (88). We accounted for the negative sigmoidal dependence of muscle stress development (*S*) on P_{i} concentration as follows:
*R*, when its rate of production is ρ, as follows:
_{Pi} is the diffusion constant for P_{i} [2.73 × 10^{−10} and 4.05 × 10^{−10} m^{2}/s, respectively, at 20 and 35°C (78)]. The results show that it is necessary to increase the simulated rates of P_{i} production by at least an order of magnitude to simulate the data shown in Fig. 13, A–D. These simulations show that, even if anoxia were to develop, the resulting accumulation of P_{i} would be insufficient to explain the observed radius dependence of stress production.

### Conclusions

From the various studies (both our own and those of others) and our mathematical simulations (using literature values for both 1-D and 2-D model parameters), we conclude that the radius-dependent decline of twitch stress production arises, at least in part, from cellular or tissue mechanisms that are independent of diffusive oxygen supply. Equivalently, we infer that the radius-dependent decline of muscle twitch heat production reflects a decline of oxygen demand and not insufficiency of oxygen supply.

## GRANTS

This work was supported by Marsden Fund (Royal Society of New Zealand) Project Grant 06-UOA-123.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## ACKNOWLEGDMENTS

The authors thank Dr. M. T. Cooling for performing the parameter sensitivity analysis.

## APPENDIX A PARAMETER VALUES FROM THE LITERATURE

Tables A1⇓–A3 show literature values for the various parameters.

## APPENDIX B PARAMETER SENSITIVITY ANALYSIS

In the Morris method, a parameter space is defined by setting upper and lower bounds on parameter values. Analysis trajectories through the parameter space are calculated to form a statistical sample of the effect of parameter changes on one or more objective functions. Here, 2 sets of 5,000 trajectories for each parameter were run independently on the 2 objective functions. This number of trajectories gave comparable results for two independent runs. The 2 sets of runs were then combined to achieve an overall population of 10,000 sensitivity measures. The average absolute effect of each parameter, within its defined range, on each objective function, was computed.

Following Cooling et al. (13), we determined the sensitivity of our model parameters (using the values appropriate for 22°C; Table 1) on two objective functions, namely, *R*_{crit} and *A*_{oxy} at *R* = 100 μm. Each of the parameters was varied within ±32% of its tabulated magnitude (Table 1) except for P, which was constrained to the range of 0.2 × 10^{5}–1 × 10^{5} Pa.

Table B1 shows the sensitivity of the parameters and ranks their effects on the first objective function (*R*_{crit}). The four most sensitive parameters are P, *m*, *D*, and σ. For example, to change *R*_{crit} by 10 μm, only a change of ∼8% of P, on average, is required. Consequently, a 10% change in P, on average, changes *R*_{crit} by ∼12.5 μm. In contrast, the parameters representing the myoglobin terms (*D*_{M}, C_{M}, P_{M50}, and *n*_{M}) have comparatively negligible effects. Similar sensitivity rankings were obtained for the second objective function (*A*_{oxy} at *R* = 100 μm).

- Copyright © 2011 the American Physiological Society

## REFERENCES

- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.↵
- 98.↵
- 99.↵
- 100.↵
- 101.↵