Loss of mitochondrial function is a fundamental determinant of cell injury and death. In heart cells under metabolic stress, we have previously described how the abrupt collapse or oscillation of the mitochondrial energy state is synchronized across the mitochondrial network by local interactions dependent upon reactive oxygen species (ROS). Here, we develop a mathematical model of ROS-induced ROS release (RIRR) based on reaction-diffusion (RD-RIRR) in one- and two-dimensional mitochondrial networks. The nodes of the RD-RIRR network are comprised of models of individual mitochondria that include a mechanism of ROS-dependent oscillation based on the interplay between ROS production, transport, and scavenging; and incorporating the tricarboxylic acid (TCA) cycle, oxidative phosphorylation, and Ca2+ handling. Local mitochondrial interaction is mediated by superoxide (O2. −) diffusion and the O2. −-dependent activation of an inner membrane anion channel (IMAC). In a 2D network composed of 500 mitochondria, model simulations reveal ΔΨm depolarization waves similar to those observed when isolated guinea pig cardiomyocytes are subjected to a localized laser-flash or antioxidant depletion. The sensitivity of the propagation rate of the depolarization wave to O2. − diffusion, production, and scavenging in the reaction-diffusion model is similar to that observed experimentally. In addition, we present novel experimental evidence, obtained in permeabilized cardiomyocytes, confirming that ΔΨm depolarization is mediated specifically by O2. −. The present work demonstrates that the observed emergent macroscopic properties of the mitochondrial network can be reproduced in a reaction-diffusion model of RIRR. Moreover, the findings have uncovered a novel aspect of the synchronization mechanism, which is that clusters of mitochondria that are oscillating can entrain mitochondria that would otherwise display stable dynamics. The work identifies the fundamental mechanisms leading from the failure of individual organelles to the whole cell, thus it has important implications for understanding cell death during the progression of heart disease.
The spatial and temporal organization of the mitochondrial network is crucial for understanding its function in cells [1], [2]. Complex spatiotemporal factors not only contribute to physiological signaling [3]–[6], but determine the fate of the cell under stress. In heart cells, we have previously studied how the lattice-like packing arrangement of the mitochondrial network lends itself to propagation of bioenergetic signals [7] and to synchronization of self-organized oscillations in ROS, redox potential, and mitochondrial inner membrane potential (ΔΨm) in response to pathological stimuli [8]–[10]. Moreover, the scaling of instability of the mitochondrial network to the whole-cell and whole-organ levels has been shown to underlie the electrophysiological and contractile dysfunction associated with cardiac disease [11]–[13]. Only recently have we begun to understand the mechanisms responsible for stable, unstable, and oscillatory modes of mitochondrial bioenergetics, by combining experiments with mathematical model development. In studies of isolated cardiomyocytes subjected to localized oxidative stress, we elucidated the role of ROS in triggering autonomous synchronized oscillations of mitochondrial energetics [8], [9]. An emergent low frequency (0. 01Hz), high amplitude, limit cycle oscillation of ΔΨm was observed that spanned the length and breadth of the cardiomyocyte, but only after a considerable delay following the initial perturbation (local laser flash). We found that the abrupt transition in ΔΨm was preceded by the gradual increase of oxidative stress in the network, more specifically, when ∼60% of the mitochondria accumulated ROS to a threshold level. We referred to the state of the mitochondrial network just before depolarization as the point of “mitochondrial criticality” [9]. In the critical state, a small perturbation anywhere in the network can lead to the propagation of a ΔΨm depolarization wave. The mechanism of this phenomenon involves a mitochondrial ROS-induced ROS release mechanism (a term originally coined by Zorov et al [14], [15]), triggering an energy dissipating inner membrane anion channel (IMAC) that is inhibited by mitochondrial benzodiazepine receptor (mBZR) ligands, but not by inhibitors of the permeability transition pore [8]–[10]. A computational model of a ROS-dependent single mitochondrion oscillator was constructed to study the dynamics of the system [16] and it could reproduce the main features observed experimentally, including bursts of ROS release to the extramitochondrial space during rapid ΔΨm depolarization. Interestingly, in addition to the slow, large amplitude oscillation mode most closely associated with the pathological state, the model displayed a wide range of frequencies and amplitudes, suggesting the possibility that mitochondrial ROS release may be operating under physiological conditions as well. A subsequent study provided experimental evidence to support the idea that the mitochondrial network is organized as a collection of weakly coupled oscillators under physiological conditions that couple more strongly under stress to produce slow, synchronized waves and oscillations [3]. Because the initial collapse of the mitochondrial network initiates a cascade of events including activation of sarcolemmal KATP channels, alteration of the electrical excitability of the cardiomyocyte, and ultimately cardiac arrhythmias [11], [17], it is important to gain a precise understanding of the mechanisms responsible for signal propagation between mitochondria in the network, in order to find ways to interrupt potentially catastrophic events. Thus, in the present work, we test the proposal that O2. − diffusion from one mitochondrion to its nearest neighbors is responsible for the propagation of ΔΨm depolarization waves and synchronization of oscillation in the mitochondrial network by building one dimensional (1D) and two dimensional (2D) reaction-diffusion models of the mitochondrial network. We show that model simulations closely resemble ΔΨm depolarization wave propagation observed in experiments in response to oxidative stress, and remarkably, we find that a few oscillating mitochondria can entrain the entire network into an oscillatory mode even if the majority of the mitochondrial lattice is not in the oscillatory parametric domain. We propose that perturbing a minimal number of mitochondria is sufficient to trigger cell-wide responses through ROS-dependent coupling of mitochondria in the network.
The present results demonstrate that the observed emergent macroscopic properties of the mitochondrial network can be reproduced in a reaction-diffusion model of RIRR. Moreover, the simulations have uncovered a novel aspect of the synchronization mechanism, which is that clusters of mitochondria that are in the oscillatory domain of the parametric space can entrain mitochondria that would otherwise display stable dynamics. Propagation of electrical and Ca2+ signals within and between cells in the heart is essential for cardiac function, but less is understood about propagation of signals between organelles [18]. Several hypothetical mechanisms have been previously suggested to explain how ΔΨm depolarization may spread throughout the mitochondrial network of the cardiac myocyte. Amchenkova et al [19] proposed that there may be direct electrical continuity between mitochondria; however, more recent data showing that individual mitochondria or groups of mitochondria can be depolarized with little effect on their immediate neighbors argues against direct coupling. Ichas et al [20] provided evidence for propagated mitochondrial Ca2+-induced Ca2+ release mediated ΔΨm depolarization that involved the activation of the PTP. A mitochondrial ROS-induced ROS release mechanism was described by Zorov et al [14] to explain how laser-induced depolarization generates a large burst of O2. − generation from the electron transport chain upon activation of the PTP, which can occur independently of Ca2+. In previous work we investigated whether or not the PTP and/or Ca2+ was involved in the mechanism of whole-cell ΔΨm oscillation ([8]; see [21] for a review). Briefly, several lines of evidence ruled out both the PTP and Ca2+ as playing a role in the mitochondrial oscillations observed in heart cells. With regard to Ca2+: (i) the myocytes were studied under quiescent, minimally Ca2+-loaded conditions, and no sarcomere shortening was evident; (ii) inhibition of the sarcoplasmic reticulum Ca2+ pump or mitochondrial Ca2+ handling did not influence flash-induced mitochondrial oscillations; (iii) extensive buffering of intracellular Ca2+ with 1 mM EGTA did not affect flash-induced oscillations. The possible contribution of the PTP was ruled out by the following evidence: (i) cyclosporin A did not block the transitions, and (ii) small (600 MW) fluorophores were not lost from the mitochondrial matrix upon depolarization. More recently, we have provided extensive evidence showing that IMAC and PTP open sequentially as a function of the glutathione redox status in permeabilized cardiomyocytes [22]. Studies from our laboratory demonstrated that in cardiomyocytes, mitochondrial ΔΨm depolarization and redox potential during metabolic stress can be highly synchronized throughout the mitochondrial network, displaying complex behaviors including wave propagation within and between cells [7] and limit cycle oscillations [7], [8], [16]. Using highly localized laser excitation of less than 1% of the cellular volume to induce O2. − release and ΔΨm depolarization, we showed that after several minutes of delay, spatiotemporally synchronized oscillations in ΔΨm, O2. −, NADH, and GSH [8], [9], [16], spanning the entire mitochondrial network of the cardiomyocyte, can occur. This phenomenon could be prevented or acutely reversed by interrupting mitochondrial O2. − generation, increasing antioxidant capacity, or blocking IMAC, and neither Ca2+ nor PTP opening were involved in this response [8]. These responses could be described by a computational model involving RIRR and the activation of IMAC [16]. Employing the RD-RIRR mitochondrial network model, the present work successfully reproduces several experimental observations, including: i) ΔΨm redox wave propagation and its spatial-dependence on O2. − diffusion, production, and scavenging, and ii) synchronization of independent mitochondrial oscillators. Interestingly, a new feature was revealed by the 2D model simulation - entrainment of mitochondrial ΔΨm oscillation in mitochondria that would otherwise show stable behavior and low O2. − production. Finally, another important achievement was the direct confirmation of a key component of the oscillatory mechanism and model, which was the experimental demonstration that O2. − itself can induce ΔΨm depolarization and mitochondrial O2. − accumulation in permeabilized cardiomyocytes. The role of O2. − was previously inferred from the effects of superoxide dismutase mimetic compounds [8]. Reaction-diffusion theory (pioneered by Turing [23]), as a basis for pattern formation in biological or chemical systems, emphasizes the importance of two components; an autocatalytic reaction producing a local product (mediator), and the transport of this product by diffusion away from the source. This process can give rise to spontaneous symmetry-breaking and the appearance of self-organized spatial patterns including waves and oscillations [24]–[30]. With respect to the present model, the reaction consists of the reduction of O2 to produce ROS (specifically O2. −) driven by mitochondrial electron donors (e. g. , NADH). The local concentration of O2. − around the mitochondrion is shaped by several others factors, including buffering by the antioxidant reactions and transport of O2. − across the mitochondrial membrane. As in many reaction-diffusion systems, the mitochondrial oscillator also displays an inhibitory, or a self-limiting, mechanism; the concentration of O2. − at the activator site on IMAC decreases during ΔΨm depolarization because i) the rate of scavenging by SOD increases as O2. − accumulates and ii) the driving force for O2. − production and transport diminishes. Diffusion of the O2. − to neighboring mitochondria is shaped by the O2. − diffusion coefficient (DO2. −i) and the amount of the O2. − scavenger enzyme, etSOD, which consequently determines the rate of propagation of ΔΨm depolarization through the network. As expected, increasing etSOD slowed down the depolarization wave. The rate of propagation of the depolarization wave in the model corresponded to 26 µm s−1 with low DO2. −i (on the order of 10−14), which compares well with the experimentally determined 22 µm s−1 [9] at 37°C. A restricted diffusion range of O2. − in cells is consistent with experimental data; however, the actual diffusion coefficient of O2. − in cells (with antioxidant systems disabled) has not been determined and is likely to be influenced by local reactions with other molecules and molecular crowding around mitochondria, which would decrease the effective volume and increase the viscosity of the medium. This assumption of restricted diffusion is represented by the low DO2. −i in the model. Scaling of local interactions in complex dynamic systems to produce global emergent behavior is common in physical, social, financial and biological networks [31]. The RD-RIRR model illustrates how local neighbor-neighbor interactions (1–2 µm distance) can lead to long distance spatiotemporal patterns. Mechanistically, in cell-wide mitochondrial oscillations, propagation is mediated by regenerative RIRR between neighboring mitochondria. Our previous work treated the network as a percolation lattice [9] and we postulated that near the percolation threshold of the system, any mitochondrion within the spanning cluster (i. e. , one which is close to the critical level of ROS accumulation in the matrix) can depolarize, producing a burst of O2. −, which diffuses to its neighbor with a particular spatial concentration profile that is a function of the rate of O2. − scavenging, to elicit a cell-wide response. A suprathreshold level of cytoplasmic O2. − must reach the neighboring mitochondrion for the response to be regenerative, eliciting the O2. −-mediated opening of IMAC. In the model, the open probability of IMAC increases as a function of O2. −, which alters the IMAC conductance versus ΔΨm relationship (see Eq. 50 in the Supplementary Materials). The initial IMAC opening leads to a regenerative increase in ROS because the increased energy dissipation accelerates respiration, in turn increasing the number of electrons shunting to O2. −, which crosses the inner membrane to further activate IMAC. It is important to note that the global response is very much dependent on both the arrangement of mitochondria in the network and also the number of mitochondria close to the threshold. Before the first ΔΨm depolarization, the balance of O2. − production to O2. − scavenging must approach the threshold of oxidative stress in approximately 60% of the mitochondria, at which point a small perturbation can cause the synchronous collapse (and/or oscillation) of ΔΨm in the mitochondrial network. We referred to this vulnerable state as mitochondrial criticality [2], [9]. A limitation of this first RD-RIRR model is that it does not entirely reproduce the collective properties of network close to the critical state. For example, when we trigger ΔΨm oscillations by laser flash, the local increase in ROS load on the network is enough to induce a gradual accumulation of oxidative stress throughout the cell, and the bulk of the network gets closer and closer to criticality. The first global wave of depolarization does not usually happen instantaneously after the trigger, but often occurs after several minutes of delay: it presumably originates from a region having a suprathreshold level of ROS to trigger IMAC opening, and neighbors that are also close to the threshold. The phase of initial spreading of the “oxidative load” after laser stimulation is a collective network property that is difficult to represent with this simple reaction-diffusion model (indeed, modeling the behavior of critical systems is a nascent science). At the edge of criticality, depolarization of even a single mitochondrion can evoke collapse of almost the entire network (Fig. 10, see also Supplemental Materials Video S4). Alternatively, we can induce a critical state of the mitochondrial network by depleting the antioxidant defenses of the cell (globally) with diamide (see [10]). In this condition, we hypothesize that a local increase of ROS above threshold in any part of the 3D network could evoke RIRR. Specifically, the fast synchronized ΔΨm depolarization-repolarization cycles and waves are the aspects of the phenomenon that are well-reproduced by the 2D RD-RIRR model. Another limitation of current 1D and 2D models is that in the actual experiments, the O2. − can, of course, diffuse both laterally and vertically to other mitochondria, that is, in all 3 dimensions. Presumably, the ROS will have a similar diffusion rate in all directions, although this might depend on the structural organization of the organelles, membranes, and myofilaments. A 3D network would introduce additional triggerable ROS sources above and below the mitochondrion, which could, theoretically, alter the rate of wave propagation. This could be explored in the future by extending the model to three dimensions; however, we do not think that the major conclusions of the study regarding the properties of the excitable system will be substantially different. The RD-RIRR simulations revealed a novel aspect of RIRR-mediated coupling. A small number of oscillating mitochondria could elicit oscillations in the entire network in both the 1D and 2D models, even though the parameters were set for stable behavior in most of the network (i. e. , the percentage of respiration leaking to ROS was below threshold for independent oscillation). This indicates that entrainment by forced oscillation might also occur in the mitochondrial network. In this case, the combined effects of subthreshold ROS release and ROS diffusion from nearby mitochondria exceed the threshold for a regenerative response. This behavior is evident as recruitment of more and more mitochondria into the oscillatory cluster after several cycles of oscillation (as shown in Fig. 8D and the corresponding movie in the Supplemental Materials, Video S3). While focusing on a specific mechanism of RIRR (i. e. , IMAC-mediated), the present findings also provide general theoretical support for mitochondrial communication via RIRR. The RD-RIRR model simulations confirm that O2. − diffusion occurring locally between neighboring mitochondria over a distance of a few microns is sufficient for propagation and synchronization of ΔΨm depolarization over a larger distance. RIRR involving PTP-dependent depolarization [15], [32] can be readily incorporated into the model in the future; however, at present it is unclear what additional factors besides ROS may be required to evoke this more pronounced, and typically irreversible, cell death event. In our own experiments, we have demonstrated that the IMAC mechanism occurs with mild-to-moderate oxidative stress (by GSH depletion or lower O2. − levels) while PTP activation occurs during more severe stress, e. g. with further GSH depletion [10], or, in the present study, at higher O2. − concentrations. Mitochondrial Ca2+ overload is also thought to be a requirement for PTP activation as well [33], [34]. Hence, it is necessary to gain a better understanding of the interplay between the ROS species and the mitochondrial Ca2+ load in the activation of the PTP in the intact myocyte before a comprehensive PTP activation model can be constructed. In summary, we show that the autocatalytic release of O2. − and diffusion between mitochondria are the essential components of propagated RIRR between mitochondria in a closely packed array like that found in the cardiomyocyte. The kinetics and emergent spatiotemporal patterns of ΔΨm depolarization in the RD-RIRR mitochondrial network are modulated by O2. − production, O2. − scavenging, and diffusion. The highly interdependent nature of the mitochondria as a network of oscillators also suggests that synchronization by forced oscillation may occur when only a part of the network is perturbed.
Cardiac cell injury and death is a key component of cardiac diseases such as heart failure or myocardial infarction, thus it is important to understand the earliest steps leading up to irreversible cell damage. Mitochondria are the organelles responsible for generating the energy required to keep the cell running, yet they are particularly vulnerable to damage by toxic byproducts of metabolism, which include reactive oxygen species (ROS). ROS wreak havoc on cellular functions by attacking proteins, lipids and DNA, so the cardiac cell has evolved sophisticated defenses to remove them. The work we present in this paper using a computer model of the mitochondrial network describes how ROS generated inside the cell can spread from one mitochondrion to the next in a positive feedback process known as ROS-induced ROS release. Understanding this process will help us to find ways to intervene in this catastrophic mechanism to prevent loss of cell function and the associated cardiac arrhythmias and contractile problems leading to sudden cardiac death.
biochemistry/experimental biophysical methods
cell biology/cellular death and stress responses
cardiovascular disorders/arrhythmias, electrophysiology, and pacing
physiology/cardiovascular physiology and circulation
biophysics/theory and simulation
pathology/cellular pathology
physiology/integrative physiology
biochemistry/theory and simulation
biophysics/membrane proteins and energy transduction
pathology/molecular pathology
cardiovascular disorders/myocardial infarction
biochemistry/membrane proteins and energy transduction
cardiovascular disorders/heart failure
A Reaction-Diffusion Model of ROS-Induced ROS Release in a Mitochondrial Network