Quantification of mitochondrial distribution in single cell cross-sections

We collected serial block face scanning electron microscopy (SBF-SEM) images of one complete cell [17, 18]. Fig 1A presents a longitudinal section from the 3D dataset, while Fig 1B, 1C and 1D show the cross sections corresponding to the positions marked by the white lines on Fig 1A. The first cross-section in Fig 1B (section 1) was located near a branching point of the cell. On the other hand, section 2 is in the proximity of a nucleus, although none of the three cross-sections contain any portions of the nucleus. As evident from these images, mitochondria were highly clustered in both lateral and longitudinal sections. In order to quantify the distribution of mitochondria in the cross sections, we used IMOD [19], an open source image processing software package, to segment the cell membrane and mitochondrial cluster boundaries [20] in the three cross sectional images (Fig 2A). Following this, we calculated the local area density of mitochondria (simply referred as mitochondrial density here on) at different points in the images, which is defined as the normalized mitochondrial area (Area_M) present inside a square local window (Area_T) centred at a given point (Fig 2B). We tested the sensitivity of this measurement to the size of the square window used in the analysis. The difference between the first and third quantile of the mitochondrial density distribution was maximum when a small local window size (<1 μm) was used (Fig 2C) and decreasing as the window size increased. Increasing the window size also leads to a decrease in the difference between median of the mitochondrial density distribution and the fractional area of the total mitochondria present in the cross section. However, increasing the window size beyond 1.6 μm did not significantly change the median of the distribution, and it also resulted in many outliers (Fig 2C). Therefore, a square window of 1.6 μm x 1.6 μm was used for our analyses (Fig 2D–2F).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Analysis of mitochondrial distribution in different cross-sections.

(A) Cross-section 1 with mitochondria segmented in green and sarcolemma in red. (B) Square windows representing the ratio of mitochondrial area vs myofibrilar area present inside these windows (C) Statistical box plots of the distribution of these ratios corresponding to different sizes of the square windows. The star in the middle of box plots represents the area fraction of total mitochondria present in section 1. (D-F): Box plot representing the distribution of mitochondrial density at various points in the cross-section 1–3. The stars represent the total mitochondrial area fraction present in each cross section. Mitochondrial area fraction is 55% in section 2 and 3 and 45% in section 1.

https://doi.org/10.1371/journal.pcbi.1006640.g002

We found a significant variation in the mitochondrial density distribution of each cross section (Fig 2D–2F). The first and third quantile of the distribution differs by a margin (IQR) of 20% - 25% in these images. The margin of this difference is lowest in section 2. The star marks inside the boxplots represent the fractional area of total mitochondria present in these sections. Section 2 and 3 have similar area fraction of mitochondria (55%), followed by section 1 (45%).

A high resolution two-dimensional spatial model of cardiac bioenergetics and its validation with experimental results

The descriptions of mitochondrial distribution were incorporated into a 2D finite element (FE) model of cardiac bioenergetics. We first assumed that the distribution of phosphagens within the cardiomyocyte cross sections will be similar across the length of several sarcomeres. We further assumed that individual mitochondria that cluster together share a common inner membrane space (IMS) (Fig 3A). We made this assumption based on a previous EM imaging study by Picard et al. [21], who reported that adjacent mitochondrial outer membranes in cardiac myocytes are connected through electron dense inner-membrane junctions (IMJ). Picard et al. also found a close alignment between the cristae belonging to adjacent mitochondria along the IMJ.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Model assumptions and FE mesh generation from EM image cross-sections.

(A): We assumed that individual mitochondria that cluster together will share a common inner membrane space. (B): Based on this assumption, each mitochondrial cluster was divided into two regions: an IMS region, depicted in green and a matrix region, depicted in blue. A third region depicted in red covers the myofibrils present in the cell. (C): Image of a FE mesh developed from the EM image of cross-section 1 with the three different regions.

https://doi.org/10.1371/journal.pcbi.1006640.g003

Assuming there is metabolic connectivity between juxtaposed mitochondria, we divided every mitochondrial cluster present in a cross-section EM image into two regions: (i) an IMS region where we implemented partial differential equations (PDE) to simulate diffusion of phosphagens; (ii) and a mitochondrial matrix region with negligible metabolite diffusivity (Fig 3B). The mitochondrial matrix reactions associated with oxidative phosphorylation (OXPHOS) were simulated using an ordinary differential equation (ODE) model developed by Beard [22]. This model had explicit biophysical descriptions of the electron transfer chain, ATP synthase and substrate transport, but a phenomenological relation for the dehydrogenase flux from the TCA cycle. The Beard et al. ODE model was solved at each computation node on IMS regions and coupled with PDEs for metabolite diffusion within the IMS region. As such, we developed a FE mesh from EM images of section 2 using this protocol (Fig 3C).

