Modeling Reveals That Dynamic Regulation of c-FLIP Levels Determines Cell-to-Cell Distribution of CD95-mediated Apoptosis*

The expression levels of caspase-8 inhibitory c-FLIP proteins play an important role in regulating death receptor-mediated apoptosis, as their concentration at the moment when the death-inducing signaling complex (DISC) is formed determines the outcome of the DISC signal. Experimental studies have shown that c-FLIP proteins are subject to dynamic turnover and that their stability and expression levels can be rapidly altered. Even though the influence of c-FLIP on the apoptotic behavior of a single cell has been captured in mathematical simulation studies, the effect of c-FLIP turnover and stability has not been investigated. In this study, a mathematical model of apoptosis was developed to analyze how the dynamic turnover and stability of the c-FLIP isoforms regulate apoptotic signaling for both individual cells and cell populations. Intercellular parameter and concentration distributions were used to describe the behavior of cell populations. Monte-Carlo simulations of cell populations showed that c-FLIP turnover is a key determinant of death receptor responses. The fact that the developed model simulates the state of whole cell populations makes it possible to validate it by comparison with empirical data. The proposed modeling approach can be used to further determine limiting factors in the DISC signaling process.

Apoptosis, or programmed cell death, is a mechanism of active controlled cellular elimination, which, together with cell growth and division, regulates key processes such as embryonic development, immune responses, and tissue homeostasis. Upon stimulation, death receptors such as Fas/CD95/Apo-1 are activated and recruit cytoplasmic proteins, forming a deathinducing signaling complex (DISC). 4 In the DISC, the initiator caspase-8 can be activated by cleavage of procaspase-8, resulting in downstream activation of effector caspases and propaga-tion of the apoptotic signal. Caspase-8 activation is regulated by the c-FLIP (cellular FLICE inhibitory protein) inhibitory proteins, which exist in three splice variants, a long form (c-FLIP L ) and two short forms (c-FLIP S and c-FLIP R ). The DISC-interacting proteins procaspase-8 and c-FLIP S all contain two death effector domains, but c-FLIP S and c-FLIP R lack the longer protease-like domain of procaspase-8 and c-FLIP L (1). The c-FLIP isoforms have been shown to primarily delay or inhibit death receptor-mediated apoptosis (2), although the long c-FLIP variant is known to have a dual function and may also act as a caspase activator (3)(4)(5)(6). Recently, Fricker et al. (7) proposed an underlying mechanism for the dual role of c-FLIP L and developed a mathematical model, which was verified experimentally.
Experiments on human K562 erythroleukemia cells and WM35 melanoma cells (8) demonstrated that c-FLIP S is more apt to degradative ubiquitylation and consequent degradation than c-FLIP L . Furthermore, the anti-apoptotic activity of c-FLIP S and its ubiquitin-mediated degradation are functionally and structurally independent of each other (8). The destabilizing sequence in c-FLIP S implies distinct degradation rates of the c-FLIP isoforms, allowing differential regulation of their concentration levels.
It has also been found that heat shock sensitizes cells to CD95-mediated apoptosis in Jurkat and HeLa cells, independently of the expression of Hsp70 (9,10). The study demonstrated that a 30-min heat shock treatment is sufficient to decrease the levels of c-FLIP L and c-FLIP S and to sensitize the cells to apoptosis mediated by death receptors. Hyperthermia also sensitizes activated primary human T lymphocytes to CD95-induced apoptosis, primarily due to induced down-regulation of c-FLIP S (10). It was also demonstrated that the hyperthermia does not affect the synthesis of c-FLIP but rather rapidly enhances the proteasomal degradation of c-FLIP L and c-FLIP S .
Mathematical modeling and in silico analysis of biological systems have become important tools for simulation-based hypothesis testing (11)(12)(13)(14). The apoptotic signaling pathway has a complex dynamic behavior involving several reactions and feedback loops, and it may be difficult to infer its overall response by examining individual parts of the system separately. Mathematical modeling and simulation can therefore provide valuable insights into the behavior of apoptotic signaling pathways. A number of mathematical models of apoptosis have been proposed, including ordinary differential equation models (15)(16)(17)(18)(19)(20)(21)(22), stochastic models (18,23), and Petri net models (24,25). One of the most detailed differential equation approaches was presented by Bentele et al. (18). The model predicted CD95 signaling to be characterized by a threshold behavior, which depends on c-FLIP concentration in the DISC and the initial concentration of the death ligand CD95L (6).
Experimental studies (8 -10) indicate that c-FLIP proteins are subject to active ubiquitylation, thereby yielding high turnover rates. Consequently, their expression levels can change rapidly when the ubiquitylation and stability of c-FLIP are affected. In this study, mathematical modeling was used to mimic experimental studies of apoptosis in which c-FLIP levels are dynamically regulated. Synthesis and degradation of the isoforms c-FLIP S and c-FLIP L were introduced within the framework of the apoptosis model of Bentele et al. (18), and it is demonstrated that the dynamic regulation of c-FLIP plays an important role for apoptosis. A problem in verifying these findings is that experimental studies of apoptosis usually report the proportion of cell deaths during a given time interval. To overcome this limitation of single-cell models, we extended the model to take into account the fact that cells are not identical copies of each other (26) but may show differences in size, intercellular structure, and protein concentrations. We studied cell populations in which the reaction rate constants and initial protein concentrations in the various cells are statistically distributed. The statistical treatment and Monte Carlo simulations made it possible to compare the outcome of the simulations with results from empirical studies, in which the proportion of apoptotic cells is reported. This approach was also employed to study the effect of down-and up-regulation of c-FLIP under conditions similar to those used in the previously reported heat shock experiments (9,10). To mimic experimental studies of apoptosis in which c-FLIP levels are dynamically regulated and in which this regulation is a key factor for the outcome (6,8,10), the model was here modified to encompass dynamic c-FLIP turnover.
The results corroborate the role of c-FLIP-targeted protein degradation as a fundamental determinant of death receptor responses. A further asset of the developed model is the compatibility with cell population data, enabling convenient communication between empirical and modeled data. This feature was employed to validate the modeling results using empirical cell population data. The model can be used to study, apart from ubiquitylation-mediated degradation, other rate-limiting steps in DISC signaling such as the recently defined isoformspecific phosphorylation processes of c-FLIP (27,28) that partly regulate ubiquitylation and partly regulate other aspects.

