Bistability Analyses of a Caspase Activation Model for Receptor-induced Apoptosis *

From the ‡Institute for Systems Theory in Engineering, University of Stuttgart, Pfaffenwaldring 9, 70550 Stuttgart, Germany, the ¶Institute for System Dynamics and Control, University of Stuttgart, Pfaffenwaldring 9, 70550 Stuttgart, Germany, the Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstr. 1, 39106 Magdeburg, Germany, and the **Institute for Cell Biology and Immunology, University of Stuttgart, Allmandring 31, 70569 Stuttgart, Germany


INTRODUCTION
Apoptosis is a genetically defined, major form of programmed cell death enabling the organism to remove unwanted cells, e.g. during embryonal development and after immune responses, to select educated immune cells and to eliminate virally infected and transformed cells (1,2). Enhanced or inhibited apoptotic cell death can be involved in severe pathological alterations, including developmental defects, autoimmune diseases, neurodegeneration or cancer. Extrinsic and intrinsic apoptotic pathways can be distinguished, although partly employing overlapping signal transduction pathways. A hallmark of the ongoing apoptotic process is the activation of a family of aspartate directed cysteine proteases, the caspases.
Caspases are produced as proenzymes and become activated upon cleavage (3). Activation of caspases finally dismantles the cells via the cleavage of important regulatory and structural proteins and enables phagocytic removal of the dying cell (4). A simplified outline of the extrinsic pathway of apoptosis induction after death receptor stimulation is depicted in Figure 1.
Mathematical modeling and systems theory can provide valuable tools to get insight into complex dynamical systems, to test hypotheses and to identify weak points (5,6).
Previous modeling approaches in apoptosis focused on the extrinsically triggered pathways, resulting in complex models (7,8). The model parameters were fitted to data derived from cell population studies showing caspase activation in the range from 30 min to several hours.
These models can describe and nicely illustrate certain aspects of the signal transduction pathway. However, more recent experimental results performed at the single cell level show that the majority of caspases is activated within a very short time interval (< 15 min) (9)(10)(11)(12). represents the main signaling route in so called type I cells (13) (Fig. 1, yellow background).
We therefore translated the current picture of the extrinsically triggered caspase cascade in a very elementary form into a mathematical model enabling a thorough investigation through the application of analytical methods. Our results show that within large parameter ranges, including values from literature, this straightforward model structure is unable to appropriately describe the expected behavior that can be deduced from experimental data. We then show a way of extending our model structure to reconcile these observed differences and present a model now able to describe key characteristics like a fast execution phase and bistability.
In addition, results from our model studies show a way to reconcile the fast kinetics of caspase 3 activation observed at the single cell level with the much slower kinetics found at the level of a cell population in terms of understanding and modeling.