We assumed a uniform distribution of various enzymes like cytosolic creatine kinase enzymes, adenylate kinase and contractile proteins over the myofibril region within the FE mesh. The equations describing ATP consumption in the myofibril was modelled as a Michaelis–Menten function previously developed by Wu et al [23]. Here, [Pi_C] denotes the cytosolic concentration of inorganic phosphate, R is a constant mass-action ratio and X_ATPase is a model parameter that can be varied to simulate steady state ATP consumption at various workloads. In order to simulate the action of enzymes such as adenylate kinase and various cytosolic (M-CK) and mitochondrial (mt-CK) creatine kinases, the model also incorporated equations describing the enzyme kinetics from a reaction-diffusion model of cardiac energy transfer described in [13].

We modelled free diffusion of Pi, PCr and Cr between the shared IMS region and the myofibrillar region in the segmented images. On the other hand, diffusivity of ATP, ADP and AMP were set to a value 1% of that in the IMS region, compared to their cytosolic values. This was done to simulate low permeability of the VDACs (voltage depended anion channels) to these large anions as suggested by previous studies [13, 24, 25]. We also assumed that none of these large molecules can diffuse through the mitochondrial matrix due to intramatrix barriers like cristae structure and protein complexes [26].

Owing to its small molecular size, O₂ was assumed to diffuse through the mitochondrial matrix with similar diffusivity as the myofibrils and IMS. Therefore, O₂ was modelled as a species diffusing equally through all the three regions of the EM images. We imposed a constant Dirchlet boundary condition of O₂ in the cell membrane to simulate the O₂ flow from a capillary into the cell. All the other diffusing species in the model were subjected to a no flux Neumann boundary condition at the cell membrane.

We further utilized a thermodynamically-consistent and energy-metabolite-sensitive model of cross-bridge cycling previously developed by Tran et al [27] to predict the force generation at different points in the myofibril area based on the concentration of the following phosphagens: ATP, ADP and Pi. The overall model was solved using the finite element method. A detailed description of the model formulation and choice of various parameters and initial/boundary conditions can be found in Supplementary Material (S1 Text).

Before using our model to predict the metabolic landscape of the cells, we validated the phosphocreatine shuttle and substrate feedback mechanism that we have implemented in our model. Previously, Vendelin et al.’s [13] compartmental model of a rat heart cell suggested that PCr/ATP ratio remains effectively constant over low and moderate workloads but declines significantly at high workloads (Fig 4A). As evident from Fig 4A, our model was able to reproduce similar behaviour for different levels of O₂ consumption. The temperature in our model is the same as that used in the study by Vendelin et al. (298 K). Vendelin et al. validated their model predictions of PCr/ATP ratio by comparing them with experimental PCr/ATP ratio values which were collected from in vivo NMR spectroscopy on several canine hearts [28]. We have also reproduced this data in Fig 4A for comparison with our model predictions. As evident from the figure, our model predictions were consistent with the experimental PCr/ATP ratios from canine hearts.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Model validation with experimental results.

(A): Comparison between our model-predicted, spatially averaged PCr/ATP ratio vs V_O2 with model predictions from Vendelin et al. [13] and experimental results from Valdur Saks et al. [28] (B): Comparison between our model predictions for spatially averaged ADP vs V_O2 with predictions from Vendelin et al.’s [13] zero dimensional model (C): Comparison between model predicted radial profiles of O₂ with intracellular PO₂ levels calculated by Takahashi et al. [29] in six different isolated cells.

https://doi.org/10.1371/journal.pcbi.1006640.g004

Our model also predicted a parabolic relationship between average cytosolic ADP concentration vs V_O2, similar to the prediction from Vendelin et al.’s computational study (Fig 4B). These results indicate accurate modelling of the mechanism of substrate feedback between the mitochondria and myofibrils. The V_O2 value of 100 μmol/min/g tissue in the plot corresponded to a high workload level of X_ATPase = 0.05 μM/sec. According to Vendelin et al.’s study [13], this V_O2 value also corresponds to around 60% of the maximal filling rate of isolated rat heart perfused according to the ‘‘working heart” protocol of Neely [3]. Therefore, we have used X_ATPase = 0.05 μM/sec to simulate high physiological levels of workload in all the subsequent simulations.

