An Analysis of the Truncated Bid- and ROS-dependent Spatial Propagation of Mitochondrial Permeabilization Waves during Apoptosis*

Apoptosis is a form of programmed cell death that is essential for the efficient elimination of surplus, damaged, and transformed cells during metazoan embryonic development and adult tissue homeostasis. Situated at the interface of apoptosis initiation and execution, mitochondrial outer membrane permeabilization (MOMP) represents one of the most fundamental processes during apoptosis signal transduction. It was shown that MOMP can spatiotemporally propagate through cells, in particular in response to extrinsic apoptosis induction. Based on apparently contradictory experimental evidence, two distinct molecular mechanisms have been proposed to underlie the propagation of MOMP signals, namely a reaction-diffusion mechanism governed by anisotropies in the production of the MOMP-inducer truncated Bid (tBid), or a process that drives the spatial propagation of MOMP by sequential bursts of reactive oxygen species. We therefore generated mathematical models for both scenarios and performed in silico simulations of spatiotemporal MOMP signaling to identify which one of the two mechanisms is capable of qualitatively and quantitatively reproducing the existing data. We found that the explanatory power of each model was limited in that only a subset of experimental findings could be replicated. However, the integration of both models into a combined mathematical description of spatiotemporal tBid and reactive oxygen species signaling accurately reproduced all available experimental data and furthermore, provided robustness to spatial MOMP propagation when mitochondria are spatially separated. Our study therefore provides a theoretical framework that is sufficient to describe and mechanistically explain the spatiotemporal propagation of one of the most fundamental processes during apoptotic cell death.

