Theoretical Biology and Medical Modelling Construction and Analysis of a Modular Model of Caspase Activation in Apoptosis

which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract Background: A key physiological mechanism employed by multicellular organisms is apoptosis, or programmed cell death. Apoptosis is triggered by the activation of caspases in response to both extracellular (extrinsic) and intracellular (intrinsic) signals. The extrinsic and intrinsic pathways are characterized by the formation of the death-inducing signaling complex (DISC) and the apoptosome, respectively; both the DISC and the apoptosome are oligomers with complex formation dynamics. Additionally, the extrinsic and intrinsic pathways are coupled through the mitochondrial apoptosis-induced channel via the Bcl-2 family of proteins.

Results: A model of caspase activation is constructed and analyzed. The apoptosis signaling network is simplified through modularization methodologies and equilibrium abstractions for three functional modules. The mathematical model is composed of a system of ordinary differential equations which is numerically solved. Multiple linear regression analysis investigates the role of each module and reduced models are constructed to identify key contributions of the extrinsic and intrinsic pathways in triggering apoptosis for different cell lines.
Conclusion: Through linear regression techniques, we identified the feedbacks, dissociation of complexes, and negative regulators as the key components in apoptosis. The analysis and reduced models for our model formulation reveal that the chosen cell lines predominately exhibit strong extrinsic caspase, typical of type I cell, behavior. Furthermore, under the simplified model framework, the selected cells lines exhibit different modes by which caspase activation may occur. Finally the proposed modularized model of apoptosis may generalize behavior for additional cells and tissues, specifically identifying and predicting components responsible for the transition from type I to type II cell behavior.