In order to make accurate predictions about spatial variation in metabolite concentrations or reaction rates, the model needs additional validation with experimental results, where spatial distributions of diffusing species in the cells have been measured using live imaging techniques (i.e. fluorescence microscopy). In a previous study, Takahashi et al [29] visualised intracellular gradients of myoglobin oxygen saturation using high spatial resolution spectrophotometry and further converted the observed oxygenation level into partial pressure of O₂ (PO₂) using linear regression of data and the Hill equation. Fig 4C represents a comparison between Takahashi et al’s data and model estimations for radial profile of PO₂ inside an isolated rat cardiomyocyte. In the experiment, individual cells were subjected to a superfusing O₂ pressure of 15.2 mmHg under a high O₂ consumption rate at room temperature. We simulated high O₂ consumption in a similar range in cross section 1. Subsequently, the model predicted image of PO₂ profile was reduced in resolution and convolved with a confocal microscope point-spread-function to generate a simulated confocal image (S1 Fig), similar to a previous study [30]. Fig 4C shows that the concentration profile in the simulated confocal image of O₂ distribution follows a similar profile to the radial PO₂ gradient experimentally measured by Takahashi et al. [29]. Both the model and experimental measurements demonstrate steep concentration gradients proximal to the sarcolemma and shallow gradients within the cell interior.

The overall analysis presented in Fig 4A–4C illustrates that our model successfully reproduces both the whole-cell averaged and spatially profiled experimental results described in literature.

Heterogeneous mitochondrial distribution leads to large gradients in metabolite concentration across the cross-section of the cell

Following model validation with the simulation results from cross-section 1, we used the same protocol to simulate cardiac energy metabolism in the three cross sections presented in Fig 1B–1D under a high workload (V_O2 = 100 μmol/min/g tissue). It was assumed that the cell is exposed to a normal oxygen pressure from capillaries, and the concentration of O₂ at the sarcolemma was maintained at 45 μM based on Mik et al.’s experimental study [31].

Profiles of cytosolic ADP/ATP ratio, Pi and ATP hydrolysis rate (J_ATPase), corresponding to these cross sections (Fig 5) show significantly large gradients (600–1000 μM) of Pi concentration across the cell. On the other hand, ATP and ADP show spatial gradients in the order of 10 μM, which is small compared to the average ADP level (approximately 100 μM) and negligible compared to the average ATP level (approximately 10 mM). The ATP hydrolysis rate, being a function of both Pi and ADP/ATP, varies by a margin larger than the ADP/ATP ratio variation but lower than the gradients of Pi. By comparing the different cross-sections (Fig 5), we also observe that section 1 has a higher average level of ADP and Pi compared to the other two sections for the same cytosolic X_ATPase activity. Section 2 and 3 also show a larger gradient in the V_ATPase. These differences can be explained by the lower mitochondrial area fraction of section 1 as shown in Fig 2D. Although section 1 is closer to section 2 than section 3 (12.5 μm vs 25 μm of longitudinal distance), section 1 is located near a branching point of the cell unlike the other two sections. Besides similar average values, cross-sections 2 and 3 also have very similar spatial profiles of Pi. The Pi concentration is higher in the upper corner of these cross sections and it gradually increases towards the lower left of the sections. These results are consistent with our initial assumption of symmetric metabolic landscape between different cross-sections of the cell.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Spatial variation in model predicted cytosolic phosphagen levels and reaction rates represented using colour spectrum.

All the figures were generated by using the same range of colour spectrum for a particular phosphagen. (A-C): EM images of sections 1–3 (D-F): Distribution of ADP and ATP in the cross-sections represented using the ADP/ATP ratio. The ADP/ATP ratio remains uniform throughout the cross-sections (G-I): Inorganic phosphate (Pi) shows large spatial variation at the order of 1 mM. (J-L): Spatial variation in the ATP hydrolysis rate in the myofibrils corresponding to X_ATPase = 0.05 μM/sec.

https://doi.org/10.1371/journal.pcbi.1006640.g005