Apoptosis is a form of programmed cell death that is essential for the efficient elimination of surplus, damaged, and transformed cells during metazoan embryonic development and adult tissue homeostasis. Situated at the interface of apoptosis initiation and execution, mitochondrial outer membrane permeabilization (MOMP) represents one of the most fundamental processes during apoptosis signal transduction. It was shown that MOMP can spatiotemporally propagate through cells, in particular in response to extrinsic apoptosis induction. Based on apparently contradictory experimental evidence, two distinct molecular mechanisms have been proposed to underlie the propagation of MOMP signals, namely a reaction-diffusion mechanism governed by anisotropies in the production of the MOMP-inducer truncated Bid (tBid), or a process that drives the spatial propagation of MOMP by sequential bursts of reactive oxygen species. We therefore generated mathematical models for both scenarios and performed in silico simulations of spatiotemporal MOMP signaling to identify which one of the two mechanisms is capable of qualitatively and quantitatively reproducing the existing data. We found that the explanatory power of each model was limited in that only a subset of experimental findings could be replicated. However, the integration of both models into a combined mathematical description of spatiotemporal tBid and reactive oxygen species signaling accurately reproduced all available experimental data and furthermore, provided robustness to spatial MOMP propagation when mitochondria are spatially separated. Our study therefore provides a theoretical framework that is sufficient to describe and mechanistically explain the spatiotemporal propagation of one of the most fundamental processes during apoptotic cell death.
Mitochondrial outer membrane permeabilization (MOMP) 3 is centrally implicated in a sequence of molecular signaling events that control cellular life-death decisions. MOMP is required for the efficient induction of the apoptosis execution phase and requires the formation of pores that are composed of homo-oligomers of activated Bcl-2 family proteins Bax and Bak (1). Consequently, Bax/Bak deficiency confers resistance to both extrinsic and intrinsic apoptosis-inducing agents and cellular stresses, and gives rise to developmental defects associated with excess cell material and causes high rates (Ͼ90%) of perinatal death in mouse models (2,3).
The activation of Bax and Bak is regulated by a complex interplay with other Bcl-2 family proteins (1). These include anti-apoptotic members such as Bcl-2, Mcl-1, and Bcl-xL, as well as an additional group of pro-apoptotic proteins, the BH3only proteins. Bcl-2, Mcl-1, and Bcl-xL antagonize activated Bax/Bak, but also neutralize BH3-only proteins. The BH3-only proteins are the main activators of Bax/Bak, and individual members or combinations of BH3-only proteins transduce distinct pro-apoptotic stresses (4). Although other BH3-only proteins are either transcriptionally induced or post-translationally activated, the BH3-only protein Bid, a central mediator of extrinsic apoptosis, is activated by proteolysis. Caspases-8 and -10, aspartate-specific proteases activated in response to ligands that bind to cell surface death receptors, cleave fulllength Bid in an exposed, unstructured loop and thereby generate truncated Bid (tBid), a Bid fragment with increased apoptotic activity (5,6). tBid inhibits anti-apoptotic Bcl-2 family members and also directly activates Bax/Bak (5-7). The high mitochondrial affinity of tBid depends on the presence of cardiolipin, an important lipid-cofactor of tBid-induced Bax pore formation, and MTCH2/MIMP, a Bid binding partner, in the outer mitochondrial membrane (8 -10). Cardiolipin is predominantly found in the inner mitochondrial membrane, but its mobility is increased by peroxidation to cardiolipin hydroperoxide, a process initiated by reactive oxygen species (ROS) (11). Cardiolipin hydroperoxide more readily integrates into the outer mitochondrial membrane and thereby sensitizes to tBid membrane insertion and MOMP (11). The interaction of tBid with the outer mitochondrial membrane further increases the rate of ROS generation, probably by inhibiting state-3 respiration and ATP synthesis, and thereby augments cardiolipin peroxidation (12,13). Cardiolipin peroxidation also mobilizes cytochrome c for release through Bax/Bak pores (14), and the mitochondrial loss of cytochrome c release causes an additional increase in ROS production (15). Therefore, tBid and ROS signaling during apoptosis are positively coupled.
The process of MOMP has been kinetically described in living cells by time-lapse imaging of mitochondrial cytochrome c release. Cells from clonal populations initiate MOMP individually at different times, with MOMP events occurring asynchronously over a period of many hours (16). Once initiated, MOMP proceeds rapidly in individual cells, with the entire pool of cytochrome c redistributing throughout the cell body within 2-5 min (16,17). Because cells with a common mitotic history initiate MOMP in a very narrow time window of 5-10 min, the onset of MOMP can be predicted and monitored at seconds resolution with minimal phototoxicity (17). This highlighted that MOMP is spatiotemporally coordinated and frequently spreads as a wave-like signal through epithelial cells. MOMP waves were also observed to travel through larger cells (Ͼ200 m) such as cardiomyocytes (18) and syncytiotrophoblasts (19). These observations are of importance because they demonstrate that a local induction of apoptosis can be propagated through the entire cell body and thereby can contribute to an efficient and irreversible commitment to apoptosis.
The understanding of the spatial anisotropies that coordinate MOMP waves remains very limited at the present time. So far, two mechanisms based on apparently contradictory experimental data have been proposed for the spatial propagation of MOMP, with reaction-diffusion signaling of either tBid or ROS being suggested as causative factors (17,18). However, the conditions and model systems used in the related experimental studies differed substantially (intact versus permeabilized cells; small versus large cells). Here, we therefore investigated whether mathematical models of tBid-or ROS-driven MOMP waves can successfully recapitulate the MOMP signaling characteristics observed in both experimental conditions. Although each model could reproduce a subset of experimental findings, we found that only an integrated mathematical model of coupled tBid and ROS signaling can accurately describe spatiotemporal MOMP signaling. Furthermore, the combination of tBid and ROS signaling was required for efficient MOMP propagation at conditions where mitochondria are spatially separated.