MATERIALS AND METHODS
Mathematical Model-The mathematical modeling of apoptosis originated from the approach of Bentele et al. (18), in which the apoptotic signaling process is described using ordinary differential equations for 41 chemical components in 32 reactions. The reaction rates are modeled by mass action and Michaelis-Menten kinetics. The model exhibits a threshold mechanism in apoptotic behavior, which depends on the stationary c-FLIP concentration in the DISC and the initial con-centration of ligand. The model parameters were determined using experiments performed on the human B lymphoblastoid cell line SKW 6.4. To fit experimental data, the model incorporates both the long and short c-FLIP isoforms to account for their different reaction mechanisms. The two short c-FLIP isoforms, c-FLIP S and c-FLIP R , have been shown to affect caspase-8 activation and apoptosis sensitivity in a similar way (8,10,27,33), and they have similar half-lives. As the c-FLIP property that is important for this study is its half-life, and as the half-lives of the short forms of c-FLIP, i.e. c-FLIP S and c-FLIP R , are the same in all the model cell lines (8,10,27,33), we can assume that c-FLIP S and c-FLIP R will behave in the same way. Therefore, no distinction has been made between the two short isoforms in this study. Instead, both short forms are included as one factor, c-FLIP S/R . Besides its anti-apoptotic role, the long c-FLIP isoform can also act as an activator of caspase-8 (4,5). Experimental evidence supported by mathematical modeling shows that c-FLIP L has a pro-apoptotic role only at high ligand concentration or upon overexpression of one of the short c-FLIP isoforms (7). The pro-apoptotic roles of c-FLIP L have been shown only under conditions of ectopic overexpression. These conditions are, however, not relevant for this study, in which the focus is to investigate apoptosis in cell populations in which the endogenous c-FLIP levels are down-regulated due to increased turnover and/or proteasomal degradation. We have therefore used the mechanism presented in Ref. 18, in which the pro-apoptotic role of c-FLIP L is excluded.
Dynamic Regulation of c-FLIP Levels-In the mathematical model presented by Bentele et al. (18), the simplifying assumption was made that, upon stimulation, only the c-FLIP initially present in the cell takes part in the reactions, whereas the effect of synthesis and degradation of c-FLIP is not taken into account. However, the reported high turnover rate of c-FLIP (8,10) implies that c-FLIP bound in the DISC is rapidly replenished with newly synthesized c-FLIP, tending to counteract changes in c-FLIP levels. Moreover, a number of experimental reports (8,10,27) indicate that regulation of the stability and consequent concentration levels of c-FLIP plays an important role in controlling apoptotic sensitivity.
To describe the experimentally found high turnover rate of c-FLIP and the fact that the levels of c-FLIP S/R and c-FLIP L are dynamically regulated by affecting degradation rates, the rate of c-FLIP degradation was modeled as a first-order reaction (Equation 1), At steady state, the synthesis rates must equal the degradation rates, so it follows that the synthesis rates r syn S/R and r syn L are related to the steady-state concentrations [c-FLIP S/R ] ss and [c-FLIP L ] ss as in Equation 2.
gives a quantitative form for the results presented by Poukkula et al. (8), i.e. if the synthesis rate is constant, the expression level of c-FLIP can be efficiently regulated by affecting the stability of c-FLIP, defined by the degradation rate constant.
Time to Apoptosis-In the study by Bentele et al. (18), poly-(ADP-ribose) polymerase cleavage indicated the occurrence of apoptosis, and the time of poly(ADP-ribose) polymerase cleavage was used to determine relative times to apoptosis. It turns out, however, that this characterization of the time to apoptosis gives consistently too short times compared with experimental results (6,8,10) in which different approaches using flow cytometry had been used to determine the occurrence of apoptosis. Harrington et al. (22) used the maximum value of activated caspase-3 as an indicator of apoptosis. To introduce a definition that would better reflect the time at which apoptosis is detected, the time to apoptosis (t ap ) was taken as the time elapsed from exposure to ligand to the instant when the concentration of activated caspase-3 had decreased to half of its maximum value. In addition, the original model (18) was modified by assuming that apoptosis does not affect ligands in the cell exterior, receptors on the cell membrane, or ligand-receptor complexes.
Modeling of Cell Populations-For a given cell type, individual cells may show large variations in apoptotic sensitivity. Experimental studies of apoptosis therefore often report the proportions of cells that undergo apoptosis during given time intervals. In contrast, mathematical models of apoptosis based on differential equations describe the behavior of a single cell with specified kinetic parameters and molecular concentrations. Therefore, simulation results from differential equation models are not directly compatible with experimental data, making both model validation and use for prediction purposes difficult.
There are two major sources of random variations in cellular reactions. The first one is intrinsic noise (29) caused by the inherent stochastic nature of chemical reactions. The intrinsic noise causes large fluctuations when the numbers of the relevant molecules are small, but its effect decreases rapidly with increasing concentrations. Stochastic models, which take into account the intrinsic noise, describe the reactions at a molecular level (30). The second important source of randomness in cellular processes stems from variations between individual cells, which may differ both in the concentrations of various molecules and in the chemical reaction rates. Such variations can be described by modeling a cell population in which the parameters in the cells are randomly distributed. The behavior of the cell population can be studied using Monte Carlo simulations, in which the concentrations and kinetic parameters are selected randomly from the associated distributions. This approach has been applied to investigate robustness properties of biochemical networks (31,32), including apoptosis models (23).
To account for differences between the behavior of different cells in a population, we introduced distributions on initial concentrations, [X i ](0), and chemical reaction rate constants, k i . Randomly distributed parameters that take only positive values can be modeled using the log-normal distribution (23), in which the logarithm of a parameter is normally distributed. Also measurements have indicated log-normal distributions of concentrations in cells (26). We thus assumed the initial concentrations and the rate constants in the cells to be independent stochastic variables given by Equations 3 and 4, where N(0,) denotes a normally distributed stochastic variable with zero mean and variance , and [X i ] 0 and k i 0 are the nominal initial concentrations and rate factors, respectively, which were here taken as those given in the study by Bentele et al. (18). The model was programmed in MATLAB (The MathWorks, Inc., Natick, MA).