In the Fig 6A we have presented a boxplot of the model predicted gradients of all the diffusing metabolite levels present in cross section 1. The differences in spatial gradients of phosphagens can be attributed to the activity of creatine kinase (CK) enzyme isoforms in the IMS and myofibrils. Based on previous work by Vendelin et al [13, 25], we assumed ADP, ATP and AMP diffusivity in the IMS region is 1% of their corresponding diffusivity in the cytosolic space. Despite restricted diffusivity, the estimated ADP and ATP gradients were less than one tenth of the gradients present in PCr and Cr (Fig 6A). These results point towards the significance of the cytosolic MCK enzyme reaction, which acts as a spatial buffer of the ADP/ATP ratio by converting excess ADP into ATP at the cost of the PCr reserve of the cell. We contrasted the PCr and Cr concentration at different regions of cross section 1 with the corresponding local mitochondrial density (Fig 6B). This model predicts that myofibril regions with lower mitochondrial densities contain a lower reserve of PCr. The model predictions are also consistent with the accepted role of the PCr shuttle, in facilitating a rapid metabolic exchange of the equivalent energy of ATP between the IMS and myofibrils.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Spatial variation in phosphagens concentration in cross section 1.

(A): Phosphagen concentrations were calculated at all the myofibril points in the mesh and subtracted from the spatially averaged concentration of each phosphagen (Phosg_ave). Box plots represent the statistical distribution of the difference between phosphagen concentration at each point and Phosg_ave value of each phosphagen. The zero value corresponds to the Phosg_ave. (B): Normalized concentration of PCr and Cr displayed as a function of local mitochondrial density at various points in the cell. PCr and Cr concentrations were normalized by a factor equivalent to the difference of PCr and Cr maximum and minimum concentrations in cross section 1.

https://doi.org/10.1371/journal.pcbi.1006640.g006

Force dynamics is insensitive to heterogeneity in concentration of phosphagens and distribution of mitochondria

To analyse the impact of the spatial distribution of energy phosphagens on mechanical force production, we have used Tran et al.’s [27] cross-bridge cycling model to simulate isometric twitches at each of the discretised nodes. The isometric twitches are activated by a prescribed Ca2+ transient. It was assumed that the average value of various metabolite concentrations remains unchanged during a twitch, and accordingly the behaviour of the cross section was modelled as a function of cytosolic ATP, ADP, Pi, MgATP, MgADP and pH. The peak twitch force (F_peak) was simulated (Fig 7A–7C) for the cross sections, as were twitch durations (Fig 7D–7F), defined as the duration that the force was above 5% of F_peak (t₉₅). F_peak shows slight spatial variation in sections 2 and 3 (~7% difference between maximum and minimum spatial values), while it is almost uniform across section 1 (<3% difference). Similarly, differences between the maximum and minimum values of t₉₅ is around 12% in section 2 and 3, while it is around 5% in cross-section 1. These results imply that the force production is relatively insensitive to the large gradients in inorganic phosphate observed in Figs 5 and 6. The results also predict that rapid diffusion of PCr and Cr, coupled with CK enzyme reactions in myofibrils and IMS, is sufficient to maintain a uniform force generation across the cell cross section despite the heterogeneous distribution of mitochondria.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Force dynamics of the cell corresponding to the model predicted phosphagens distribution.

(A-C): F_peak is uniform across section 1 and exhibits slight spatial variation in section 2 and 3. (D-F): t₉₅ shows similar spatial gradients across the cross sections.

https://doi.org/10.1371/journal.pcbi.1006640.g007

Diffusion of PCr & Cr is not sufficient to maintain uniform force dynamics under hypoxic conditions

According to previous studies [32, 33] on isolated rat cardiomyocytes, cardiac energetics is impacted when the extracellular O₂ diminishes below a margin of 20 μM, as the mitochondrial oxidative phosphorylation system is limited. In an in-vivo non-invasive study on wistar rat hearts [31], Mik et al. found an average PO₂ of 35 mm Hg (roughly equivalent to 45 μM concentration assuming an O₂ solubility coefficient 1.35×10−6 M/mmHg [34]) present over the surfaces of cardiac mitochondria, which is well above this margin. However, Mik et al. also found that a significant percentage (25%) of in-vivo mitochondria was exposed to PO₂ in the range of 10–20 mm Hg under normal O₂ supply. Reducing the O₂ supply to the heart resulted in an even larger fraction of cardiac mitochondria (40%) with low surface PO₂ in the range of 0–10 mm Hg. Mik et al.’s findings were consistent with predictions from previous modelling studies by Beard et al. [34] and Groebe et al. [35]. Both modelling studies predicted that cardiac myocytes in-vivo can be exposed to a variety of PO₂ depending on location of the cells relative to capillary networks.

Therefore, we investigated the effect of low oxygen on phosphagens distribution in the 2D cross-sections. As shown in the previous results, we first applied a Dirichlet boundary condition with O₂ concentration of 45 μM along the cell membrane, which represents the majority of cells present in heart. In the second simulation we imposed a 15 μM oxygen concentration along the cell membrane, representing the cells exposed to lower levels of oxygen supply (hypoxic condition). Both the simulations corresponded to the same high workload levels.