Modification and Extension of a tBid Reaction-Diffusion
Model-The model code for spatiotemporal tBid reactiondiffusion signaling was modified and extended from a code that was documented as part of a prior study on spatiotemporal MOMP signaling (17). The model represents cells as a one-dimensional slab and was adapted to allow for the simulations of variable cell sizes as well as for simulations that reflect intact and permeabilized cells. Cell sizes were implemented with each modeled mesh point reflecting 1 m distance. For the simulations of intact cells, boundary conditions that prevent flux out of the modeled system were implemented.
These boundary conditions were lifted when modeling permeabilized cells. tBid was added into the system at mesh posi-tion 1 using a sigmoid input function that reflects experimentally determined Bid cleavage kinetics (20,21). The input function was initiated at time 0, and the local concentrations of tBid were normalized to a percentage scale, with 100% reflecting the maximum concentration expected for evenly distributed tBid. A detailed description of the normalization procedure was provided previously (17).
Implementation of a ROS Reaction-Diffusion Model-To develop a reaction-diffusion model for ROS-induced MOMP, space, time, and boundary parameterization was implemented as in the tBid reaction-diffusion model. The initial ROS input was positioned at mesh point 1, and ROS was added at a constant input rate of 1 arbitrary unit/s. ROS absorbed at neighboring mesh points was considered to initiate MOMP once a threshold amount of locally absorbed ROS was reached. If MOMP occurred, these mesh points contributed as additional ROS input sites. ROS inputs lasted for 2 min, taking into account that cells actively counteract ROS generation and remove damaged mitochondria. Depending on the modeled cell size, between 30 and 200 ROS input functions therefore contributed to the overall ROS distribution that was calculated according to Equation 3. The amount of ROS required to induce MOMP at a given mesh point was adjusted for size and boundary conditions for permeabilized cardiomyocytes. This allowed us to re-model the experimentally determined end-toend velocity of 0.19 m/s (18).

Implementation of a Combined tBid/ROS Reaction-Diffusion
Model-For the integration of the above models into one combined model, the tBid input was defined as the near-end master trigger for MOMP. MOMP induction was then coupled with the ROS input function. From there on, both tBid and ROS were considered to contribute toward inducing MOMP in more distant mesh points. For the modeling of discontinuous mitochondrial patterns (Fig. 5), mesh points were defined as devoid of mitochondria. These mesh points consequently did not contribute to tBid adsorption and neither served as a source for ROS generation. With 1 and 0 representing presence or absence of mitochondria at a mesh point, the periodic patterns used for the simulations were 1 (continuous presence of mitochondria), 10, 100, 1000, and 10000. For comparison, these mitochondrial patterns were also implemented in simulations with the tBid and ROS models, respectively.
Mathematical Simulations-All models were implemented as MATLAB (version 7.12.0 (2011a); The Mathworks, UK) code for numerical analysis. Partial differential equations were solved using the PDE-PE subroutine, which uses an adaptive step Runge-Kutta ODE solver (Gear74). Similar to our previous partial differential equation modeling studies (17,22), modeling was performed in one spatial dimension, as solvers for partial differential equations for a four-dimensional analysis of three spatial and a temporal dimension and anisotropic inputs are not available. Where needed, noise was reduced in ROS or tBid/ROS models by calculating moving averages. The MATLAB code for all model variants is available as supplemental Information 1.

The Comparison of Spatiotemporal MOMP Signaling in Distinct Experimental Settings Provides a Rationale for a Systems
Level Mathematical Analysis-To develop and evaluate mathematical models that can describe the subcellular dynamics of MOMP signaling, we compared the available data from two previous studies that quantitatively investigated the spatiotemporal propagation of MOMP, and also compared the associated experimental conditions and model systems used to gather these data (17,18). In our own prior work, we measured MOMP waves in living HeLa cervical cancer cells. We suggested that a transient and localized production of tBid, together with its Bcl-2 family protein interactions and tBid membrane adsorption could explain experimentally measured velocity profiles of MOMP (17). The MOMP waves in HeLa cells progressed swiftly through the cell body, and wave velocities increased with distance. In contrast, Garcia-Perez et al. (18) showed in large permeabilized and mitochondria-rich cardiomyocytes that tBid-dependent MOMP signals may be propagated by bursts of ROS, and that MOMP waves propagate with a constant velocity. The experimental model systems employed in these studies therefore varied substantially in factors that potentially influence the propagation of signals by reaction-diffusion mechanisms (cell sizes and distances traveled by the MOMP waves; open versus closed cell systems; Fig. 1A). We therefore developed mathematical models for both mechanisms of MOMP waves to evaluate these for their compatibility with the available experimental data. In particular, models were analyzed for their capability to reproduce the signaling dynamics observed in both experimental conditions, and transferability or non-transferability was used for experiment-based model validation or rejection (Fig. 1B).

Characteristics of Spatiotemporal MOMP Propagation by tBid-dependent Reaction-Diffusion
Signaling-To analyze tBid-driven MOMP propagation, we extended a previous partial differential equation model of tBid reaction-diffusion signaling (17). At its core, this model builds on the application of the Kolmogorov-Petrovsky-Piskounov equation (23) to the biological scenario described for HeLa cervical cancer cells.
Here, u(x,t) represents the distribution of free tBid in space and time. This distribution calculates from the diffusion of tBid according to Fick's second law of diffusion (24), with D tBid representing the diffusion coefficient of tBid (17), and from an additional reaction term u MOMP . This reaction term reflects the loss in free tBid due to interactions with other Bcl-2 family members and/or mitochondrial membrane interactions. u MOMP (x,t) was therefore considered to indicate the amount of tBid that directly or indirectly contributes to MOMP. In the modeled scenario, tBid is entered into the system at the near end of the modeled cell ( Fig. 2A), and the tBid influx was based on experimentally determined Bid cleavage kinetics (20,21). In the original reaction-diffusion model, the simulations covered a distance of 30 m to resemble typical HeLa cell sizes, and boundary conditions were implemented that prevented tBid from leaving the cell system (17). In the extended model used here, we introduced cell size as a variable and also allowed to select between calculations for open and closed cell systems. These extensions allowed us to upscale the model toward sizes of cardiomyocytes used in the experiments performed by Garcia-Perez et al. (18) and to investigate the consequences of cell permeabilization on spatiotemporal signal propagation.
As a control for correct model implementation, we first performed simulations for HeLa cell conditions. As only a fraction of the cellular tBid pool is required to induce MOMP (6,20) we plotted thresholds that reflect the mitochondrial adsorption of 5, 10, or 20% tBid to visualize the spatiotemporal profiles for MOMP propagation as well as the associated velocity profiles. As expected, our extended model successfully replicated the accelerating MOMP waves previously described for HeLa cells (Fig. 2B). We next investigated whether this model can successfully be transferred to re-model the experimental findings made in cardiomyocytes. To this end, we performed calculations for stepwise increases in cell size (30, 60, 120, and 200 m) for both intact and permeabilized cell systems (Fig. 2, C and D). Simulations for intact cells demonstrated that at cell sizes Ն120 m the critical amount of tBid to induce MOMP is no longer reached at the far end (Fig. 2C). When conducting these simulations for permeabilized cells, the spatial spread of MOMP was attenuated even further (Fig. 2D). We next calculated the corresponding velocity profiles for the MOMP waves (Fig. 2, E and F). As mentioned before, the MOMP waves accelerated for HeLa cell conditions (Fig. 2E), closely resembling velocities measured experimentally (17). However, velocities were found to drop with increases in cell size and eventually the far end was no longer reached (Fig. 2E). In permeabilized cells, MOMP waves never accelerated with distance traveled FIGURE 1. A, overview of experimental systems previously used to investigate the spatiotemporal spread of MOMP. Systems differed in size and boundary conditions. B, design of a systems modeling approach testing the plausibility and compatibility for tBid-or ROS-dependent propagation of MOMP waves. Systems models for tBid-or ROS-dependent MOMP propagation were developed for the described experimental model systems (HeLa cells or cardiomyocytes, respectively). Models were then transferred between conditions and tested for their validity using the reported experimental data.