RESULTS
We used the modified apoptosis model to study how c-FLIP turnover and intercellular variability affect death receptor-mediated apoptosis. First, the effect of the synthesis and degradation rates of c-FLIP on the apoptotic process was analyzed. Then, intercellular variability was introduced to describe the apoptotic behavior of cell populations. Finally, the cell population model with dynamic regulation of c-FLIP levels was used to reproduce reported experimental results (10), studying recovery of a cell population that has been sensitized to apoptosis by a heat shock, followed by exposure to ligands after various time periods of recovery.
Effect of c-FLIP Turnover and Degradation Rate-The effect of c-FLIP turnover in the nominal (single-cell) model is shown in Fig. 1A, which depicts the concentrations of some key components in the apoptotic process for different turnover rates. The initial (steady-state) concentration of c-FLIP was taken to be the same as the initial concentration in the previously reported model (18) Fig. 1A. Although the c-FLIP levels are now better maintained by synthesis of new molecules, c-FLIP turnover has only a minor effect on the apoptotic behavior or on the time to apoptosis (t ap ϭ 144 min). To confirm that the apoptotic signaling depends mainly on the initial (steady-state) concentration of c-FLIP and not on the replacement of consumed c-FLIP by synthesis, the responses for infinite turnover rates (r syn S/R and r syn L 3 ∞) were simulated and are illustrated in the lower panels of Fig. 1A. These results show that, although the dynamic responses of the c-FLIP concentrations depend strongly on the turnover rates, they have little effect on the concentrations of caspases associated with the DISC or on the apoptotic process. The results are similar for other initial ligand and c-FLIP concentrations (data not shown). Hence, the apoptotic behavior of the model can be concluded to depend primarily on the concentration of c-FLIP during DISC formation, whereas the turnover rate has a much smaller impact.
For a constant synthesis rate, the steady-state concentration of c-FLIP (and hence, apoptotic sensitivity) is affected by the degradation rate of c-FLIP (cf. Equation 2). This is illustrated in Fig. 1B, which depicts the results of simulations in which the death receptors were activated by death ligands at 5 ng/ml (ligand/receptor ϭ 0.005) at t ϭ 200 min. In the upper panels of Fig. 1B, the c-FLIP turnover rates and levels are the ones reported by Meinander et al. (10) (Fig. 1A, middle panels), and no apoptosis is triggered. The middle panels of Fig. 1B show the effect of increased degradation rates of c-FLIP, e.g. caused by heat shock (10).  Fig. 1B show the response obtained when DISCbound c-FLIP was degraded 100 times faster than free c-FLIP. As expected, an increase in apoptotic sensitivity was obtained, and apoptosis now occurred 228 min after exposure to ligand. This and other simulation studies (data not shown) indicate that the fraction (DISC-bound versus free) of c-FLIP that is degraded does not have a major effect on apoptotic sensitivity but that the overall degradation is the determinant.
Thus, Fig. 1 indicates that the apoptotic behavior is determined primarily by the c-FLIP concentration at ligand exposure and that this concentration can be regulated by the degradation rate of c-FLIP. This is supported by a study (10) that concluded that the apoptotic sensitivity of the cells depends mainly on the c-FLIP levels at the time of receptor ligation. This conclusion was based on experiments in which no sensitization to apoptosis was observed when cells were exposed to CD95 ligand prior to heat shock. In this situation, the c-FLIP concentrations are at normal levels at the time of ligand stimulation, in contrast to the situation in which cells have been subjected to heat shock prior to ligand stimulation, when c-FLIP levels have been down-regulated due to proteasomal degradation (10). These findings (10) could be confirmed by simulations indicating that, when the heat shock is initiated after ligand stimulation, its effect on the apoptotic process is negligible.
In the model used here, the long and short c-FLIP forms affect the apoptotic process in similar ways, as the nominal reaction rates in the DISC are the same for both isoforms. Fig. 2 summarizes how the apoptotic sensitivity in the nominal model depends on the total initial concentration of c-FLIP, showing that the threshold initial ligand concentration, [L] thr (0), for triggering apoptosis increases almost linearly with the total initial c-FLIP concentration.
Apoptotic Response of Cell Populations-Monte Carlo simulations were performed to study the apoptotic response of a cell population with intercellular variations. The number of cells used in the simulations was 1000, and all initial concentrations and kinetic parameters in the simulation model were assumed to have the same standard deviation () in the log-normal distribution (Equations 3 and 4). Fig. 3 shows the simulated distributions of the time to apoptosis after exposure to ligands obtained for different initial ligand concentrations and various values. The simulation time was limited to 1000 min, and simulations in which apopto-sis had not occurred within this time frame were considered non-apoptotic. This approximation was motivated by the fact that, in extensive tests with a longer time frame, there were extremely few simulated apoptotic events with t ap Ͼ 1000 min. Fig. 3 shows that, for moderate intercellular parameter variations (upper panels), the time to apoptosis is quite well defined but becomes smeared out when the parameter variations are large. It is also interesting to note that, despite the statistical parameter distributions, a quite clear threshold behavior for apoptosis is observed with respect to initial ligand concentrations. Fig. 4 illustrates the proportion of simulated non-apoptotic events as a function of the initial ligand concentration in the simulation runs. An apoptotic threshold at [L](0) Ϸ 7.2 ng/ml is observed for small values. This is consistent with the nominal model (Fig. 2). Note that, at the nominal threshold ligand concentration for apoptosis, only the direction of parameter variation, and not its magnitude, determines whether the cell will undergo apoptosis or not, and therefore, the proportion of simulated non-apoptotic events is independent of the standard deviation of the parameter distributions at this point (where the curves intersect in Fig. 4).      0)) when using the previous definition of t ap (in minutes) and expressing ligand concentration in nanograms/ml. Observe that this relationship could not be predicted using a single-cell model with nominal parameters ( ϭ 0), which results in a completely different relationship between ligand concentration and mean time to apoptosis (Fig. 5).
In experimental studies of apoptosis, the results of standard measurement methods such as different types of flow cytometry analysis can be summarized in the number of apoptotic cells in a cell population, and the apoptotic process is then often illustrated using bar graphs showing the proportion of apoptotic cells obtained at given times after exposure to ligands, for example (6,10). The results of the stochastic simulations can be represented in a form compliant with those of the experimental studies by depicting the proportion of apoptotic events as a function of time. Fig. 6 shows a magnification of the panel of  Fig. 3) are 4.2, 19.2, 26.6, and 34.7%, respectively. Thus, the degree of stochastics in the initial protein concentrations and reaction rate parameters has a substantial effect on the proportion of cells that are predicted to undergo apoptosis within a given time interval. Example of Modeling Recovery after Heat Shock-To more specifically compare modeling results with experimental data, we investigated a specific case of sensitization to apoptosis caused by c-FLIP destabilization. Meinander et al. (10) studied the effect of heat shock treatment on human T lymphocytes and found that the degradation rate of c-FLIP increased considerably, whereas synthesis rates were unaffected. The half-lives of the two isoforms decreased from 40 and 120 min to 15 and 30 min, respectively. The investigators also studied how the downregulation of c-FLIP during a 30-min heat shock treatment affected the sensitivity to apoptosis when the cells were exposed to ligands after various periods of recovery from the heat shock.
The single-cell simulation in Fig. 1B indicated that downregulation of c-FLIP levels made the cell more sensitive to apoptosis when subjected to ligands. To mimic the experimental conditions reported by Meinander et al. (10), the following program was simulated. A 30-min period of steady-state conditions, in which the synthesis of both c-FLIP S/R and c-FLIP L counterbalances their net degradations, was followed by a 30-min period in which the degradation rates were kept high (using the same parameters as in the simulation shown in the middle panels of Fig. 1B). Next, the original degradation rates were restored, resulting in a recovery of the c-FLIP concentrations to their original levels. This dynamic behavior is illustrated in the upper panel of Fig. 8, which depicts the responses of the c-FLIP concentrations of the original single-cell model (18) complemented by synthesis and degradation of c-FLIP. The curves clearly show the more rapid response of c-FLIP S/R and also illustrate that the c-FLIP levels did not reach their new steady-state values during the 30-min heat shock period. The implications of these changes on sensitivity to apoptosis for the nominal single-cell model are illustrated in the middle and lower panels of Fig. 8, which show how the threshold ligand  Fig. 9B. The simulations give a somewhat faster recovery after heat shock and a slightly higher proportion of apoptotic cells than experimentally observed. These discrepancies are actually what might be expected, considering the fact that the dynamics of the stabilization and destabilization of c-FLIP have been ignored, and the real down-and up-regulation rates of c-FLIP would therefore be somewhat slower than those shown in Fig. 8. It is also worthwhile to observe that, although the previously compiled model (18), which was further developed in this study, was developed for a human B lymphoblastoid cell line (SKW 6.4), its simulations fitted well with data obtained from other cells, indicating that the model has universal applicability.
Note that the results shown in Fig. 9 cannot be directly deduced from the single-cell simulations in Fig. 8 alone. However, the time to apoptosis in the deterministic model (Fig. 8, lower panel) is slightly above 2 h throughout the heat shock experiment. Thus, the single-cell model predicts that, for a cell population, the proportion of apoptotic cells obtained 2 h after ligand stimulation is strongly correlated with the apoptotic sensitivity during the experiment.

