Apoptosis is a highly regulated cell death mechanism involved in many physiological processes. A key component of extrinsically activated apoptosis is the death receptor Fas which, on binding to its cognate ligand FasL, oligomerize to form the death-inducing signaling complex. Motivated by recent experimental data, we propose a mathematical model of death ligand-receptor dynamics where FasL acts as a clustering agent for Fas, which form locally stable signaling platforms through proximity-induced receptor interactions. Significantly, the model exhibits hysteresis, providing an upstream mechanism for bistability and robustness. At low receptor concentrations, the bistability is contingent on the trimerism of FasL. Moreover, irreversible bistability, representing a committed cell death decision, emerges at high concentrations which may be achieved through receptor pre-association or localization onto membrane lipid rafts. Thus, our model provides a novel theory for these observed biological phenomena within the unified context of bistability. Importantly, as Fas interactions initiate the extrinsic apoptotic pathway, our model also suggests a mechanism by which cells may function as bistable life/death switches independently of any such dynamics in their downstream components. Our results highlight the role of death receptors in deciding cell fate and add to the signal processing capabilities attributed to receptor clustering.
Many prominent diseases, most notably cancer, arise from an imbalance between the rates of cell growth and death in the body. This is often due to mutations that disrupt a cell death program called apoptosis. Here, we focus on the extrinsic pathway of apoptotic activation which is initiated upon detection of an external death signal, encoded by a death ligand, by its corresponding death receptor. Through the tools of mathematical analysis, we find that a novel model of death ligand-receptor interactions based on recent experimental data possesses the capacity for bistability. Consequently, the model supports threshold-like switching between unambiguous life and death states; intuitively, the defining characteristic of an effective cell death mechanism. We thus highlight the role of death receptors, the first component along the apoptotic pathway, in deciding cell fate. Furthermore, the model suggests an explanation for various biologically observed phenomena, including the trimeric character of the death ligand and the tendency for death receptors to colocalize, in terms of bistability. Our work hence informs the molecular basis of the apoptotic point-of-no-return, and may influence future drug therapies against cancer and other diseases.
Apoptosis is typically viewed as a bistable system, with a sharp all-or-none switch between attracting life and death states. This bistability is important for conferring robustness . Consequently, researchers have used computational models to identify and study potential sources of bistability in apoptosis, including positive caspase feedback , inhibition of DISC by cFLIP , cooperativity in apoptosome formation , double-negative caspase feedback through XIAP , and double-negative feedback in Bcl-2 protein interactions . In this work, we propose that bistability may be induced upstream by the death receptors themselves.
The current model of death ligand-receptor dynamics assumes that FasL activates Fas by direct crosslinking, producing a DISC concentration that varies smoothly with the ligand input . However, recent structural data  suggests a different view. In particular, Fas was found in both closed and open forms, only the latter of which allowed FADD binding and hence transduction of the apoptotic signal. Moreover, open Fas were observed to pair-stabilize through stem helix interactions. This affords a mechanism for bistability, similar to the Ising model in ferromagnetism , where open Fas, presumably disfavored relative to their native closed forms , are able to sustain their conformations even after removal of the initial stimulus promoting receptor opening, past a certain critical density of open Fas. This induces hysteresis in the concentration of active, signaling receptors and therefore in apoptosis.
We studied this proposed mechanism by formulating and analyzing a mathematical model. The essential interpretation is that FasL acts as a clustering platform for Fas, which establish contacts with other Fas through pairwise and higher-order interactions to form units capable of hysteresis (Figure 1). At low receptor concentrations, the model exhibits bistability provided that the number of receptors that each ligand can coordinate is at least three. This hence gives a theory for the trimeric character of FasL. Furthermore, at high concentrations, for example, through receptor pre-association – or localization onto lipid rafts , irreversible bistability is achieved, implementing a permanent cell death decision. Thus, our model suggests a primary role for death receptors in deciding cell fate. Moreover, our results offer novel functional interpretations of ligand trimerism and receptor pre-association and localization within the unified context of bistability.
The transmembrane death receptor Fas natively adopts a closed conformation, but can open to allow the binding of FADD, an adaptor molecule that facilitates apoptotic signal transduction. Open Fas can self-stabilize via stem helix and globular interactions, which is enhanced by receptor clustering through association with the ligand FasL.
Constructing a mathematical model of Fas dynamics is not entirely straightforward as receptors can form highly oligomeric clusters , . A standard dynamical systems description would therefore require an exponentially large number of state variables to account for all combinatorial configurations. To circumvent this, we considered the problem at the level of individual clusters. Each cluster can be represented by a tuple denoting the numbers of its molecular constituents, the cluster association being implicit, so only these molecule numbers need be tracked.
In our model, a cluster is indexed by a tuple , where represents FasL and , , and are three posited forms of Fas, denoting closed, open and unstable, and open and stable, i.e., active and signaling, receptors, respectively. Within a cluster, we assumed a complete interaction graph and defined the reactions(1a)(1b)(1c)(1d)The first reaction describes spontaneous receptor opening and closing; the second, constitutive destabilization of open Fas; the third, ligand-independent receptor cluster-stabilization; and the fourth, ligand-dependent receptor cluster-stabilization (Figure 2). The orders of the cluster-stabilization events are limited by the parameters and , which capture the effects of receptor density and Fas coordination by FasL, respectively. Although only pair-stabilization () has been observed experimentally , higher-order analogues, for example, as facilitated by globular interactions, are not unreasonable.
Examples of ligand-independent cluster-stabilization reactions involving unstable () and stable () open receptors of molecularities two (A), three (B), and four (C). Higher-order reactions follow the same pattern. Ligand-dependent reactions are identical except that FasL () must be added to each reacting state.
Formally, these reactions are to be interpreted as state transitions on the space of cluster tuples. However, the reaction notation is suggestive, highlighting the contribution of each elementary event, which we modeled using constant reaction rates (for simplicity, we set uniform rate constants and for all ligand-independent and -dependent cluster-stabilization reactions of molecularity , respectively). Then on making a continuum approximation, we reinterpreted the molecule numbers as local concentrations and applied the law of mass action to produce a dynamical system for each cluster in the concentrations of . Validity of the model requires that the molecular concentrations are not too low and that the timescale of receptor conformational change is short compared to that of cluster dissociation.
To study the long-term behavior of the model, we solved the system at steady state (denoted by the subscript ). Introducing the nondimensionalizations(2a)(2b)(2c)(2d)(2e)where is a characteristic concentration and is time, and(3a)(3b)(3c)(3d)this is(4a)(4b)where is the nondimensional total receptor density, and is given by considering(5)and solving with , a polynomial in of degree . Clearly, the model is bistable only if (two stable nodes must be separated by an unstable node as the model is effectively one-dimensional in ).
We used as a measure of the apoptotic activation of a cluster. In principle, all open receptors contribute to apoptotic signaling, but is small, at least at steady state (since due to the assumed prevalence of the closed form ), and so can be neglected.
Bistability and receptor clustering
While measures the coordination capacity of FasL and hence may be equated with its oligomeric order (e.g., in the biological context), an appropriate value for , relating to the total receptor concentration, is somewhat more elusive. Therefore, we began our analysis by performing a simple receptor density estimate. Approximating the cell as a cube of linear dimension m, the associated volume of pL implies the correspondence nM molecules molecules/nm on restricting to the membrane, i.e., by averaging over the surface area of m. Thus, for a conservative receptor concentration estimate of nM , , , , the number of Fas molecules in the neighborhood of each receptor is only , assuming a charateristic size of nm. We hence found that receptors may be very sparsely distributed. In this low density mode, high-order Fas interactions in the absence of ligand can be neglected (). Therefore, in this context, bistability is possible only if , and the trimerism of FasL thus demonstrates the lowest-order complexity required for bistability.
From the form of , this bistability is reversible as a function of the FasL concentration since the governing polynomial for is of degree only at . This suggests that at the cluster level, the cell death decision can be reversed, which may have adverse effects on cellular and genomic integrity. However, irreversible bistability at higher receptor densities may also be achieved. Researchers have observed tendencies for death receptors both to pre-associate as dimers or trimers – and to selectively localize onto membrane lipid rafts . The result of either of these processes may be to increase the local receptor concentration. In this high density mode, we set , as the preceeding approximation is no longer valid. Irreversible bistability then becomes attainable, representing a committed cell death decision.
For the remainder of the study, we incorporated both the low and high receptor density regimes into a single model by setting , using as a continuous transition parameter. Furthermore, we set to correspond to observed biology.
Characterization of the steady-state surface
Calculation of the steady-state activation curves showed that the model indeed exhibits bistability (Figure 3) for reasonable parameter choices (Methods). Thus, we established the possiblity of a novel bistability mechanism in extrinsic apoptosis. The associated hysteresis enables threshold switching between well-separated low and high activation states. Biologically, these define local signals of life and death, which are integrated at the cell level to compute the overall apoptotic response.
The steady-state active Fas concentration shows bistability and hysteresis as a function of the FasL concentration (stable, solid lines; unstable, dashed lines). At low receptor concentrations , the bistability is reversible, but irreversibility emerges for sufficiently high, representing a committed cell death decision. All parameters set at baseline values unless otherwise noted.
As per the previous analysis, reversibility of the bistability is dependent on , with irreversibility emerging for sufficiently high. This suggests a bivariate parameterization of , producing a multivalued steady-state surface over -space (Figure 4). The result is a cusp catastrophe, an elementary object of catastrophe theory, which studies how small perturbations in certain parameters can lead to large and sudden changes in the behavior of a nonlinear system . A more instructive view of the dependence of the model's qualitative structure on and is shown in Figure 5.
The steady-state surface for the active Fas concentration as a function of the FasL and total Fas concentrations and , respectively, is folded, indicating the existence of singularities, across which the system's steady-state behavior switches between monostability and bistability (stable, blue; unstable, red). All parameters set at baseline values unless otherwise noted.
Steady state diagram identifying the regions of parameter space supporting monostability (colored) or bistability (gray) as a function of the FasL and total Fas concentrations and , respectively. The monostable region is colored as a heat map corresponding to the steady-state active Fas concentration . Irreversible bistability is indicated by the extension of the bistable region to the axis .
Sensitivity and robustness analyses
We then focused on the activation and deactivation thresholds , respectively, defining the bistable regime. These are the points at which the steady state switches discontinuously from one branch to the other, and are given by the values of at which the hysteresis curve turns, i.e., at (Figure 6). We performed a sensitivity analysis of by measuring the effects of perturbing the model parameters about baseline values (Methods). For each threshold-parameter pair, we computed a normalized sensitivity by linear regression.
The activation (red) and deactivation (blue) thresholds characterizing the bistable regime (green) are defined as the concentrations of FasL at which the steady-state active Fas concentration (black) switches discontinuously from one branch to the other (stable, solid line; unstable, dashed line).
Strong effects of , , and were observed (Figure 7); for the corresponding Fas thresholds at , respectively, the parameters , , , and were emphasized. Thus, the bistability thresholds do not appear particularly robust. However, the data reveal that essentially all parameter sets sampled were bistable. This suggests a weaker form of robustness, namely, robustness of bistability, which nevertheless supports life and death decisions over a wide operating range.
The robustness of the bistability thresholds is investigated by measuring the effects of perturbating the model parameters about baseline values. For each threshold-parameter pair, a normalized sensitivity is computed by linear regression. Top, sensitivities for the FasL thresholds ; bottom, sensitivities for the corresponding Fas thresholds at FasL concentrations , respectively.
To probe this further, we sampled parameters with increasing spread about baseline values and computed the fraction of parameter sets that remained bistable (Methods). The results show that has an exponential form (Figure 8). Extrapolating to , the data suggest an asymptotic bistable fraction of . Hence, robustness of bistability remains substantial even under significant parameter variation.
The fraction of parameter sets that exhibit bistability as a function of the sampling variability follows the exponential form , where is the asymptotic bistable fraction. The fitted value of suggests that this robustness remains substantial even as .
Cell-level cluster integration
Thus far, we have considered only the apoptotic activation of an individual cluster. To obtain the more biologically relevant cell-level activation, we must integrate over all clusters. In principle, this integration should account for intercluster transport as well as any intrinsic differences between clusters, e.g., as due to spatial inhomogeneities. Here, however, we provide as demonstration only a very simple integration scheme. Specifically, we assumed that clusters are identical (apart from their parameter values, which are drawn randomly) and independent, and that FasL is homogeneous over the cell membrane. Then we can express the normalized cell activation as(6)where the subscript denotes reference to cluster .
A characteristic cell-level hysteresis curve is shown in Figure 9. As is immediately evident, such integration is a smoothing operator, averaging over the sharp thresholds of each cluster. Thus, the cell-level signal may be graded even though its constituents are not. Note, however, that the lack of a sudden switch from low to high Fas signaling does not necessarily imply the same at the level of the caspases which ultimately govern cell death, as downstream components may possess switching behaviors , , , , .
The apoptotic signals of all Fas clusters are integrated to produce a normalized cell activation . The resulting hysteresis curve on as a function of the FasL concentration is graded due to the heterogeneity of the bistability thresholds across the clusters (top). Despite this variability, a strong linear dependence persists between (bottom; the valid region is shown in green).
Finally, we sought to outline protocols to experimentally discriminate our model against the prevailing crosslinking model , which we considered in a slightly simplified form . To be precise, the crosslinking model that we used has the reactions(7a)(7b)(7c)where is FasL, is Fas, and is the complex FasL∶Fas for . With(8a)(8b)(8c)(8d)(8e)(continuing the notational convention that lowercase letters denote the concentrations of their uppercase counterparts), the steady-state solution under mass-action dynamics is(9a)(9b)(9c)where(10a)(10b)and(11a)(11b)are the total ligand and receptor concentrations, respectively. In analogy with our proposed model, hereafter called the cluster model, we used(12)as a measure of the apoptotic signal.
Clearly, the crosslinking model has only one steady state, while the cluster model is capable of bistability. This hence provides a ready discrimination criterion. Although tracing out the associated hysteresis curve may be problematic, we can nevertheless probe for bistability by using hyperactive mutants, e.g., the mutation of Ile 313 to Asp in Fas, which stabilizes the open conformation and enhances apoptotic activity .
Specifically, we considered an experimental setup in which the concentrations of FasL and Fas, both wildtype and mutant, can be controlled, and in which the apoptotic signal can be measured, e.g., through the degree of FADD binding or of caspase activation. Hence we can map out the response curves at various levels of mutant penetrance. Denoting mutant Fas by (we assumed that mutant Fas cannot close, so there is no distinction between the stable and unstable open forms), we characterized the mutant penetrance by the mutant population fraction , where is the nondimensional mutant Fas concentration and is the total receptor concentration, composed of contributions from both wildtype () and mutant () forms. We assumed no other functional differences between wildtype and mutant Fas.
Proceeding first for the crosslinking model, at fixed , the amount of Fas bound by FasL is determined only by . Hence we assumed that a fraction of all receptors is bound. Since wildtype and mutant Fas are functionally identical by assumption, the fraction bound for each of the wildtype and mutant populations is also by independence of the recruitment process. Therefore, under the crosslinking model, varying at fixed yields an invariant response curve for the active wildtype Fas fraction .
In contrast, for the cluster model, we expected mutant receptor cluster-interactions to affect the wildtype response. Accordingly, the reactions (1c) and (1d) were amended for interaction with by replacing with(13a)(13b)respectively. This gives the analogue(14)of (5).
As seen in Figure 10, receptor interactions indeed cause the apoptotic signal to increase with even after accounting for the effect of mutants. This is because mutants can activate wildtype receptors by pushing the cluster past its switching threshold. Furthermore, the convergence to the active cluster state at high provides evidence for bistability. Thus, the variance of the -response curve at various can be used for model discrimination.
The wildtype response curve, giving the steady-state active wildtype Fas fraction as a function of the FasL concentration (stable, solid lines; unstable, dashed lines), of the cluster model varies with the mutant population fraction , reflecting receptor interactions absent in the crosslinking model. The total receptor concentration is fixed at . All parameters set at baseline values unless otherwise noted.
Alternatively, if working with mutants should prove difficult, we provide also a discrimination test based on steady-state invariants, i.e., functions that vanish at steady state. Clearly, for each model, provides a steady-state invariant since necessarily at steady state. However, the difficulty lies in expressing solely in terms of variables that are experimentally accessible. For example, current technology may not allow the concentrations to be measured accurately, if at all. Therefore, all such variables must be eliminated. Rate constants were considered parameters and so were not subject to this rule.
We assumed the same experimental setup as above and hence expressed each model invariant in terms of , , and , giving functions of the form , where encompasses all model parameters (Methods). The task then is to use , with provided by experiment, to assess the fit of a model. However, the parameters remain unknown, so this assessment cannot proceed directly. Instead, we considered the best possible fit over all parameters. A high value of indicates a poor best-case fit and hence that a model is unlikely to be correct. Clearly, prior knowledge of can be used to guide the invariant to biologically plausible fits.
To demonstrate that model discrimination using steady-state invariants is practical, we generated synthetic data from each model, calculating the accessible concentrations for each parameter set. This gives two sets of model-generated data. For each data set, we computed the best-fit invariant error for each model, where is the root mean square operator. The results suggest that this test can correctly identify the model from the data (Figure 11).
Steady-state invariants are fit to synthetic data generated from each model. For each model-data pair, the invariant error is minimized over the model parameter space. The results suggest that invariant minimization can correctly identify the model from the data.
The systems that we have presently considered are simple enough that experimentally inaccessible variables can be eliminated by hand. For more complicated systems, the tools of computational algebraic geometry, notably Gröbner bases, may prove useful; for such an application, see .
In this work, we showed through analysis of a mathematical model that receptor clustering can support bistability and hysteresis in apoptosis through a higher-order analogue of biologically observed Fas pair-stabilization . Hence we add to the signal processing activities in which receptor clustering has been suggested to participate –. This bistability plays an important functional role by enabling robust threshold switching between life and death states. Significantly, our results indicate potential key roles for ligand trimerism  and receptor pre-association – and localization onto membrane lipid rafts . Thus, we provide novel interpretations for these phenomena within the unified context of bistability.
Our model suggests an additional cell death decision, supplementing those that have been studied previously , , , , . Critically, the proposed decision is implemented upstream at the very death receptors that initially detect the death signal encoded by FasL. This decision is therefore apical in that it precedes all others in the system. Consequently, it operates independently of all intracellular components and so offers a general mechanism for bistability, even in cell lines with, for example, only feedforward caspase-activation networks , , , . Thus, receptor cluster-activation may explain how an effective apoptotic decision is implemented in such cells. Moreover, this suggests a novel target for induced cell termination in the treatment of disease .
We believe that our model provides an attractive theory for the observed biology. Although unlikely to be correct in mechanistic detail, the model may nevertheless reflect reality at a qualitative level. The significance of our work hence lies in its capacity to guide future research. We therefore readily invite experiment, which can reveal the true nature of the molecular mechanisms involved. Given their structural and functional homology, similar investigations on other members of the tumor necrosis factor receptor family may also prove fruitful. Such work serves to further our understanding of the formation and mode of action of complex signaling platforms such as the DISC, which in this view may be considered the macromolecular aggregates of active Fas.
The rationale for the choices is presented in the text; here, we further defend these by noting that no new behaviors are introduced with or . The remaining parameter values were guided by the following considerations. Specifically, we required and due to the assumed stabilities of the receptor species; all other parameters were assumed to be close to . Within these constraints, parameters were selected to ensure that is of the correct order of magnitude , , , . The baseline parameter values used were , , , , , , and .
To analyze the effects of variability in the model parameters, parameter values were sampled from a log-normal distribution, characterized by a variation coefficient , defined as the ratio of the standard deviation to the median of the distribution. For the sensitivity analysis, parameter sets were drawn at ; for the robustness analysis, parameter sets were drawn over ; and for the cell-level integration, parameter sets were drawn at . All parameters were drawn about baseline median values.
For each threshold-parameter pair, linear regression was performed on the threshold data against the parameter data, each normalized by reference values. For parameters, the reference is the baseline (median) value; for thresholds, the reference is the threshold computed at baseline parameters. The normalized sensitivity was defined as the slope of the linear regression.
The cluster invariant was derived by considering at steady state, i.e., with , and identifying . Similarly, the crosslinking invariant was derived by considering at steady state and identifying . For the full forms of the invariants, see Protocol S1.
For the model discrimination computation, parameter sets were drawn from log-normal distributions with median at . The active Fas concentration was calculated for each parameter set for each model at baseline parameters ( for the crosslinking model); for the cluster model, if bistability was observed, one of the stable values of was chosen at random. The invariant error was minimized using SLSQP with a lower bound of for all parameters.
All calculations were performed with Sage 4.5 , using NumPy/SciPy  for numerical computation and matplotlib  for data visualization. The Sage worksheet containing all computations is provided in the Supporting Information (Protocol S1) and can also be downloaded from  or http://www.courant.nyu.edu/~ho/.
(5.44 MB GZ)
We thank Leslie Greengard for useful discussions and for facilitating our research. We also thank the anonymous reviewers for their very helpful comments and suggestions.
Conceived and designed the experiments: KLH HAH. Performed the experiments: KLH HAH. Analyzed the data: KLH HAH. Wrote the paper: KLH HAH.
- 1. Thompson CB (1995) Apoptosis in the pathogenesis and treatment of disease. Science 267: 1456–1462.
- 2. Raff M (1998) Cell suicide for beginners. Nature 396: 119–122.
- 3. Meier P, Finch A, Evan G (2000) Apoptosis in development. Nature 407: 796–801.
- 4. Fulda S, Debatin KM (2006) Extrinsic versus intrinsic apoptosis pathways in anticancer chemotherapy. Oncogene 25: 4798–4811.
- 5. Taylor RC, Cullen SP, Martin SJ (2008) Apoptosis: controlled demolition at the cellular level. Nat Rev Mol Cell Biol 9: 231–241.
- 6. Fussenegger M, Bailey JE, Varner J (2000) A mathematical model of caspase function in apoptosis. Nat Biotechnol 18: 768–774.
- 7. Bentele M, Lavrik I, Ulrich M, Stößer S, Heermann D, et al. (2004) Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis. J Cell Biol 166: 839–851.
- 8. Eißing T, Conzelmann H, Gilles ED, Allgöwer F, Bullinger E, et al. (2004) Bistability analyses of a caspase activation model for receptor-induced apoptosis. J Biol Chem 279: 36892–36897.
- 9. 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.
- 10. Bagci EZ, Vodovotz Y, Billiar TR, Ermentrout GB, Bahar I (2006) Bistability in apoptosis: Roles of Bax, Bcl-2, and mitochondrial permeability transition pores. Biophys J 90: 1546–1559.
- 11. Legewie S, Blüthgen N, Herzel H (2006) Mathematical modeling identifies inhibitors of apoptosis as mediators of positive feedback and bistability. PLoS Comput Biol 2: e120.
- 12. 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.
- 13. 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: e299.
- 14. Okazaki N, Asano R, Kinoshita T, Chuman H (2008) Simple computational models of type I/type II cells in Fas signaling-induced apoptosis. J Theor Biol 250: 621–633.
- 15. Budihardjo I, Oliver H, Lutter M, Luo X, Wang X (1999) Biochemical pathways of caspase activation during apoptosis. Annu Rev Cell Dev Biol 15: 269–290.
- 16. Danial NN, Korsmeyer SJ (2004) Cell death: Critical control points. Cell 116: 205–219.
- 17. Ashkenazi A, Dixit VM (1998) Death receptors: Signaling and modulation. Science 281: 1305–1308.
- 18. Peter ME, Krammer PH (1998) Mechanisms of CD95 (APO-1/Fas)-mediated apoptosis. Curr Opin Immunol 10: 545–551.
- 20. Nicholson DW, Thornberry NA (1997) Caspases: killer proteases. Trends Biochem Sci 22: 299–306.
- 21. Nuñez G, Benedict MA, Hu Y, Inohara N (1998) Caspases: the proteases of the apoptotic pathway. Oncogene 17: 3237–3245.
- 22. Thornberry NA, Lazebnik Y (1998) Caspases: Enemies within. Science 281: 1312–1316.
- 23. Nicholson DW (1999) Caspase structure, proteolytic substrates, and function during apoptotic cell death. Cell Death Differ 6: 1028–1042.
- 24. Kitano H (2004) Biological robustness. Nat Rev Genet 5: 826–837.
- 25. 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.
- 26. Lai R, Jackson TL (2004) A mathematical model of receptor-mediated apoptosis: dying to know why FasL is a trimer. Math Biosci Eng 1: 325–328.
- 27. Scott FL, Stec B, Pop C, Dobaczewska MK, Lee JJ, et al. (2009) The Fas-FADD death domain complex structure unravels signalling by receptor clustering. Nature 457: 1019–1022.
- 28. Ising E (1925) Beitrag zur theorie des ferromagnetismus. Z Phys 31: 253–258.
- 29. Huang B, Eberstadt M, Olejniczak ET, Meadows RP, Fesik SW (1996) NMR structure and mutagenesis of the Fas (APO-1/CD95) death domain. Nature 384: 638–641.
- 30. Chan FKM, Chun HJ, Zheng L, Siegel RM, Bui KL, et al. (2000) A domain in TNF receptors that mediates ligand-independent receptor assembly and signaling. Science 288: 2351–2354.
- 31. Siegel RM, Frederiksen JK, Zacharias DA, Chan FKM, Johnson M, et al. (2000) Fas preassociation required for apoptosis signaling and dominant inhibition by pathogenic mutations. Science 288: 2354–2357.
- 32. Chan FKM (2007) Three is better than one: Pre-ligand receptor assembly in the regulation of TNF receptor signaling. Cytokine 37: 101–107.
- 33. Muppidi JR, Siegel RM (2004) Ligand-independent redistribution of Fas (CD95) into lipid rafts mediates clonotypic T cell death. Nat Immunol 5: 182–189.
- 34. Arnol'd VI (1992) Catastrophe Theory. Berlin, Germany: Springer-Verlag, 3rd edition.
- 35. Harrington HA, Ho KL, Ghosh S, Tung K (2008) Construction and analysis of a modular model of caspase activation in apoptosis. Theor Biol Med Model 5: 26.
- 36. Manrai AK, Gunawardena J (2008) The geometry of multisite phosphorylation. Biophys J 95: 5533–5543.
- 37. Bray D, Levin MD, Morton-Firth CJ (1998) Receptor clustering as a cellular mechanism to control sensitivity. Nature 393: 85–88.
- 38. Sourjik V (2004) Receptor clustering and signal processing in E. coli chemotaxis. Trends Microbiol 12: 569–576.
- 39. Endres RG, Oleksiuk O, Hansen CH, Meir Y, Sourjik V, et al. (2008) Variable sizes of Escherichia coli chemoreceptor signaling teams. Mol Syst Biol 4: 211.
- 40. Stein WA (2008) Can we create a viable free open source alternative to Magma, Maple, Mathematica and Matlab? pp. 5–6. doi:10.1145/1390768.1390771. Available: .
- 41. Oliphant TE (2007) Python for scientific computing. Comput Sci Eng 9: 10–20.