tBid and ROS Dependence of MOMP Waves
and in large cells also did not reach the far end (Fig. 2F). Modifying the input kinetic of tBid to remodel different amounts of tBid-producing caspase-8 activity affected wave velocities but not the distance traveled (not shown). Of note, for large cells (200 m) the behavior of the models for intact or permeabilized cells converged (Fig. 2, C-F). These findings demonstrate that a model assuming a localized input of tBid is not sufficient to replicate the spatial propagation of MOMP in larger cell systems.
Characteristics of Spatiotemporal MOMP Propagation by ROS-dependent Reaction-Diffusion Signaling-We next developed a mathematical model in which sequential bursts of ROS are assumed to propagate the MOMP waves, reflecting conditions for large cells such as cardiomyocytes (18). Furthermore, to model an open system, boundary conditions were not defined. To initiate the MOMP wave, a local burst of ROS was used as a trigger at the near end of the simulated cell (Fig. 3A). The spatiotemporal spread of ROS was calculated by the following equation.
The diffusion coefficient D ROS for ROS species was set to 2.752 ϫ 10 3 m 2 /s, based on calculations that applied the Wilke-Chang equation for conditions of oxygen diffusion in aqueous environments at 37°C (25). The reaction term u ROS was parameterized based on the reactivity of ROS species by using a half-time t1 ⁄ 2 ϭ 0.0001 s that was averaged from values previously reported for eukaryotic cells (26,27). ROS spreading through the cell was considered to damage other mitochondria upon reacting. If sufficiently damaged, these mitochondria were simulated to respond with MOMP and ROS bursts identical to the initial near-end trigger, thereby propagating the MOMP wave through the cell (Fig. 3A). The mitochondrial threshold for ROS responsiveness was then chosen so that MOMP waves propagated with the experimentally reported end-to-end velocity of ϳ0.19 m/s (18). The resulting model outputs demonstrated that the MOMP wave could reach the far end of the cell, and also that overall the wave velocity remained largely constant (Fig. 3B), a finding in line with the reported experimental data (18). Analogous to the tBid reaction-diffusion model, we tested whether the ROS model can be transferred from large permeabilized cardiomyocytes to remodel MOMP waves in small intact cells. Our simulations demonstrated that reducing the cell size to HeLa cell conditions did not notably affect the velocity or propagation of the waves (Fig.  3C). Similar findings were observed when modeling large intact cells or small permeabilized cells (not shown). Changing the amount of ROS required for MOMP induction affected wave velocities, with lower velocities calculated for higher MOMP thresholds. Again, wave velocities were not affected by dis-tance traveled (not shown). Therefore, whereas the ROS model could reproduce wave propagation and velocity profiles described for permeabilized cardiomyocytes, it failed to reproduce the acceleration of the waves reported for intact living HeLa cells.
Combining tBid-and ROS-dependent Reaction-Diffusion Signaling Is Sufficient to Accurately Recapitulate the Spatiotemporal Propagation of MOMP-Both the tBid-and ROS-dependent MOMP models failed to fully reproduce the published experimental evidence. Given that tBid and ROS signaling are tightly interlinked biologically (11)(12)(13)15), we hypothesized that integrating both models into a combined tBid-ROS reaction-diffusion model might allow a more accurate reproduction of the experimental data. To develop a combined model, we considered that the reaction terms for both tBid and ROS contribute toward reaching the MOMP initiation threshold. At the near end, we considered tBid as the initial inducer of MOMP, with the responding mitochondria in turn producing ROS. The produced ROS then diffuses and reacts in parallel to the amount of tBid that spreads from the near-end input (Fig. 4A). We then conducted simulations for conditions of intact HeLa cells as well as for permeabilized cardiomyocytes and plotted the propagation of MOMP waves and the associated wave velocities. Corresponding to the published data, MOMP waves in small intact cells accelerated with distance (Fig. 4B). In contrast, the MOMP waves in large permeabilized cells spread with overall largely constant velocities over longer distances (Fig. 4C) and, as expected, similar results were obtained when modeling a large intact cell (Fig. 5E, left  panel). Therefore, whereas the tBid-and ROS-dependent models had limited explanatory power, the combined model correctly described the key characteristics of tBid-and ROSdependent MOMP waves.
Combined tBid-and ROS-dependent Reaction-Diffusion Signaling Confers Robustness to Spatial MOMP Propagation-In our previous simulations we assumed that mitochondria are present continuously along the entire modeled distance. However, mitochondrial networks and distributions are known to be highly plastic and inhomogeneous (28). We therefore analyzed how discontinuities in mitochondrial networks may influence the propagation of MOMP in the tBid-, ROS-, and combined tBid/ROS models. To this end, we implemented periodic gaps of increasing size between mitochondria in the models (Fig. 5A). Positions at which mitochondria were absent did not contribute to tBid membrane insertion and neither could serve as ROS sources. For the tBid model, we observed that MOMP failed to reach the far end for continuous mitochondrial networks, whereas gaps between mitochondria promoted the propagation of the MOMP wave to the far end (Fig. 5B). This can be explained by reduced mitochondrial tBid adsorption along the modeled distance, allowing higher amounts of tBid to FIGURE 2. Boundary conditions and cell size significantly affect tBID-dependent profiles of MOMP propagation. A, schematic overview of tBid-dependent reaction-diffusion modeling of MOMP propagation, reflecting conditions previously investigated experimentally in intact epithelial cells (17). tBid (green dots) is produced at the near end and, over time, spreads by diffusion through the cell body. Although diffusing, tBid is adsorbed at mitochondria, contributing to MOMP. B, spatiotemporal profiles of tBid amounts and velocity/distance profiles of MOMP waves in intact epithelial cells. Profiles are largely independent of the amounts of tBid required to induce MOMP. C, spatiotemporal profiles of tBid amounts calculated for intact cells with increasing sizes. D, spatiotemporal profiles of tBid amounts calculated for permeabilized cells with increasing sizes. E, velocity/distance profiles of MOMP waves calculated for intact cells with increasing cell sizes. F, velocity/distance profiles of MOMP waves calculated for permeabilized cells with increasing cell sizes. FEBRUARY 26, 2016 • VOLUME 291 • NUMBER 9 reach distant mitochondria. As a reference condition for subsequent integration with ROS-dependent signaling, we selected the scenario of a 200-m cell and a 20% tBid threshold (Fig. 5C). A behavior opposite to the tBid model was observed when using the model for ROS-induced MOMP. Here, increases in the distance between mitochondria reduced wave velocities and eventually prevented the signal to reach the far end of the cell (Fig.  5D). This can be attributed to the very short half-times of ROS, which limit its mean travel distance and reduce the strength of the MOMP-inducing signal at more distant mitochondria. When repeating these simulations with the combined tBid/ ROS model, the MOMP wave reached the far end of the mod-eled cell in all scenarios, with the characteristics of ROS-dependent MOMP propagation dominating for high mitochondrial densities and velocity patterns for tBid-dependent MOMP propagation dominating for discontinuous mitochondrial networks (Fig. 5E). These findings therefore demonstrate that the combined action of tBid-and ROS-dependent MOMP induction confers robustness against inhomogeneous mitochondrial distributions and ensures that MOMP can be induced at distant sites. The finding that high mitochondrial densities promote ROS dependence of MOMP waves is in line with measurements in cardiomyocytes (18), because these cells contain high amounts of mitochondria to meet the ATP demands of the  (18). ROS (red dots) is produced at the near end and spreads by diffusion through the cell body, triggering MOMP in neighboring mitochondria once a critical ROS threshold is reached. Mitochondria undergoing MOMP serve as an additional source of ROS, contributing to the propagation of the signal. B, spatiotemporal profile for the time required to reach the ROS threshold for MOMP in permeabilized cardiomyocytes, as well as the associated velocity/distance profile. C, spatiotemporal profile for the time required to reach the ROS threshold for MOMP in intact epithelial cells, as well as the associated velocity/distance profiles. cardiac muscle. Taken together, our findings therefore allow us to visualize the dependence of MOMP waves on tBid and ROS signaling for different cell sizes and mitochondrial densities (Fig. 6).