EXPERIMENTAL PROCEDURES
The mathematical model. For each reaction considered in the Results one reaction rate can be deduced (v 1 -v 13 ). The cleavage reactions (1, 2 and 4) are treated as being irreversible and it is assumed that the intermediary cleavage products ("enzyme-substrate complexes") only achieve very low levels and can thus be eliminated, reasonable estimations that have been confirmed by simulation experiments (data not shown). The reaction rate equations are deduced according to the law of mass action, which we here consider prior to other kinetic approaches, like Michaelis-Menten kinetics, although theoretical considerations show that the results would be very similar in our case. From this, molecular balances can be derived for each considered molecular species resulting in a system of ordinary differential equations (ODEs; Eq. 1-8). Two different model structures are considered in greater detail during this study. The basic model includes the reactions 1-10 and translates into the equations 1-6 (not including v 11 The models were implemented in both Matlab (for simulation experiments) and Mathematica (for analytical analysis).  (14,15). Estimating a cell volume of 1 pl yields that 600 molecules per cell equals 1 nM. Accordingly, these values are roughly in the same order of magnitude and were used as initial concentrations. The other compounds were considered not to be present in the absence of a stimulus. In the extended model, the concentration of the newly introduced molecule BAR was assumed to be 40,000 molecules/cell. We consider the unit molecules/cell more illustrative for cellular concentrations than the unit M, but on the other hand we prefer and use units such as M -1 s -1 for the K m /k cat ratios. Table 1 lists the parameters as used in the "single set" simulations (unless indicated otherwise). The respective values are also provided in more common units (in brackets). For the reactions 3 and 5-10 the parameter values were taken from literature as stated in the text.
The respective references are summarized in Table 2 The Hurwitz criterion provides conditions for stability based on the coefficients of the characteristic polynomial. The most restrictive for the basic model is that the coefficient c (below) is positive, which was also used to construct the red area shown in Figure 3 (transcritical bifurcation manifold). The stability of the steady states other than the life steady state were evaluated numerically.
Deriving an input distribution. We assume a population behavior as depicted in Figure 5A, which can be interpreted as a probability distribution. We further assume that 100 % of caspase activation corresponds to 100 % cell death, i.e. caspase 3 becomes significantly activated in every cell within the population. Further, based on the simulations of the deterministic single cell model described in Figure 4, we can describe the maximal caspase 3 activation as a function of C8* input. Hereby we assume the maximal caspase activation to define the time point of cell death. This correlates the stochastic time point of cell death to a stochastic input signal for the single cells within a population. From the original distribution of Figure 5A we thus obtain a distribution of cell death probability as a function of input activation. The corresponding probability density function can be derived by differentiation as shown in Figure 5B.

RESULTS
The biology of the model system. The type I cell-like (13) model (called basic model in the following) was constructed with the purpose of being as uncomplicated as possible without neglecting essential steps concerning our analyses (see below), i.e. simplifications represent conservative estimations. As a model input we use a pulse of activated caspase 8, which is produced by the death inducing signaling complex (DISC) formed at the membrane after death receptor stimulation (16) (although the initial steps seem to be more complex in the  (17)). The model is outlined in Figure 1 (yellow background) and contains the following reactions: C3* + IAP ↔ iC3*~IAP [3] C3* + IAP → C3* [4] degradation & turnover reactions [5][6][7][8][9][10]  "noise") (28). Also, caspases are known to possess zymogenicity (18) and partial activation is observed in some physiological processes (29). However, if the apoptotic initiation signal is beyond a certain threshold, the cell must irreversibly enter the pathway to develop apoptosis.
In the following, we use this information in a reverse engineering manner (6) and take bistability as a "conditio sine qua non" to evaluate possible models with respect to this expected behavior.  Figure 2). Figure 3 shows the bifurcation point (red area) in dependence of three parameter classes with fixed ratios in each class. The parameter combination that can be deduced from literature (Fig. 3, yellow dot) is far away from those combinations enabling a stable life steady state (below the red area in Fig. 3 Bistability within a small parameter domain. The two additional steady states besides the life steady state provide further information. Theoretically, one additional steady state within the positive concentration range is sufficient to achieve bistability (because other phenomena than an unstable steady state could separate the two stable steady states). However, in our model this configuration is only possible above the red area (Fig. 3, see Experimental Procedures) and for that setting the life steady state is unstable. Accordingly, bistability is only possible if both additional steady states are within the positive concentration range (above the blue area in Fig. 3). Another restriction is imposed by our biological system as the solutions must contain real numbers (above the green area in Fig. 3). Thus, bistability is only possible in a very restricted parameter area far away from values reported in literature (Fig. 3, yellow dot). In accordance with general considerations on this topic (25,30) the feedback from caspase 3 onto caspase 8 is necessary for bistability (data not shown).
IAPs and their cleavage. Interestingly, by and large the stability of the life steady state seems to be independent of the IAP cleavage reaction (Fig. 3). This is not expected and Further analysis of the model reveals that apoptosis can only proceed after the IAP pool is exhausted, as otherwise most of the active caspase 3 molecules become neutralized (data not shown). As the binding of active caspase 3 to IAPs is a reversible reaction, a slower degradation of the complexes would elevate the levels of free active caspase 3 and therefore promote apoptosis. Thus, our results also argue for the view of IAPs as altruistic proteins, sacrificing themselves to prevent cell death (31). In our investigations we did not change the initial concentrations in order to restrict the number of free parameters. However, in the mathematical formulas it can be seen that changing a parameter value has similar effects as changing an initial concentration, and so we have indirectly evaluated these influences as well.
Bistability through model extension. An inherent problem of the described basic model of caspase activation is that relatively fast activation kinetics must be realized in order to be consistent with parameter values from literature and observations in various experimental setups (9-12,32,33). On the other hand, only if the kinetics are slow, the system displays the desired bistability. There are several possibilities how a model extension could reconcile those facts as described in the Discussion. One possibility is a mechanism to inhibit active caspase 8. In fact, recent reports describe that activated caspase 8 can be functionally inactivated in mitochondrial membranes. There, the molecule BAR (bifunctional apoptosis regulator) has been proposed to bind activated caspase 8 via its pseudo death effector domain leading to effective neutralization of its proteolytic activity (34)(35)(36)(37). We therefore introduced the molecule BAR into our model, and assumed it to bind to activated caspase 8 with an affinity similar to XIAP binding to activated caspase 3 (23,38). In addition, we introduced a normal turnover rate for the new compounds (see Experimental Procedures).
The  Table 1). We performed simulations with input signals of different strength (Fig. 4), using a set of parameters where the system displays a bistable behavior (Table 1) activation is typically observed within a range of a few hours (13,39,40). Figure 5A shows an example of a slow time course of effector caspase activation. As described above, the strength of the apoptotic stimulus mainly translates into the delay between the apoptotic trigger and significant effector caspase activation. Thus, different delay times due to the stochastic nature of biological systems enable small fast single-cell segments of caspase activation to integrate, forming an overall slow caspase activation at the macroscopic level. Figure 5B shows the distribution of input signals into the extended BAR model (as described above) required to produce a population behavior as depicted in Figure 5A