Profiles of O₂ concentration in the cardiomyocyte, corresponding to these two simulations, are shown in Fig 8A. In both the cases, the O₂ distributions show steep gradients from the sarcolemma. While in the first case O₂ concentration remains above 20 μM throughout the cell, our model predicts an anoxic cell core under hypoxic conditions. This prediction is consistent with observations of Takahashi et al. [29, 36] on isolated cardiomyocytes as well as Mik et al.’s study [31] on mitochondria in the in-vivo heart. To better understand the effects of an oxygen deprived cell core on mitochondrial electron transfer system (ETS), we have plotted the mitochondrial complex IV reaction rate (V_C4) as a function of the IMS space (Fig 8B). Complex IV is the main site of oxygen consumption in the ETS. As expected in normoxia, activity of V_C4 is nearly uniform across the cell cross section, with small variations due to mitochondrial cluster sizes. However, in the hypoxic simulation, a negligible amount of O₂ passes to the interior of the cell, which in turn inhibits the J_C4.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Spatial variation in model predicted mitochondrial reaction rates and O₂ profiles for cross section 1.

All the figures were generated by using same colour spectrum range. (A): Steady state O₂ profiles showing steep O₂ gradients in the vicinity of the sarcolemma. In the hypoxic case, it also leads to the formation of an anoxic core with negligible O₂ level (B): The reduction of oxygen at complex IV is inhibited in hypoxia, while it is uniform across the section in normoxia. (C): Hypoxia leads to a steep gradient in the model predictions for mitochondrial membrane potential (D): Under normoxia Complex V (F₁-F₀ ATP synthase) reaction rate is nearly uniform across the cell with slight variation near the larger mitochondrial clusters. However, ATP synthesis at Complex V is inhibited in the anoxic core under hypoxia.

https://doi.org/10.1371/journal.pcbi.1006640.g008

We also investigated the influence of intracellular O₂ profiles on mitochondrial membrane potential (Δψ) for both the normoxic and hypoxic simulations (Fig 8C). In the first case, Δψ was nearly uniform across the cell. This uniform distribution of Δψ is clearly in accordance with the uniform profile of V_C4 observed in Fig 8B. However, in the second case, Δψ is significantly depressed in the core of the cell. This spatial depression of Δψ can be attributed to the decreased ETS flux in the centrally located IMF mitochondria. The lower row in Fig 8D demonstrates the effect of this depression in Δψ on reaction flux of complex V. The rate of ATP synthesis at complex V follows the same spatial profiles as Δψ. As a result, ATP synthesis declines in the centre of the cell in hypoxia.

The diminished rate of ATP synthesis in the mitochondria at the core of the cell also had a significant impact on transportation of ATP/ADP though the phosphocreatine shuttle (Fig 9). For instance, the reaction rate of mtCK in the forward direction (ATP + Cr → PCr + ADP) was around 15% lower in the central parts of the cell compared to the periphery under hypoxia (Fig 9A). As a result, the model predicts a strong radial gradient in the ADP/ATP ratio (Fig 10A) in the hypoxic simulations. This contrasts with the simulations of normoxic situation where both the mtCK reaction rates and the ADP/ATP ratio shows negligible spatial variation (Figs 9A and 10A). Fig 10B similarly demonstrates a radial profile in the ATP hydrolysis rate for the second case. The abundance of ADP in the anoxic cell core also leads to an elevated cytosolic MCK reverse reaction rate (PCr + ADP → ATP + Cr) which is substantially higher in the anoxic core compared to the periphery of the cell (Fig 9B). Together these results suggest that activity of mtCK and MCK is not sufficient for maintaining uniform ADP/ATP under limited O₂ supply.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Spatial variation in rates of CK enzymatic reactions involved in phosphocreatine shuttle under normoxia and hypoxia.

All the figures were generated by using same colour spectrum range. (A): Steady state profiles of reaction rates of mitochondrial mtCK in forward direction (ATP + Cr → PCr + ADP) (B): Steady state profiles showing reaction rates of myofibrilar MCK in reverse direction (PCr + ADP → ATP + Cr).

https://doi.org/10.1371/journal.pcbi.1006640.g009

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Spatial variation in cytosolic metabolite concertation and force dynamics under normoxia and hypoxia.