Open Access
Background Apoptosis, or programmed cell death, is a highly regulated cell death mechanism involved in many physiological processes including development, elimination of damaged cells, and immune response [1][2][3][4][5][6][7][8][9]. Dysregulation of apoptosis is associated with pathological conditions such as developmental defects, neurodegenerative disorders, autoimmune disorders, and tumorigenesis [10][11][12][13][14][15][16]. The apoptotic pathway is characterized by complex interactions of a large number of molecular components which are involved in the induction and execution of apoptosis. Although scientists do not fully understand the entire pathway, key characteristics have been identified which motivates further study of this cellular process.
The extrinsic initiator caspase (caspase-8) couples the two pathways by initiating the mitochondrial apoptosisinduced channel (MAC), leading to the activation of the intrinsic pathway [24]. The subsequent cell death for either pathway is executed through a cascade activation of effector caspases (e.g., caspase-3) by initiator caspases (e.g., caspase-8 and -9) and the amplification of death signals implemented by several positive feedback loops and inhibitors in the network [4,15,16,[25][26][27][28].
Both the extrinsic and intrinsic pathways may converge at the destruction of the mitochondrial membrane. The extrinsic pathway may activate the intrinsic pathway through a mitochondrial apoptosis-induced channel (MAC) of intracellular signals involving the Bcl-2 protein family, which includes both pro-apoptotic (e.g., Bid, tBid, Bax, Bad, Bcl-xs) and anti-apoptotic (e.g., Bcl-2, Bcl-xL) members [51,52].
Mathematical models have been employed recently to gain further insights on the complex regulation of caspase activation in apoptosis [57,[64][65][66][67][68][69][70][71]. Most of these models focus on specific components of the full apoptotic machinery. Models by Eissing et al. [65] and Legewie et al. [66] emphasized only either the extrinsic or intrinsic pathways, respectively. The model of Fussenegger et al. [67] implemented both pathways but did not consider the coupling between them; however, Bagci et al. [57], Albeck et al. [72] and Cui et al. [73] modeled the mitochondrial apoptosis-induced channel. Stucki et al. [68] modeled only the caspase-3 activation and degradation but none of the aforementioned models closely track the upstream formation dynamics of the DISC and the apoptosome, which have since been modeled in detail by Lai and Jackson [74], and by Nakabayashi and Sasaki [75], respectively. Hua et al. [69,70] formulated complete system models that incorporate the differences in type I and II signaling as well as include more species, such as Smac; however not all dynamics (e.g. feedbacks) are included from previous component models [65,66,74,75]. More recently, Okazaki et al. [71] formulated a model based on Hua et al. of the phenotypic switch from type I and type II apoptotic death, but their model does not incorporate protein synthesis or degradation.
The primary focus of this work is to construct the simplest model of caspase-3 activation featuring the oligomerization kinetics of the DISC, mitochondrial apoptosis-induced channel (MAC) and the apoptosome; the dynamics of the extrinsic and intrinsic caspase subnetworks, as well as the coupling between the extrinsic and intrinsic pathways. To accomplish this, we constructed three independent functional modules [76][77][78][79]. These are implemented for the abstraction of oligomerization kinetics that simplify the full system. Analysis of the system generates predictions of key system components; furthermore, reduced models are constructed to validate the analysis for different cell types.

Model formulation
The full reaction network of the model is built from three component subnetworks (see Figure 1): the extrinsic, coupling, and intrinsic subnetworks; and three oligomerization modules (represented by gray areas in Figure 1): the DISC, MAC, and apoptosome modules. Each subnetwork captures a vital part of the full apoptotic reaction network and borrows heavily from previous work [57,65,66,70,71], while each module abstracts the oligomerization kinetics of an apoptotic complex to give a simplified net synthesis function using steady-state results [74,75].
The extrinsic subnetwork follows Eissing et al. [65] and captures the dynamics of the extrinsic pathway. The subnetwork contains the species FasL, FasR, DISC, procaspase-8 (Casp8), caspase-8 (Casp8*), procaspase-3 (Casp3), caspase-3 (Casp3*), XIAP, and BAR. The subnetwork is driven by DISC, whose formation dynamics from FasL and FasR are encapsulated by the DISC module using the results of Lai and Jackson [74]. DISC induces the cleavage of Casp8 to Casp8*, which then activates Casp3 to produce Casp3*. Positive feedback between Casp8* and Casp3* is provided by the activation of Casp8 by Casp3*. XIAP and BAR act as regulators by binding to Casp3* and Casp8*, respectively. Furthermore, degradation of XIAP is enhanced by Casp3*.
The extrinsic subnetwork can drive the intrinsic pathway through the coupling subnetwork, which describes the role of Casp8* in inducing mitochondrial membrane permeabilization and triggering the release of cytochrome c and Smac. The coupling subnetwork takes after a combination of Bagci et al., Hua et al., and Okazaki et al. [57,70,71], and contains the additional species Bid, tBid, Bax, cytochrome c (mitochondrial, Cytc; cytosolic Cytc*), and Smac (mitochondrial, Smac; cytosolic, Smac*). The subnetwork receives input from Casp8*, which cleaves Bid to produce tBid. Bax then dimerizes with tBid to form tBid-Bax 2 , which is taken as a representation of the MAC that controls the release of Cytc and Smac from the mitochondria to produce Cytc* and Smac*, respectively; the formation dynamics of tBid-Bax 2 are abstracted in the MAC module using similar methods as for the DISC module. Morever, Smac* acts as a regulator by binding to XIAP.
The intrinsic subnetwork follows the intrinsic pathway from the assembly of the apoptosome to the resulting caspase interactions. The oligomerization of the apoptosome is abstracted in the apoptosome module using the results of Nakabayashi and Sasaki [75], while the remainder of the subnetwork is simplified from Legewie et al. [66]. Additional species contained in the subnetwork include Apaf-1 (Apaf), apoptosome (Apop), procaspase-9 (Casp9), and caspase-9 (Casp9*). The subnetwork is driven by Cytc*, which binds to Apaf; activated Apaf then oligomerizes to form Apop, which cleaves Casp9 to produce Casp9*. As in the extrinsic subnetwork, positive feedback exists between Casp9* and Casp3*. Furthermore, Casp9* binds XIAP.
Constitutive synthesis and degradation rates are assumed for all appropriate species.

Steady-state abstraction of oligomerization kinetics
The oligomerization kinetics of the DISC, MAC, and the apoptosome are abstracted using steady-state results; this abstraction is a demonstration of a simple technique for modularization and model reduction. For an oligomer X with intermediate structures The abstractions for each of the DISC, MAC, and apoptosome modules are described below, where the notation is understood to apply only within each module.

DISC module
The DISC oligomerization kinetics are simplified from the crosslinking model [80][81][82] of Lai and Jackson [74] and follow the reactions where K D = k r /k f . Apply the conservation relations where r ss is given by solving , which has at most one positive root. Assume now that FADD is in excess (see, e.g., [70,71]) to obtain where it is assumed that both FasL-FasR 2 and FasL-FasR 3 can propagate the death signal [74]. Externally, in the full reaction network, the oligomerization rate function will be called as f DISC ([FasL] 0 , [FasR] 0 ; K DISC ). This abstraction reduces the order of the system by four.

MAC module
The oligomerization kinetics of the MAC module are assumed to follow a similar crosslinking model and therefore obey the reactions With the analogous notation l ≡ which again has at most one positive root. Therefore, , where the dynamical concentration of tBid is used as input. The abstraction reduces the order of the system by three.

Apoptosome module
The oligomerization kinetics of the apoptosome follow the model of Nakabayashi and Sasaki [75] with no dissociation, which considers bimolecular interactions of the form where Apop ≡ (Apaf-Cytc*) 7 . With the nondimensionalizations where τ = aa 0 t, l = k 2 /k 1 , and δ is the Kronecker delta. Integration of this system until steady state over a range of c 0 generates a curve for x 7 that may be accurately fit with a piecewise exponential function

Remarks on modularization
The steady-state profiles of the oligomerization kinetics (as shown in Figure 2) are supported by the models that motivated this simplification [74,75] and experimentally for tBid inducing a switch [49]. The abstraction enables these module simplifications to operate as inputs into the full dynamical system of apoptosis.

Model dynamical equations
The model species and reactions are summarized in Tables 1 and 2. Reaction kinetics are described by mass action, with the corresponding ordinary differential equation (ODE) system given in Table 3. Initial conditions to solve the ODEs for HeLa cells (from [65]) and Jurkat T cells (based on [70,71]), as well as steady-state abstraction parameters, are given in Table 4, where in particular the baseline value of [FasL] 0 = 2 nM corresponds to a dose which has been used to experimentally induce apoptosis (see [70]). Table 5 summarizes all model parameters (forward and reverse reactions, synthesis and degradation rates and parameters for the steady-state abstractions). Additionally, a variant of the Jurkat T cell, denoted Jurkat T*, is considered, which has the the same parameter values as Jurkat T but with k 2 = k 5 = k 12 = 0 following Hua et al. and Okazaki et al. [70,71].
The model ODEs are implemented in MATLAB R2007a (The MathWorks, Inc., Natick, Mass., USA) and solved using ode15s.

Regression analysis and model reduction
Integration of the model ODEs at baseline parameter values (Table 5) gives the [Casp3*] time courses shown in Figure 3. Both the HeLa and Jurkat T cells (the Jurkat T* case will be addressed in the results) demonstrate a characteristic behavior, whereby [Casp3*] stays low initially, then quickly switches to a high state at some threshold time.
Two quantitative descriptors are used to capture the form of these time courses: the peak activation, the maximum value of [Casp3*] attained over the time course; and the activation time, the time at which this peak is achieved. To determine the most significant aspects of the model within a given parameter regime, sensitivity analysis is performed with respect to these descriptors according to the following procedure: For a given set of baseline parameter values, we generate normally distributed random parameters about the baseline with standard deviation 5% of the baseline values. Then we simulate the model at these parameters, compute the descriptors and repeat this 100 times (the model has 54 parameters) to collect a set of synthetic data.
Since only local parameter perturbations have been considered, linear relationship y = (1X)b is assumed between the standardized descriptors y (y being one of [Casp3*] max and τ in standardized form) and the standardized random parameters X, where each row of X is a concatenation of the 54 model parameters in the order given by Table 5. The relation b is solved by multiple linear regression and large regression coefficients are taken to indicate essential components of the network. This information is used to guide the formulation of reduced models.

Results and discussion
Regression analyses and reduced models for FasL induction Regression analysis as described previously is performed for baseline HeLa parameter values. Regression coefficients for each of the descriptors show isolated peaks, indicating that only a small subset of the network is responsible for the system behavior. Particularly, the coefficients for the peak activation (r 2 = 0.9991) show strong components only at the synthesis and degradation rates a Casp3 and μ Casp3 , which together control the initial concentration [Casp3] 0 ; evidently, this turns out to largely be the case for all parameter sets considered (not shown), so the peak activation will not be generally further discussed. More interesting is the result for the activation time (r 2 = 0.9958; see Figure 4a), which, notably, shows that only the reactions of the extrinsic subnetwork appear to be essential. Accordingly, a reduced model (Figure 5a) consisting only of the extrinsic subnetwork is formulated, and validation of the reduction is given by comparison of the [Casp3*] time courses between the full and reduced models.
Note that this result should be expected since the HeLa cell was used in Eissing et al. [65] to study type I apoptosis. Surprising, a similar analysis of the Jurkat T cell, whose initial concentration parameters were used to study type II apoptosis by Hua et al. and Okazaki et al. [70,71], leads to a similar reduction. The regression coefficients (for the activation time; r 2 = 0.9903) are shown in Figure 4b, with reduction shown in Figure 5b, which is just that for the HeLa case but with XIAP omitted. It should be noted that the regression analysis does not show a strong component at k 2 , perhaps due to the corresponding reaction occurring at saturation; therefore not sensitive to small perturbations.  . This is to be compared against the low FasL cases for the HeLa and Jurkat T cells, which do not exhibit such a transition (not shown). The activation time regression coefficients for the Jurkat T* cell induced by low FasL case are shown in Figure  4d (r 2 = 0.9569), which in particular has strong components at k 7 and k 8 , which describe Bid truncation Model species and description are given. In the model, synthesis and degradation rates are given for the model system and labeled a and μ, respectively. Each reaction described highlights whether the reaction is a forward or reversible reaction by the arrows. The rates are provided from previous work. Reaction are illustrated in Figure 1. and the release of Cytc. Moreover, the peak activation regression coefficients (r 2 = 0.9972, not shown) exhibit a strong contribution by a Smac . The reduced model (Figure 5c) is correspondingly dominated by the intrinsic pathway; indeed, there is no direct interaction between Casp8 and Casp3 at all. Furthermore, as implicated by the synthesis rate of its inactive form, Smac*, and correspondingly its target XIAP, plays a vital role in achieving the correct activation level, which in particular illustrates the critical role of the shared-inhibitor motif in apoptosis as discussed by Legewie et al. [66].

Regression analyses and reduced models for mitochondrial apoptosis
The behavior of the system pathways under mitochondrial apoptosis can also be studied. Cell stressors that cause the depolarization and permeabilization of the mitochondrial membrane are functionally represented in the model by an input [tBid] 0 = 25 nM (now [FasL] 0 = 0). As for the FasL case, peak activation regression coefficients for the cases considered below are dominated by a Casp3 and μ Casp3 ; therefore, will not be further discussed.
Performing the regression analysis on the HeLa cell induced by tBid produces the activation time regression coefficients shown in Figure 4e (r 2 = 0.9705). Strong components corresponding to the reactions of the intrinsic subnetwork are observed; interestingly, the system behavior is sensitive to several extrinsic reactions as well. The model reduction is shown in Figure 6a, which demonstrates that the extrinsic caspase feedback Ordinary differential equations for the full system are given in the left hand column. Corresponding reaction velocities use mass-action kinetics are found in the right hand column. between Casp8 and Casp3 is essential to capturing the correct dynamics (compare the time course with k 2 = 0). Thus, the HeLa cell displays an apoptotic mechanism that involves the intrinsic pathway triggering the extrinsic pathway. Furthermore, the role of Smac* as an indirect activator of Casp3 through the sequestration of XIAP is recovered. Although Casp9* possesses a similar sequestration ability, the analysis reveals that the primary role of Casp9* is through direct activation of Casp3.
Analysis of the Jurkat T cell induced by tBid gives similar results (Figure 4f, r 2 = 0.9879; reduced model not shown), though the magnitude of the regression coefficient of k 13 , which describes the activation of Casp3 by Casp9*, is larger than in the HeLa case, suggesting a stronger role for the intrinsic caspase. For completeness, the Jurkat T* cell is induced by tBid is also considered. The activation time regression coefficients are shown in Figure 4g. In this case, the fit is relatively poor (r 2 = 0.8873) and some parameters are selected in error (e.g., k 1 , which has no effect on the system by construction; also note the larger number of significant components). Nevertheless, the regression serves to guide the model reduction, which in this case required manual correction. The reduced model ( Figure 6b) reveals a purely intrinsic mechanism of caspase activation. Similarly to the HeLa and Jurkat T cells, the sequestration of XIAP by Smac* is essential, while that by Casp9* may be neglected.
Although the peak activation for each of the HeLa, Jurkat T, and Jurkat T* cells is essentially identical to that obtained under FasL induction, the activation time shows a significant increase (factor increase of 2.1457, HeLa; 1.3003, Jurkat T; 1.9920, Jurkat T*). This is in general agreement with experimental evidence that caspase activation through the intrinsic pathway is delayed relative to that through the extrinsic pathway [62]. The counter on the left hand columns totals the 54 model rates and parameters for the full system. Each subscript for k, a and μ corresponds to its reaction number. The final column are the parameters used in the abstraction of oligomerization kinetics for the three modules. Type II apoptosis prediction In the preceding cases considered, type II apoptosis was observed only for the Jurkat T* cell under low FasL induction. This may be unsatisfactory since the Jurkat T* cell omits caspase feedback interactions which suggest potentially questionable biological relevance. Thus, a natural idea is to determine whether parameters leading to type II apoptosis may be predicted for the full reaction network rather than resorting to the Jurkat T* formulation.
An attempt to use the regression analysis for this task was made based on the idea of performing regression with respect to differences in the peak activation and in the activation times between a given parameter set and the corresponding set with k 7 = 0 (no Bid truncation, i.e., no extrinsic-intrinsic coupling). The intuition in this   analysis is that strong regression coefficients (assuming the modified descriptors are taken with the appropriate sign) now select parameters whose increase may effect a transition from type I to type II behavior. Furthermore, the parameters randomly perturbed are now restricted to only the synthesis and degradation rates, and to [FasL] 0 , [FasR] 0 , [Bax] 0 , and [Apaf] 0 , i.e., the parameters that control only the initial concentration, referred to as cellspecific parameters, as these are presumably the only parameters which may vary between different cell types.
Unfortunately, the regression coefficients for this analysis give poor fits (0.3884 ≤ r 2 ≤ 0.7714 for the activation time difference) for the cases considered, so the method fails. However, progress may nevertheless be made by considering the result from the case of Jurkat T* induced by low FasL. The strategy is to transform the conditions of that case into equivalent cell-specific parameter conditions. For example, the Jurkat T* cell mutes the reactions involving the action of Casp3* on other molecules. This effect may be achieved in principle by increasing [XIAP] 0 and hence the inhibition of Casp3*, which turns out to be insufficient as a result of the strong positive feedback between Casp8 and Casp3. Therefore, it is further necessary to decrease the rate at which Casp8 is activated. This may be controlled at the DISC module, so accordingly decrease [FasR] 0 (for the dependence, see Figure 2a).
At the assumed rate parameters, however, the changes in [XIAP] 0 and [FasR] 0 required to achieve type II apoptosis are rather dramatic. Note though that the dissociation rates k -4 and k -6 of Casp3*-XIAP and Casp8*-BAR, respectively, as estimated from Eissing et al. [65] are suspiciously large; if the estimate k -4 = k -6 = 10 -3 s -1 is taken instead, more consistently with, e.g., [66,70,71], then the changes required are no more than an order of magnitude. Specifically, starting with Jurkat T parameters, increasing [XIAP] 0 from 20 to 200 nM, and decreasing [FasR] 0 from 10 to 1 nM gives a cell type for which the intrinsic pathway is significant even under high FasL induction. The sensitivity regression analysis for this cell is shown in Figure 4h (r 2 = 0.9264), which displays significant components corresponding to the intrinsic subnetwork, notably at k 13 . The influence of the intrinsic pathway demonstrated by the comparison of time courses in Figure 7 shows a significant delay of caspase activation upon disabling the pathway coupling through tBid. In comparison, control results for the HeLa and Jurkat T cells show no such dependence (not shown).
Perhaps in light of this result, an alternative interpretation of the fact that the modified regressions produced poor fits occurs due to type II transitions requiring large changes that the local character of the linear regression cannot capture. This is consistent with the changes that Hua et al. and Okazaki et al. [70,71] report to effect transitions in their models (without caspase feedback), where, effectively, [Casp8] 0 was modified by a similar amount. As a final note, changes of this magnitude are likely reasonable given the inherent variability experimentally observed (see, e.g., [83]).

Activation thresholds and stability
It should be noted that the model in its present formulation is unstable, even to transient signals, as Figure 7 suggests. However, some notion of stability may nevertheless be achieved by considering activation times. This is shown for HeLa, Jurkat T, and Jurkat T* cells in Figure 8.
Consider first the case of receptor-mediated apoptosis, i. e., by FasL induction. For HeLa and Jurkat T cells, the peak activation is essentially constant (Figure 8a  signal no longer affects the response time. Analogously, the initial drop appears to define a transition region, wherein the cell switches from slow to fast apoptotic dynamics over a narrow range of the death signal input; this is indicative of some threshold-like behavior. Although this is not bistability, a sense of the existence of both low and high apoptotic states is nevertheless furnished, which, furthermore, may be made precise by introducing an artificial cutoff on the activation time to discount activations which take too long to occur. The case of the Jurkat T* cell is similar, though now the peak activation does show nontrivial variation with [FasL] 0 . However, the peak activation remains uniformly rather high which questions biological significance. Corresponding data for mitochondrial apoptosis (variable [tBid] 0 with [FasL] 0 = 0) are shown in Figure 8c, d.
The cases for the HeLa and Jurkat T cells are similar as the receptor-mediated case; however, the Jurkat T* cell appears to exhibit bistability (Figure 8c). For low FasL (approximately [FasL] 0 < 10 -2 nM), the peak activation stays low (near zero), whereas for high FasL ([FasL] 0 > 1 nM), the peak activation reaches a high state around 145 nM. Intermediate concentrations define a transition region where the cell may be interpreted to switch from life to death.
The present data was computed with a constant input for the receptor-mediated case and an exponentially decaying (i.e., transient) input for the mitochondrial case (since tBid has a constitutive degradation rate in the model). Interestingly, instituting a transient FasL signal, with estimated degradation rate μ FasL = 10 -5 s -1 , gives no discernable change to the receptor-mediated data, while setting μ tBid = 0 degrades the quality of the bistability result of the Jurkat T* case for mitochondrial apoptosis (not shown). This affords some insight into why the noted bistability is observed: by virtue of the delay of apoptosis incurred by the intrinsic pathway through the necessary activation of the mitochondrial apoptogenic factors and the assembly of the apoptosome (compare Figures 8b and 8d), the intrinsic pathway is able to better filter out transient signals.
With regard to stability, perhaps this implies that further models of apoptosis should be careful to include potentially important regulators such as cFLIP, which inhibits DISC and hence imposes a delay on Casp8 activation [70,71]. Moreover, it may likewise be prudent to expand in full any series of activations occuring sequentially; an example might be the interactions between Bid and Bax to form the MAC, which currently is not mechanistically understood and consequently may be inappropriately abstracted; although in defense of the abstraction, it is experimentally suggested that component levels (such as Apaf-1 or Casp9) may determine how quickly some cells die [83]. Finally, in this view, the steady-state abstractions presented are actually quite unsuitable for the purposes of model stability; however, modulation of the given abstracted dynamics by appropriate time-dependent functions (e.g., by an appropriate Heaviside function) may suffice.

Conclusion
This study has presented a methodological construction of a straightforward and informative mathematical model of apoptosis. This was done by combining both the extrinsic and intrinsic pathways through the implementation of functional modules and subnetworks motivated by previous models and findings [65-67, 69-71, 74, 75]. The subnetworks, responsible for the activation of Casp3 and ultimately apoptosis, included descriptions of both the extrinsic and intrinsic pathways as well as the coupling between them. Modularization of the oligomerization kinetics of the DISC, MAC, and apoptosome were achieved through the implementation of steady-state abstraction techniques.
Sensitivity analysis by linear regression was used to identify key components of the apoptotic network under various cell conditions. This allowed for the formulation  of reduced models to capture only the essential dynamics of the system. Importantly, these reductions allowed the extraction of biological insight and helped clarify the roles of specific molecular components. For example, the model predicts for the parameter regimes considered that Casp9* contributes to the activation of Casp3 by direct catalytic activation rather than through sequestration of their common inhibitor XIAP. Furthermore, the reduced models validated many previous findings, including the critical role of XIAP and the shared-inhibitor motif in mediating apoptosis [66,[84][85][86][87], as well as the transition from type I to type II apoptosis as the induction of the extrinsic pathway is decreased [70,71]. Finally, the analysis revealed the variety of modes through which caspase activation can be achieved. In the cases considered, caspase activation was observed to occur 1) solely through the extrinsic pathway 2) solely through the intrinsic pathway 3) through the extrinsic triggering the intrinsic pathway and 4) through the intrinsic triggering the extrinsic pathway.
Whether cells employ all of these modes is an interesting experimental question, with possibly profound biological significance.
The results of the regression analyses were also used to predict cell parameters (i.e., initial concentrations) that would elicit type II apoptosis, even under high FasL induction, without having to use the Jurkat T* model of Hua et al. and Okazaki et al. [70,71], which omits important caspase feedback interactions [65,66]. This adheres to the notion of highly conserving the apoptosis pathway [1,6,15], and in principle, achieving both type I and type II apoptosis using the same network. Naturally, the type II cell prediction invites experimental investigation.
Furthermore, remarks on caspase activation thresholds and stability were given. The critical element in achieving bistability in the system (at least to transient signals) appears to be related to whether sufficient delays are included. In particular, this implies the importance of modeling regulators, especially inhibitors, of the system, as well as the correct dynamical description of complex formation. Specifically, for this latter point, the present formulation neglects the time dependence of the oligomerization rate and assumes that the formation of a given final complex proceeds without delay. This, however, does not reflect actual dynamics; for example, the model of apoptosome assembly by Nakabayashi and Sasaki [75] at the parameter values considered in this study exhibits a characteristic time delay on the order of 100 min. A simple improvement is the delayed initiation of the present approximation by an appropriate time. A general theory of oligomerization that gives such a time would be particularly useful. Finally, of special interest is whether the incorporation of such delays can recover the expected type II behavior of the Jurkat T cell while maintaining the type I behavior of the HeLa cell.
Future directions for model refinement include more sophisticated treatment of oligomerization kinetics as described. A more comprehensive procedure for model reduction would also be helpful. The current method of sensitivity analysis is unable to eliminate reactions near saturation; however these cases should intuitively be treatable analytically. Moreover, the implementation of a faithful model exhibiting bistability is of primary biological interest as this would allow the formal definition and investigation of a point of no return in apoptosis. Furthermore, it may be profitable to adapt and apply the model to other cell types, e.g., mature neurons, which have repressed Apaf-1 expression and hence apoptosome formation [88,89]. Extending the presented work to model apoptosis at a cell population level may predict key mechanisms; and perhaps, prove fruitful for understanding drug sensitivity in various cell lines.
The model thus presented serves as a guide for future theoretical and experimental work in analyzing apoptosis and achieves progress toward a full model of this important biological process.