DISCUSSION
In the past, systems theory has been mainly used to understand dynamic processes of technical structures, but recently it has been increasingly applied in biology (5,6). Many articles focus on the mathematical modeling of aspects of the metabolism, where (quasi) steady state approximations can be used and the mere structure provides valuable information (41). The initiation and development of apoptosis, however, is a highly dynamic process where steady state approximations cannot be easily applied. We employ a reduced model in order to evaluate the principle behavior, but also take dynamic information into account using data from various sources trying to fuse them together with the help of our model. The class of biological systems exemplified here has previously not been evaluated with respect to bistability, although bistability can be shown for special classes of larger systems (30). Our results show that (partly) analytic approaches for additional classes are feasible for larger systems than widely believed.
The present study analyses the basic core reactions of caspase activation as they are widely accepted to occur after death receptor activation (2,3,42). Because the molecular events occurring at the level of the receptor induced signaling complex are only in part understood (13,17), the input signal of the model is represented by a given number of molecules of activated initiator caspase (Fig. 1), excluding the pathways how this initial caspase activation is mediated. This model, however, reveals no stable steady state in the life status, mainly caused by the strong forward reaction of the caspase cleavage reactions as such, containing a positive feedback loop. We therefore reduced further analyses to the reactions depicted in Figure 1 on the yellow background, including the caspase activation reactions. Bifurcation analyses of this basic model reveal that it is capable to display the desired bistable behavior only within a narrow range of kinetic parameter values that are far away from values reported in literature. We thus considered possible modifications in the model structure, capable to reconcile the observed fast kinetics and tolerance to subthreshold stimuli. A common way in cell biology to achieve such a behavior is cooperativity (25), which could also be assumed to exist for caspases, as they are known to act as lag phases in agreement with experimental data (9). Although additional explanations for these long lag phases have been proposed recently (at least for TNF induced apoptosis (17)), our studies demonstrate that caspase 8 capture/inactivation is a likely mechanistic explanation for the observed lag period followed by fast effector caspase activation. Whereas this lag phase might be largely controlled at the level of the mitochondria in type II cells (9,10), a similar behavior was also observed for type I cells, although initial activation of caspase 8 takes place immediately after stimulation (13,16).
In fact, in a very recent publication a new class of molecules has been described, termed caspase 8 and 10 associated RING proteins (CARPs), that are proposed to be novel caspase 8 and 10 inhibitors, acting similar to the IAP proteins (44). In the light of the results from our studies, these experimental data underscore the necessity for an effective regulation at the level of initiator caspases in death receptor mediated apoptosis. Regarding In experimental studies the activation of effector caspases both in type I as well as type II cell populations occurs gradually within a few hours (13,39,40), as schematically depicted in Figure 5A. However, at the single cell level, rapid caspase 3 activation has been observed after a lag phase in the very same experimental setup (9)(10)(11)(12). Using our extended model, we can easily reconcile the observed fast kinetics at the single cell level with the slower rates observed in cell populations. We show here that the variation of the input signal strength within single cells, as depicted in Figure 5B, results in slow kinetics at the population level as shown in Figure 5A.  parameter can only be estimated, see Table 2). Above the red area no stable life steady state exists, below the blue area no second stable steady state within the positive concentration range can exist and below the green area the solutions are complex numbers.
Thus, bistability requires parameter combinations below the red AND above the blue AND above the green area.