Bistability in Apoptosis by Receptor Clustering

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.


Introduction
Apoptosis is a coordinated cell death program employed by multicellular organisms that plays a central role in many physiological processes. Normal function of apoptosis is critical for development, tissue homeostasis, cell termination, and immune response, and its disruption is associated with pathological conditions such as developmental defects, neurodegenerative disorders, autoimmune disorders, and tumorigenesis [1][2][3][4][5]. Due to its biological significance, much effort has been devoted to uncovering the pathways governing apoptosis. Indeed, recent progress has enabled the proliferation of mathematical models, both mechanistic and integrative [e.g., [6][7][8][9][10][11][12][13][14], which together have offered profound insights into the underlying molecular interactions. The current work takes a similarly mathematical approach and hence inherits from this legacy.
There are two main pathways of apoptotic activation: the extrinsic (receptor-mediated) pathway and the intrinsic (mitochondrial) pathway, both of which are highly regulated [15,16]. In this study, we focus on the core machinery of the extrinsic pathway, which is initiated upon detection of an extracellular death signal, e.g., FasL, a homotrimeric ligand that binds to its cognate transmembrane death receptor, Fas (CD95/Apo-1), in a 1:3 ratio. This clusters the intracellular receptor death domains and promotes the ligation of FADD, forming the deathinducing signaling complex (DISC) [17][18][19]. The DISC catalyzes the activation of initiator caspases, e.g., caspase-8, through death effector domain interactions. Initiator caspases then activate effector caspases, e.g., caspase-3, which ultimately execute cell death by direct cleavage of cellular targets [20][21][22][23].
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 [24]. Consequently, researchers have used computational models to identify and study potential sources of bistability in apoptosis, including positive caspase feedback [8], inhibition of DISC by cFLIP [7], cooperativity in apoptosome formation [10], double-negative caspase feedback through XIAP [11], and double-negative feedback in Bcl-2 protein interactions [25]. 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 [26]. However, recent structural data [27] 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 pairstabilize through stem helix interactions. This affords a mechanism for bistability, similar to the Ising model in ferromagnetism [28], where open Fas, presumably disfavored relative to their native closed forms [29], 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 [30][31][32] or localization onto lipid rafts [33], 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.

Model formulation
Constructing a mathematical model of Fas dynamics is not entirely straightforward as receptors can form highly oligomeric clusters [27,33]. 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 (L, X , Y , Z), where L represents FasL and X , Y , and Z 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 LzjY 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 m and n, which capture the effects of receptor density and Fas coordination by FasL, respectively. Although only pairstabilization (m~n~2) has been observed experimentally [27], higher-order analogues, for example, as facilitated by globular interactions, are not unreasonable.

Author Summary
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 thresholdlike 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. (1d) 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 k (i) s and k (i) l for all ligand-independent and -dependent cluster-stabilization reactions of molecularity i, 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 (l, x, y, z) of (L, X , Y , Z). 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 where s is a characteristic concentration and t is time, and this is where s~jzgzf is the nondimensional total receptor density, and f ? is given by considering and solving df=dt~0 with (j, g, f).(j ? , g ? , f ? ), a polynomial in f ? of degree maxfm, ng. Clearly, the model is bistable only if maxfm, ng §3 (two stable nodes must be separated by an unstable node as the model is effectively one-dimensional in f).
We used f as a measure of the apoptotic activation of a cluster. In principle, all open receptors contribute to apoptotic signaling, but g is small, at least at steady state (since k o %1 due to the assumed prevalence of the closed form [29]), and so can be neglected.

Bistability and receptor clustering
While n measures the coordination capacity of FasL and hence may be equated with its oligomeric order (e.g., n~3 in the biological context), an appropriate value for m, 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 *10 mm, the associated volume of *1 pL implies the correspondence 1 nM *600 molecules *10 {6 molecules/nm 2 on restricting to the membrane, i.e., by averaging over the surface area of *600 mm 2 . Thus, for a conservative receptor concentration estimate of 100 nM [7,9,12,13], the number of Fas molecules in the neighborhood of each receptor is only *1, assuming a charateristic size of 100 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 (m~2). Therefore, in this context, bistability is possible only if n §3, and the trimerism of FasL thus demonstrates the lowest-order complexity required for bistability.
From the form of df=dt, this bistability is reversible as a function of the FasL concentration l since the governing polynomial for f ? is of degree only m~2 at l~0. 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 [30][31][32] and to selectively localize onto membrane lipid rafts [33]. The result of either of these processes may be to increase the local receptor concentration. In this high density mode, we set m §3, 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 m~3, using s as a continuous transition parameter. Furthermore, we set n~3 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 wellseparated 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.
As per the previous analysis, reversibility of the bistability is dependent on s, with irreversibility emerging for s sufficiently high. This suggests a bivariate parameterization of f ? , producing a multivalued steady-state surface over (l, s)-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 [34]. A more instructive view of the dependence of the model's qualitative structure on l and s is shown in Figure 5.

Sensitivity and robustness analyses
We then focused on the activation and deactivation thresholds l + , 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 l at which the hysteresis curve turns, i.e., at Ll=Lf ?~0 (Figure 6). We performed a sensitivity analysis of l + by measuring the effects of perturbing the model parameters about baseline values (Methods). For each threshold-parameter pair, we computed a normalized sensitivity S by linear regression.
Strong effects of s, k o , and k u were observed ( Figure 7); for the corresponding Fas thresholds f + at l~l + , respectively, the parameters s, k o , k (2) l , and k (3) l 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.
To probe this further, we sampled parameters with increasing spread D about baseline values and computed the fraction f of parameter sets that remained bistable (Methods). The results show that f has an exponential form (Figure 8). Extrapolating to D??, the data suggest an asymptotic bistable fraction of f ? &0:4. Hence, robustness of bistability remains substantial even under significant parameter variation.

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 where the subscript i denotes reference to cluster i. 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 [7,8,10,11,25].

Model discrimination
Finally, we sought to outline protocols to experimentally discriminate our model against the prevailing crosslinking model [26], which we considered in a slightly simplified form [35]. To be precise, the crosslinking model that we used has the reactions LzR 3k z k{ C 1 , ð7aÞ where L is FasL, R is Fas, and C i is the complex FasL:Fas i for i~1, 2, 3. With (continuing the notational convention that lowercase letters denote the concentrations of their uppercase counterparts), the steadystate solution under mass-action dynamics is where and L~lzc 1 zc 2 zc 3 , ð11aÞ are the total ligand and receptor concentrations, respectively. In analogy with our proposed model, hereafter called the cluster model, we used as a measure of the apoptotic signal. Hyperactive mutants. 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 [27].
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 Z D (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 D~f D = s s, where f D~zD =s is the nondimensional mutant Fas concentration and s s~szf D is the total receptor concentration, composed of contributions from both wildtype (s) and mutant (f D ) forms. We assumed no other functional differences between wildtype and mutant Fas.
Proceeding first for the crosslinking model, at fixed s s, the amount of Fas bound by FasL is determined only by L. Hence we assumed that a fraction~(L) 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 (L) by independence of the recruitment process. Therefore, under the crosslinking model, varying D at fixed s s yields an invariant response curve for the active wildtype Fas fraction~f ? =s.
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 Z D by replacing with i~2, . . . ,n, j~1, . . . ,i, k~0, . . . ,i{j, k'~1, . . . ,j, respectively. This gives the analogue of (5). As seen in Figure 10, receptor interactions indeed cause the apoptotic signal to increase with D 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 l provides evidence for bistability. Thus, the variance of theresponse curve at various D can be used for model discrimination.
Steady-state invariants. 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, _ f f:df=dt provides a steadystate invariant since _ f f~0 necessarily at steady state. However, the difficulty lies in expressing _ f f solely in terms of variables that are experimentally accessible. For example, current technology may not allow the concentrations c i 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 L, s, and f ? , giving functions of the form _ f f(L, s, f ? ; a), where a encompasses all model parameters (Methods). The task then is to use D _ f fD, with (L, s, f ? ) provided by experiment, to assess the fit of a model. However, the parameters a remain unknown, so this assessment cannot proceed directly. Instead, we considered the best possible fit min a D _ f fD over all parameters. A high value of min a D _ f fD indicates a poor best-case fit and hence that a model is unlikely to be correct. Clearly, prior knowledge of a 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 (L, s, f ? ) for each parameter set. This gives two sets of model-generated data. For each data set, we computed the best-fit invariant error E~min a rmsD _ f fD for each model, where rms is the root mean square operator. The results suggest that this test can correctly identify the model from the data (Figure 11).
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 [36].

Discussion
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 [27]. Hence we add to the signal processing activities in which receptor clustering has been suggested to participate [37][38][39]. 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 [17] and receptor preassociation [30][31][32] and localization onto membrane lipid rafts [33]. 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 [7,8,10,11,25]. 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 [7,9,13,14]. Thus, receptor clusteractivation 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 [1].
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.

Parameter selection
The rationale for the choices m~n~3 is presented in the text; here, we further defend these by noting that no new behaviors are introduced with m or nw3. The remaining parameter values were guided by the following considerations. Specifically, we required k o and k u %1 due to the assumed stabilities of the receptor species; all other parameters were assumed to be close to O(1). Within these constraints, parameters were selected to ensure that l z is of the correct order of magnitude [7,9,12,13]. The baseline parameter values used were s~1, k o~2 |10 {3 , k u~1 0 {3 , k (2) s~0 :1, k (3) s~0 :5, k (2) l~1 , and k (3) l~5 .

Parameter sampling
To analyze the effects of variability in the model parameters, parameter values were sampled from a log-normal distribution, characterized by a variation coefficient D, defined as the ratio of the standard deviation to the median of the distribution. For the sensitivity analysis, 500 parameter sets were drawn at D~0:25; for the robustness analysis, 250 parameter sets were drawn over 0ƒDƒ5; and for the cell-level integration, 100 parameter sets were drawn at D~0:25. All parameters were drawn about baseline median values.

Sensitivity analysis
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 S was defined as the slope of the linear regression.

Steady-state invariants
The cluster invariant was derived by considering _ f f at steady state, i.e., with (j, g).(j ? , g ? ), and identifying l.L. Similarly, the crosslinking invariant was derived by considering _ f f at steady state and identifying r ? .s{f ? . For the full forms of the invariants, see Protocol S1.
For the model discrimination computation, 100 parameter sets (L,s) were drawn from log-normal distributions with median (0:2,1) at D~(0:2,0:5). The active Fas concentration f ? was calculated for each parameter set for each model at baseline parameters (k~0:1 for the crosslinking model); for the cluster model, if bistability was observed, one of the stable values of f ? was chosen at random. The invariant error E was minimized using SLSQP with a lower bound of 10 {3 for all parameters.