(A): Steady state profiles of ADP/ATP ratio. ADP/ATP is significantly higher in the anoxic cell core under low O₂ supply (B): ATP hydrolysis rate is similarly lower in the cell core under hypoxic condition (C): Peak twitch force (F_peak) is significantly higher in the cell core under hypoxia (D): Increased F_peak is complemented by increased twitch duration in the cell core, which can lead to loss of energy in the form of shear strain.

https://doi.org/10.1371/journal.pcbi.1006640.g010

For the second simulation, F_peak is also higher in the anoxic core of the cell (Fig 10C) and the twitch durations are significantly longer (Fig 10D). The difference between the maximum and minimum spatial values of t₉₅ is around 23% in the second simulation compared to only 5% in the first simulation. This is a result of lower ATP, and higher ADP and Pi in the centre of the cell—which slows down the cycling rate, but at the same time, the cross bridges spend more time in the force-producing states (due to lower ATP), leading to higher force development. This result is consistent with several experimental studies [37, 38] as well as a recent modelling study by Tewari et al. [39] which suggests that force production increases with increasing ADP and decreasing ATP levels.

A consequence of having longer twitch durations is that at high stimulus frequencies, the twitch does not have time to completely relax resulting in a higher diastolic force. Thus, according to the model predictions, different regions of the same cell cross section will be experiencing different levels of force in both systole and diastole. This may lead to shear strain inside the cell and wastage of metabolic energy through production of heat–which might not be a favourable situation during a high cardiac workload. The effect of contraction on the cellular ultrastructure was not considered in the model. Contraction of the cell can change the cell dimensions, particularly by increasing the cell cross section area due to conservation of volume. If the hypoxic cell core is contracted longer due to a reduction in ATP, it is likely to lead to an increase in the average diffusion distance between the cell core and the cell periphery. Increases in diffusion distance can lead to a further decrease in the supply of ATP in the hypoxic cell core–thus creating a negative feedback loop that might contribute to irreversible contraction of the cell.

These results imply that unlike under normoxic conditions, diffusion of Cr and PCr and associated CK enzyme activities is not sufficient to maintain a uniform ATP/ADP concentration and force dynamics under hypoxic conditions.

Distribution of phosphagens in cell cross-sections

We have developed a computational model of cardiac bioenergetics based on EM images to study the impacts of the heterogeneous distribution of mitochondria on the energy metabolic landscape of the cell. Note that this differs markedly from the earlier held concept that cardiomyocyte mitochondria are arranged within rigid crystalline lattice networks, which is thought to enhance diffusion and exchange of compounds such as the phosphagens. When less structured mitochondrial arrangements are incorporated into models, our results reveal that significant gradients in the phosphagen concentrations are likely. While concentration gradients are currently difficult to observe directly in experimental studies, previous computational studies, defined over idealised crystal lattice-like geometries, predicted smaller concentration gradients between adjacent mitochondria and myofibril units. Our current study indicates that significant concentration gradients may occur over the entire cross-sectional area of the cell, that this is dependent on local mitochondrial densities at a given point in the cell, and that these gradients appear acutely sensitive to hypoxia, such as would occur in infarcts.

Distribution of phosphagens along the longitudinal axis of cell

Our model was based on 2D cross sectional images, but we expect similar results if a 3rd dimension along the longitudinal axis of the cell is added to the model. To elaborate, Figs 5 and 6 shows that spatial profiles of various metabolites and quantities like ATP hydrolysis rate and peak twitch force are influenced by two different structural factors–(i) total area fraction of mitochondria which determines the average concentration of phosphagens in cross sections and (ii) density distribution of mitochondria which determines the spatial distribution of phosphagens. The model also predicts nearly uniform ATP and ADP levels across the cell cross sections in normoxia due to spatial buffering by PCr and Cr. In a 3D simulation, we can similarly expect a nearly uniform ATP and ADP concentration across the cell volume irrespective of the axial distribution of mitochondrial columns. The average concentrations of the phosphagens will be determined by total volume fraction of mitochondria in the cell.

S2 Fig illustrates that different cross sections along a continuous length of the cell have similar mitochondrial area fraction as well as similar spatial variation in the mitochondrial density distribution. Since spatial distribution of metabolites like Pi, PCr and Cr depend on the distribution of mitochondria, similar mitochondrial distribution in the cross sections should also lead to similar metabolic landscapes in these cross sections. As a result, there will be negligible concentration gradient along the longitudinal axis except near areas where the myofibrils and mitochondria branch out. This also implies that a 3D simulation, which includes a 3rd dimension along the longitudinal axis of the cell, will not result in a metabolic landscape that is different from the current 2D results.

Significance of the phosphocreatine shuttle

In our model we have employed a previously postulated formulation of phosphocreatine shuttle where ATP, ADP and AMP exchange between mitochondria and myofibril is severely restricted [7, 13, 25]. As a result, more than 90% of the ATP/ADP transfer between the mitochondria and myofibril occurs via facilitated diffusion by Cr molecules instead of direct diffusion of ATP and ADP. This reliance on rapid diffusion of PCr and Cr helps the cell to maintain a uniform level of ATP and ADP in the cardiomyocytes. In the absence of the PCr and Cr facilitated diffusion, large gradients of ATP and ADP will be required to supply ATP to myofibril spaces which are far from mitochondrial clusters–thus giving rise to a non-uniform distribution of ATP and ADP. However, the cytosolic ATP pool not only acts as a source of ATP for the cross-bridge cycle, but it supports ion channel activities (e.g. SERCA pumps) and can also act as divalent metal ion chelator [40], crucial to the viability of the cell. According to a previous modelling study [41] SERCA pumps, which are intrinsic to cellular Ca²⁺ management and EC-coupling, appear to be highly sensitive to the ATP/ADP ratio in the cell. Therefore, it is vital for cardiomyocytes to maintain a uniform ATP and ADP concentration, in spite of the heterogeneous distribution of mitochondria. Our model thus suggests that phosphocreatine shuttle might be a necessary aspect of energy metabolism in cardiomyocytes with heterogeneously distributed mitochondria.

Effects of hypoxia

We used an O₂ diffusivity value (D_O2) of 2.41x10^-5 cm²/s in our current simulation based on previous modelling work of Beard et al. [34]. This diffusivity value corresponds to the overall oxygen diffusion through mammalian striated muscles [42], and can be considered as an average diffusivity through all the organelles present inside a cardiomyocyte. Therefore, we did not separately model O₂ transport through the transverse tubules which are near uniformly distributed along the Z disks [43]. We also neglected the effect of myoglobin facilitated oxygen diffusion since diffusivity of myoglobin is reported to be only a small fraction (<1/10_th) of the free O₂ diffusivity in cardiomyocytes [44]. Despite these assumptions, our model prediction of radial profiles in oxygen is consistent with previous experimental and theoretical works [29, 35]. Our model indicates that in a normoxic state, profiles of O₂ have no effect on the metabolism, and that the phosphagen distribution is dictated by the mitochondrial distribution.

However, in hypoxia, our model predicts a significant impact of O₂ profile on the phosphagen distributions and the effect of the heterogeneous mitochondrial distribution is minimal. The level of mitochondrial membrane potential in the core of the cell was substantially lower than the periphery of the cell. The ATP/ADP in the cell core was similarly depressed due to a lack of O₂. However, these predictions are inconsistent with findings from Takahashi et al.’s 2008 work [36] exploring isolated cardiomyocytes. Takahashi et al. showed that isolated cells sustain depleted, but near homogenous levels of ATP and ΔΨ, in spite of formation of an anoxic core, when the respiration rate is elevated under a hypoxic condition. Takahashi et al. also suggested that CK mediated diffusion of PCr and Cr from the sub-sarcolemmal space supports ATP supply and ΔΨ levels in the anoxic core for a substantial time interval before necrotic cell death. Our results contrast with this view, as they indicate that CK reaction and diffusion of PCr and Cr is not sufficient to generate a homogenous metabolic landscape under hypoxia. This inconsistency between our model predictions and experimental data suggest that other mechanisms may be at play to ensure uniform ATP supply. Other mechanisms could include glycolysis, creatine phosphate stores or ΔΨ conductivity between sub-sarcolemmal and inter-myofibrilar mitochondria. Glancy et al recently proposed that ΔΨ can be transmitted across the cell through mitochondrial networks in both skeletal [45] and cardiac muscle cells [11]. Transmission of ΔΨ might be facilitated through IMJs present between adjacent mitochondria. Glancy et al. also proposed that membrane potential conduction via mitochondrial reticulae is a dominant pathway for skeletal muscle energy distribution.

In the current study, we assumed that there is no conduction of membrane potential between adjacent mitochondria, and accordingly modelled ΔΨ as a spatially distributed non-diffusing variable without any direct interaction between adjacent mitochondrial points. With this modelling approach, we did not find any significant difference in the mitochondrial Δψ values across the cross section of the cell in normoxia. These results are consistent with Glancy et al.’s [11] observation that mitochondrial Δψ conductivity is much weaker along the cell cross section compared to the longitudinal axis of the cardiomyocytes. We propose that, in normoxia_, the rapid diffusion of cytosolic PCr and Cr is sufficient to maintain a uniform Δψ as well as myofibrillar ATP distribution across the cell cross section without the need for Δψ conduction. However, our results also imply that in a hypoxic cardiomyocyte, the phosphocreatine shuttle is not sufficient to maintain a uniform distribution of ATP. Δψ conduction might play a role in maintaining uniform levels of intracellular ATP across the cell cross-section. In the future studies, it needs to be investigated whether hypoxia can lead to formation of larger number of IMJs along the cell cross section which can subsequently lead to higher Δψ conductivity compared to normoxic condition.

Potential implications for diabetic cardiomyopathy

The total cellular PCr pool and activities of mitochondrial CK enzyme are severely impacted in several metabolic diseases, such as diabetic cardiomyopathy. According to previous studies [46–48], mtCK and MCK activities are suppressed by a substantial margin in diabetes (35%-50%). Without normal CK activities, homogenous ATP distribution within cardiomyocytes can be severely impacted due to heterogeneous mitochondrial distribution. Our previous studies also revealed that mitochondrial arrangement in type 1 diabetes is significantly irregular than control cells due to mitochondrial fission and clustering [9], which can also affect the distribution of cellular phosphagens. In future work, we will study these potential changes in metabolite distribution and their effect on force dynamics using a 3D geometry based FE model of cardiac bioenergetics. Future studies will also require a modification in the currently used model of mitochondrial membrane potential to incorporate the conduction of membrane potential.

Ethics statement

All animal procedures followed guidelines approved by the University of Auckland Animal Ethics Committee (for animal procedures conducted in Auckland, Application Number R826).

Tissue sample preparation and serial block face imaging

The tissue sample used for SBF-SEM imaging was collected from the left ventricular wall of a sixteen-week-old male Sprague-Dawley rat at the University of Auckland, New Zealand. We followed the guidelines approved by the University of Auckland Animal Ethics Committee (for animal procedures conducted in Auckland, Application Number R826) for the process.

Following euthanasia of the Sprague-Dawley rat, the heart was excised and quickly cannulated and connected to a Langendorff apparatus operating at a hydraulic pressure of 90 cm. We perfused the heart with a Tyrode solution including 20 mM 2,3-Butanedione Monoxime for 2–3 minutes. This was followed by perfusion with 2.5% glutaraldehyde, 2% paraformaldehyde and 50 mM CaCl₂ in 0.15 M sodium cacodylate buffer. We dissected a tissue block from the left ventricular free wall and stored in the same fixative that was precooled on ice for 2 hours. The fixative was subsequently replaced with 2% osmium tetroxide and .8% potassium ferrocyanide in 0.15 M sodium cacodylate and left overnight. Afterwards, we stained the block with ice-cold 2% uranyl acetate for 60–120 minutes. Following this, the block was washed of excess uranyl acetate and punched into a 1.5 mm diameter sample. The sample was then progressively dehydrated with ethanol followed by transition to room temperature in acetone. The sample was finally embedded in epoxy resin (Durcupan ACM resin from EM sciences).

We trimmed the tissue sample to a square block face of 1 mm² and 300 μm deep with fibers running perpendicular to the future cutting face. The sample was subsequently mounted on a Teneo VolumeScope (FEI,Hillsborough, USA) and imaged in low vacuum mode (50 mbar) with a 10 nm pixel size. We acquired a total of four thousand and nineteen serial sections of 50 nm in thickness. The sections were then aligned using IMOD. We rotated the sections by 90 degrees around the Y axis and binned the data to an isotropic voxel size of 50 x 50 x 50 nm. A detailed description of the animal procedures, sample preparation and imaging techniques is available at our previous publications [17, 18].

Model development and computational methods

As discussed earlier, the results section provides a brief description of the model development, while the supplementary document S1 Text contains description of all the PDEs, ODEs and algebraic equation used in the model. The finite element implementation of the model was developed using reaction-diffusion module of OpenCMISS [49], an open source software package for FE simulation. Besides PDEs, the OpenCMISS model also contains ODEs and algebraic rate equations that have been described in CellML [50] models–which were solved using the Strang operator splitting method. We offer the OPECMISS implementation of the FE model for research use at https://github.com/CellSMB/cardiac_bioenergetics.

The FE models were solved on an IBM iDataplex x86 system. Individual FE models contained 300,000 FE nodes on average and required 20 hours of runtime for simulating every 1000 milliseconds of heartbeat. The results in this paper correspond to steady state solutions at 3 seconds of heartbeat time.