Hydrogen peroxide diffusion and scavenging shapes mitochondrial network instability and failure by sensitizing ROS-induced ROS release

The mitochondrial network of cardiac cells is finely tuned for ATP delivery to sites of energy demand; however, emergent phenomena, such as mitochondrial transmembrane potential oscillations or propagating waves of depolarization have been observed under metabolic stress. While regenerative signaling by reactive oxygen species (ROS)-induced ROS release (RIRR) has been suggested as a potential trigger, it is unknown how it could lead to widespread responses. Here, we present a novel computational model of RIRR transmission that explains the mechanisms of this phenomenon. The results reveal that superoxide mediates neighbor-neighbor activation of energy-dissipating ion channels, while hydrogen peroxide distributes oxidative stress to sensitize the network to mitochondrial criticality. The findings demonstrate the feasibility of RIRR as a synchronizing factor across the dimensions of the adult heart cell and illustrate how a cascade of failures at the organellar level can scale to impact cell and organ level functions of the heart.

. Mitochondrial arrangement in a myocyte. (a) Example confocal microscopy image of canine cardiomyocyte with green stained sarcolemma and t-tubules, and red stained mitochondria. Adult canine cardiomyocyte image was obtained using a 60 × lens on a two-photon laser-scanning microscope with a Ti-sapphire mode-locked laser at 780 nm. Cell was stained with tetrarhodamine methyl ester (TMRM) to visualize polarized mitochondria (red image), as previously described 9 , and surface membranes were stained with Thioflavin S (green image) to visualize t-tubules. A bandpass filters at 605 ± 25 nm was used to collect the TMRM emission and a cutoff filter was used to collect light < 490 nm for Thioflavine S. (b) Zoomed in schematic of mitochondrial arrangement with myofilaments colored yellow.
Scientific RepoRtS | (2020) 10:15758 | https://doi.org/10.1038/s41598-020-71308-z www.nature.com/scientificreports/ tion in m ), while the stimulated mitochondria at the edge of the region (the immediate neighbors of the central mitochondrion) only partially depolarized to 115 mV after 10 s, and only a bit further to 105 mV after 100 s (Fig. 3b). The remaining mitochondria in the network did not depolarize for the simulation duration of 300 s (Fig. 3c). Figure 3c shows the evolution of the membrane potential of each mitochondrion in the network. Figure 4 presents the simulations with H 2 O 2 diffusion. Like before, the stimulation protocol involved 3 mitochondria in the center of the network with elevated SO stim (Fig. 4a). Unlike the simulations without H 2 O 2 diffusion, multiple mitochondria depolarized (Fig. 4b). First, at time 0, all mitochondria were polarized. Next, the 3 stimulated mitochondria depolarized after 4 s. The middle mitochondrion, which had the largest SO generation rate, permanently depolarized. The behavior of this mitochondrion matched the behavior found in experiments, where laser flash stimulated mitochondria irreversibly depolarized after a few seconds 2,9 . In contrast to the central mitochondrion, the surrounding two mitochondria were stimulated to a lesser degree and exhibited sustained depolarizations only 201 s after the laser flash. At 4 s in the simulation, these two mitochondria first underwent a transient depolarization lasting 2.5 s. Afterwards, the two mitochondria displayed partial depolarizations. Unlike the results of the simulations without H 2 O 2 diffusion, these mitochondria continued to gradually depolarize and eventually fully depolarized after 201 s. Thus, the response of the 3 stimulated mitochondria is represented by the thick vertical black line in Fig. 4c. Collectively, these simulation results reproduced the experimental observations where mitochondria gradually irreversibly depolarized within the laser flash stimulated region 2,9 . The black horizontal extrusions starting at 147 s in Fig. 4c represent the mitochondrial depolarizations that propagated outward from the stimulated mitochondria. This delay in the response of the bulk of the network to oxidative stress has been consistently observed in both laser flash stimulated mitochondrial networks and chemically stressed mitochondrial networks (The paper by Zhou et al. demonstrated such behavior in their supplemental video (S3)) 27 . In our simulation results, the extent of the propagating depolarizations ( Fig. 4c inset length where color represents the m value. The horizontal axis represents the spatial extent of the 1D mitochondrial network, and the vertical axis represents time. The white space in the plot represents the mitochondria at spans of time where they are polarized and have a normal membrane potential around 150 mV. Mitochondria that depolarized had a reduced membrane potential and are represented by the different colors depending on depolarization magnitude. The central mitochondrion, which is permanently depolarized, is represented as a black vertical line in this plot.  www.nature.com/scientificreports/ x r and x d ) gradually increased after each consecutive depolarization wave (Fig. 4c). This network behavior is consistent with behavior in the stimulation experiments where depolarization propagations did not, at first, extend globally throughout the network. The depolarization waves in our simulations repeated roughly every 30 s, which was within the range of normal frequencies found in experiments. Each of these outward propagating depolarization waves were then followed by inwardly-directed repolarization waves. This last-in/first-out depolarization/repolarization pattern exhibited a compressed diamond like shape (Fig. 4c). Our simulation results reproduced the bias in faster depolarization propagation velocity versus the repolarization propagation velocity ( Fig. 4c inset ( x d / t d ):( x r / t r )). The total time frame of the depolarization wave and repolarization wave took 3 s at the first oscillation, and increased to 5 s for the last observed oscillation. This trend of increasing the average depolarization duration of the mitochondria in these depolarization waves was also observed in experiments 2,9,27 . Mechanisms driving mitochondrial network failure patterns. From the simulation with H 2 O 2 diffusion, Fig. 5a displays the evolution of m alongside the evolution of ROS accumulation. We found that the mitochondrial network accumulated H 2 O 2 around the stimulus location ( Fig. 5a H 2 O 2 ). The H 2 O 2 concentration monotonically increased in the stimulus region over time for the entire simulation duration (300 s). The center mitochondrion had an average H 2 O 2 accumulation rate of 1.7µM/s. The immediate neighboring mitochondria had a similar average H 2 O 2 accumulation rate but only after a delay of 3-5 s, and remained lower in concentration by 25 µM compared with the stimulus region. This pattern of delay and reduced H 2 O 2 concentration was repeated for successive neighboring mitochondria with roughly the same step decrease of 25 µM per mitochondrion. This pattern ended when the next successive neighboring H 2 O 2 concentration reached the baseline (< ~ 10 nM) and thus marked the boundary of the extent of H 2 O 2 accumulation. A roughly linear H 2 O 2 gradient formed, extending from the stimulus location to this boundary. Altogether, both the size of the H 2 O 2 accumulation region and the H 2 O 2 of the mitochondria within the region increased over time.
In contrast, the patterns of SO accumulation did not form a simple linear gradient over time and instead displayed complex fluctuating patterns both along the string of mitochondria (i.e. horizontal direction in Fig. 5a) and over time (vertically in Fig. 5a), both in the matrix and the inter-mitochondrial space (SO m & SO i ). The laser flash stimulated mitochondria begins generating ROS at the start of the simulation and as a result, the matrix SO of the central mitochondrion quickly increased to above 150uM. At this concentration, the SO leakage into the intermembrane space triggered the central mitochondrion's depolarization at 4 s. The central mitochondrion then released more SO into the intermitochondrial space, reducing its matrix SO to 13.6uM. This SO concentration remained relatively stable, and only increased at a rate of 113 nM/s as long as the oxidative stress stimulus was maintained, reaching an observed maximum of 45 µM after 300 s. Interestingly, besides the initial release of SO at the onset of the central mitochondrion's depolarization, the mitochondria directly adjacent to the central mitochondrion exhibited low matrix SO (< 1 µM) compared to its outer neighbors for the first 221 s of simulation. These mitochondria displayed partial depolarizations, where only a slight reduction of m (from 151 to 148 mV) occurred. These partial depolarizations were caused by the partial conductance increase of IMAC, which was caused by SO stimulation that diffused from the central mitochondrion. SO generated in these partially stimulated mitochondria continuously escaped and prevented matrix SO accumulation. The generated SO instead transferred to the intermitochondrial space, raising the concentration of SO in that space to ~ 50 nM, which was substantially larger than the < 1 nM of the more distant polarized neighbors, but not large enough to induce positive feedback of RIRR bursting. The remaining mitochondria in the network displayed matrix SO accumulation patterns similar to their H 2 O 2 accumulation patterns, but only prior to their next depolarization event. Also, during this period before the next depolarization, intermitochondrial SO did not accumulate. Only the region of mitochondria with enough accumulated matrix SO (> 1 nM) depolarized. Sufficient accumulation of matrix SO was required such that a triggered release of matrix SO into the intermitochondrial space would induce further self-release, which then needed to be large enough to trigger depolarizations in neighboring mitochondria. This transport process resulted in a rapid transient increase in intermitochondrial SO depending on the amount of accumulated stored SO, and also a rapid decrease in matrix SO. After sufficient SO scavenging, the depolarized mitochondria repolarized once the intermitochondrial SO returned to baseline levels. These patterns of matrix SO accumulation and release continued after each subsequent depolarization propagation. Further, the length of the extent of the accumulation region increased with each cycle where the additional number of mitochondria recruited per depolarization ranged between 2 and 4. Interestingly, the extent of matrix SO accumulation region was contained within the extent of the H 2 O 2 accumulation region and the extent of intermitochondrial SO accumulation was contained within the depolarized region (Fig. 5b).
The H 2 O 2 scavenging system component GSH formed a depletion zone (Fig. 6a GSH) that enveloped the H 2 O 2 accumulation region (Fig. 6b GSH|H 2 O 2 ). NADPH, the energetic supplier to that system, also formed a depletion zone but was instead enveloped by the GSH depletion region (Fig. 6b GSH | NADPH). Interestingly, the growth of the GSH depletion zone underwent two phases. First, a depletion zone of ~ 12 um formed at the center over the first 40 s while NADPH was elevated. After the NADPH of that region depleted, the second phase of GSH depletion occurred and the GSH and NADPH depletion region gradually increased in extent over time. This change in extent growth is denoted by the dashed red line in Fig. 6a. Thus, depletion of the cascaded energetic drivers of the H 2 O 2 scavenging system was required for initiating H 2 O 2 accumulation and for further increasing the extent of this H 2 O 2 accumulation region. Figure 7 presents the results of simulations where the HP stim and SO stim rates were varied independently to compare different ratios of H 2 O 2 and SO generation rates. Like before, the HP stim was only assigned to the central laser flash stimulated mitochondrion. The different SO stim rates were assigned to the immediate neighboring stimulated mitochondria Scientific RepoRtS | (2020) 10:15758 | https://doi.org/10.1038/s41598-020-71308-z www.nature.com/scientificreports/ www.nature.com/scientificreports/ since the central mitochondrion's SO stim was defined by depolarization status protocol like described before. The case shown in Fig. 7 second row & second column involved the same stimulation parameters as in Fig. 4 where SO stim was 1.4 × 10 -5 and HP stim was 2 × 10 -3 . Reducing HP stim to 1.5 × 10 -3 resulted in the m patterns presented in Fig. 7 second row & first column. Reducing HP stim from 2 × 10 -3 to 1.5 × 10 -3 increased the delay of the first propagating depolarization by ~ 77 s, decreased the number of observed propagating depolarizations by over half, increased the time between depolarizations by ~ 30 s, and decreased the size of the extent of depolarizing mitochondria by ~ 30%. The start times for all propagating depolarization waves in each plot were different for each ROS combination (Table 1). No propagating depolarization waves were observed when HP stim was reduced to 1 × 10 -3 . Increasing SO stim decreased the delay in the first propagating depolarization wave by ~ 8 s. This decrease was explained by the increased rate of matrix SO accumulation of the stimulated mitochondria, which causes increased leak to the intermitochondrial space sooner, which causes intermitochondrial SO to reach the IMAC opening trigger threshold sooner. This additional flux into the intermitochondrial space can also reduce the matrix SO accumulation of neighbors while their neighbor's matrix SO levels are still low. The additional intermitochondrial SO slightly increases their neighbor's IMAC conductance that leaks out some matrix SO. This results in an increased time between subsequent depolarization waves by ~ 10 s, thus reducing the frequency of oscillations. In contrast, increasing HP stim from 1.5 × 10 -3 to 2 × 10 -3 decreased the time between depolarizations by about ~ 30 s, and therefore increased the propagating depolarization wave frequency. Since H 2 O 2 diffused more readily to the mitochondria near the H 2 O 2 source, those mitochondria accumulated H 2 O 2 and thus also matrix SO more rapidly.