Discussion
Here, we performed a mathematical analysis of spatiotemporal MOMP signaling driven by tBid-or ROS-dependent reaction diffusion processes. We found that neither process alone is sufficient to replicate the entirety of the available experimental data. However, an integration of both mechanisms into a combined reaction-diffusion model readily allowed the replication of all available measurement data and conferred robustness against discontinuities in mitochondrial distributions. Our study therefore suggests that a coupled signaling mechanism in which tBid and ROS act

. A merged model of tBid-and ROS-dependent MOMP propagation replicates experimental findings in both epithelial cells and cardiomyocytes.
A, schematic overview of tBid-and ROS-dependent reaction-diffusion modeling of MOMP propagation. Both tBid (green dots) and ROS (red dots) contribute to signal propagation. B, spatiotemporal profile for the time required to reach the MOMP threshold in intact epithelial cells for the combined tBid/ROS model as well as the associated velocity/distance profile. Note that the wave velocity increases while approaching the far end. C, spatiotemporal profile for the time required to reach the MOMP threshold in permeabilized cardiomyocytes for the combined tBid/ROS model as well as the associated velocity/distance profile. Note that the wave velocity remains largely constant along the modeled distance. FEBRUARY 26, 2016 • VOLUME 291 • NUMBER 9

JOURNAL OF BIOLOGICAL CHEMISTRY 4609
cooperatively is propagating MOMP signals through the cell body.
The mathematical theory for diffusion and reaction-diffusion processes is particular well developed, with seminal contributions from the fields of physics, chemistry, and engineering disciplines having established a framework for a highly accurate treatment of these processes at well controlled in vitro conditions (23)(24)(25)29). Based on this background, we were able to implement multiple mathematical models, and to test whether the results of our simulations were compatible with data that were collected at different experimental conditions (17,18). The possibility to apply our mathematical models to different experimental conditions (such as intact and permeabilized cells as well as to cells of different sizes) allowed us to determine that only the combined action of tBid and ROS signaling is compatible with the experimental evidence. In contrast, a well controlled experimental investigation of the spatiotemporal propagation of MOMP beyond the currently available data, at sufficient spatial and temporal resolution and, in large cells, also across sufficient distance would be a very challenging task. For example, trade-offs need to be made between the observable field of view (number of cells and observable distance), magnification, and resolution in space and time, while minimizing bleaching and phototoxicity of MOMP indicators at the same time. Because MOMP is initiated at apparently random time points, which can stretch over many hours between individual cells within cell populations treated with common apoptosis inducers, MOMP used to be investigated at a temporal resolution of minutes (16). However, a resolution of seconds is required to capture the characteristics of MOMP and MOMP waves in epithelial cells. This can be achieved in proliferating cells, where sibling cells initiate MOMP highly synchronously. MOMP in one sibling cell can therefore serve as a predictor of MOMP initiation in the other sibling cell, which then can be imaged for short periods at high spatiotemporal resolution and in the absence of detectable phototoxicity or fluorophore bleaching (17,22). However, this approach cannot be applied to postmitotic cells, such as neurons or cardiomyocytes. Instead, permeabilization of the cell membrane is routinely used to externally prime such cells for MOMP induction within a narrow, predictable time window (18). This comes at the cost of eliminating the cell membrane as a diffusion barrier, which significantly affects the characteristics of spatial signal propagation. To experimentally investigate the robustness of MOMP propagation identified in Fig. 5, it would be necessary to establish conditions in which tBid and ROS signaling are uncoupled. Although it could be assumed that this can be achieved by adding external ROS scavengers, we noted in control experiments that the absence of ROS also impairs the activation of caspase-8 and therefore the production of tBid. In these experiments we used mouse embryonic fibroblasts deficient in Bax and Bak expression and expressing an IETD Förster resonance energy transfer probe (30,31) and induced caspase-8 activation by 100 ng/ml of tumor necrosis factor ␣/1 g/ml of cycloheximide in the presence or absence of the ROS scavenger Tempol (5 mM). FRET probe cleavage, indicative of caspase-8 activation, was reduced by ϳ50% within the first 14 h of treatment (not shown). Therefore, due to the challenges and limitations associated with the experimental approaches, the application of a mathematical modeling strategy to investigate spatiotemporal MOMP signaling is helpful to mechanistically understand this process and to assess the impact different conditions may have on signal propagation.
However, our mathematical analysis also has limitations. First, our model is limited to one spatial dimension (distance). Even though for symmetrical systems the simplification to investigate reaction-diffusion processes in one spatial dimension usually provides good approximations, our model cannot capture effects arising from variability in cellular three-dimensional shapes. In the future, more advanced spatiotemporal modeling approaches may therefore be required to investigate the impact of cell shapes and also the impact of other features, such as intracellular diffusion barriers (for example, other organelles and their membrane barriers that need to be circumvented). Approaches by which the intracellular space is discretized into voxels and that employ rules-based logic modeling appear to be attractive strategies to achieve this (32,33). Such modeling approaches may also be suited to take into account complex changes in mitochondrial connectivity and trafficking that go beyond the simple scenarios of mitochondrial discontinuities analyzed in our study. For example, it has been described that apoptotic stress conditions can cause a subcellular re-distribution of mitochondria and may promote mitochondrial fission (34,35), both of which may differentially affect tBid-and ROS-dependent MOMP wave propagation. For tBiddependent waves, mitochondrial fission could possibly favor greater propagation. In contrast, fusion may contribute toward mitochondria escaping MOMP. However, this would be counteracted by ROS signaling in tBid/ROS co-dependent MOMP waves. Fusion and fission also cause additional effects through altering mitochondrial membrane curvature. Changes in cur-   FEBRUARY 26, 2016 • VOLUME 291 • NUMBER 9 vature affect membrane insertion of Bcl-2 proteins and mitochondrial MOMP responsiveness, with fragmented mitochondrial networks being less susceptible to Bax-induced pore formation (36,37). Changes in mitochondrial shapes likewise could be taken into account using more complex modeling approaches than employed here. Consequences on the diffusive spread of reactants by molecular crowding have been suggested to considerably influence intracellular diffusion and signal spread for many years. However, deviations from Brownian diffusion behavior arising from molecular crowding within cells appear to be far less pronounced than previously thought (38).

tBid and ROS Dependence of MOMP Waves
Interestingly, studies by Gottlieb and colleagues (39,40) demonstrated that upon induction of extrinsic apoptosis, caspase-8 can redistribute from death-inducing signaling complexes to mitochondria to locally produce membrane-associated tBid and drive MOMP. Our major conclusions would be valid also for such a scenario, because it could be modeled by the spatiotemporal spread and mitochondrial adsorption of caspase-8, with mitochondrial caspase-8 then acting as a local sink where Bid is converted to tBid. However, because caspase-8 is first activated at death-inducing signaling complexes and in living cells is capable of cleaving its substrates there prior to being active at the mitochondria (41), it is currently unclear whether mitochondrial caspase-8 activities contribute to MOMP waves. Furthermore, it is not yet known how mitochondrial caspase-8 activities quantitatively compare with caspase-8 activities at the death-inducing signaling complexes. Mitochondrial caspase-8 activities may contribute at later stages to ensure the irreversibility of tBid-dependent MOMP decisions, thereby preventing cell recovery. A similar role could be attributed to caspase-3, an effector caspase that can cleave Bid but which is activated several minutes after MOMP has been triggered (42)(43)(44). Due to the delay between MOMP initiation and caspase-3 activation, the MOMP wave reaches the far end of small epithelial cells before effector caspases are activated (17,22). A feedback of effector caspases onto Bid activation therefore does not need to be considered in this scenario. For the situation of large postmitotic cells, such as the cardiomyocytes considered in our study, effector caspases could in principle be activated before the wave reaches the far end. However, differentiated cardiomyocytes and other long lived postmitotic cells such as skeletal muscle cells and neurons, deregulate their balance of XIAP to Apaf-1 expression so that effector caspases can no longer be activated (45)(46)(47). Therefore, also for the setting of large postmitotic cells effector caspases cannot contribute to the propagation of MOMP waves.
Taken together, our study mechanistically explains the observed discrepancies in MOMP propagation in different experimental systems and highlights the requirement for cooperative tBid and ROS signaling to efficiently and robustly propagate MOMP waves over long distances. Building on this, in the future our model could be coupled with other mathematical models that more explicitly describe the complex interplay of MOMP-regulating Bcl-2 family proteins as well as the biochemical processes involved in the generation of mitochondrial ROS (48 -52). This could provide further insight into the regu-lation of one of the most fundamental processes of apoptotic cell death signaling.
Author Contributions-S. F. J. performed the programming and mathematical data analyses, data interpretation, and contributed to study design and manuscript writing. M. L. W. contributed to data interpretation and study design. M. E. D. performed the experiments and contributed to manuscript writing. M. R. conceived, coordinated and designed the study, contributed to data interpretation, and wrote the manuscript. All authors approved the final version of the manuscript.