A simplified Bcl-2 network model reveals quantitative determinants of cell-to-cell variation in sensitivity to anti-mitotic chemotherapeutics

Anti-mitotic drugs constitute a major class of cytotoxic chemotherapeutics used in the clinic, killing cancer cells by inducing prolonged mitotic arrest that activates intrinsic apoptosis. Anti-mitotics-induced apoptosis is known to involve degradation of anti-apoptotic Bcl-2 proteins during mitotic arrest; however, it remains unclear how this mechanism accounts for significant heterogeneity observed in the cell death responses both within and between cancer cell types. To unravel quantitative determinants underlying variability in anti-mitotic drug response, we constructed a single-cell dynamical Bcl-2 network model describing cell death control during mitotic arrest, and constrained the model using experimental data from four representative cancer cell lines. The modeling analysis revealed that, given a variable, slowly accumulating pro-apoptotic signal arising from anti-apoptotic protein degradation, generation of a switch-like apoptotic response requires formation of pro-apoptotic Bak complexes with hundreds of subunits, suggesting a crucial role for high-order cooperativity. Moreover, we found that cell-type variation in susceptibility to drug-induced mitotic death arises primarily from differential expression of the anti-apoptotic proteins Bcl-xL and Mcl-1 relative to Bak. The dependence of anti-mitotic drug response on Bcl-xL and Mcl-1 that we derived from the modeling analysis provides a quantitative measure to predict sensitivity of distinct cancer cells to anti-mitotic drug treatment.

the slowly accumulating pro-apoptotic signal in prolonged mitotic arrest and identified depletion of Mcl-1, due to transcriptional silence, was one key pro-apoptotic trigger to activate mitotic death 9 . Moreover, by imaging a live-cell fluorescent reporter of mitochondrial outer membrane permeabilization (MOMP) 10 , the committed step of intrinsic apoptosis, we have shown that MOMP preceded nearly all cell death activated during mitotic arrest, and was rapid and switch-like, completing within minutes. MOMP is known to be regulated by Bcl-2 family proteins, such as Mcl-1; however, it is unresolved how a long, gradual pro-apoptotic signal from Mcl-1 depletion, which decays exponentially in the time scale of hours, may give rise to a sharp, all-or-none induction of apoptosis within minutes. In this study, we will perform both analytical and numerical analysis of the dynamics of a simplified Bcl-2 network to elucidate the quantitative mechanism that links a gradual, exponential signal to MOMP and the rapid MOMP induction across distinct timescales.
The other key question that we will address in this computational study is the quantitative origins of cell-to-cell variation in both sensitivity and kinetics to apoptosis during anti-mitotics-induced mitotic arrest. We chose to focus on analyzing mitotic death control, but not death after slippage, as it is the most variable point in the response to anti-mitotic drugs. Mcl-1 is known to be depleted to similar final levels in both apoptosis-sensitive and -resistant cell lines, thus loss of Mcl-1 alone cannot account for the apoptosis regulation during mitotic arrest. Based on results from gene knockdown by RNA interference (RNAi), we previously pinpointed Bcl-xL, but not Bcl-2, Bcl-w or pro-apoptotic BH3 proteins, as the other key regulator of apoptosis in mitotic arrest 9 . Variation in expression levels of Mcl-1 and Bcl-xL largely determine variability in sensitivity to mitotic death induced by anti-mitotics, such as paxlitaxel and Kinesin-5 inhibitor, across different cultured cancer cell lines. That is, the threshold for triggering cell death during mitotic arrest is mainly determined by basal expression levels of Mcl-1 and Bcl-xL. However, in order to employ Mcl-1 and Bcl-xL as diagnostic markers to predict anti-mitotic drug response in patients with distinct cancer types and heterogeneous tumor mass, we need to establish the quantitative, beyond qualitative, dependence of anti-mitotic drug response on Mcl-1 and Bcl-xL expression levels and their depletion kinetics, as well as determine to what extent the variation in the above dynamic parameters impacts the degree of variability in drug response both between cancer cell types (i.e., inter-cell line variability) and within a cancer type (intra-cell line variability). Therefore, in this study we conducted computational simulation of the simplified Bcl-2 network model for mitotic death control to analyze cell-to-cell variation at the single cell level, profiled the parameter space of levels and kinetics of Mcl-1 and Bcl-xL, and then derived the quantitative dependence of individual cell response to mitotic death induced by anti-mitotic drugs.