Discussion
In this study we tested the hypothesis that H 2 O 2 diffusion underlies the progression of mitochondrial network failure. To do so, we used a computational model of a mitochondrial network. The model included a detailed representation of the ROS scavenging system, ROS production from both intrinsic origins and externally induced alterations, and ROS diffusion between mitochondria. We regionally stimulated the network with ROS and observed the network behavior with and without H 2 O 2 diffusion. We found that H 2  In ventricular cells, where there is less space between mitochondria versus neuronal cells, the dominant key messenger molecule in the RIRR phenomenon was found to be SO 8,9 . Propagating RIRR, where ROS release of a mitochondrion induces ROS release of a neighboring mitochondrion, occurs in ventricular cells when SO emission from a mitochondrion triggers opening of the IMAC of a neighboring mitochondrion 9,34 . The resulting expulsions of stored SO and other anions also depolarize the mitochondria's m , which is a mitochondrial failure. Thus, SO transport was believed to play the dominant role in producing mitochondrial network failure 2,9,34 . www.nature.com/scientificreports/ However, propagation of this RIRR process relies on SO transport between neighbors that is large enough to trigger another IMAC opening 27 . Any initial SO release dissipates in magnitude over distance according to the laws of diffusion, and thus continuation of this propagation requires regenerative SO release across the network. A mitochondrial network where enough oxidative stress has accumulated, such that any triggering of SO release can cause a cascade of SO releases across the network, has been termed "mitochondrial criticality" 34 . However, a chicken-and-egg problem arises when assuming SO is responsible for bringing the network to criticality because propagation of SO releases depends on prior SO accumulation. A different communicator of oxidative stress was required to bring the network to criticality prior to mitochondrial depolarization. In agreement with the observation that widespread CM-DCF oxidation (primarily an indicator of H 2 O 2 ) precedes the first cell-wide depolarization 2 , the present results show that H 2 O 2 can function as this alternative communicator of oxidative stress. Further, since H 2 O 2 diffusion was necessary to observe mitochondrial depolarizations across the network, H 2 O 2 communication underlies SO accumulation. As H 2 O 2 diffused and accumulated throughout the network, scavenging of SO via SOD was inhibited 35,36 , which enabled SO to accumulate within the mitochondria matrices to produce a "critical" network. Thus, H 2 O 2 still plays an important role in cell types where SO was believed to be the most important ROS in governing mitochondrial network failure.
Delay in network response to oxidative stress is from scavenging capacity depletion time course. While H 2 O 2 diffusion is the mechanism for how oxidative stress can communicate and spread throughout the network prior to mitochondrial depolarizations, it is unknown how this form of oxidative stress determines the progression of the network failure behavior. Laser stimulus-induced oxidative stress in a small region of the mitochondrial network consistently shows a delay of more than a minute before a cell-wide depolarization is observed 9,27 . Previous mitochondrial network models have been unable to reproduce this delayed behavior. Some type of "counter" mechanism must exist for any time-invariant system to display a delay effect, where the counter is some pool that is depleted or accumulated over time. Consistent with the behavior of a relaxation oscillator, when the counter crosses a threshold, downstream effects are triggered. Given a steady input of ROS to the mitochondrial network, we looked for patterns of system quantities gradually accumulating or depleting.
The mitochondrial network within ventricular cells contains a complex scavenging system that has a cascaded topology 37 . SO, which is an intrinsic by-product of the processes of the electron transport chain and the production of ATP, is continuously reduced by SOD to produce H 2 O 2 25 . H 2 O 2 is then scavenged by GSH peroxidase and thioredoxin-dependent peroxiredoxins to prevent H 2 O 2 accumulation 25 . Regeneration of reduced GSH and thioredoxin is accomplished through NADPH-driven reductases, while NADPH redox potential ultimately depends on mitochondrial intermediates provided by the Krebs cycle. Modeling helps to reveal the dynamics of these redox-coupled reactions, and our results demonstrated that, indeed, GSH and NADPH are gradually depleted over time and across the network, prior to observing mitochondrial depolarizations. Further, our simulations also displayed the delay in response of the network to oxidative stress. SO accumulation and therefore the sensitivity of the network to perturbation, are potentiated by H 2 O 2 accumulation, which occurs prior to GSH and NADPH depletion. The encasement of the region of matrix SO accumulation within the regions of H 2 O 2 accumulation (Fig. 6b) and the NADPH and GSH depletion regions (Fig. 5b), demonstrates the sequential cascaded topology of the scavenging system. Thus, the delay in response of the network is caused by the time it takes to deplete the NADPH and GSH reserves of the network, and simultaneously, to accumulate H 2 O 2 and SO throughout the network. Our model is able to represent the evolution of a mitochondrial network gradually overwhelmed by oxidative stress.
Modeling variations in mitochondrial network failure behavior. There are variations in the protocols for how oxidative stress is introduced to mitochondrial networks. The common methods of experiments include bathing the network in exogenous ROS solution 27 , providing a laser stimulus that locally produces ROS 9 , metabolically stressing the system via ischemia/reperfusion protocols 38,39 , or by reducing scavenging capacity of the system via targeted inhibitory chemicals 17,40 . Each method affects the distribution of oxidative stress and the overall response of the network. Models have been developed to try study these systems and reproduce some of the behavior of these experiments. A previous model by Yang 28 examined how the short range effects of SO and the longer range effects of H 2 O 2 could participate to independently activate IMAC and PTP, respectively; however, the impact of H 2 O 2 accumulation (with consequent NADPH and GSH depletion) on the SO triggering mechanism was not explored. In a subsequent model from Nivala et al. 3 , SO-mediated IMAC activation was modeled as a stochastic process that caused local propagations that could, depending on SO levels, sum to produce propagating waves of ∆Ψ m depolarization, but the effect of H 2 O 2 was not included. Also, models from our lab 27 have been used to study propagation mechanisms during oxidatively stressed conditions. These models, however, do not capture the delay in response of the network to oxidative stress or demonstrate the relationship between H 2 O 2 diffusion, the SO accumulation, or the ultrasensitive nature of the ∆Ψ m oscillations. By varying the ratio of H 2 O 2 and SO, our model shows that altering the H 2 O 2 stimulus modulates the growth rate of the oscillatory region, whereas altering the SO stimulus alters the frequency of oscillations. Full validation of the model will require accurate measurements of the different distributions of ROS across the network of a cell and examination of the impact of network topology on propagations and synchronization. Nevertheless, we believe that the specific dependencies of H 2 O 2 and SO on long-and short-range coherent behavior, respectively, should remain.
The present model supports the feasibility of long range ROS-dependent communication in the mitochondrial network, but does not exclude the possibility of direct physical communication between mitochondria, for example, through nanotunnel structures 41 . However, such structures have only been shown to span distances Scientific RepoRtS | (2020) 10:15758 | https://doi.org/10.1038/s41598-020-71308-z www.nature.com/scientificreports/ of ~ 1-5 µm (roughly 1-2 sarcomeres) within cardiomyocytes, which could not account for synchronization of oscillations over the entire length of the adult heart cell (> 100 µm). Glancy and Balaban have proposed connections over even larger distances 11 , however our model is based on experimental observations that mitochondria in adult myocytes predominantly behave independently 8,20,42 . Assuming mitochondria are connected in a reticulum implies strong electrical coupling that would alter the sensitivity of the network and change the dynamics/spatial patterns of the depolarization waves. A further modeling study would be needed to properly characterize these changes. Further study is also needed to see how the effects of H 2 O 2 communication are important in experiments with ischemia/reperfusion protocols, and how this network behavior scales up in tissue and in the whole heart. We hypothesize that the high mobility of H 2 O 2 is important in distributing oxidative stress at the tissue and whole heart levels and may be important for producing the heterogenous distributions of network oscillation phases that may play a role in arrhythmogenesis. Incorporating our network model into a cardiac tissue model can be used to investigate these unknowns.

Methods
Computational model. The model developed here uses as its basis the mitochondrial network model by Zhou et al. 27 , which describes self-organized m oscillations mediated by Ca 2+ -independent RIRR via IMAC. Currently, there is no precise identification of the molecular structure of IMAC, so like previous models 25,27 , IMAC is a model assumption. The model does not incorporate Ca 2+ -dependent activation of PTP, which can be sensitized to open by higher levels of sustained oxidative stress 40 . The model represents the mitochondrial network as a collection of closely apposed nodes (mitochondria), with inter-connections represented by SO diffusion. Each individual mitochondrion is represented by the model of Cortassa et al. 25 , where ordinary differential equations describe the evolution of the mitochondrial matrix and intermitochondrial space. The model includes descriptions of the TCA cycle, oxidative phosphorylation, Ca 2+ handling, SO production, SO transport via IMAC from matrix to intermitochondrial space, and both SO and H 2 O 2 scavenging.
We added new components to the single mitochondrion model (Fig. 2a) to enable us to better represent ROS scavenging. First, we added a representation of the mitochondrial matrix scavenger superoxide dismutase (SOD) to address its possible role in H 2 O 2 dynamics, as matrix SO is converted to H 2 O 2 . Second, the H 2 O 2 scavenging system was further developed to include the thioredoxin system, as described by Kembro et al 37 , since it supplements the glutathione (GSH) system in scavenging H 2 O 2 and thus affects H 2 O 2 dynamics. Third, we incorporated parameterized formulations of nicotinamide adenine dinucleotide phosphate (NADPH) dynamics. NADPH drives the production of GSH and thioredoxin and thus plays a role in maintaining H 2 O 2 scavenging capacity. We parameterized the size of the NADPH pool and its replenishment rate to set limits on the transient and steady state capacity for maintaining a functioning H 2 O 2 scavenging system. Finally, we removed the representation of catalase to highlight the effects of the energetically driven H 2 O 2 scavenging systems, because the catalase content in the cardiac ventricular myocyte is low 43 , and because local catalase concentrations within the different cell compartments are not known.
In addition to further developing ROS scavenging, we redefined some and added other ROS production components to the single mitochondrion model. In their model, Cortassa et al. represented the basal SO produced from the electron transport chain as a function of oxygen consumed. This model mechanism is a positive feedback loop where SO production increases as SO accumulates within the matrix. However, to control for the contributions of SO production rate, we redefined the mitochondrial intrinsic SO production rate to be parameterized as SO base (7.5 × 10 −6 mM ms −1 ). By doing so, we excluded electron transport chain driven oscillatory pressure and instead highlighted the pressure due to scavenging dynamics. We also added the two parameters for ROS stimulation, HP stim and SO stim , which are the respective H 2 O 2 and SO production rates during any given stimulation protocol.
Completing model development, we added H 2 O 2 diffusion between mitochondria to enable H 2 O 2 to act as a communicator of oxidative stress across the network (Fig. 2b). We used a higher diffusion coefficient for H 2 O 2 , as compared to SO, to represent the higher diffusibility of H 2 O 2 over SO between mitochondria 44,45 . Similar to Zhou et al. 27 , mitochondria in our network model were arranged in a line and spaced 1um apart, the average distance measured experimentally in cells 46 . A total of 100 mitochondria were represented in the model. This number was determined iteratively to be sufficiently large to observe all the morphological changes in the depolarization patterns of the network while simultaneously minimizing computational expense. Similar to Zhou et al. 27 , we assumed a no flux boundary condition for SO diffusion and reapplied this assumption for H 2 O 2 diffusion. Previous mitochondrial network models had assumed the same boundary condition for H 2 O 2 diffusion 28 .
Oxidative stress stimulation protocol. To test our hypothesis, we simulated an oxidatively stressed mitochondrial network for 300 s. This is a typical time span to observe network failure progression in oxidative stress experiments 9,27 . We introduced oxidative stress to the network by stimulating a small group of three neighboring mitochondria in the network (Fig. 3a). Stimulation of these mitochondria was represented by assigning a value to their SO and H 2 O 2 production rate (SO stim and HP stim respectively). The value of the stimulus was maintained constant during the entire simulation duration (300 s). To be able to compare simulation results to experimental, our stimulation protocol mimicked the stimulation protocol used in laser flash experiments, where a brief (~ 500 ms) localized laser flash induces sustained ROS production in a localized region of the network 9 . The central mitochondrion in the group of three was assigned the largest SO stim of 5 × 10 −5 mM ms −1 (Fig. 3a), which was iteratively determined to be the smallest value needed to induce a permanent depolarization in that mitochondrion. The remaining 2 stimulated mitochondria were assigned a smaller SO stim of 1.4 × 10 −5 mM ms −1 to represent the reduced laser intensity from being located at the edge of the laser flash region. Production of Scientific RepoRtS | (2020) 10:15758 | https://doi.org/10.1038/s41598-020-71308-z www.nature.com/scientificreports/ H 2 O 2 was represented by assigning the center mitochondrion's HP stim equal to 2 × 10 −3 mM ms −1 . To determine the effect of H 2 O 2 transport on the progression of mitochondrial network failure, we performed this stimulation protocol with and without H 2 O 2 diffusion. Further, to explore the effect of different ROS stimulus makeups, we repeated these simulations with different HP stim rates of 1 × 10 −3 mM ms −1 , 1.5 × 10 −3 mM ms −1 , and 2 × 10 −3 mM ms −1 , and different SO stim rates of 12 × 10 −6 mM ms −1 , 14 × 10 −6 mM ms −1 , and 16 × 10 −6 mM ms −1 . These alterations to the previous stimulus rates were iteratively determined to be large enough to demonstrate a change in network behavior.

