text:
When exposed to tumor necrosis factor (TNF) or TNF-related apoptosis-inducing ligand (TRAIL), a closely related death ligand and investigational therapeutic, cells enter a protracted period of variable duration in which only upstream initiator caspases are active. A subsequent and sudden transition marks activation of the downstream effector caspases that rapidly dismantle the cell. Thus, extrinsic apoptosis is controlled by an unusual variable-delay, snap-action switch that enforces an unambiguous choice between life and death. To understand how the extrinsic apoptosis switch functions in quantitative terms, we constructed a mathematical model based on a mass-action representation of known reaction pathways. The model was trained against experimental data obtained by live-cell imaging, flow cytometry, and immunoblotting of cells perturbed by protein depletion and overexpression. The trained model accurately reproduces the behavior of normal and perturbed cells exposed to TRAIL, making it possible to study switching mechanisms in detail. Model analysis shows, and experiments confirm, that the duration of the delay prior to effector caspase activation is determined by initiator caspase-8 activity and the rates of other reactions lying immediately downstream of the TRAIL receptor. Sudden activation of effector caspases is achieved downstream by reactions involved in permeabilization of the mitochondrial membrane and relocalization of proteins such as Smac. We find that the pattern of interactions among Bcl-2 family members, the partitioning of Smac from its binding partner XIAP, and the mechanics of pore assembly are all critical for snap-action control.
Author Summary
In higher eukaryotes, tissue development and homeostasis involves a subtle balance between rates of cell birth and death. Cell death (apoptosis) is triggered by activation of caspases, specialized enzymes that digest essential cellular constituents and trigger degradation of genomic DNA. Under normal circumstances receptor-dependent cell death is very tightly repressed, but it is irreversibly induced upon receipt of an appropriate signal. Mutations that interfere with this all-or-none control contribute to developmental abnormalities, autoimmune disease, and cancer. The biochemical properties of most apoptotic proteins are quite well understood, but it is unclear how these proteins work together. By combining live-cell microscopy, genetic perturbation, and mathematical modeling, we seek quantitative insight into cell death with a focus on network dynamics and control. We find that cells vary dramatically in the time between receipt of an apoptotic signal and the commitment to death. This variability arises from cell-to-cell differences in the activities of receptor-proximal biochemical reactions. Rapid all-or-none progress from commitment to actual death is achieved downstream by pro-apoptotic proteins found in the mitochondrial membrane. Our work provides a quantitative picture of apoptosis that advances understanding of oncogenic mechanisms and should eventually assist in the development of pro-apoptotic cancer therapies.
Citation: Albeck JG, Burke JM, Spencer SL, Lauffenburger DA, Sorger PK (2008) Modeling a Snap-Action, Variable-Delay Switch Controlling Extrinsic Cell Death. PLoS Biol 6(12): e299. https://doi.org/10.1371/journal.pbio.0060299
Academic Editor: Andre Levchenko, Johns Hopkins University, United States of America
Received: October 16, 2007; Accepted: October 20, 2008; Published: December 2, 2008
Copyright: © 2008 Albeck et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was funded by National Institutes of Health (NIH) grant GM68762 to PKS.
Competing interests: PKS is a member of the board of directors of Applied Precision Inc., manufacturer of the microscopes used in this study.
Abbreviations: C3, caspase-3; C6, caspase-6; C7, caspase-7; C8, caspase-8; C9, caspase-9; CyC, cytochrome c; DISC, death-inducing signaling complex; EARM, extrinsic apoptosis reaction model; EC-RP, effector caspase reporter protein; FRET, Förster resonance energy transfer; GFP, green fluorescent protein; IAP, inhibitor of apoptosis protein; IC-RP, initiator caspase reporter protein; IMS-RP, mitochondrial intermembrane space reporter protein; MOMP, mitochondrial outer membrane permeabilization; ODE, ordinary differential equation; RNAi, RNA interference; siRNA, small interfering RNA; TNF, tumor necrosis factor; TRAIL, TNF-related apoptosis-inducing ligand
Apoptosis is essential for the development of multicellular organisms but is misregulated in diseases as diverse as cancer and autoimmunity [1,2]. Activation of the potent effector caspases (caspases-3 and −7; hereafter C3 and C7), the hallmark of apoptosis, is triggered via the intrinsic cell death pathway by intracellular events such as DNA damage and oxidative stress, and via the extrinsic cell death pathway by extracellular stimuli such as TNF (tumor necrosis factor) and TRAIL (TNF-related apoptosis-inducing ligand) [3]. C3 and C7 directly degrade the proteome and, by activating DNAses, also dismantle the chromosomes of cells committed to die [4]. Caspase activation represents an irreversible change in cell fate and is consequently regulated at multiple levels, including assembly of complexes on death receptors [5], binding of pro- and anti-apoptotic members of the Bcl-2 family to each other in cytosolic and mitochondrial compartments [6,7], mitochondria-to-cytosol translocation of Smac and cytochrome c (CyC) [8–10], and direct repression of caspases by inhibitor of apoptosis proteins (IAPs) [11]. Studies of extrinsic apoptosis at the single-cell level reveal a long and variable delay prior to effector caspase activation but rapid and sudden progression to substrate cleavage once activation has begun [12,13], a behavior that we term "variable-delay, snap-action" switching. Failure in snap-action switching generates an indeterminate physiological state and sublethal cellular damage that may predispose cells to genomic instability [14,15].
A variety of kinetic models of apoptotic cell death have been published to date [16–26], most of which are based on coupled systems of differential equations. These models incorporate different mechanisms for achieving all-or-none caspase activation, including positive feedback via caspase-8 (C8) [17,20] or caspase-9 (C9) [22,23], and ultra-sensitivity in C9 activation [20]. Most models also focus on subsets of reactions, such as those stimulated directly by death receptors [17] or those downstream of mitochondrial outer membrane permeabilization (MOMP) [23] rather than on the interplay between upstream and downstream reactions. Moreover, most studies rely on previously published data rather than cycles of model-based hypothesis generation and experimental test. In this paper we attempt to overcome these limitations by including upstream and downstream reactions in a single model and by tightly coupling modeling and experimentation. The ordinary differential equation (ODE)-based model of C3 regulation we describe is based on mass action kinetics and has been trained against population-based and single-cell data obtained from cells perturbed by RNA interference (RNAi) and protein overexpression. We establish a role for C8 and its substrates in pre-MOMP delay and explore how competition among pro- and anti-apoptotic Bcl-2 family members determines when MOMP occurs. We also examine a series of alternative topologies for reactions involving Bcl-2-like proteins and uncover a subtle interplay between protein compartmentalization, translocation, and multimerization in the regulation of snap-action switching.
Modeling Extrinsic Cell Death Pathways
A mathematical model of proteins known to regulate C3 during extrinsic apoptosis was constructed on the basis of mass-action kinetics, with elementary reactions represented as ODEs (Figure 1; Protocol S1; Tables S1–S6). All biochemical transformations were represented as unimolecular or bimolecular reactions, and rate laws were therefore expressed as r = k × [A], for a reaction involving one copy of protein A, r = k × [A] × [B], for a bimolecular reaction of A and B, or r = k × [A] × [A], for dimerization of A (Tables S2 and S6). Transport between cellular compartments was also modeled as an elementary unimolecular reaction, and the assembly of multiprotein complexes as a series of bimolecular reactions. Because no complex algebraic forms such as Hill functions were used in our model, ultrasensitivity and other nonlinear behaviors arise from interactions among simple elementary reactions rather than the properties of higher-order equations. Where possible, estimates for model parameters (rates and initial protein concentrations) were obtained from the literature (Tables S4 and S5). In the absence of such information, kinetic rate constants and initial conditions were set to intermediate values within a physically plausible range and then fitted so as to optimize model performance [27].
(A) The convention of Kitano et al. [94] is followed. The major features of the network are highlighted by color: gray, receptor module; blue, direct caspase cascade; green, positive feedback loop; yellow, mitochondrial feed-forward loop.
(B) A condensed alternate representation of the network.
Well-substantiated biochemical reactions comprising four interacting cell death subcircuits were included in the "extrinsic apoptosis reaction model" described here (EARM v1.0), but some regulatory processes were simplified or omitted (Figure 1). EARM v1.0 contains 58 species corresponding to 18 gene products having nonzero initial conditions and 40 additional species representing complexed, cleaved, or differentially localized forms of the initial species, which interact via 28 reactions described by 70 nonzero rate constants (including forward, reverse, and kcat rates for each reaction; Tables S1 and S2). The four subcircuits in EARM v1.0 comprise (i) a lumped-parameter representation of receptor binding by TNF or TRAIL and the subsequent activation of pro-C8 by receptor-bound death-inducing signaling complexes (DISC) to form C8* (Figure 1, gray); (ii) an enzyme cascade in which C8* directly cleaves C3 [28] to form active C3*, which can cleave effector caspase substrates (a process represented in our model by cleavage of PARP to form cPARP [29]) but not when bound to XIAP (X-linked IAP [11] Figure 1, blue); (iii) a mitochondrial feed-forward pathway in which C8* cleaves Bid (into tBid) [30] to activate Bax (to Bax*) and promote formation of pores in the mitochondrial membrane through which CyC [10] and Smac [8,9] translocate into the cytosol following MOMP; cytosolic CyC then binds Apaf-1 and C9 to form the apoptosome (which also cleaves pro-C3 [31,32]), and Smac neutralizes XIAP [33–36], thereby de-inhibiting C3* (Figure 1, yellow); (iv) a positive feedback loop in which pro-caspase-6 (pro-C6) is cleaved by C3* to form C6* [37], which then activates additional pro-C8 (Figure 1, green) [38,39].
EARM v1.0 aims to be reasonably complete with respect to biochemical mechanism, but three simplifications were made to reduce the number of species and free parameters. First, the details of DISC and apoptosome assembly, both of which involve multiple copies of several protein species [40,41], were omitted in favor of simplified "lumped parameter" representations. Second, protein synthesis was omitted because all experiments were performed in the presence of cycloheximide (which is commonly used to sensitize cells to the action of TNF [42], but which, in our experiments, also simplifies modeling by eliminating source terms). Third, proteins with similar biochemical activities were represented by a single species: C8 and caspase-10 (C10) by C8 alone; C3 and C7 by C3 alone; and the Bcl-2-like family of proteins by three prototypical examples: Bid, a pro-apoptotic "activator," Bcl-2, an apoptosis inhibitor, and Bax, a pore-forming protein. We are aware of controversy regarding the precise mechanism by which Bcl-2-like proteins regulate MOMP and have implemented the simplest form of "direct activation" [43]. Further research will be required to determine whether kinetic modeling can help to distinguish this scheme from alternative "indirect activation" hypotheses [25,44].
Experimental Determination of C3 Activation Dynamics
To gather data for model training, we first sought to establish precise dose-response relationships for death ligands in HeLa cells, which are widely used in cell death studies and a robust biological setting in which to combine RNAi and single-cell imaging [45–48]. HeLa cells were exposed to TRAIL over a 500-fold range of concentrations spanning roughly physiological to saturating. Cell death was monitored by live-cell microscopy using either of two Förster resonance energy transfer (FRET)-based reporter proteins whose fluorescence changes upon cleavage (effector or initiator caspase reporter proteins, EC-RP and IC-RP) and a reporter for mitochondrial outer membrane permeabilization (mitochondrial intermembrane space reporter protein [IMS-RP]) whose cytosolic translocation mimics that of Smac and CyC [15].
Following TRAIL exposure, C3* activity remained low, as measured by EC-RP fluorescence, for a period of time that varied widely from cell to cell before rising rapidly to a plateau, at which point cells died (Figure 2A and 2B) [15]. Prehn and colleagues [12] have shown that these dynamics have the general form:
(A) Time-lapse images of HeLa cells expressing EC-RP treated with 50 ng/ml TRAIL. Imaging was performed at 465 nm and 535 nm; cleavage of EC-RP caused a shift in the 465/535 ratio that was pseudo-colored as a shift from blue to yellow. Inset numbers indicate the time in hours after treatment with TRAIL.
(B) Time courses of EC-RP cleavage in individual HeLa cells treated with 250, 50, or 10 ng/ml TRAIL; five cells are shown for each dose. Cell-to-cell variation in dynamics for a single dose is observed reproducibly and is not a consequence of imprecision in the assay.
(C) Idealized single-cell time course for EC-RP cleavage, showing relationships among Td (delay time), Ts (switch time), and f (fraction substrate cleaved).
(D, E) Frequency distributions for Td (D) and Ts (E) determined by live-cell microscopy in EC-RP-expressing HeLa cells treated with varying concentrations of death ligand (n > 100 for each condition); binning intervals were 30 min (E), or 5 min (F).
(F) Flow cytometry of PARP cleavage, as assayed with an antibody selective for PARP cleaved at Asp214 by C3 and C7, in HeLa cells treated with 10 to 250 ng/ml TRAIL for 1–9 h as indicated. Colored regions underlying the histogram represent intervals used to discretize data into PARP cleavage levels of <5% (green), 5%–25% (yellow), or >25% (gray; percentages are relative to the median fluorescence intensity of the fully positive population).
(G) Discretized flow cytometry data. The fraction of cells having low, intermediate, or high PARP cleavage is color-coded to match the intervals in (F).
(H) Immunoblot analysis of PARP cleavage in HeLa cells treated with 10 to 250 ng/ml TRAIL.
(I) Quantitation of the cleaved 89 kDa form of PARP for blots shown in (H).
where c(t) is the amount of substrate cleaved at time t, f is the fraction cleaved at the end of the reaction, Td is the delay period between TRAIL addition and half-maximal substrate cleavage (c(t)= 0.5 f), and Ts is the switching time between initial and complete effector substrate cleavage (the reciprocal of the slope at t = Td; Figure 2C). A fourth parameter, the time constant of the cleavage reaction Tc = Ts × f, provides a measure of switching time that is independent of the final amount of substrate cleaved. When live-cell data from ∼150 TRAIL-treated cells were parameterized using Equation 1, Td varied from 1 to 15 h depending on ligand dose (Figure 2D and Table 1). TNF treatment, even at saturating concentrations of ligand, elicited a response that was even more heterogeneous with respect to Td. In contrast, Ts and f varied little from cell to cell regardless of ligand dose or identity (with mean values ranging from 19 to 27 min for Ts; Figure 2E and Table 1). Thus, in accordance with previous studies [12], we find the extrinsic pathway of cell death to involve dose- and ligand-dependent variation in Td concomitant with maintenance of Ts at a constant value.
Values for Switching Parameters Determined by Live-Cell Microscopy
A caveat of using fluorescent reporter proteins to study apoptosis is that cleavage of synthetic substrates is not necessarily representative of endogenous substrate cleavage. Cleavage of the endogenous C3* substrate PARP was therefore monitored by flow cytometry of cells stained with antibodies specific for the cleaved form (cPARP) [49]. Over a range of TRAIL concentrations, a bimodal distribution in fluorescence intensities was observed in which antibody-nonreactive cells became less abundant with time, while antibody-reactive cells became more abundant (Figure 2F). Similar data were also obtained for the C3* substrate cytokeratin (Figure S1). To quantify these changes, fluorescence signals were discretized into levels corresponding to no detectable PARP cleavage (Figure 2G, green bars), maximal cleavage (grey bars), and cleavage of a subset of the PARP in each cell (yellow bars; as opposed to complete PARP cleavage in a subset of cells). Regardless of ligand dose or identity, <5% of cells scored as having partially cleaved PARP, as expected for cells that undergo a sudden transition to from life to death [50]. Thus, flow cytometry shows cleavage of endogenous C3* substrates to be rapid and complete at the single-cell level. To ascertain whether maximal fluorescence as measured by flow cytometry or FRET actually corresponds to complete cleavage, we monitored PARP levels and molecular weight by immunoblotting. Conversion of >80% of 112 kDa full-length PARP into 89 kDa cPARP was observed by 9 h at all doses of TRAIL (Figure 2H and 2I), and we therefore conclude that f ≈ 1.0 under our experimental conditions.
Merging Data from Single-Cell and Population-Based Measurements
To merge data obtained by imaging, flow cytometry, and immunoblotting and thereby quantify the mean and variance of Ts, Td, and f, we simulated the experimental procedures involved (Protocol S2). Synthetic live-cell data were generated for 10,000 idealized TRAIL-treated cells assuming snap-action, variable-delay caspase activation as per Equation 1 with f = 1.0 and normal distributions for Ts = 20 ± 10 min and Td = 180 ± 40 min (values were based on live-cell data for cells treated with 250 ng/ml; ranges represent standard deviations). 10,000 synthetic live-cell trajectories were then averaged at discrete points in time to generate synthetic immunoblot data (Figure 3A). Synthetic flow cytometry data were generated from single-cell trajectories by computing distributions of caspase substrate cleavage at fixed points in time. Empirically derived values for background fluorescence and measurement noise were added, the data were discretized into 1,024 bins (corresponding to 10-bit detection), and results were plotted on a log-linear scale (Figure 3A). Varying the mean values of Ts, Td, or f in the synthetic trajectories revealed that: (i) when Td alone increases, cPARP levels remain bimodal when monitored by flow cytometry (Figure 3B); (ii) when Ts is increased with f constant, bimodality is lost and cells having intermediate levels of cPARP accumulate at short but not long times (because f → 1.0 as t → ∞) (Figure 3C, yellow bars); (iii) when f decreases with or without changes in Td and Ts, many cells exhibit partial PARP cleavage even at long times (Figure 3D). Analogous inverse procedures make it possible to compare live-cell imaging, flow cytometry, and immunoblotting data quantitatively. In addition, they emphasize that the three methods differ in their ability to estimate mean values and distributions for Ts, Td, or f (Table 2). Flow cytometry, for example, yields indirect measures of the mean and variance of Ts, Td, and relative values for f at discrete points in time, while immunoblotting provides a good estimate of the average absolute value of f. Moreover, one measurement—suitably processed—usually adds information to another; determination of f by immunoblotting, for example, enables calibration of flow cytometry and live-cell microscopy values.
(A, B, C, D) Simulated live-cell (top, twenty individual cells shown as blue lines), immunoblot (top, red lines), raw flow cytometry (middle, computed from 10,000 simulated live-cell profiles), and discretized flow cytometry data (bottom) for four different sets of switching parameters: (A) Baseline, corresponding to 250 ng/ml TRAIL; (B) an increase in mean Td to 240 min. (C) An increase in Ts to 100 min (D) or a decrease in f to 0.1.
(E, F, G) Concordance of data on EC substrate cleavage as obtained from immunoblot, flow cytometry, and live-cell assays in HeLa cells treated with 250 ng/ml (E), 50 ng/ml (F), or 10 ng/ml (G) TRAIL. Values for immunoblots are derived from quantitating the cleaved 89 kDa form of PARP by immunoblot; flow cytometry values represent the percentage of cPARP-positive cells computed as described in the legend of Figure 2, and live cell data were computed as the cumulative percentages of dead cells for populations of >150 cells. To facilitate comparisons, all three signals were normalized to their final value.
Experimental Determinability of Switching Parameters
When experimental live cell, flow cytometry, and immunoblot data from HeLa cells exposed to TRAIL at three different doses were merged, excellent quantitative agreement in death dynamics was observed (Figure 3E–3G). This confirms that EC-RP and IC-RP are effective reporters of endogenous caspase substrate cleavage and that an increase in TRAIL dose from 2 ng/ml to 1,000 ng/ml causes Td to vary 4- to 5-fold even as Ts and f remain constant at ∼20 ± 10 min and 1.0, respectively. Quantitative models of apoptosis must account for the overall efficiency of this process, the dose-dependence and length of Td, the rapidity of Ts, and the independence of Ts and Td.
Linking Models and Experiment Via Perturbation
Parameters in EARM v1.0 were manually adjusted to minimize the difference between simulated trajectories and experimental data including: (i) mean values for Ts, Td, and f in cells exposed to a range of TRAIL concentrations (Figure 4A); (ii) composite live-cell time-courses of initiator and effector caspase activity and MOMP (Figure 4B); (iii) estimates of Ts, Td, and f obtained by flow cytometry of TRAIL-treated cells perturbed by small interfering RNA (siRNA) or protein overexpression (Figures 5 and S3). To obtain composite time courses, ∼100 live-cell measurements of individual cells expressing IMS-RP and either IC-RP or EC-RP were aligned by the time of MOMP. As we have recently reported [15], gradual cleavage of IC-RP but not EC-RP was apparent during the pre-MOMP delay. MOMP occurred once the IC-RP signal rose to ∼30%–50% of its maximum value and EC-RP was then cleaved to near completion within 20–30 min (Figure 4B).
(A) Simulation of Td (left) or Ts and f (right) as a function of TRAIL dose (lines) alongside corresponding experimental values (points with error bars indicating standard deviations). For predicted values of Td, an envelope of constant coefficient of variation is shown, as estimated from experimental data (CV ≈ 20%); the source of variation is not known.
(B) Composite plot of IC-RP and EC-RP cleavage for >50 cells treated with 50 ng/ml TRAIL in the presence of CHX and aligned by the average time of MOMP (left panel) and model-based simulation of the corresponding species (right panel). Data in the left panel were originally reported elsewhere [15].
First column: perturbation values measured by immunoblot (see Figure S2). Second column: comparison of EARM v1.0-simulated and experimentally derived values for Ts, Td, f, and Tc. Simulated values were computed by EARM v1.0 with the perturbation conditions indicated in the first column. Experimental values were derived by fitting the simFACS data model to discretized flow cytometry data (See Materials and Methods and Protocol S2 for details). Third column: comparison of EARM v1.0-simulated cPARP cleavage (red) to time courses derived from flow cytometry using a data model (blue). Blue curves were produced by Equation (1) parameterized with the experimental Ts, Td, and f values shown in the second column. The dashed line in the control condition shows simulation under non-siRNA conditions; all other simulations were performed under siRNA conditions (see "Linking Models and Experiment Via Perturbation" for details). Fourth column: comparison of experimental (blue) and predicted (red) flow cytometry plots at the indicated time points (see Figure S3 for comparisons of all time points). Predicted flow cytometry data were produced by simFACS data simulation using the EARM v1.0-simulated values for Ts, Td, and f shown in the second column. Raw data for several of the conditions shown were originally reported elsewhere [15].
Data on perturbed cells were obtained by siRNA-mediated protein depletion (typically, three independent siRNA oligos were validated for each gene; Figures 5 and S2) or by protein overproduction in cells stably expressing cDNAs under the control of viral promoters. Simulating a reduction in the abundance of a species 5- to 20-fold often had a different outcome from simply eliminating the species. Thus, RNAi was modeled by adjusting the initial value of a species ([X]0 for species X) to the measured extent of protein depletion as determined by semiquantitative western blotting . Uniformity of depletion was established using flow cytometry (Figure S2). A subtlety that arose when linking models to data was a ∼100–150 min reduction in mean Td following transfection of cells with nontargeting siRNA oligos. The reduction probably reflects siRNA-mediated induction of interferon-stimulated genes [51], which include the DR5 TRAIL receptor and C8 [52–54]. To account for this effect, cells transfected with targeting and nontargeting RNAi were modeled as having more TRAIL receptors (1 × 105 versus 2 × 102 per cell; Figure 5 and Table S5). Experimentally, targeting oligos were always compared in parallel to nontargeting controls.
Perturbation of the mitochondrial pathway (Figure 1, yellow circuit) by Bid depletion or Bcl-2 overexpression blocks cell death in response to TRAIL (Figure 5). To monitor the dynamics of the cytosolic reactions in the absence of the full mitochondrial feed-forward pathway, it is necessary to combine Bcl-2 overexpression with XIAP depletion [15]. Under these conditions, only a fraction of cells died in response to TRAIL, with Ts 3- to 5-fold longer than in unperturbed cells and C3 substrate cleavage incomplete (f = 0.15–0.3; Figure 5). The phenotype of Smac depletion was very similar, and both represent a highly undesirable state of "partial" cell death in which effector caspases only achieve sublethal levels [15]. These data demonstrate that efficient snap-action cleavage of effector caspase substrates absolutely requires the mitochondrial feed-forward pathway: the direct C8* → C3* cascade and C3* → C6* → C8* feedback loop are insufficient.
Reconciling dose-response, composite time-course, and perturbation data placed significant demands on EARM v1.0. Accommodating the conflicting requirements that XIAP be sufficiently abundant to fully block C3* in the pre-MOMP interval and that XIAP be efficiently sequestered by cytosolic Smac following MOMP involved careful adjustment of species concentrations and reaction rates, a fragility that implies the existence of additional as-yet unknown regulatory processes [15]. Nonetheless, simulating TRAIL treatment over a range of concentrations resulted in nearly constant Ts of ∼30 min, f = 1.0, and dose-dependent variation in Td, all of which are in good agreement with experimental data (Figure 4A). Simulated time-courses for Bid cleavage (representing cumulative C8* activity), Smac translocation (MOMP), and cPARP levels (representing cumulative C3* activity) also matched experimental observations closely (Figure 4B). Efforts to improve the goodness of fit between model and experiment, assess model identifiability, and quantify parametric uncertainty are ongoing but are complicated by cell-to-cell variability and the consequent necessity of finding best-fit distributions of initial conditions and rate parameters rather than single values. To address this challenge, new methods will be required, but recent singular perturbation analysis of a reduced version of the EARM v1.0 model nonetheless suggests that parameter values reported here (Tables S3 and S5) represent reasonable order-of-magnitude estimates (J.M. Burke and P.K. Sorger, unpublished data).
Transition from a Graded to a Snap-Action Signal
Differing dynamics of initiator caspase and effector caspase substrate cleavage raise the question of where in the extrinsic apoptosis pathway a steady and gradual increase in C8* activity is converted into a snap-action downstream signal. To address this issue, simulation was used to monitor each step in a "typical" cell (Figures 6A and S4), something that cannot be done experimentally. The binding of TRAIL to DR4/5 receptors results in gradual C8* accumulation, steady processing of Bid into tBid, and gradual accumulation of active Bax* (Figure 6B). In contrast, MOMP and the consequent release of Smac and CyC (which are identical in their release kinetics) from mitochondria via Bax*-containing pores is sudden and rapid, beginning ∼3 h post-TRAIL addition and reaching completion within ∼15 min.
For simplicity, positive feedback was omitted in all simulations by setting [C6]0 = 0; results with positive feedback produce highly similar conclusions and can be found in Figure S4.
(A) Overlay of multiple species. To accommodate the wide range of concentrations, some species are scaled according to the values in parentheses in the caption.
(B) The same data as in (A) except that each vertical axis is scaled independently to better depict the full dynamic range for each species. The pink vertical line denotes the duration of MOMP. The estimated LD50 for caspase activity (corresponding to ∼10% PARP cleavage) [15] is denoted by a dashed line in the cPARP plot; the gray shaded region denotes points in time after this lethal dose has been reached and cells are already destined to die.
(C) Simulated time courses for C8*, tBid, Bax*, and mitochondrial pores as in (A) shown on a truncated vertical axis to display the times at which each species achieves a concentration of one molecule/cell. The full time course for mitochondrial Smac/CyC is shown as a dashed line on a full vertical axis (right) to show MOMP.
(D) Top panel, simulation of the total number of MOMP pores as in (A) (green) in comparison to the subset of pores bound to Smac or CyC during the process of release (red). Bottom panel, simulation of released Smac/CyC as in (A) (pink) and the discrete-time derivative of the release reaction (red).
(E) Live-cell imaging of IMS-RP (red) and GFP-Bax (green) in the same cell. Frames correspond to 30-s intervals, with the first frame of MOMP denoted by a pink box. Partial mitochondrial localization of IMS-RP following MOMP is an artifact of the overexpressed reporter; endogenous CyC is fully cytosolic at this point, as determined by immunofluorescence (not shown).
(F) Quantitation of Bax pores (top) and rate of IMS-RP release (bottom) for two cells (orange curves correspond to the cell shown in [E]); puncta and rate of release were quantified as described in Materials and Methods.
(G) Live-cell images for the cell in (E) shown at higher zoom and with pseudocoloring to highlight the temporal and spatial relationship between Bax puncta and IMS-RP release. Frames correspond to 30-s intervals; the magnified region (yellow box) encompasses the area of the cell in which both GFP-Bax puncta and IMS-RP release are first visible.
Two complementary effects appear to account for the transition from graded C8* activation to snap-action Smac/CyC release: maintenance of a very tight "off" state during the pre-MOMP delay followed by creation of a low-impedance "on" state. When we examined the time at which each species reaches a threshold of one molecule/cell (an arbitrary level that is useful for illustration), C8*, tBid, and mitochondrial Bax* all exceeded this threshold within minutes of TRAIL addition, but Bax-containing pores (M*) remained below the threshold for hours until seconds before Smac/CyC translocation began (Figure 6C). Thus, pores are essentially absent during the long pre-MOMP delay. Simulation also shows that only about ∼100 pores (which require <15 min to form once assembly is initiated) are required to translocate >105 molecules of Smac/CyC into the cytosol, because movement down a very steep concentration gradient is involved (Figure 6D). (This corresponds to the low-impedance "on" state.) Thus, the presence of very few active pores is massively amplified during Smac/CyC translocation, resulting in rapid depletion of protein from the intermembrane space and a switch-like response. Simulation also reveals that pore formation can continue for up to ∼2 h (assuming the cell does not lyse first; Figure 6D) so that it vastly exceeds what is required for Smac translocation. Ongoing pore formation appears to be an example of extreme overshoot: to induce efficient early translocation, the number of pores continues to rise well beyond what is required for Smac/CyC transport. Indeed, increasing the concentration of Smac in the model from 105 to 109 molecules/cell (an unrealistically large number) confirms that rapid release depends on excess pore capacity relative to the pool of proteins to be released (Figure S5).
Predictions I and II: MOMP Is Complete by the Time Dying Cells Have Assembled Relatively Few Pores
Two testable predictions arise from the simulation of MOMP described above: (i) Smac release should begin nearly simultaneously with the formation of the first Bax-containing pores, and (ii) pore formation should continue long after Smac release is complete. To test these predictions, cells expressing green fluorescent protein (GFP)-Bax and IMS-RP were treated with TRAIL (50 ng/ml) and imaged at 60× resolution at 30-s intervals for 1 h before and after MOMP. The rate of release of IMS-RP was estimated using an edge-detection algorithm sensitive to the transition from clustered IMS-RP signals diagnostic of mitochondrial localization to diffuse signals diagnostic of cytosolic localization. In agreement with previous studies [55] and the assumptions in EARM v1.0, GFP-Bax had a diffuse cytosolic localization prior to MOMP but, in dying cells, formed bright puncta that colocalized with mitochondria (as marked by IMS-RP; Figure 6E). Moreover, GFP-Bax puncta appear to be identical to puncta detected by immunofluorescence microscopy (not shown), suggesting that GFP-Bax is representative of endogenous Bax. Appearance of the first GFP-Bax puncta coincided with IMS-RP translocation, which reached a maximal level within 1–2 frames (∼1 min), after which puncta continued to form for 20–30 min more, typically rising to >100 per cell by the time cells began to fragment (at which point further observation was unreliable; Figure 6F and 6G). Live-cell studies by others using similar methods have also demonstrated a close temporal link between initial formation of Bax puncta and MOMP (D.R. Green, personal correspondence). We often observed that the timing of MOMP varied with location in a cell, such that IMS-RP release from some mitochondria preceded release from other mitochondria by ∼1 frame (30 s; Figure 6G). In these cases, the first observable Bax punctum was associated with the earliest-releasing subset of mitochondria. A similar relationship was observed between aggregation of GFP-Bak (a second pore-forming protein involved in MOMP) and IMS-RP translocation, with the exception that GFP-Bak was found on the mitochondrial membrane prior to MOMP [56]. The first appearance of GFP-Bak puncta was coincident with IMS-RP translocation and puncta once again continued to form for >20 min thereafter (Figure S6). In contrast, GFP-Bcl-2 and GFP-Bcl-XL exhibited diffuse mitochondrial localization throughout, with no apparent changes during MOMP (unpublished data). Bid-GFP remained in the cytosol before, during, and after MOMP, suggesting that tBid dissociates from Bax and Bak prior to pore formation (some Bid-GFP aggregates were visible once membrane blebbing and cell shrinkage had begun but only long after MOMP was complete; unpublished data).
The precise relationship between functional translocation pores and visible GFP-Bak/Bax puncta is not known, and puncta visible by live-cell imaging certainly contain more than the 4–8 Bax or Bak subunits thought to comprise functional pores [56]. However, the close temporal and spatial association between puncta and IMS-RP release (particularly at different locations in a single cell) implies that puncta may be clusters of pores. Regardless, these experiments clearly confirm our two model based predictions (i) that cytosolic translocation of mitochondrial intermembrane proteins is complete by the time a relatively small number of pores have formed and (ii) that pores continue to form long afterward. Modeling also provides a possible explanation for this latter phenomenon: overshoot in pore forming reactions guarantees that the release of intermembrane proteins is sudden and complete regardless of variations in the rate of initial pore formation.
Prediction III: Snap-Action Switching Does Not Require Feedback
It is an open question whether feedback from processes downstream of Smac/CyC translocation is necessary to ensure rapid all-or-none induction of MOMP. The existence of feedback has been proposed previously [37,39] and is likely to account for the biphasic cleavage of IC-RP following exposure of cells to TRAIL (Figure 4B). Previous analysis shows that slow pre-MOMP cleavage is entirely C8* dependent [15], but the rapid post-MOMP phase probably involves one or more feedback processes including (i) cleavage by C3* itself, (ii) cleavage by C8*, whose levels are expected to rise rapidly as a consequence of the C3* → C6* → C8* feedback loop, and (iii) cleavage by C9, an initiator caspase activated by CyC translocation. In simulation, inhibition of these feedback loops eliminates the second, rapid phase of IC-RP cleavage but has little effect on the dynamics of MOMP (Figure 7A). Modeling predicts that one simple way to attenuate feedback experimentally is depletion of Smac, since failure to inactivate XIAP blocks all three feedback loops (Figure 7A). We observed Smac-targeting siRNA oligonucleotides to eliminate the rapid phase of IC-RP cleavage, confirming the prediction that feedback from downstream to upstream processes had been inhibited (Figure 7B). Nonetheless, the kinetics of IMS-RP release were unaltered over a wide range of TRAIL concentrations. Most strikingly, when Smac-depleted cells were exposed to very low concentrations of TRAIL (2 ng/ml), IC-RP was processed very slowly over a 10-h period, but IMS-RP translocation was as rapid as in unperturbed cells in which IC-RP cleavage was ∼5-fold faster (Figure 7B, bottom panel). Thus, in agreement with simulation, rapid induction of MOMP does not appear to require any of the feedback loops that impact initiator caspase substrate cleavage. These include feedback-mediated C8 activation [39] and cleavage of Bid by C3* [57]. The data do not exclude a role for feedback acting downstream of tBid (such as self-activation of Bax [58]), but on the basis of model analysis, we speculate that feedback—if it exists—is not actually necessary for MOMP to achieve its snap-action character in HeLa cells (Figure 7C).
(A) Simulated values for the average rate of C8 substrate cleavage during the early (pre-MOMP) and late (post-MOMP) phases (left axis) and Ts of MOMP (right axis). Simulations were performed under conditions corresponding to stimulation by 2 ng/ml TRAIL.
(B) Live-cell measurements of IC-RP cleavage and IMS-RP release in control or Smac-depleted cells treated with 10 or 2 ng/ml TRAIL. IC-RP and IMS-RP signals were normalized individually to allow comparison on the same axes.
(C) Schematic diagram of potential feedback pathways regulating MOMP. The data shown in (B) are inconsistent with snap-action behavior resulting from feedback loops that act upstream of Bid cleavage (red arrows) but do not rule out the possibility that snap-action may arise from feedback loops that act downstream of Bid (green arrows). Some of the depicted loops are hypothetical and are shown for logical completeness.
Prediction IV: Dose-Dependent Activation of C8* Controls the Duration of Td
Modeling suggests that the dose-dependence of Td is determined by the time required to saturate Bcl-2 with Bax* and thereby generate Bax* active in pore-formation. This, in turn, depends on the rate at which tBid is generated by C8* (Figure 8A) (and by other factors as well). Thus, we predict that Td is controlled by dose-dependent changes in C8* activity. Consistent with this hypothesis, we observed a ∼2-fold increase in the rate of IC-RP cleavage at 10 ng/ml TRAIL relative to 2 ng/ml and a further ∼1.5-fold increase at 250 ng/ml TRAIL (Figure 8B). To test directly the role of C8* in delay duration, we asked whether the short delay characteristic of cells exposed to 250 ng/ml TRAIL could be converted into a long delay by modulating C8* activity with a small-molecule inhibitor (Z-IETD-FMK; Figure 8C and 8D). At 10 μM, Z-IETD-FMK blocked C8 activation and MOMP in >95% of cells, as previously reported [59]. However, at 2 μM Z-IETD-FMK, C8* was only partially inhibited, and the rate of IC-RP processing in cells exposed to 250 ng/ml TRAIL was reduced to that of cells exposed to 50 ng/ml TRAIL. MOMP occurred under this condition, with Td increased from 2.1 to 3.3 h (Figure 8D), thereby showing that dose-dependent changes in C8* activity are sufficient to alter Td.
(A) Simulation of C8 substrate cleavage (tBid, top) and free mitochondrial Bcl-2 (bottom) in comparison to Smac/CyC release at 4 TRAIL concentrations (spanning approximately 0.2 to 200 ng/ml). Red circles indicate the level of tBid or Bcl-2 when Smac/CyC release is 50% complete.
(B) Live-cell measurement of IC-RP cleavage and IMS-RP release in individual cells stimulated with 250, 50, or 10 ng/ml TRAIL. IC-RP and IMS-RP signals were normalized individually to allow comparison on the same axes. Circles outlined in red indicate the IC-RP signal at the time of IMS-RP release.
(C) Averaged IC-RP cleavage for ten to 50 HeLa cells treated with 0 (black line), 50 (blue), or 250 (red) ng/ml TRAIL (as indicated) in the presence of cycloheximide. At each dose of TRAIL, time-courses were aligned by the average time of MOMP (indicated by symbols with red outlines).
(D) Average IC-RP cleavage in 50 Hela cells treated with 250 ng/ml TRAIL in the presence of 10 or 2 μM C8 inhibitor. Time courses for cells treated with 2 μM inhibitor were aligned at the average time of MOMP. Average time courses from (C) are shown for comparison in light blue, green, and gray. In (C) and (D), MOMP did not occur in cells not treated with TRAIL or treated with TRAIL in the presence of 10 μM inhibitor, and time courses were aligned by the time of TRAIL treatment
Prediction V: Ratiometric Control of Snap-Action Transmission
Two sets of reactions "transmit" the death signal generated by MOMP to C3*: (i) binding of cytosolic Smac to XIAP and (ii) assembly of the apoptosome. Smac binds and neutralizes XIAP, making it critical to C3* activation. The role of the apoptosome in extrinsic apoptosis is less clear, but simulation suggests that Apaf-1 and C9 should vary in importance depending on Smac levels. When Smac is present at high levels (5 × 105 molecules/cell), removal of either Apaf-1 or C9 has little effect on f, Ts, and Td (32 min versus 23 min; Figure 9A–9C). However, at lower levels of Smac, depletion of apoptosome components abolishes snap-action C3* activation (f = 0.5, Ts > 100 min; in all cases, MOMP was unaffected by changes in Apaf-1/C9 levels). We hypothesize that the key quantity controlling signal transmission from MOMP to C3* is the ratio of XIAP to the sum of the Smac molecules and apoptosome complexes in the cytosol (Figure 9D, dotted line). In support of this idea, depletion of Apaf-1 by siRNA resulted in a loss of snap-action switching that could be reversed by simultaneous depletion of XIAP and intensified by simultaneous depletion of Smac (Figure S7). Unfortunately, we were not able to confirm this result with multiple Apaf-1 oligos and our prediction therefore remains only weakly substantiated. Nonetheless, model analysis shows that, under conditions of limiting Smac levels, the phenotype associated with a C9* variant unable to bind XIAP is similar to that of C9 or Apaf-1 deletion (Figure 9E), whereas C9* lacking catalytic activity does not substantially alter the kinetics of PARP cleavage (Figure 9F). Thus, we predict that the XIAP-binding properties of the apoptosome are more critical in extrinsic cell death than C3* activation and that other XIAP-binding proteins such as ARTS [60] and Omi/HtrA2 [61] may also be important. A proper exploration of "MOMP transmission" will therefore require measuring and manipulating the levels of multiple overlapping XIAP binding proteins.
(A–C) Simulation of the time course of C3 (EC) substrate cleavage in control cells (A), or cells lacking C9 (B), or Apaf-1 (C). For each condition, simulations were performed using baseline conditions for Smac (black lines), or 5-fold overexpression of Smac (red lines).
(D) Simulation of f as a function of initial concentrations of Smac and the apoptosome; for simplicity, the apoptosome components C9 and Apaf-1 were assumed to be present at equal concentrations. Orange dotted line indicates the region where [XIAP]0 = [Smac]0 + [apoptosome]0; numbered circles and cyan arrows indicate the positions of the indicated conditions. Point 1a corresponds to the baseline model, in which either Smac or Apaf-1 knockdown leads to a decrease in f, while point 1b corresponds to the baseline model with 5-fold higher Smac, in which Apaf-1 depletion does not reduce f.
(E, F) Simulation of the time course of PARP cleavage in cells in which apoptosome function has been altered so that it cannot process pro-C3 (E) or so that it cannot bind XIAP (F).
Prediction VI: C6*-Mediated Feedback Controls Td but Not Ts
Depletion of C6 10-fold by RNAi did not significantly alter the dynamics of effector caspase substrate cleavage in HeLa cells treated with TNF or TRAIL (Figure 5, relative to control depletions; Td ≈ 80 min; Ts ≈ 15 min; f ≈ 1.0) and the lack of dependence on C6 was also captured by modeling (Figure 5). Remarkably, however, simulation showed C6* to play a major role in determining overall levels of C8* in dying cells. If the generation of C8* is separated into its two sources, DISC and C6*, it is apparent that a dramatic acceleration in C8* formation is induced upon C6* activation (Figure 10A). A priori, it might be assumed that this acceleration would constitute a key feature of snap-action switching. However, C8* generated by C6* lags temporally behind C8* generated by DISC and is present only after rapid cleavage of PARP has begun, making it inconsequential to the dynamics of TRAIL-mediated cell death under most circumstances (Figure 10B). As the dose of TRAIL is reduced, levels of DISC-generated C8* are also reduced, resulting in a larger relative contribution of C6* → C8 during the pre-MOMP period (Figure 10C and 10D). Thus, modeling predicts that C6 should have an impact on Td at very low TRAIL concentrations. In addition, in cases in which [C6]0, or the rate of C8 cleavage by C6* is high, the dynamics of C3* substrate cleavage should depend on the C6* feedback loop. This prediction appears to be supported by data showing C6 overexpression to enhance the sensitivity of cells to various apoptotic stimuli [62]. We conclude that positive feedback mediated by C6 (or by other topologically analogous processes) plays a minor role in snap-action C3* activation in Hela cells but may be more important under other conditions, including type-I apoptosis, in which MOMP is not required [63].
(A) Simulation showing the effect on Td of increasing the levels of [C6]0 10- (green line) or 100-fold (blue lines) above baseline values (black line) or of eliminating C6 altogether (red line).
(B–D) Simulations showing concentrations of total C8* (purple lines) and the subset of C8* generated by DISC (red lines) or by C6* (blue lines) at three different doses of TRAIL. The bottom panels magnify the period immediately before and after snap-action switching; bar plots to the right depict levels of DISC- or C6*-generated C8* (normalized to a total C8* value of 1.0) at the time of MOMP onset (dotted vertical lines).
Building a Threshold-Sensitive Snap-Action Switch
Modeling strongly suggests that both the biochemistry and physical partitioning of apoptotic regulators play roles in MOMP. To undertsand these roles, we generated a series of models of increasing complexity that describe reactions linking C8* activation to cytosolic translocation of Smac/CyC; as described above, these are the reactions shown experimentally to transform a graded initiator caspase input into switch-like effector caspase output. The input-output responses of variant models were analyzed over a range of C8* input concentrations for Td(MOMP), Ts(MOMP), and Ton(MOMP). The later metric represents the time at which 1% of Smac/CyC is translocated and is designed to capture the efficiency of repression in the pre-MOMP interval. In cells, Ts(MOMP) is short and invariant over a wide range of C8* levels (that is, "input-independence"), while Td(MOMP) and Ton(MOMP) should be virtually identical and increase in a dose-dependent manner (Figure 11A, inset).
Models representing varying topologies of the MOMP module were analyzed for snap-action behavior. For each model, the input species (C8*) was introduced at values ranging from 1 to 103 molecules/cell, and the resulting release behavior of Smac determined by simulation (middle column). The kinetic characteristics Ts, Td, and Ton (defined as the time at which Smac release reached 1% of its final value) for Smac release were determined from these simulations as a function of input strength (right column). Note that f = 1 in all cases and therefore Tc = Ts. Binding of Bid to Bcl-2 was omitted for simplicity, although this interaction is included in the full EARM v1.0.
(A) Basic motif model. Simulations are shown for parameter sets representing physiologically realistic rate constants (gray) and for irreversible, faster-than-diffusion binding (orange) for Bax-Bcl-2 association, and Ts, Td, and Ton curves shown only for physiologically realistic values. The inset plot in the right column shows the expected behavior for an idealized variable-delay snap-action switch.
(B) Model of known topology, including Bid cleavage.
(C) Model including Bax oligomerization.
(D) Model including separate mitochondrial reaction compartment.
(E) Model corresponding to the MOMP module in EARM v1.0. Note that the cytosolic pool of Bcl-2 anti-apoptotic proteins, which binds to tBid in EARM v1.0, has been omitted for simplicity.
(F) Model as in (F), but with rate constants for Bax oligomerization and insertion adjusted to represent positive cooperativity.
The simplest topology for a threshold-sensitive switch, and one implicit in the oft-cited concept of a Bax-Bcl-2 "rheostat" [64], involves an active species (Bax*) that is generated by an input (C8*) and then antagonized by binding to an inhibitor (Bcl-2). As Bax* increases in concentration, the pool of Bcl-2 is exhausted until [Bax*] > [Bcl-2] and the output (Smac*, representing cytosolic transport of Smac through Bax-containing pores) is induced. This system functions as a threshold-sensitive switch if Bax and Bcl-2 are assumed to associate at a rate faster than diffusion and to bind irreversibly (Figure 11A, orange curves). However, when realistic biochemical constants are used (e.g., a diffusion-limiting on-rate of ∼106 M−1s−1 and Kds in the measured range of 1–10 nM, [7]), the circuit is ineffective as a switch (Figure 11A): Ts(MOMP) is highly input-dependent while Ton(MOMP) is always short. Addition of Bid-dependent Bax activation (C8* → tBid → Bax*; Figure 11B) and a further requirement for Bax multimerization during pore formation (Bax* → Bax2* → Bax4* → Smac*; Figure 11C) result in successive improvements in performance as judged by input-independence of Ts(MOMP) and input-dependence of Ton(MOMP). Further increases in the input-independence of Ts(MOMP) are achieved by introducing a mitochondrial membrane compartment in which Bax and Bcl-2 interact (Figure 11D). With Bax and Bcl-2 confined to a membrane compartment that is ∼7% of the cytosolic volume, protein concentrations and rates of association rise 14-fold. As a consequence, the Bax*-Bcl-2 association/dissociation reaction remains close to equilibrium even as free Bcl-2 falls to very low levels. A final improvement in performance is achieved if Bax4* is required to undergo a conformational change and membrane insertion reaction during pore formation (Figure 11E); motivation for this step comes from studies showing that Bax-mediated pores form only at specific sites in the mitochondrial outer membrane (an insertion step is likely to be even more importantfor Bak, which localizes to mitochondria even in normally growing cells) [65,66]. Overall, the steady improvement in performance observed with successive models is primarily attributable to imposition of tighter repression on pore formation during the delay phase (Figure S8); as noted above, very tight negative regulation of MOMP is critical to sustaining long Td.
Model E, which corresponds to the topology in EARM v1.0, generates Td(MOMP), Ts(MOMP), and Ton(MOMP) that are largely consistent with the data, although Ts(MOMP) remains somewhat input-dependent at low TRAIL concentrations in the model but not in experiments. In model F, further input-independence in Ts(MOMP) can be achieved by inclusion of positive cooperativity in Bax oligomerization (this is modeled by making the rate constants of Bax dimerization, Bax tetramerization, and Bax insertion successively more rapid; Figure 11F). The impact of such cooperativity on the full ODE model remains to be determined, as does a demonstration that it exists in vivo. Nonetheless, it is noteworthy that robust all-or-none switching can be achieved in our simulations by a network that does not include feedback. It has been suggested that binding of Bcl-2 to both Bid and Bax multimers creates an "implicit positive feedback loop" [22,24], because the same pool of Bcl-2 inhibits sequential steps in the tBid → Bax* → Bax*2 → Bax*4 reaction series. However, implicit feedback does not appear to be essential for rapid switching, since model performance was not substantially impaired when separate, non-communicating, pools of Bcl-2 were implemented for all Bax-binding reactions (Figure S8). Instead, our findings suggest that snap-action control over MOMP arises from interplay between the biochemistry of protein-protein interaction and translocation between physical compartments of different volumes.
Variable-delay, snap-action regulation of effector caspases has been observed in a number of studies of extrinsic apoptosis, but the mechanisms responsible for this behavior have not, to our knowledge, been identified. Here, we develop and test an ODE-based model of pathways linking TRAIL-receptor binding to an extended pre-MOMP delay of variable duration followed by rapid snap-action cleavage of effector caspase substrates. Modeling cellular biochemistry inevitably involves a tradeoff between tractability, implying fewer species, and detail or scope, implying more species. Previously published models of caspase regulation by death receptors have focused on subsets of the reactions in EARM v1.0 including the C3/C6/C8 feedback loop [17], XIAP [22,23], C3 degradation [19], the apoptosome [20,67], or mitochondrial permeabilization [20,24–26]. However, our work suggests that models spanning reactions upstream and downstream of MOMP yield insight that cannot be obtained from models having less scope. In light of the the trade-off between scope and detail, it is interesting to note that significant mechanistic detail (e.g., multimerization, physical compartmentalization, etc.) was necessary to reproduce the features of TRAIL-induced MOMP identified experimentally. We therefore expect further insight as we add to EARM v1.0 details of DISC assembly [68], multiple initiator and effector caspases [3], and the 20 or so members of the Bcl-2-protein family [6]. Expanded models will be harder to manipulate and train but will make it possible to distinguish among alternate proposals for the reactions driving MOMP during intrinsic apoptosis and to explore the mechanisms of action of Bcl-2-binding drugs such as ABT-737 (see the EARM documentation page at http://www.cdpcenter.org for updates) [69].
Several hypotheses regarding the origins of snap-action control of apoptosis have emerged from previous studies, but relatively few have been explored experimentally at the single-cell level. Our primary guide has been single-cell data collected from normally growing and RNAi-depleted cells rather than mathematical analysis of EARM v1.0 equations (for which, after all, multiple formulations are possible). The importance of integrating modeling and experimentation is often emphasized [70] but is difficult to achieve in practice. One challenge is obtaining sufficiently accurate and quantitative experimental data, particularly for a process such as apoptosis that varies significantly from cell to cell. Here we use data fusion and simulation to integrate measurements made by live-cell imaging, flow cytometry, and immunoblotting with the result that mean values and variances could be obtained for key descriptors of C3* dynamics in normally growing and perturbed cells across a range of ligand concentrations. The use of synthetic data as a means to troubleshoot and analyze experiments is common in the measurement sciences, and the apparent effectiveness of this approach here and elsewhere [71] hints at its potential value in biology. A more fundamental challenge in linking simulations to experimental data is the difficulty of estimating unknown parameters. Given available data, EARM v1.0 is nonidentifiable and in the current work we rely on order-of-magnitude estimates. Singular perturbation analysis suggests that many model-based predictions are fairly robust to parametric uncertainty (J.M. Burke and P.K. Sorger, unpublished data), but further work is clearly required on this topic.
Determinants of Snap-Action Activation
We find that MOMP is the point in the extrinsic apoptosis network at which a graded TRAIL → → C8* → → Bax* signal is transformed into an all-or-none snap-action signal. This manifests itself as sudden formation of GFP-Bax and GFP-Bak puncta, which are thought to correspond to clusters of pores [56], concomitant with rapid CyC and Smac translocation [47,48]. In contrast, cleavage of an IC-RP reporter (a proxy for the production of the Bax activator tBid) is gradual over many hours prior to MOMP. Our conclusion that snap-action behavior originates at the level of MOMP is consistent with previous single-cell [47,48,72] and biochemical [73,74] studies.
Reactions linking Smac and CyC translocation to effector caspase activation are also essential for execution of an all-or-none death decision. Data and modeling agree that depletion of Smac or overexpression of XIAP [23] creates a mismatch between XIAP levels and those of its binding partners, thereby preventing XIAP from being fully sequestered from C3* by proteins released during MOMP. Reactions involved in apoptosome assembly can potentially generate bistable or ultrasensitive behavior [20,22], but this capacity does not appear to be essential for snap-action C3 activation in our simulations: any such bistability manifests itself after the live-die decision is made.
To better understand how snap-action switching emerges from the reactions controlling MOMP, we have analyzed a series of models in which processes linking activation of C8* and release of Smac/CyC are represented with increasing complexity and realism. The models were scored for their ability to maintain a durable "off" state, a delay whose duration was determined by input strength (corresponding to TRAIL dose), and input-independence in the rate of transition from "off" to "on." While a number of studies have characterized biochemical switches from the perspective of steady-state behavior and Hill functions [75,76], our metrics attempt to capture the temporal dynamics of switching, a topic that has received less attention [77]. An important but relatively obvious finding is that simple competition between an activator (Bax*) and an inhibitor (Bcl-2) makes for a poor switch under any set of realistic biophysical assumptions. Switch-like behavior relies on additional components, physical partitioning of reactants into compartments, and multimerization of Bax into functional complexes. One unusual aspect of MOMP is that its regulators are found in both the cytosol, where they diffuse in three dimensions, and the mitochdonrial membrane, where they diffuse in two dimensions. This transition could potentially restrict reaction volumes and further facilitate formation of the Bcl-2:Bax complexes critical for the inhibition of pore formation during the delay phase. When Bax- and Bak-containing pores finally form, rapid release of Smac/CyC is ensured by the favorable kinetics of moving concentrated protein stores down a steep concentration gradient (from the mitochondrial intermembrane space to the cytosol). Overshoot in the number of pores and rapid translocation kinetics do not guarantee complete input-independence of Ts(MOMP) in our models but do produce values in reasonable agreement with experimental data (5–10 min). However, at very low input levels, Ts(MOMP) in the model exceeds the range of values observed experimentally (increasing to ∼40 min). One potential way to overcome this is cooperativity in the Bax1* → Baxn* assembly reactions, but other mechanisms involving additional Bcl-2-like proteins also exist and require further exploration. Moreover, the remarkable overshoot in the number of pores formed relative to the number needed to translocate Smac/CyC raises the question of whether pores have another function or whether overshoot is an inevitable byproduct of rapid pore assembly. Finally, while the models described in this paper are effective at capturing the kinetics of apoptosis over relatively short time periods (6–24 h), protein synthesis and survival signaling must be added to capture long-term outcomes. Analysis of more complete models of MOMP, coupled with live-cell imaging of additional reactants, should help to further refine our understanding of how cells reconcile the competing demands of a long pre-MOMP delay and fast dose-independent switching.
Feedback and Network Bistability
Biological switches involving changes in cell fate are often analyzed in terms of positive feedback and the potential for bistability [17,20,24]. However, we find that at least one large class of positive feedback mechanisms can be ruled out by our data: snap-action behavior at the level of MOMP occurs independently of caspase-dependent feedback. This finding, based on simultaneous monitoring of initiator caspases and MOMP in individual cells, is in agreement with data from others who have found that the kinetics of CyC release are unaffected by pan-caspase inhibition [48,72]. Moreover, we find no effect of depleting C6, a putative mediator of feedback from C3* to C8. These results do not rule out a role for local feedback in extrinsic cell death (involving, for example, the proposed conversion of Bax to Bax* by already active Bax* [24,58]). However, we find that even positive feedback loops that improve the dose-independence of Ts do so at the cost of reducing agreement between simulation and other experimental data, such as the dose-dependence of Td and the kinetics of Bax cluster formation (unpublished data). Since putative switching mechanisms must be assessed in the context of all available data rather than their impact on Ts alone, our preliminary conclusion is that positive feedback in the induction of MOMP cannot be justified by existing data or simulations.
Bistability is an attractive concept in a network that decides between two fates, but does not appear to be important in our model. Protein synthesis and caspase-mediated feedback, necessary for bistability in many published models, are not required for snap-action behavior in our system. Moreover, even at the lowest TRAIL concentrations (in the presence of cycloheximide), all cells eventually die, implying the absence of a bifurcation in cell fate (this is not true for TNF, probably because of additional complexities in DISC formation). In agreement with this finding, singular perturbation analysis of EARM reveals the existence of a monostable transcritical bifurcation in the Bax-Bcl-2 binding reactions controlling MOMP rather than the bistable cusp bifurcation observed in pathways with strong feedback (J.M. Burke and P.K. Sorger, unpublished data). In rationalizing the absence of bistability, it is important to note that cell fate cannot be equated with the steady-state level of C3; under many conditions, transient effector caspase activation seals the fate of a cell long before steady-state is reached, with the result that a quasi-steady state having low C3* does not necessarily correspond to sustained survival [78]. It is clear however, that caspase networks play a role in cellular processes other than apoptosis [79,80], and it is possible that bistability arises in these networks in other physiological settings.
Determinants of Variable Delay
What features of the pathways regulating extrinsic apoptosis control the pre-MOMP delay? One very stringent requirement is that no pores form during the delay phase since experiments and simulation suggest that a very small numbers of pores can kill a cell. During the delay, the C8* substrate tBid generates Bax*, which binds immediately to inhibitory Bcl-2 (in EARM v1.0 this takes place in the mitochondrial membrane compartment); when the active forms of Bax* exceed the number of Bcl-2, a threshold is reached, pores form, and cells proceed rapidly to death. We show by simulation and experiment that varying C8* activity is sufficient to explain the dose-dependence of Td. Model analysis indicates that the delay corresponds to the time necessary to generate enough Bax* to exhaust the pool of mitochondrial Bcl-2, in agreement with the long-standing notion that MOMP is controlled by the Bax:Bcl-2 ratio [81–83]. Although the available pool of mitochondrial Bcl-2 ultimately sets the threshold for pore formation, analysis predicts that changes in the levels of any protein upstream of Bax* can also affect delay. Increases in the levels of protective proteins such as FLIP, Bar, or cytosolic Bcl-2 should lengthen Td, while increases in pro-apoptotic proteins such as TRAIL receptor, C8, Bid, and Bax should shorten Td. Because caspase-mediated feedback can regulate delay at very low TRAIL concentrations (even if it does not play a role in snap-action switching), proteins downstream of MOMP such as C6, C3, and XIAP also have an impact on Td under some circumstances. Td also increases overall in cells that have not been treated with cycloheximide, largely as a consequence of a rise in the threshold at which MOMP occurs [15].
An oft-stated assumption about cellular networks is that they are robust to changes in protein concentrations and other parameter values [84]. Our results suggest that while Ts and f are robust to ligand dose and identity, Td is sensitive to the concentrations of many reactants in the extrinsic cell death pathway and to external factors. The sensitivity of Td can easily be rationalized as being important for a network that must adjust its responsiveness to different extra-cellular cues, depending on physiological context. Conversely, the robustness of Ts to natural variation in internal and external parameters is likely to be important in preventing accumulation of "half-dead" states. Ascertaining precisely how tumor cells exploit sensitive and robust features of cell death networks so as to escape normal control over proliferation promises to provide a new window into the role of apoptosis in oncogenic transformation.
Physiological Importance of Variable-Delay Snap-Action Switching
At least three adaptive advantages might accrue to a cell in which apoptosis is under the control of a variable-delay, snap-action switch. First, the rapidity of switching should prevent cells from initiating but failing to complete the cell death program. Otherwise, a "half-dead" state might result, running the risk that chromosomal rearrangements and genomic instability typical of leukemia and other cancers might occur [14,85]. Second, dose-dependent variation in Td should make it possible to regulate the strength of a death signal independent of an all-or-none response at the single-cell level. Longer delays prior to caspase activation are likely to correspond to lower overall sensitivity to cell death (under physiological conditions), since caspase activation pathways do not operate in isolation; induction of intracellular survival pathways can override TRAIL-induced death signals. Indeed, the number of cells committing apoptosis in response to death ligands can be modulated by treatment with EGF or insulin, which stimulate pro-survival pathways [86]. Because the duration of the pre-MOMP delay in TRAIL-treated cells is similar to the time required for induction of new protein synthesis, a race is run between pro- and anti-apoptotic processes, and longer delays in the pro-death cascade favor survival. Third, cell-to-cell variability in Td within a population, a phenomenon observed even in the absence of new protein synthesis, prevents cells from dying en masse following exposure to a death stimulus. Recent work from our lab suggests that naturally occurring cell-to-cell variation in the concentrations of apoptotic regulators is responsible, in large part, for variability in Td (S.L. Spencer and P.K. Sorger, unpublished data). We hypothesize that variability in Td makes it possible for an all-or-none response at the single-cell level to be graded at the population level, such that the concentration of ligand controls the fraction of cells responding to a stimulus.
Materials and Methods
Reagents and cell culture.
HeLa cells stably expressing pMIG-Bcl-2 or pMIG (empty vector control) were gifts of Fei Hua and Michael Cardone; HeLa cells stably expressing EC-RP, IC-RP, and IMS-RP have been described elsewhere [15]. Plasmids containing the coding sequences of Bax and Bak were obtained from Richard Youle and the Harvard Institute of Proteomics. GFP-Bax and GFP-Bak were derived by PCR amplification of Bax and Bak coding sequences and ligation in frame with EGFP in pExchange-1 (Stratagene). HeLa cells stably expressing IMS-RP in combination with GFP-Bax and GFP-Bak were derived by standard transfection and selection techniques. Multiple stable clones were examined for each construct.
Antibodies were obtained as follows: cPARP, cleaved C3, XIAP, Bid, and Apaf-1 (BD Biosciences); C6, C8, and PARP (Cell Signaling); Smac (Calbiochem); Bax (Santa Cruz); cleaved cytokeratin (Roche). Secondary antibodies for flow cytometry and immunoblotting were from Molecular Probes (Invitrogen) and Rockland Immunochemicals. SuperKiller TRAIL was from Alexis Biochemicals, and cycloheximide (CHX) was from Sigma. siRNA oligos were purchased from Dharmacon; sequences are listed in Table S7.
Live cell microscopy.
HeLa cells stably expressing fluorescent reporters were grown in chambered coverglass slides (LabTek) for 24–48 h. Immediately prior to imaging, medium was replaced with phenol red-free CO2-independent medium (Gibco) containing the appropriate concentrations of SuperKiller TRAIL and CHX. Cells were then imaged at 37 °C on a Deltavision Spectris AX71 microscope equipped with a temperature control chamber using 10×, 20×, or 60× objectives. CFP, YFP, GFP, or RFP images were acquired at 3-min intervals for 12–18 h. Cleavage of the reporters, which results in a separation of the CFP and YFP proteins and loss of intramolecular FRET, was monitored by the change in CFP:YFP ratio. The changes in CFP:YFP ratio over time for individual cells were measured using ImageJ software with a custom-made plugin, which is available from the authors, and Td and Ts were determined by fitting cells to Equation (1) in Matlab (for Ts, time courses were fit using only data prior to cell shrinkage, since data after this point were unreliable). Mitochondrial localization of IMS-RP, GFP-Bax, and GFP-Bak was assessed by quantitating the total intensity for individual cells after images were processed using the "Find Edges" algorithm in ImageJ. Due to photobleaching, analysis of the IMS-RP signal in each cell was limited to the period beginning 1–2 h prior to MOMP and ending with the appearance of apoptotic morphology; no significant changes in IMS-RP localization were visible outside of this period. Rate of release of IMS-RP was computed as the difference in the edge-detected signal for successive frames. For alignment of multiple cells, Td(MOMP) for each single cell time course was set to the population average value for Td(MOMP), and the average FRET cleavage ratio was then calculated for each time point.
siRNA transfection.
siRNA transfections were performed as described [87], with a final concentration of 3–5 nM in the growth medium for each oligo. Medium was changed one day after transfection, and TRAIL treatment was performed on the second day after transfection. Oligo sequences are shown in Table S7.
Immunoblots were carried out as described [88], detected using an Odyssey Scanner (Li-Cor Instruments), and quantitated using ImageJ.
Inclusion of cycloheximide in death assays.
All TRAIL or TNF treatments were carried out with 2.5 μg/ml CHX present to inhibit protein synthesis. Cells treated with 2.5 μg/ml CHX alone showed no evidence of PARP cleavage or cell death for at least 24 h.
Flow cytometry assays.
Fixation, staining, and analysis procedures for flow cytometry were performed as previously described [15]. Antibody concentrations were 1:250 for anti-cleaved PARP and 1:500 for anti-cleaved C3 and anti-cleaved cytokeratin. Analysis was performed using the FlowJo software package (Treestar).
Simulation of flow cytometry data.
Simulated live-cell, flow cytometry, and immunoblot data were produced using the Matlab program simFACS, which is included as supplementary Protocol S2. Estimation of Td, Ts, and f was performed by manually adjusting these parameters to fit experimentally observed frequencies of negative, intermediate, and positive cells.
Mathematical model.
Biochemical reaction equations were derived from the canonical death receptor network (see [89] for review). For each given binary reaction i the biochemical equation is represented by one of following general mass-action paradigms:
where an enzyme or other protein E reacts with its substrate or binding partner S forming complex E:S, and depending on the specific reaction, i forming product P such that ki, k−i, and κi are the forward, backward, and catalytic rates, respectively. Tables S1 and S2 list the species that are used in the biochemical equations.
Cytosolic and mitochondrial compartments (see below) are each assumed to be well mixed. Let x = [X] represent the number of molecules of species X in each compartment. Let the rate at which the number of molecules per compartment of species X changes with respect to time t, in seconds, be represented by x = dx/dt. Using the rates found in Table S3, and the state variables and initial conditions in Table S5, the biochemical reactions in Table S2 are modeled by the system of ODEs found in Table S6.
Simulation and visualization.
Matlab stiff solver ODE15s was used to numerically solve the system of ODEs; the full model and Matlab routines for producing all simulations are contained in Protocol S1. Graphs and visualizations were made using Matlab and Axum software packages. Experimental overexpression and knockdown perturbations are simulated in the model by adjusting the baseline initial conditions appropriately in accordance with immunoblot data (Figure S2).
Cytosolic and mitochondrial compartments.
Since many of the biochemical processes of the mitochondrial pathway occur on or in the outer membrane of the mitochondria, these processes are effectively physically separated from processes occurring in the cytosol. To reflect this segregation, the model is comprised of two compartments, cytosolic and mitochondrial. Each compartment is a separate, well-mixed volume; specific molecules, such as CyC, may move from one compartment to another at a defined rate. Based on measurements of the mitochondrial and cytosolic volume of HeLa cells [90], the volume of the mitochondrial compartment is estimated to be ∼7% of the volume of the cytosolic compartment. The rates of second-order (or higher-order) reactions, in which two molecules must collide, increase as volume decreases, due to the increased likelihood of a collision. To reflect this, the reaction rates of intermolecular reactions occurring in the mitochondrial compartment are divided by the ratio of the mitochondrial compartment volume to the cytosolic compartment volume (v = 0.07).
Rate constants.
We assume a cellular volume of 10−12 liters (1 pL), and diffusion-limited bimolecular reaction rate constants ki. Therefore the forward reactions are bound by ki ≤ 4.2 × 106 M−1 s−1 and 4.2 × 104 M−1 s−1 in the cytosolic and mitochondrial compartments, respectively, for all reactions i [91,92]. All dissociation rates are approximately 10−3 s−1, resulting in dissociation constants in the range of 1–10 nM in the cytosolic compartment (with one exception: for activation of receptor the dissociation rate is 10−5 s−1, since this is a lumped parameter fit to match data). All catalytic rates are 0, 1, 0.1, or 10 s−1. All translocation reactions (between cytosol and mitochondrial compartments) have a rate of 10−2 s−1.
Initial protein concentrations.
Relatively few references for initial concentrations of proteins are available (Table S5; note also that values used in the model do not always correspond to those in the provided references due to fitting). However, many signaling proteins have been found to be present in the range of 1nM to 1μM (∼103–106 copies per cell) [93]. Therefore, we restricted initial protein concentrations to fall within these bounds (with the exception of c-Flip, which does not figure significantly in the model and which was estimated to be present at lower levels).
Protein synthesis.
To reflect the use of cycloheximide in all experiments, the mathematical model contains no source terms.
Protein degradation.
Protein degradation is assumed to be slow compared to signaling dynamics, so with one exception, the model contains no sink terms. The one exception is the case of C3* degradation by XIAP, which was experimentally determined to occur at significant levels on the timescale of the experiment [15]. In this case, degradation is approximated by assuming that ubiquitinated C3 is catalytically inactive.
MOMP variant models.
Models were constructed corresponding to the diagrams shown in Figure 11. Rate constants and initial conditions for Bid, Bax, Bcl-2, and mitochondrial binding sites were identical to the corresponding parameters in EARM v1.0 unless otherwise indicated. For simplicity, MOMP was represented by the release of Smac alone, with an initial concentration of 106 molecules/cell.
Supporting Information
Figure S1. Observation of Snap-Action Behavior in Multiple Endogenous Effector Caspase Substrates
HeLa cells were treated with 50 ng/ml TRAIL + 2.5 μg/ml CHX for the indicated times, and then analyzed by flow cytometry for anti-cleaved PARP or anti-cleaved cytokeratin antibodies. The data were first gated to exclude debris with low FSC signal (top panels, blue boxes) and then plotted as histograms (lower panels). In the case of both cleaved PARP and cleaved cytokeratin, >95% of cells display either positive or negative staining at all times, with <5% of cells falling into the intermediate range. (Note however, that the position and width of the positive population is different for the two antigens, probably as a result of differences in antibody affinity and subcellular distribution of the antigens.)
Figure S2. Measurement of Perturbation Magnitudes
(A) Quantitation of protein depletion and overexpression by immunoblot. For each siRNA oligo used, HeLa cells were transfected as described in the Methods section. After ∼48 h, cells were lysed in sample buffer and separated by SDS-PAGE. Blots were then probed with an antibody against the targeted protein (top image in each panel) and an antibody against a nontargeted protein as a loading control (bottom image in each panel). The level of depletion induced by each oligo (relative to nontargeting control oligo, and corrected for loading where possible) was quantitated using ImageJ software and is shown beneath each lane.
(B,C) Assessment of knockdown homogeneity by flow cytometry. HeLa cells were transfected with XIAP siRNA as in (A) and then trypsinized, fixed, and analyzed for XIAP levels by flow cytometry. Data are shown either as histograms (with data collected in linear amplification mode [B]), or scatter plots (with data collected in logarithmic amplification mode [C]; blue region denotes the untransfected population). Analysis of the fraction of untransfected cells (percentages listed above gated regions) indicated that knockdown occurred in 90%–95% of cells under typical conditions.
Figure S3. Agreement between Simulated and Experimental Flow Cytometry Data
Complete timecourses for each of the perturbations shown in Figure 5; simulated distributions were computed using simulated values for Ts, Td, and f (from EARM v1.0) as input into the flow cytometry data simulation (see Protocol S2). Although the overall qualitative agreement between simulation and experiment is high, several points of discrepancy can be observed. First, the cell-to-cell variability parameters assumed by the flow cytometry simulation are not appropriate under all perturbed conditions. For example, under Smac depletion conditions, the assumed variance in Td is too great at early times, causing a shift in the unstimulated (0 h) population that is not observed experimentally, while the assumed variance in f at late times is too small, leading to an equilibrium population (4.5 h) that is narrower than observed experimentally. Furthermore, while measurement noise is assumed to be the same for all cells in the flow cytometry simulation, experimental data suggest that there is considerably more variation in the measured values for cPARP-positive cells than for cPARP-negative cells (compare negative and positive peaks for control cells), which is likely explained by the increased morphological heterogeneity of cells that are being actively degraded in the later stages of apoptosis. Second, the shift in cPARP fluorescence in MOMP-blocked cells (Bid-depleted or Bcl-2-overexpressing) is larger in EARM v1.0 than seen in experiment; this shift is due to the difficulty of maintaining an effective C3* "off" state using reported values for cellular XIAP and C3 concentrations, which we have discussed in detail elsewhere [15]. Third, in Bcl-2-overexpressing, XIAP-depleted cells, model and experiment display minor differences in f (0.32 versus 0.15 for 87% XIAP depletion, and 0.57 versus 0.3 for 94% depletion) that appear exaggerated on a log scale, despite the fact that the model accurately captures the key qualitative behavior of the system, switching with f < 1. Thus, the model is in reasonable quantitative agreement with experimental observations, but additional calibration of both EARM v1.0 and the data simulation model may eventually yield closer matches.
Figure S4. Multiple-Species Simulation in the Presence of Feedback
The simulations shown in Figure 6A were also run in the presence of C6-mediated feedback. (Feedback was omitted in Figure 6A to make clear that the observed kinetics were feedback-independent). Very similar kinetics are observed whether C6-mediated feedback is present or absent. The largest differences are seen in the kinetics of C8* and tBid, which show a second, rapid phase of activation following MOMP in the presence of feedback.
Figure S5. Dependence of Smac Release Kinetics on the Size of the Mitochondrial Smac Pool
Simulations were performed under identical conditions to those in Figure 6D, except that the initial concentration of Smac in the mitochondria was increased to 109 molecules/cell and CyC was set to 0 (so that the total pool of molecules to be released was 109). Under these conditions, Smac release takes much longer (∼2 h, compare to Figure 6D), and many more pores are involved in the release event. This result confirms that the rapid, switch-like kinetics of Smac release seen with baseline Smac concentration (Figure 6D) require that the capacity of mitochondrial pores be in large excess relative to the size of the Smac pool to be released.
Figure S6. Formation of GFP-Bak Puncta In Synchrony with MOMP
HeLa cells expressing GFP-Bak and IMS-RP were treated with TRAIL and imaged at 30-s intervals with a 60× objective. The localization of GFP-Bak and IMS-RP are shown as raw images and as processed images after the application of an edge-detection filter. The first frame in which IMS-RP is detected is highlighted by a yellow box; in the same frame, the first GFP-Bak puncta appear (arrow heads). Thus, as expected for functionally redundant proteins involved in pore formation, GFP-Bax and GFP-Bak both form oligomeric clusters simultaneous with the initiation of MOMP and continue to form clusters long after MOMP has been completed.
Figure S7. Effects of Apaf-1 Depletion
(A–C) HeLa cells were transfected with siRNA oligos targeting Apaf-1 (A), Apaf-1 and XIAP (B), or Apaf-1 and Smac (C), treated with 50 ng/ml TRAIL, and assayed for cPARP by flow cytometry.
(D) Quantitation of Apaf-1 knockdown by immunoblot. For one Apaf-1 oligo (oligo number 1 in [D]), Apaf-1 produces similar results to Smac depletion (Figure 5), resulting in a persistent population of intermediate-cPARP cells that indicates a decrease in f (A). Co-depletion of XIAP and Apaf-1 restores efficient PARP cleavage (B), while codepletion of Smac and Apaf-1 has an additive inhibitory effect on PARP cleavage (more cPARP-negative cells persist at longer times [C]). These results favor a model in which both Smac and the apoptosome cooperate to antagonize XIAP, so that depletion of either Smac or the apoptosome component Apaf-1 results in a decrease in f (point 1 in Figure 9D). However, other Apaf-1 oligos, which appeared to deplete Apaf-1 with similar efficiency (D), did not produce a significant change in f, and did not act in an additive manner with Smac codepletion (not shown). Thus, the role of the apoptosome in snap-action switching cannot be conclusively determined from these results. Nonetheless, EARM v1.0 can be reconciled with both a scenario where the apoptosome plays a significant role, and one in which it does not (Figure 9D).
Figure S8. Role of Bcl-2/Bax Configuration in Maintaining the Off-State of the Snap-Action Switch
Using the current MOMP module (Model E in Figure 11) as a starting point, several variants were examined.
(A) Inhibition of Bax oligomeric forms by distinct pools of Bcl-2. Although it is more likely that a common pool of anti-apoptotic Bcl-2-like proteins bind to all oligomeric forms of Bax, we examined the behavior of a model in which a separate pool of Bcl-2 proteins mediates inhibition of each successive oligomeric form of Bax. Although this model is less plausible, it removes the potential for "implicit feedback," in which upstream steps of Bax oligomerization are accelerated by the sequestration of Bcl-2 by downstream oligomeric forms. Switching behavior is only slightly impaired in this model, indicating that implicit feedback is not a central component of the switching seen in Model E. Moreover, the slight loss of snap-action behavior is attributable, at least in part, to factors other than implicit feedback, such as the fact that the separate pools of Bcl-2 do not reach the point of depletion simultaneously.
(B) Effect of selective inhibition of Bax oligomeric forms. The six models shown are identical to that in the top panel of (A), with the exception that Bcl-2 was only permitted to bind to the oligomeric forms of Bax shown in the corresponding diagrams. In nearly all cases, snap-action behavior is impaired relative to a model in which Bcl-2 binds to all oligomeric forms of Bax (A, top), with models in which Bcl-2 inhibits only a single oligomeric form of Bax* performing weakest. Therefore, the ability of Bcl-2 to bind to multiple oligomeric forms of Bax may contribute substantially to snap-action behavior.
Protocol S1. A Matlab File Containing the Full EARM v1.0 Model and Routines for Producing the Simulations Shown in Figures 4–10
Instructions for choosing which simulations to run are included within the file.
Protocol S2. A Matlab Program, simFACS, for Producing Simulated Single Cell, Flow Cytometry, and Immunoblot Data
Running the program produces a graphical interface with user-editable fields for adjusting all switching parameters; pressing the "Plot" button updates the visualizations.
Table S2. EARM v1.0 Biochemical Reactions
Table S3. EARM v1.0 Rate Constant Values
Table S4. References for Rate Constants
Table S5. Variable Definitions and Initial Concentrations for EARM v1.0
Table S7. Sequences of siRNA Oligonucleotides
We thank F. Hua and M. Cardone for Bcl-2-overexpressing and control HeLa cells, Richard Youle for pECFP-Bax, and Doug Green for sharing unpublished data. We thank S. Gaudet and B. Schoeberl for help with modeling.
Author Contributions
JGA, JMB, SLS, DAL, and PKS worked closely together to design, perform and interpret the experiments; JGA and PKS wrote the paper.
- 1. Green DR, Evan GI (2002) A matter of life and death. Cancer Cell 1: 19–30.
- 2. Wang J, Zheng L, Lobito A, Chan FK, Dale J, et al. (1999) Inherited human Caspase 10 mutations underlie defective lymphocyte and dendritic cell apoptosis in autoimmune lymphoproliferative syndrome type II. Cell 98: 47–58.
- 3. Fuentes-Prior P, Salvesen GS (2004) The protein structures that shape caspase activity, specificity, activation and inhibition. Biochem J 384: 201–232.
- 4. Fischer U, Janicke RU, Schulze-Osthoff K (2003) Many cuts to ruin: a comprehensive update of caspase substrates. Cell Death Differ 10: 76–100.
- 5. Kischkel FC, Hellbardt S, Behrmann I, Germer M, Pawlita M, et al. (1995) Cytotoxicity-dependent APO-1 (Fas/CD95)-associated proteins form a death-inducing signaling complex (DISC) with the receptor. Embo J 14: 5579–5588.
- 6. Kim H, Rafiuddin-Shah M, Tu HC, Jeffers JR, Zambetti GP, et al. (2006) Hierarchical regulation of mitochondrion-dependent apoptosis by BCL-2 subfamilies. Nat Cell Biol 8: 1348–1358.
- 7. Letai A, Bassik M, Walensky L, Sorcinelli M, Weiler S, et al. (2002) Distinct BH3 domains either sensitize or activate mitochondrial apoptosis, serving as prototype cancer therapeutics. Cancer Cell 2: 183.
- 8. Verhagen AM, Ekert PG, Pakusch M, Silke J, Connolly LM, et al. (2000) Identification of DIABLO, a mammalian protein that promotes apoptosis by binding to and antagonizing IAP proteins. Cell 102: 43–53.
- 9. Du C, Fang M, Li Y, Li L, Wang X (2000) Smac, a mitochondrial protein that promotes cytochrome c-dependent caspase activation by eliminating IAP inhibition. Cell 102: 33–42.
- 10. Liu X, Kim CN, Yang J, Jemmerson R, Wang X (1996) Induction of apoptotic program in cell-free extracts: requirement for dATP and cytochrome c. Cell 86: 147–157.
- 11. Deveraux QL, Roy N, Stennicke HR, Van Arsdale T, Zhou Q, et al. (1998) IAPs block apoptotic events induced by caspase-8 and cytochrome c by direct inhibition of distinct caspases. Embo J 17: 2215–2223.
- 12. Rehm M, Dussmann H, Janicke RU, Tavare JM, Kogel D, et al. (2002) Single-cell fluorescence resonance energy transfer analysis demonstrates that caspase activation during apoptosis is a rapid process. Role of caspase-3. J Biol Chem 277: 24506–24514.
- 13. Tyas L, Brophy VA, Pope A, Rivett AJ, Tavare JM (2000) Rapid caspase-3 activation during apoptosis revealed using fluorescence-resonance energy transfer. EMBO Rep 1: 266–270.
- 14. Vaughan AT, Betti CJ, Villalobos MJ (2002) Surviving apoptosis. Apoptosis 7: 173–177.
- 15. Albeck JG, Burke JM, Aldridge BB, Zhang M, Lauffenburger DA, et al. (2008) Quantitative analysis of pathways controlling extrinsic apoptosis in single cells. Mol Cell 30: 11–25.
- 16. Fussenegger M, Bailey JE, Varner J (2000) A mathematical model of caspase function in apoptosis. Nat Biotechnol 18: 768–774.
- 17.
Eissing T, Conzelmann H, Gilles ED, Allgower F, Bullinger E, et al. (2004) Bistability analyses of a caspase activation model for receptor induced apoptosis. J Biol Chem.
- 18. Bentele M, Lavrik I, Ulrich M, Stosser S, Heermann DW, et al. (2004) Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis. J Cell Biol 166: 839–851.
- 19. Stucki JW, Simon HU (2005) Mathematical modeling of the regulation of caspase-3 activation and degradation. J Theor Biol 234: 123–131.
- 20.
Bagci EZ, Vodovotz Y, Billiar TR, Ermentrout GB, Bahar I (2005) Bistability in Apoptosis: Roles of Bax, Bcl-2 and Mitochondrial Permeability Transition Pores. Biophys J.
- 21. Hua F, Cornejo MG, Cardone MH, Stokes CL, Lauffenburger DA (2005) Effects of Bcl-2 levels on Fas signaling-induced caspase-3 activation: molecular genetic tests of computational model predictions. J Immunol 175: 985–995.
- 22. Legewie S, Bluthgen N, Herzel H (2006) Mathematical modeling identifies inhibitors of apoptosis as mediators of positive feedback and bistability. PLoS Comput Biol 2: e120.
- 23. Rehm M, Huber HJ, Dussmann H, Prehn JH (2006) Systems analysis of effector caspase activation and its control by X-linked inhibitor of apoptosis protein. Embo J 25: 4338–4349.
- 24. Cui J, Chen C, Lu H, Sun T, Shen P (2008) Two independent positive feedbacks and bistability in the Bcl-2 apoptotic switch. PLoS ONE 3: e1469.
- 25. Chen C, Cui J, Zhang W, Shen P (2007) Robustness analysis identifies the plausible model of the Bcl-2 apoptotic switch. FEBS Lett 581: 5143–5150.
- 26. Chen C, Cui J, Lu H, Wang R, Zhang S, et al. (2007) Modeling of the role of a Bax-activation switch in the mitochondrial apoptosis decision. Biophys J 92: 4304–4315.
- 27. Aldridge BB, Burke JM, Lauffenburger DA, Sorger PK (2006) Physicochemical modelling of cell signalling pathways. Nat Cell Biol 8: 1195–1203.
- 28. Stennicke HR, Jurgensmeier JM, Shin H, Deveraux Q, Wolf BB, et al. (1998) Pro-caspase-3 is a major physiologic target of caspase-8. J Biol Chem 273: 27084–27090.
- 29. Tewari M, Quan LT, O'Rourke K, Desnoyers S, Zeng Z, et al. (1995) Yama/CPP32 beta, a mammalian homolog of CED-3, is a CrmA-inhibitable protease that cleaves the death substrate poly(ADP-ribose) polymerase. Cell 81: 801–809.
- 30. Luo X, Budihardjo I, Zou H, Slaughter C, Wang X (1998) Bid, a Bcl2 interacting protein, mediates cytochrome c release from mitochondria in response to activation of cell surface death receptors. Cell 94: 481–490.
- 31. Zou H, Henzel WJ, Liu X, Lutschg A, Wang X (1997) Apaf-1, a human protein homologous to C. elegans CED-4, participates in cytochrome c-dependent activation of caspase-3. Cell 90: 405–413.
- 32. Li P, Nijhawan D, Budihardjo I, Srinivasula SM, Ahmad M, et al. (1997) Cytochrome c and dATP-dependent formation of Apaf-1/caspase-9 complex initiates an apoptotic protease cascade. Cell 91: 479–489.
- 33. Zhang XD, Zhang XY, Gray CP, Nguyen T, Hersey P (2001) Tumor necrosis factor-related apoptosis-inducing ligand-induced apoptosis of human melanoma is regulated by smac/DIABLO release from mitochondria. Cancer Res 61: 7339–7348.
- 34. Sun XM, Bratton SB, Butterworth M, MacFarlane M, Cohen GM (2002) Bcl-2 and Bcl-xL inhibit CD95-mediated apoptosis by preventing mitochondrial release of Smac/DIABLO and subsequent inactivation of X-linked inhibitor-of-apoptosis protein. J Biol Chem 277: 11345–11351.
- 35. Li S, Zhao Y, He X, Kim TH, Kuharsky DK, et al. (2002) Relief of extrinsic pathway inhibition by the Bid-dependent mitochondrial release of Smac in Fas-mediated hepatocyte apoptosis. J Biol Chem 277: 26912–26920.
- 36. Deng Y, Lin Y, Wu X (2002) TRAIL-induced apoptosis requires Bax-dependent mitochondrial release of Smac/DIABLO. Genes Dev 16: 33–45.
- 37. Slee EA, Harte MT, Kluck RM, Wolf BB, Casiano CA, et al. (1999) Ordering the cytochrome c-initiated caspase cascade: hierarchical activation of caspases-2, −3, −6, −7, −8, and −10 in a caspase-9-dependent manner. J Cell Biol 144: 281–292.
- 38. Sohn D, Schulze-Osthoff K, Janicke RU (2005) Caspase-8 can be activated by interchain proteolysis without receptor-triggered dimerization during drug-induced apoptosis. J Biol Chem 280: 5267–5273.
- 39. Cowling V, Downward J (2002) Caspase-6 is the direct activator of caspase-8 in the cytochrome c-induced apoptosis pathway: absolute requirement for removal of caspase-6 prodomain. Cell Death Differ 9: 1046–1056.
- 40. Yu X, Acehan D, Menetret JF, Booth CR, Ludtke SJ, et al. (2005) A structure of the human apoptosome at 12.8 A resolution provides insights into this cell death platform. Structure (Camb) 13: 1725–1735.
- 41. Yang JK, Wang L, Zheng L, Wan F, Ahmed M, et al. (2005) Crystal structure of MC159 reveals molecular mechanism of DISC assembly and FLIP inhibition. Mol Cell 20: 939–949.
- 42. Wallach D (1984) Preparations of lymphotoxin induce resistance to their own cytotoxic effect. J Immunol 132: 2464–2469.
- 43. Chipuk JE, Green DR (2008) How do BCL-2 proteins induce mitochondrial outer membrane permeabilization. Trends Cell Biol 18: 157–164.
- 44. Willis SN, Fletcher JI, Kaufmann T, van Delft MF, Chen L, et al. (2007) Apoptosis initiated when BH3 ligands engage multiple Bcl-2 homologs, not Bax or Bak. Science 315: 856–859.
- 45. Luo KQ, Yu VC, Pu Y, Chang DC (2003) Measuring dynamics of caspase-8 activation in a single living HeLa cell during TNFalpha-induced apoptosis. Biochem Biophys Res Commun 304: 217–222.
- 46. Takemoto K, Nagai T, Miyawaki A, Miura M (2003) Spatio-temporal activation of caspase revealed by indicator that is insensitive to environmental effects. J Cell Biol 160: 235–243.
- 47. Rehm M, Dussmann H, Prehn JH (2003) Real-time single cell analysis of Smac/DIABLO release during apoptosis. J Cell Biol 162: 1031–1043.
- 48. Goldstein JC, Waterhouse NJ, Juin P, Evan GI, Green DR (2000) The coordinate release of cytochrome c during apoptosis is rapid, complete and kinetically invariant. Nat Cell Biol 2: 156–162.
- 49. Li X, Darzynkiewicz Z (2000) Cleavage of Poly(ADP-ribose) polymerase measured in situ in individual cells: relationship to DNA fragmentation and cell cycle position during apoptosis. Exp Cell Res 255: 125–132.
- 50. Hafeman DG, McConnell HM, Gray JW, Dean PN (1982) Neutrophil activation monitored by flow cytometry: stimulation by phorbol diester is an all-or-none event. Science 215: 673–675.
- 51. Sledz CA, Holko M, de Veer MJ, Silverman RH, Williams BR (2003) Activation of the interferon system by short-interfering RNAs. Nat Cell Biol 5: 834–839.
- 52. Fulda S, Debatin KM (2002) IFNgamma sensitizes for apoptosis by upregulating caspase-8 expression through the Stat1 pathway. Oncogene 21: 2295–2308.
- 53. Shigeno M, Nakao K, Ichikawa T, Suzuki K, Kawakami A, et al. (2003) Interferon-alpha sensitizes human hepatoma cells to TRAIL-induced apoptosis through DR5 upregulation and NF-kappa B inactivation. Oncogene 22: 1653–1662.
- 54. Siegmund D, Wicovsky A, Schmitz I, Schulze-Osthoff K, Kreuz S, et al. (2005) Death receptor-induced signaling pathways are differentially regulated by gamma interferon upstream of caspase 8 processing. Mol Cell Biol 25: 6363–6379.
- 55. Wolter KG, Hsu YT, Smith CL, Nechushtan A, Xi XG, et al. (1997) Movement of Bax from the cytosol to mitochondria during apoptosis. J Cell Biol 139: 1281–1292.
- 56. Nechushtan A, Smith CL, Lamensdorf I, Yoon SH, Youle RJ (2001) Bax and Bak coalesce into novel mitochondria-associated clusters during apoptosis. J Cell Biol 153: 1265–1276.
- 57. Slee EA, Keogh SA, Martin SJ (2000) Cleavage of BID during cytotoxic drug and UV radiation-induced apoptosis occurs downstream of the point of Bcl-2 action and is catalysed by caspase-3: a potential feedback loop for amplification of apoptosis-associated mitochondrial cytochrome c release. Cell Death Differ 7: 556–565.
- 58. Tan C, Dlugosz PJ, Peng J, Zhang Z, Lapolla SM, et al. (2006) Auto-activation of the apoptosis protein Bax increases mitochondrial membrane permeability and is inhibited by Bcl-2. J Biol Chem 281: 14764–14775.
- 59. Seol DW, Li J, Seol MH, Park SY, Talanian RV, et al. (2001) Signaling events triggered by tumor necrosis factor-related apoptosis-inducing ligand (TRAIL): caspase-8 is required for TRAIL-induced apoptosis. Cancer Res 61: 1138–1143.
- 60. Gottfried Y, Rotem A, Lotan R, Steller H, Larisch S (2004) The mitochondrial ARTS protein promotes apoptosis through targeting XIAP. Embo J 23: 1627–1635.
- 61. Suzuki Y, Imai Y, Nakayama H, Takahashi K, Takio K, et al. (2001) A serine protease, HtrA2, is released from the mitochondria and interacts with XIAP, inducing cell death. Mol Cell 8: 613–621.
- 62. MacLachlan TK, El-Deiry WS (2002) Apoptotic threshold is lowered by p53 transactivation of caspase-6. Proc Natl Acad Sci U S A 99: 9492–9497.
- 63. Scaffidi C, Fulda S, Srinivasan A, Friesen C, Li F, et al. (1998) Two CD95 (APO-1/Fas) signaling pathways. Embo J 17: 1675–1687.
- 64. Korsmeyer SJ, Shutter JR, Veis DJ, Merry DE, Oltvai ZN (1993) Bcl-2/Bax: a rheostat that regulates an anti-oxidant pathway and cell death. Semin Cancer Biol 4: 327–332.
- 65. Crompton M, Barksby E, Johnson N, Capano M (2002) Mitochondrial intermembrane junctional complexes and their involvement in cell death. Biochimie 84: 143–152.
- 66. Karbowski M, Lee YJ, Gaume B, Jeong SY, Frank S, et al. (2002) Spatial and temporal association of Bax with mitochondrial fission sites, Drp1, and Mfn2 during apoptosis. J Cell Biol 159: 931–938.
- 67. Nakabayashi J, Sasaki A (2006) A mathematical model for apoptosome assembly: the optimal cytochrome c/Apaf-1 ratio. J Theor Biol 242: 280–287.
- 68. Holler N, Tardivel A, Kovacsovics-Bankowski M, Hertig S, Gaide O, et al. (2003) Two adjacent trimeric Fas ligands are required for Fas signaling and formation of a death-inducing signaling complex. Mol Cell Biol 23: 1428–1440.
- 69. Oltersdorf T, Elmore SW, Shoemaker AR, Armstrong RC, Augeri DJ, et al. (2005) An inhibitor of Bcl-2 family proteins induces regression of solid tumours. Nature 435: 677–681.
- 70. Fisher J, Henzinger TA (2007) Executable cell biology. Nat Biotechnol 25: 1239–1249.
- 71. Dorn JF, Jaqaman K, Rines DR, Jelson GS, Sorger PK, et al. (2005) Yeast kinetochore microtubule dynamics analyzed by high-resolution three-dimensional microscopy. Biophys J 89: 2835–2854.
- 72. Goldstein JC, Munoz-Pinedo C, Ricci JE, Adams SR, Kelekar A, et al. (2005) Cytochrome c is released in a single step during apoptosis. Cell Death Differ 12: 453–462.
- 73. von Ahsen O, Renken C, Perkins G, Kluck RM, Bossy-Wetzel E, et al. (2000) Preservation of mitochondrial structure and function after Bid- or Bax-mediated cytochrome c release. J Cell Biol 150: 1027–1036.
- 74. Madesh M, Antonsson B, Srinivasula SM, Alnemri ES, Hajnoczky G (2002) Rapid kinetics of tBid-induced cytochrome c and Smac/DIABLO release and mitochondrial depolarization. J Biol Chem 277: 5651–5659.
- 75. Ferrell JE Jr., Machleder EM (1998) The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science 280: 895–898.
- 76. Ferrell JE, Xiong W (2001) Bistability in cell signaling: How to make continuous processes discontinuous, and reversible processes irreversible. Chaos 11: 227–236.
- 77. Legewie S, Dienst D, Wilde A, Herzel H, Axmann IM (2008) Small RNAs establish delays and temporal thresholds in gene expression. Biophys J 95: 3232–3238.
- 78. Aldridge BB, Haller G, Sorger PK, Lauffenburger DA (2006) Direct Lyapunov exponent analysis enables parametric study of transient signalling governing cell behaviour. Syst Biol (Stevenage) 153: 425–432.
- 79. Shimizu S, Kanaseki T, Mizushima N, Mizuta T, Arakawa-Kobayashi S, et al. (2004) Role of Bcl-2 family proteins in a non-apoptotic programmed cell death dependent on autophagy genes. Nat Cell Biol 6: 1221–1228.
- 80. Rosado JA, Lopez JJ, Gomez-Arteta E, Redondo PC, Salido GM, et al. (2006) Early caspase-3 activation independent of apoptosis is required for cellular function. J Cell Physiol 209: 142–152.
- 81. Lu L, Chaudhury P, Osmond DG (1999) Regulation of cell survival during B lymphopoiesis: apoptosis and Bcl-2/Bax content of precursor B cells in bone marrow of mice with altered expression of IL-7 and recombinase-activating gene-2. J Immunol 162: 1931–1940.
- 82. Bannerji R, Kitada S, Flinn IW, Pearson M, Young D, et al. (2003) Apoptotic-regulatory and complement-protecting protein expression in chronic lymphocytic leukemia: relationship to in vivo rituximab resistance. J Clin Oncol 21: 1466–1471.
- 83. Raisova M, Hossini AM, Eberle J, Riebeling C, Wieder T, et al. (2001) The Bax/Bcl-2 ratio determines the susceptibility of human melanoma cells to CD95/Fas-mediated apoptosis. J Invest Dermatol 117: 333–340.
- 84. Stelling J, Sauer U, Szallasi Z, Doyle FJ 3rd, Doyle J (2004) Robustness of cellular functions. Cell 118: 675–685.
- 85. Vaughan AT, Betti CJ, Villalobos MJ, Premkumar K, Cline E, et al. (2005) Surviving apoptosis: a possible mechanism of benzene-induced leukemia. Chem Biol Interact 153–154: 179–185.
- 86. Janes KA, Albeck JG, Gaudet S, Sorger PK, Lauffenburger DA, et al. (2005) A systems model of signaling identifies a molecular basis set for cytokine-induced apoptosis. Science 310: 1646–1653.
- 87. Meraldi P, Draviam VM, Sorger PK (2004) Timing and checkpoints in the regulation of mitotic progression. Dev Cell 7: 45–60.
- 88. Gaudet S, Janes KA, Albeck JG, Pace EA, Lauffenburger DA, et al. (2005) A compendium of signals and responses triggered by prodeath and prosurvival cytokines. Mol Cell Proteomics 4: 1569–1590.
- 89. Bouralexis S, Findlay DM, Evdokiou A (2005) Death to the bad guys: targeting cancer via Apo2L/TRAIL. Apoptosis 10: 35–51.
- 90. Posakony JW, England JM, Attardi G (1977) Mitochondrial growth and division during the cell cycle in HeLa cells. J Cell Biol 74: 468–491.
- 91. Zhou HX (1993) Brownian dynamics study of the influences of electrostatic interaction and diffusion on protein-protein association kinetics. Biophys J 64: 1711–1726.
- 92. Berg OG (1985) Orientation constraints in diffusion-limited macromolecular association. The role of surface diffusion as a rate-enhancing mechanism. Biophys J 47: 1–14.
- 93. Wu JQ, Pollard TD (2005) Counting cytokinesis proteins globally and locally in fission yeast. Science 310: 310–314.
- 94. Kitano H, Funahashi A, Matsuoka Y, Oda K (2005) Using process diagrams for the graphical representation of biological networks. Nat Biotechnol 23: 961–966.