DISCUSSION
We used a novel approach based on mathematical modeling and simulation of cell populations to study the role of dynamic regulation of the expression levels of c-FLIP proteins in controlling death receptor-mediated apoptosis. A completely novel feature of our model, compared with previous approaches (7,18), is inclusions of equations to account for synthesis and degradation of c-FLIP, with rate parameters derived from experimental observations. Simulations with the model confirmed that the concentrations of c-FLIP at the moment when the cell receptors are exposed to ligands largely determine whether the apoptotic machinery is triggered, whereas the high turnover rates of c-FLIP allow for efficient regulation of c-FLIP concentrations by affecting its stability.
Experimental studies of apoptosis often determine the proportions of cells in a cell population that undergo apoptosis after a given time, obtained as an integral measure of results determined for a large bulk of cells (e.g. in flow cytometry). To  (10). Degradation rates were increased at t ϭ 30 min and were restored to their initial values at t ϭ 60 min. Upper panel, concentrations of c-FLIP isoforms during heat shock and subsequent recovery period when there with no ligand exposure. Middle panel, threshold ligand concentration for apoptosis when the cells were exposed to ligands at various times during heat shock experiment. For 0, 1, 2, 3, and 4 h and infinite recovery periods after heat shock, [L] thr (indicated by circles) ϭ 4.0, 5.6, 6.3, 6.7, 6.9, and 7.2 ng/ml, respectively. Lower panel, t ap after ligand stimulation when cells were exposed to ligands (200 ng/ml) at various times during heat shock experiment. For 0, 1, 2, 3, and 4 h and infinite recovery periods after heat shock, t ap ϭ 129, 135, 138, 140, 141, and 143 min, respectively. be able to use such averaged quantities for model validation and comparison with simulation results, we modeled cell populations with intercellular variations by introducing log-normally distributed protein concentrations and reaction parameters. In this way, the proportion of apoptotic cells within a given time frame could be estimated using Monte Carlo simulations. As an example, the approach was applied to simulate heat shock experiments, in which down-regulation of c-FLIP levels was followed by recovery periods of various lengths before exposure to ligands. The proportions of apoptotic cells predicted by the simulation model were in qualitative agreement with experimental findings. These results demonstrate that a model developed for one cell line (B cell line SKW 6.4) could be applied to qualitatively predict experimental results obtained from other cells.
A recent model proposed significant non-genetic variability in cell populations subjected to TRAIL-induced apoptosis (26). Interestingly, our model predicts that variations of this kind occur among any given cell population, largely as a result of intercell variability in the fine-tuning of c-FLIP levels.
The results of this study and the developed model provide insight into how small system perturbations may have pronounced effects on cell sensitivity. In this way, the system characteristics and especially limiting factors among the key molecules governing death receptor responses may be determined under different conditions. For example, the different roles of c-FLIP proteins (7,8,34), reflected in their different chemical reactions and rate parameters as well as concentration levels in individual cell lines, may now be considered as a function of the factors affecting the stability of c-FLIP proteins. This could even be an outset for modeling different therapies targeted at c-FLIP levels and/or DISC signaling. Parameter estimation in population models based on experimental findings will also yield possibilities to estimate rate-limiting steps among the processes that determine the regulation of c-FLIP isoforms and thereby the processing of DISC signals. For example, as mentioned in the Introduction, c-FLIP isoforms are subjected to active isoform-specific phosphorylation (27,28) and ubiquitylation and degradation (9,10). All the post-translational modifications seem to have distinct roles and features not only in regulating the stability of the proteins but also in determining their functions. Further applications of the developed model will enable prediction of both the molecular functions and the relative importance of these individual processes as determinants of DISC signaling.