Detail of mitochondrial model modifications.
The following describes the details of the model additions and changes from previous models. Listed variable names correspond to the variable names in the model equations file. We made changes to the Zhou model 27 to better represent ROS scavenging. First, we added a representation of the mitochondrial matrix SOD. We reused the same formulation of the scavenging rate via intermitochondrial SOD (VSOD_i) from the Cortassa model 25 for the matrix SOD scavenging rate (VSOD_m). The concentration of matrix SOD, compared to intermitochondrial SOD, was reduced from 1.43 uM to 300 nM to represent the reduced availability of SOD in the matrix due to the differences in volume of the two spaces. Second, we updated the parameters of the glutathione based H 2 O 2 scavenging system. The concentration for glutathione peroxidase (E__GPX_T), the concentration for glutathione reductase (E__GR_T), the K constants for GSSG (K_GSSG_M) and NADPH (K__NADPH), and total glutathione pool (G_tot) was taken from the 2013 Gauthier paper 47 , Phi_2, and k_GR were taken from the other 2013 Gauthier paper 48 because the first paper's 47 supplement values appeared to have a unit conversion error. Third, we included the thioredoxin/peroxiredoxin based H 2 O 2 scavenging system. The equations for this component was taken from the first 2013 Gauthier paper's supplement 47 . Similar to the glutathione changes, the values for Phi were taken from the other 2013 Gauthier supplement 48 . Finally, we parameterized the size of the NADPH pool, using a starting capacity value of 12 mM and a replenishment rate of 7.5 nM/ms (NADPH_gen).
In addition to SO diffusion, we added H 2 O 2 diffusion. The resulting equations for ROS flux were: SOi_(left|right) was the SO concentration of the neighbors. HP_(left|right) was the H 2 O 2 concentration of the neighboring mitochondria. The diffusion coefficient for SO (C_SO) was unchanged from the Zhou model 27 and 4 mM/ms/um was used for the diffusion coefficient for H 2 O 2 (C_HP).
With the changes just described, the resulting differential equations for the affected state variables were: VGR is the rate of production of glutathione from glutathione reductase using NADPH and the precursor GSSG. Similarly, VTRXR is for replenishment of thioredoxin. VGPX and VPRX are the H 2 O 2 scavenging rates from the scavenging systems governed by GSH and thioredoxin respectively. The constants used for VGR in the Zhou et al. model 27 were updated to the more recent values in Gauthier et al. 47,48 with the changes for the constants described above. VTRXR and VPRX was also added from the Gauthier models.
Computational methods. The nonlinear system of partial differential equations describing ROS communication between nodes was spatially discretized using the finite difference method. The aggregated ODEs of all the nodes were numerically integrated using the solver CVODE, a multistep stiff ODE solver that uses a banded backward differentiation formula method and a direct linear solver to implement newton iteration. A maximum time step of 0.1 ms was used to stably simulate the model on a desktop computer. Model source code is freely available for download on github: https ://githu b.com/bmill are/RIRR-propa gatio n-sensi tivit y-model