Results
Defining Bcl-2 network components for mitotic death control. We had previously identified the key Bcl-2 family proteins responsible for mitotic death control by studying four representative cancer cell lines: HeLa, U-2 OS, OVCAR-5 and A549 9 . These lines were chosen, as they cover a wide spectrum of sensitivity to mitotic death induced by anti-mitotic drugs (e.g., paclitaxel and Kinesin-5 inhibitor), based on drug response profiling experiments 5 . By knocking down candidate Bcl-2 family proteins using siRNA treatment and then determining the resulting effects on mitotic death, we found that knockdown of the anti-apoptotic proteins Bcl-xL and, to a lesser extent, Mcl-1, enhanced cell death during drug-induced mitotic arrest (Fig. 1A) in all four cell lines, in particular the three resistant cell lines (U-2 OS, OVCAR-5 and A549), albeit to varying degrees. In contrast, knockdown of Bcl-2 or Bcl-w showed mostly minimal effect across the cell lines, suggesting that they play largely negligible roles in regulating mitotic death. Together with data showing that activator BH3 proteins, such as Bim and tBid, or up-regulation of sensitizer BH3-only proteins are not required for mitotic cell death, we concluded that Mcl-1 and Bcl-xL are the key negative regulators of cell death acting during prolonged mitotic arrest 5,9 . Mcl-1 and Bcl-xL are known to inhibit cell death by sequestering the pro-apoptotic proteins Bak and Bax, which oligomerize to form pores in mitochondrial membrane to trigger MOMP. To determine whether Bak and/or Bax are involved in mitotic death control, we knocked down both proteins by RNAi in HeLa (the cell line that is highly sensitive to mitotic death), and determined the resultant effects on cell death by measuring Parp1 cleavage (Fig. 1B). We found that loss of Bak, but not Bax, significantly attenuated the extent of mitotic death (Fig. 1C), suggesting that Mcl-1 and Bcl-xL protect cells from mitotic death primarily through inhibitory interactions with Bak during mitotic arrest.
Transcription and translation are attenuated during prolonged mitotic arrest [11][12][13] , and we proposed that such mitotic silencing of gene expression selectively depletes unstable anti-apoptotic proteins to trigger cell death. Consistent with this idea, we found that Mcl-1 protein levels decreased steadily upon mitotic arrest in all four cell lines, with half-lives shorter than the observed timescales of mitotic cell death (τ d ~ 3-8 hrs, Fig. 1C). In further agreement with this model, Bcl-xL has a measured protein half-life comparable to the timescales of mitotic death induction (τ d ~ 10 hours 14 ), whereas Bak is significantly more stable (τ d ~ 170 hours 15 ). Taken together, these results point to a mechanism, where degradation of Mcl-1 and Bcl-xL protein during mitotic arrest relieves Bak from inhibition, allowing it to form pores that permeabilize the outer mitochondrial membrane and trigger cell death.
Modeling of the simplified Bcl-2 network identifies a requirement for high-order cooperativity in Bak pore formation. After cells enter mitotic arrest upon anti-mitotic drug treatment, they typically persist for many hours in a live mitotic state before undergoing apoptosis ( Fig. 2A, see also 16 ). The switch from a live to dead state with a permeable outer mitochondrial membrane occurs rapidly, reaching completion in minutes ( Fig. 2A and Supplementary Movie SM1). Here, we first determined whether the 'inhibitor decay mechanism' described above accounts for kinetic properties of such a cell death switch. Specifically, we used mathematical modeling to determine what Bcl-2 network reaction schemes can give rise to: 1) long time delays preceding apoptosis induction, typically 10 hours or more; and 2) rapid MOMP execution, within 20 minutes or less. Our analysis combines numerical simulations of ordinary differential simulations (Fig. 2B,C) with analytical approaches (see Materials and Methods), which generate insights that hold, independent of the exact parameter choices or specific network architecture.
Based on existing biochemical evidence 17 , we first constructed a model, where Bak monomer either binds Mcl-1 to form an inactive complex, or undergoes two sequential dimerization reactions to form tetramer pores on the mitochondrial membrane, resulting in cytochrome C (CytC) translocation from mitochondria to cytoplasm (Model I, Fig. 2B). The total Bak protein level remains constant due to its high stability, whereas Mcl-1 level decreases over time with first-order kinetics. From numerical simulations (see Materials and Methods, and Table 1), we found that this model can give rise to delayed apoptosis induction after Mcl-1 decay, as intuitively expected (Fig. 2B). After the onset of mitotic arrest (t = 0 hrs), the concentration of free Mcl-1 decreases over time, causing an eventual increase in the concentration of Bak tetramer and CytC in the cytoplasm. However, while the simulation results account for the delayed timing of apoptosis switch (~10 hours after mitotic arrest), they were not able to recapitulate the sharp switch itself. With the simulations, we found that the switch timing ΔT, taken to be the time for the fraction of CytC in the cytoplasm to increase from 0.1 to 0.5, was 3.8 hours, much longer than the switching time of < 20 minutes observed experimentally ( Fig. 2A). This suggests that this Bak tetramer model, while intuitively appealing, is not sufficient to account for the switch-like properties of the apoptotic response in mitotically arrested cells.
In signaling systems, switch-like behavior can emerge from induced clustering of membrane-associated receptor proteins into large complexes [18][19][20] . Motivated by these examples, we examined an alternative possibility, where Bak oligomerizes to form large pores containing hundreds of subunits (Model II). We analyzed an example reaction scheme, where Bak monomer undergoes multiple sequential dimerization reactions to generate an active pore consisting of 256 subunits. Nonetheless, the modeling results, as shown below, do not depend on the exact number of subunits of the active pore, or on the exact reaction scheme for generating this active pore. As before, Mcl-1 forms an inactive complex with Bak, and degrades with first-order kinetics. From numerical simulations (see Materials and Methods, and Table 1), we see that this model not only gave a delayed induction of Bak pores upon Mcl-1 protein decay (Fig. 2C, top and middle), but also generated a much sharper MOMP reaction compared to the tetramer pore model (ΔT = 0.44 hours), in closer agreement with experimental observations ( Fig. 2A). These results suggest that a switch-like apoptotic response requires high-order cooperativity from the formation of mitochondrial membrane pores consisting of hundreds of Bak/Bax subunits. To establish the generality of this result, we solved the models analytically to obtain T c , the duration of the time delay before apoptosis, and ΔT, the induction time of apoptosis. For both models, we find that:  Quantitatively, in order to generate a sharp switch of MOMP after a long delay, we must fulfill the requirement that Δ T ≪ T C , which requires: The logarithm of the fold change in free Mcl-1 levels on the right hand side of this expression is of order unity; therefore, consistent with our simulation results, this analytical result shows that sharp switching can occur only when the number of subunits in the active Bak pore, A, is much greater than one. Indeed, recent imaging studies have found the Bax and Bak structures on the mitochondrial membranes of apoptotic cells are estimated to contain hundreds or even thousands of subunits [21][22][23] , consistent with this notion. Our modeling results, taken together with these experimental studies, argue that sharp apoptotic switching cannot emerge simply through formation of small pores with < 10 subunits, but requires highly cooperative assembly of large mitochondrial membrane pores with hundreds or even thousands of Bak/Bax subunits.
Bcl-2 network modeling recapitulates and quantifies single cell variability in apoptotic responses for different cancer cell lines during mitotic arrest. Given that the simplified Bcl-2 network model (i.e., model II) accounts for the basic kinetic properties of mitotic death response, we next investigate whether it can be employed to examine the variable single-cell death responses in different cancer cell types during mitotic arrest, e.g., the four representative cell lines HeLa, U-2 OS, OVCAR-5 and A549 (Fig. 1A). Previous studies have implicated cell-to-cell differences in protein concentrations as a cause of non-genetic variability in apoptosis timing in human cell lines 24 . Therefore, we combined our kinetic model with probability distributions of apoptotic protein abundances in single cells to obtain predicted mitotic survival curves, which give the fraction of surviving mitotic cells as a function of the duration of mitotic arrest (Fig. 3A, center; see also Materials and Methods). In these single cell simulations, we explicitly included Bcl-xL level, in line with experimental results (Fig. 1A, see Materials and Methods), and also included an implicit requirement for switch-like apoptosis induction by large Bak oligomers, in line with insights from kinetic modeling (Fig. 2). We then used least-squares fitting to fit these survival curves to those acquired experimentally for the four cell lines, both under control conditions and upon knockdown of either Mcl-1 or Bcl-xL. Fits were constrained using measured Mcl-1 protein half-life data (Fig. 1C) and quantitative protein level measurements from western blotting (means and standard deviations shown in Fig. 3C as crosses), which constrained the means of the single-cell protein level distributions.
From model fitting and analysis, we found that this simple model for apoptosis induction recapitulates the following key properties of the variable apoptotic response in the four cell lines that we previously studied (  (Fig. 1D). In addition to accounting for these experimental measurements, our model also recapitulates the published decay rates for Bcl-xL 14 (τ ~ 10-20 hrs, Fig. 3C), which were not used to constrain our model. While small differences in the shapes of the mitotic survival curves from experimental data were observed (Fig. 3B), the fairly close agreement of the models with data on multiple levels allowed us to derive quantitative insights and generate hypotheses beyond the experimental findings. Specifically, we made two major observations and predictions as follow: Firstly, variability in the levels of anti-apoptotic regulators underlies heterogeneity in death responses within cell lines (i.e., intra-cell line variation). In all four cell lines examined, there was considerable cell-to-cell variability in the timing of apoptosis induction within the cell line, with some cells dying shortly after mitotic arrest, and others persisting for over ten hours before dying (Fig. 3B). This intra-cell line heterogeneity could in principle arise from non-genetic variability in the levels of the anti-apoptotic regulators Mcl-1 or Bcl-xL, or in the level of the pro-apoptosis regulator Bak. Our modeling fits showed considerable cell-to-cell variation in the expression of anti-apoptotic regulators, with cells showing over three-fold variation in the levels of Mcl-1 (U-2 OS), Bcl-xL (Ovcar-5 and A549), or both proteins (HeLa). In contrast, the variation in the levels of the pro-apoptotic protein Bak was considerably smaller in all of the cell lines examined. These observations suggest that cells of the same genetic background may tightly regulate their expression levels of pro-apoptotic proteins, and generate variable death responses primarily through heterogeneous expression of anti-apoptotic proteins.
Secondly, different average levels of apoptotic regulatory proteins underlie the differential susceptibility of cell lines to mitotic death (i.e., inter-cell line variation). The inter-cell line variability is obviously much stronger than that within the same cell line. The four cell lines also showed considerably differing responses to knockdown of anti-apoptotic proteins. Mcl-1 knockdown greatly accelerated apoptotic death in HeLa cells, but had negligible effects in A549 cells; conversely, removal of Bcl-xL increased both the extent and kinetics of death in U-2 OS, OVCAR-5 and A549 cells, but showed only moderate effects in HeLa cells. These differential susceptibilities could arise from cell-line specific differences in protein expression levels, or differences in their degradation rates during mitotic arrest. Our modeling fits showed that HeLa cells, which showed the greatest sensitivity to Mcl-1 knockdown, had the highest Mcl-1/Bak ratio but the lowest Bcl-xL/Bak ratio out of the cell lines examined ( Fig. 3B-D); conversely, U-2 OS, OVCAR-5 and A549 cells, which were more susceptible to Bcl-xL knockdown, but not Mcl-1 knockdown, had comparatively lower Mcl-1/Bak ratios but progressively increasing Bcl-xL/ Bak ratios. In contrast, we found no correlation between death sensitivity and the half-lives of Mcl-1 or Bcl-xL (Fig. 3D), even though they showed considerable variation between different cell lines. Instead, these parameters appeared to affect the total duration of mitotic arrest, which was longest in cells with the highest stability of these proteins (i.e., OVCAR-5, Fig. 3B-D). Together, our data illustrate that inter-cell line variability depends on the highly expressed anti-apoptotic proteins, and suggest that ratio of the levels of an anti-apoptotic protein to its pro-apoptotic partner may be used as a useful quantitative predictor of the susceptibility of distinct cancer cell types to mitotic death activated by anti-mitotic drugs.

Discussion
Building on previous experimental data, we constructed a simplified Bcl-2 network model consisting of three key components, Mcl-1, Bcl-xL and Bak, to elucidate the dynamic control of mitotic cell death induced by anti-mitotic drugs and the associated cell-to-cell variability. Computational analysis of this simple kinetic model not only revealed critical dynamical features of the network, e.g., high-order cooperativity arising from massive Bak oligomerization, but also allowed us to examine and pinpoint specific network parameters that provide a metric to predict the susceptibility of different cancer cells to mitotic death induced by anti-mitotic drugs. We found heterogeneity in mitotic death response within a cancer type can be attributed to mainly variation in Mcl-1 and Bcl-xL expression levels in individual cells, but not Bak, with a 2 to 3-fold difference in expression being sufficient to generate significant variation in drug-induced cell fate, i.e., dead vs. live, in a clonal population with identical genetic background. Response variation between different cancer types is also mainly attributed to variable expression, in this case ratio of Bcl-xL/Bak and Mcl-1/Bak. The extent of mitotic death is particularly dependent on the ratio of Bcl-xL/Bak, as a 2-fold difference in this ratio can already distinguish cancer lines that are resistant (e.g., U-2 OS) and sensitive (e.g., HeLa) to mitotic death. In other words, our results suggest that tumor cells with a Bcl-xL/Bak ratio 2-fold larger than that of HeLa are resistant to mitotic death. Moreover, although there is considerable variation in Mcl-1 and Bcl-xL degradation kinetics between different cancer types, our dynamic modeling analysis showed that such variability mostly affects the timing to mitotic death, but not the extent.
Although the simplified Bcl-2 model for mitotic death control was constructed based on cell culture data, there are in vivo data that also support its validity. In a study of mouse xenograft models of non-small lung cancer, Tan et al. found that Bcl-xL and Mcl-1 are the key MOMP regulators with respect to paclitaxel responses in vivo 25 . Their mouse model data are clearly consistent with our results with cultured cell lines 9 , regarding the key regulatory roles of Bcl-xL and Mcl-1 in regulating anti-mitotics-induced cell death. Therefore, our Bcl-2 network model for cell death control upon anti-mitotic drug treatment is likely applicable for in vivo situation. Nonetheless, we note there are data that point to potential difference in the mechanism by which cell death is activated in vivo vs. in vitro by anti-mitotic drugs. While cell death seen in culture was all preceded by prolonged mitotic arrest, Orth et al. found that many tumor cells in a mouse xenograpt model of Fibrosarcoma did not enter or progress through paclitaxel-induced mitotic arrest, even though tumor significantly regressed 26 . The question of whether mitotic arrest is required for anti-mitotics-induced cell death in vivo, or paclitaxel can induce cell death in interphase without entry into mitotic arrest, clearly needs further mechanistic investigation. In addition, we also note that a previous study by Topam et al. showed that although anti-mitotics-induced cell death in cultured cell line was mainly repressed by Bcl-xL, it also involved up-regulation of pro-apoptotic BH-3 proteins 27 . Their result of BH-3 protein involvement in regulating mitotic death was different from our previous experimental data on a cancer cell line panel 5,9 . Such discrepancy again requires further study to elucidate and clarify its mechanistic origin.
Both numerical and analytical results from the simplified Bcl-2 network model identified a critical requirement for high-order cooperativity in the generation of a sharp switch of MOMP. The major pro-apoptotic signals that activate mitotic cell death are Mcl-1 and Bcl-xL degradation, which occur in the time scale of hours following mitotic arrest. We find that this long, exponentially decaying signal can give rise to a sharp mitotic death response that completes in minutes, only if Bak pore complexes on the mitochondrial membrane are very large (with > 10 2 subunits). In cells, assembly of such large protein complexes is typically highly cooperative, as it involves a phase transition process, occurring spontaneously when a critical parameter -in this case free Bak monomer concentration -reaches a threshold value. Such phase transition behavior is known to govern assembly in various biological systems, including cytoskeletal polymers 28 , amyloid fibers 29 , and membrane signaling clusters 30 , and our results indicate that it may also be utilized by pro-apoptotic proteins to generate all-or-none death responses. Indeed, it is been long suggested that the active pores that execute MOMP on apoptotic cells contain hundreds or even thousands of subunits [21][22][23] , consistent with this idea. Further testing of this hypothesis will require a combination of biochemical approaches, in conjunction with genetic and modeling studies.
In this study, we focused on analyzing the quantitative determinants of mitotic cell death induced by anti-mitotic drugs, as it is the drug mechanism that activates the most rapid cancer cell death and is also the most variable point in anti-mitotic drug response. However, anti-mitotic drugs are known to trigger cell death not only during mitotic arrest but also after mitotic slippage into an abnormal G1 state 5,31 . Molecular regulators that control cell death after slippage are distinct from those during mitotic arrest, mainly involving proteins associated with DNA damage response, e.g., the p53 pathway [32][33][34] . For simulation and fitting of the Bcl-2 network model in Scientific RepoRts | 6:36585 | DOI: 10.1038/srep36585 this study, we explicitly excluded the slippage process and cell population that have exited the mitotic arrest, so as to examine only dynamic features of cell death control during mitotic arrest (refer to Materials and Methods). Given that cell fate after mitotic slippage is mainly regulated by DNA damage response, the quantitative predictors of death response after slippage can be investigated in further study, using various mathematical models of DNA damage response that are already available in the literature, such as the p53 pathway models.

Methods
Cell culture. All cell lines were purchased from American Type Culture Collection (ATCC, USA) and cultured under 37 °C and 5% CO 2 in appropriate medium supplemented with 10% Fetal Calf Serum (FCS), 100 U/ml penicillin and 100 μ g/ml streptomycin. HeLa was maintained in DMEM; U-2 OS was maintained in McCoy's; OVCAR-5 was maintained in RPMI; and A549 was maintained in F-12K. The anti-mitotic drug, Kinesin-5 inhibitor (EMD534085), was provided by Merck-Serono.
Time-lapse microscopy. Cells were plated in 35 mm imaging dish (μ -dish, ibidi, Germany) and cultured in phenol red-free CO 2 -independent medium (Invitrogen) supplemented with 10% FCS, 100 U/ml penicillin and 100 μ l streptomycin. Cell images were acquired with the Nikon TE2000-PFS inverted microscope enclosed in a humidified chamber maintained at 37 °C. Cells were imaged every 10 minutes using a motorized stage and a 20X objective (NA = 0.95). Images were viewed and analyzed using the MetaMorph software (Molecular Dynamics).

Western blot analysis.
Cell lysates were obtained using LDS sample buffer (NuPAGE, Invitrogen). Proteins were resolved on 10% or 12% Tris-glycine gels and transferred onto PVDF membranes. Blots were probed with commercial primary antibodies and chemiluminescent detection using ECL-plus (Amersham). Primary antibodies: PARP1 (#9542), Bak (#3814) and Bcl-xL (#2762) were purchased from Cell Signaling; Bax (#sc-493) and Mcl-1 (#sc-819) from Santa Cruz; Anti-actin (#A5316) from Sigma was used as a loading control. For western blot analysis of Mcl-1 half-lives in synchronized mitotic cells, we grew large volume of cells in 25 cm dishes to 90% confluency, and then treated the cells with 1 μ M Kinesin-5 inhibitor to induce mitotic arrest. After 3 hours of drug treatment, the mitotic fraction of cells was collected at the indicated time points and lysed using LDS sample buffer for western blot analysis of Mcl-1.
A simplified Bcl-2 network model of Bak pore formation and apoptosis induction in mitotic arrest. Model description. We present here two ordinary differential equation (ODE) models for apoptosis induction, explicitly considering degradation of Mcl-1 protein during mitotic arrest, sequestration of Bak monomers by the anti-apoptotic protein Mcl-1, Bak oligomerization and resultant pore formation, and transport of Cytochrome C across the outer mitochondrial membrane. In Model I (Fig. 2B), the Bak tetramer forms the active mitochondrial membrane pore. This model is described by the following equations:  Here, M represents Mcl-1; B n represents the Bak n-mer formed through the indicated oligomerization reactions; α and β give the rate constants of association and dissociation for the indicated Bak n-mers; σ M gives the basal Mcl-1 synthesis rate during mitosis; and δ M gives the first-order rate constant of Mcl-1 degradation during mitotic arrest. Here, we assume that Bak shows negligible degradation over the timescales of mitotic arrest, which is consistent with previous stability measurements of this protein in isotope switching proteomic experiments 15 . For Model II (Fig. 2C), we assume the active membrane pore is Bak oligomer consisting of 256 subunits, formed through successive dimerization of smaller oligomers, as follows: Here, Mcl-1 dynamics, Mcl-1/Bak monomer interactions and the Bak dimerization reaction all remain the same, and are given by (Eqs 4-6). While this model assumes that the active pore has 256 subunits, our subsequent analysis will show that our conclusions will neither depend on the exact size of the active pore complex, nor on the detailed reaction scheme underlying the generation of the massively oligomeric pore. In both models, the total concentration of Bak is a constant: where the summation is performed over all Bak oligomeric species. Finally, for both models, the translocation of CytC from the mitochondria to the cytoplasm is given by the following equations: Numerical simulations. Simulations for both Models (Fig. 2B,C) were performed using numerical integration with a stiff ODE solver using the MATLAB SimBiology toolbox (Mathworks, Natick, MA). Initial conditions and parameter values for the simulations are shown in Table 1, and the simulation code is available as project files upon request.
Analytical solutions. To derive an analytical solution of the dynamics of the system, we first make the assumption that the timescales of Mcl-1 degradation are slow compared to that of Bak oligomerization and CytC translocation. This allows us to set the time derivatives of all differential equations, except (Eq. 4), to zero, and solve for the pseudo-equilibrium concentrations of all other species. Doing so, we find a single first-order equation that governs the dynamics of this system: M M which has the solution:  where the dissociation constants are 1 2 1 2 , with higher order dissociation constants and equilibrium expressions being similarly defined. If we further assume that only a small percentage of Bak is oligomerized during apoptosis induction, as observed experimentally 35 , such that most Bak is either bound to Mcl-1, or unbound in monomeric form, we get that: T 1 By combining this equation with (Eq. 20), we get the following relationship between free Mcl-1 level and Bak monomer level: Next, by solving for the fraction of CytC in the cytoplasm using (Eqs 16-17), we can show that: where A is the number of subunits in the active pore complex (A = 4 for Model I, and A = 256 for Model II), and K C is the critical value of [B 1 ], at which the level of cytoplasmic CytC is half its maximal value. It can be shown that this critical value is a function of the dissociation constants of the oligomerization reaction (α and β), and also the forward and back rate constants for CytC translocation (γ c and γ m ). Based on these approximations, we now derive analytical expressions for: 1) time to the MOMP transition T c , which we define to be the time at which there is half-maximal CytC translocation to the cytoplasm, and 2) sharpness of the MOMP transition ΔT, which we define to be the time between 1/10 and 1/2 maximal CytC translocation. To derive these quantities, we first use (Eqs 24-25) to derive the concentrations of free Mcl-1, where CytC translocation is 1/10 and 1/2 maximal: M A T C 1/10 1/ Next, by substituting (Eq. 26) in (Eq. 19), we find that the time delay until the MOMP transition is given by: where C = In(9) ≈ 2.2. Note that lower bound for the MOMP switch timing is independent of the initial Mcl-1 level M 0 , as to be expected, and scales inversely with the size of the Bak pore. It is also independent of the detailed rate constants of the dimerization reactions, suggesting that this result does not depend on the exact reaction scheme, through which active pores are assembled.
Computational derivation of mitotic survival curves based on the Bcl-2 network model. Using the above kinetic model (Eqs 8-22), we now generate predicted curves for the fraction of surviving mitotic cells over time, which will then be used to fit curves from experimental data. To better account for the cell death kinetics in the four cell lines studied, we first incorporate interactions between Bak and Bcl-xL by adding to the model the following rate equation: where [X] is the concentration of Bcl-xL. By solving for the steady state of this system using (Eq. 30) and the above equations, we obtain the following equation for the time evolution of free Bak monomer: T M X

1
As the rest of the equations describing Bak oligomerization remain unchanged, the relationship between the concentration of Bak monomer and the cytoplasmic CytC fraction remain the same: Taking into account our results (Fig. 2C), we now take the size of the Bak apoptotic pore to be very large, such that A→ ∞ . This allows calculation of the time of apoptotic induction t c , which occurs when [B 1 ] = K c . By incorporating this equation into (Eq. 31), we now get that: