A mathematical model for apoptotic switch in Drosophila

Apoptosis is an evolutionarily-conserved process of autonomous cell death. The molecular switch mechanism underlying the fate decision of apoptosis in mammalian cells has been intensively studied by mathematical modeling. In contrast, the apoptotic switch in invertebrates, with highly conserved signaling proteins and pathway, remains poorly understood mechanistically and calls for theoretical elucidation. In this study, we develop a mathematical model of the apoptosis pathway in Drosophila and compare the switch mechanism to that in mammals. Enumeration of the elementary reactions for the model demonstrates that the molecular interactions among the signaling components are considerably different from their mammalian counterparts. A notable distinction in network organization is that the direct positive feedback from the effector caspase (EC) to the initiator caspase in mammalian pathway is replaced by a double-negative regulation in Drosophila. The model is calibrated by experimental input–output relationship and the simulated trajectories exhibit all-or-none bimodal behavior. Bifurcation diagrams confirm that the model of Drosophila apoptotic switch possesses bistability, a well-recognized feature for an apoptosis system. Since the apoptotic protease activating factor-1 (APAF1) induced irreversible activation of caspase is an essential and beneficial property for the mammalian apoptotic switch, we perform analysis of the bistable caspase activation with respect to the input of DARK protein, the Drosophila homolog of APAF1. Interestingly, this bistable behavior in Drosophila is predicted to be reversible. Further analysis suggests that the mechanism underlying the systems property of reversibility is the double-negative feedback from the EC to the initiator caspase. Using theoretical modeling, our study proposes plausible evolution of the switch mechanism for apoptosis between organisms.


Introduction
Apoptosis, a conserved process of programmed cell death, is used by multicellular organisms to remove unnecessary or damaged cells during development and adult life [1]. Abnormal regulation of apoptosis machinery is associated with a wide variety of human diseases [2]. For instance, increased apoptosis pathologically exacerbates conditions such as autoimmune diseases and Alzheimer's disease. Conversely, inadequate apoptosis leads to the development of cancer. Deciphering the mechanisms underlying apoptosis provides insights into the complex self-organization in multicellular organisms, and may enable development of effective therapy for apoptosis-based diseases.
In mammals, the apoptotic cellular signaling is mediated by two interconnected pathways, the intrinsic and extrinsic apoptosis pathways [3] ( figure 1(a)). The intrinsic pathway responds to intracellular stress induced by mitochondrial signals. The extrinsic pathway is activated by extracellular death signals, sensed by death receptors. Both pathways rely on the cascade of caspase proteins, a family of cysteine proteases expressed in an inactive proenzyme (also termed as zymogen) [4]. The caspase proteins are categorized as initiator or effector caspases (EC) based on their temporal sequence in the apoptosis pathway, whereby the initiator caspases (IC) catalyze the EC when they receive upstream apoptotic stimuli. Once the mammalian intrinsic pathway is triggered by mitochondrial signals, cytochrome c and the second mitochondriaderived activator of caspases (SMAC) are released into cellular cytosol. Cytochrome c binds to the apoptotic protease activating factor-1 (APAF1), forming a multimeric protein complex called the apoptosome. The apoptosome catalyzes the activation of caspase 9 (CASP9), a mammalian IC. Active CASP9 then catalyzes the activation of caspase 3 (CASP3), a mammalian EC. In turn, CASP3 can activate CASP9, forming a positive feedback loop [5,6]. Both CASP3 and CASP9 are repressed by the X-linked inhibitor of apoptosis (XIAP), which itself is antagonized by SMAC. The apoptotic signaling through the extrinsic pathway is of much less complexity. When the death receptor (FasR) is stimulated by death ligand (FasL), a transmembrane complex called FADD is formed. FADD directly activates the IC caspase 8 (CASP8), which subsequently catalyzes the activation of CASP3. As the converging point for mammalian intrinsic and extrinsic apoptosis pathways, activated CASP3 triggers the apoptotic response irreversibly [7].
Experiments suggest that the activation of EC occurs in an all-or-none fashion, highlighting the function of the apoptosis pathway as a molecular switch [8][9][10]. Such switch-like decision of cellular fate is not uncommon and has been observed in other cellular activities, such as cell cycle progression and developmental cell differentiation [11,12]. Notably, for mammalian apoptosis pathway, the switch properties at the single-cell level, including the switching time of EC response (the rising phase from inactive state to active state of CASP3) and the characteristic delay time (the time lapse between the input stimulation and the activation of CASP3) have been experimentally characterized [13,14].
Using the approach of mathematical modeling and computer simulation, researchers have integrated the sophisticated molecular interactions within the mammalian apoptosis pathways, and extensively investigated their dynamical and systemic behavior [3]. Previous theoretical studies from the viewpoint of systems theory suggest that the mechanistic property of bistability can achieve the switch-like behavior of apoptosis. There has been a consensus in the community that the models of apoptosis networks are necessarily bistable, with one discrete stable steady state (inactive EC) corresponding to cell survival, and the other (active EC) corresponding to cell death [15][16][17][18].
In contrast to the intense theoretical modeling work on mammalian apoptosis pathways, the apoptotic signaling mechanisms in the fruit fly, Drosophila, have not been investigated theoretically, to our best knowledge. Drosophila is an excellent model organism for the study of apoptosis due to the advantages including easy and inexpensive culturing, and amenability to genetics and molecular biology techniques.
representing respectively a process of activation, inhibition, or synthesis, and the subscripts 'P1' and 'P2' represent respectively the name of the regulator and the target for the regulation. These regulations are further explained in the modeling results as well as in figure 2. The homologous proteins (e.g. APAF1 and DARK) are depicted via the same geometric shapes in the two pathways. In (b), we decompose the Drosophila apoptosis pathway into two functional modules, each encircled by a gray dashed box. The top module is 'Input Transduction Module' and the bottom one is 'Core Activation Module'.
Most importantly, it has become well recognized that the apoptotic machinery is in general conserved in fruit flies and vertebrates [19,20]. Indeed, all of the central signaling components in Drosophila apoptosis pathway have genetic homologs in mammalian cells (table 1). Despite of the significant genetic similarity, we find in the present study that the apoptotic regulations in Drosophila do exhibit difference from their mammalian counterparts. For instance, the organization of the apoptotic signaling presents a hybrid intrinsic-plus-extrinsic pathway in Drosophila, where the intrinsic-like pathway is stimulated directly by external signals. To probe how the distinct regulatory organizations affect the systemic behavior of the apoptotic switch in Drosophila, we develop a mass-action model. We seek to employ our mathematical model of Drosophila apoptosis pathway to explore the mechanisms underlying the switch-like execution of apoptosis in fruit flies, in comparison to the counterparts in mammalian cells, mainly based on the previous theoretical work by Legewie et al [18]. The Legewie model of mammalian intrinsic pathway theoretically identifies a positive feedback loop embedded in the network that gives rise to irreversible bistable response of EC. We ask that, for the conserved function of apoptosis, whether the cellular regulatory system in Drosophila behaves the same as that in mammals or not from the viewpoint of theoretical modeling. In addition, what are the mechanistic reasons behind the similarity or distinction of their systems features, such as feedback organization and bistability? Addressing these questions through mathematical modeling and analysis could provide insightful predictions on the design principles of an evolutionarily conserved cellular function across species.
By enumerating and integrating molecular interactions in the Drosophila apoptosis pathway, we find that there exists a positive feedback loop between IC and EC in the model. This closed loop takes a special format of double-negative feedback, which is different than the format in mammalian cells. Computer simulations and bifurcation analysis show that our model of Drosophila apoptosis pathway can produce the desired all-or-none switch response. However, the systems property of bistability with respect to the input DARK is predicted to be reversible, rather than irreversible for its mammalian counterpart. We show that the mechanistic reason for the reversibility is the organization of the double-negative feedback loop in Drosophila network. Our study theoretically predicts distinctions in apoptotic switching features between fruit flies and mammals, indicating plausible systemslevel evolution of a conserved cellular function.

Methods
Model construction A mathematical model accounting for mass-action kinetics is developed. There are four types of elementary reactions in the model: reversible binding, protein synthesis, protein degradation and catalysis of caspase proteins. The biological meaning of all the reactions, except for synthesis and degradation processes, will be given in the results and discussion section. The law of mass action is applied to the mathematical modeling of the elementary reactions, following the same formula used for the model of mammalian apoptosis pathway [18]. The developed model is 19-dimensional ordinary differential equations. The biochemical species and corresponding variable names used in the mathematical model are listed in table 2 in the results and discussion section. The complete list of the elementary reactions, the ordinary differential equations, and the nominal values of the parameters for the model are described in the supporting material. The simulations for the time trajectories of variables are calculated using MATLAB's numerical ODE45 integrator.

Parameter calibration
The nominal parameter set is determined through minimization of the objective function, defined as the squared error between the simulated and the experimentally determined delay time at varying input signal. To implement the input-output relationship as experimentally reported, we assume that the dependence of the synthesis rates of DARK and RHG on the input signal is represented by Hill functions (supporting material), the well-accepted formula for lumped gene expression and protein production, as the expression of both DARK and RHG can be upregulated by extracellular proapoptotic stimuli. That is, the model as described above is mathematically connected to the external signal via two Hill functions during parameter calibration (supporting material). The parameter values are constrained to be within two orders of magnitude above and below the values for the same type of parameters used by the Legewie model of mammalian intrinsic apoptosis pathway, which are also the initial values used for the parameter search. For parameter search, first, stochastic search is performed using the genetic algorithm (GA) implemented by the MATLAB Global Optimization Toolbox. If the optimization via GA does not converge within a preset time limit of 2 h, further minimization will be performed using deterministic nonlinear programming implemented by the 'fmincon' function of the MATLAB Optimization Toolbox, starting from selected best results in the final population of GA search. If the deterministic search eventually converges to different local optima, the value with the best performance is chosen as the final optimum.

Bifurcation and bistability analysis
The bifurcation analysis of steady-state response of active DRICE versus a single parameter is performed using XPPAUT (http://math.pitt.edu/~bard/xpp/xpp. html). The two-parameter bistable region in figure 5 of the results and discussion section is determined based on the single-parameter bistability analysis at sampled parameter points in the two-dimensional parametric plane. For figure 7 in the results and discussion section, the bistable diagram with respect to the concentration of DARK is further determined as irreversible or reversible, at sampled parameter pair of k 25 and k 26 .

Results and discussion
Modeling the Drosophila apoptosis pathway To develop the model of Drosophila apoptosis pathway, we first identify from literature the essential signaling components, which all have mammalian homologs (table 1). Specifically, the IC and EC in the fruit fly are called DRONC and DRICE, which are homologs of CASP9 and CASP3 respectively [21]. The apoptosome scaffolding protein in Drosophila is DARK, equivalent to APAF1 in vertebrates [22]. The inhibitor of apoptosis in Drosophila is called DIAP [21]. Finally, there is a family of four pro-apoptotic proteins in Drosophila, Hid, Grim, Reaper and Sickle, which are collectively known as RHG proteins [23]. The RHG family serves the same function as SMAC, its counterpart in mammals, to inhibit the DIAP protein [24]. This study seeks to highlight the distinction, if any, between the Drosophila and mammalian apoptosis mechanisms. To this end, we classify the elementary reactions within the Drosophila signaling pathway into two groups, one is the same as that in mammals, and the other is different. Such enumeration allows immediate qualitative comparison of the Drosophila apoptotic regulations to those in mammals. Interestingly, we find that although the schematic regulations along the Drosophila pathway (e.g. the regulations DARK act IC , DIAP inh IC , RHG inh DIAP , DIAP inh EC , IC act EC in figure 1(b)) generally resemble their counterparts of the mammalian intrinsic pathway (figure 1(a)) [25], the underlying molecular reactions have substantial degree of distinction. In our model, the only elementary reactions strictly conserved between Drosophila and mammals are the synthesis and degradation processes of the five signaling proteins (DARK, RHG, DIAP, IC, and EC), as well as the reversible binding reactions of DARK, RHG and DIAP to their targets, for catalysis (i.e. DARK) or inhibition (i.e. RHG and DIAP). Other elementary biochemical reactions proceed differently in Drosophila than in mammals.
The difference is caused by two major reasons. First, the expression of the most-upstream signaling proteins DARK and RHG are upregulated by extrinsic apoptotic stimuli [26,27]. Therefore we refer to the Drosophila apoptosis pathway as a hybrid type, combining the characteristics of both extrinsic and Table 2. List of the molecular species and corresponding variable names in the model.

SPECIES
DESCRIPTION VARIABLE intrinsic pathways in mammals, where the resemblance to the mammalian extrinsic apoptosis pathway is the stimulation by extracellular input. The second and most important reason is that the regulations of the initiator caspase DRONC are especially complicated, due to the combinatorial cleavage of two functional protein domains on DRONC. Specifically, the DIAP-binding domain of intact DRONC protein can be cleaved by the Drosophila active EC, denoted as DRICE*, resulting in the DRONC 0 molecular species [28][29][30]. On the other hand, when two molecules of intact DRONC protein bind with DARK, they undergo autocleavage of the so-called 'small subunit' domain and the resulting molecular species is denoted as DRONC* [29,31]. The DRONC* and DRONC 0 molecules can then be converted to DRONC 0 ⁎ through further cleavage by DRICE* and autocleavage, respectively. As a consequence of these cleavage processes, the DRONC 0 and DRONC 0 ⁎ molecules are refractory to the repression by DIAP through complexation, due to the absence of the DIAP-binding domain [28,29], while the homo-and hetero-dimers formed among the DRONC* and DRONC 0 ⁎ species, namely the (DRONC*) 2 , (DRONC 0 ⁎ ) 2 and (DRONC*·DRONC ) 0 ⁎ species, are capable of catalyzing the inactive EC, DRICE, into the active form DRICE*. Note that the active form of DRONC cannot association with DARK [32]. Therefore we assume that the DRONC* species does not form a complex with DARK, whereas in mammals active CASP9 remains associated to the APAF1-adpating apoptosome [32]. These diverse species of the DRONC protein either in monomer (DRONC, DRONC 0 , DRONC* and DRONC ) 0 ⁎ or complex format (formed among the DRONC species, or with DARK or DIAP), result in a considerable population of reacting species (table 2). Moreover, the above qualitative analysis of molecular interactions suggests a notable distinction in the Drosophila network organization. That is, there exists an inhibitory feedback from the active DRICE to the activity of DIAP, through its cleavage of the DIAP-binding domain of the DRONC protein ( figure 1(b), regulation EC inh DIAP ). Together with the inhibition of DRONC by DIAP ( figure 1(b), regulation DIAP inh IC ), it forms a double-negative feedback from EC to IC, namely EC─|DIAP─|IC, where '─|' denotes inhibition. Whereas in the mammalian apoptosis pathway, EC is positively fed back to IC directly, as active CASP 3 can enhance the activity of CASP9 via an enzymatic cleavage process [5,6,18].
Modularity is an essential systems property that has been theoretically explored for mammalian apoptosis pathways [33,34]. Specifically, a complex biomolecular network can be decomposed to several functional modules, where each module contains a cluster of components performing a relatively independent task. Based on the signaling scheme depicted in figure 1(b), we can qualitatively identify two functional modules in the Drosophila apoptosis pathway, which are similar to those in the mammalian apoptosis models [33,34]. The top module is termed as the 'Input Transduction Module', consisting of DARK and RHG ( figure 1(b)). This module relays the external apoptotic stimuli to the downstream processes. The bottom module contains IC, EC and DIAP that modulate the activation of apoptotic response ( figure 1(b)), and thus is denoted as the 'Core Activation Module'. Previous studies have indicated that the feedback loops among the caspases are critical for mammalian apoptosis pathways. We seek to probe if, analogouslyly, the Core Activation Module plays a central role for the system behavior of Drosophila apoptosis network.
The diagram of the biochemical reactions in our model of Drosophila apoptosis pathway is plotted in figure 2. We now explain the biochemical meaning of the pathway regulations enumerated in figure 1(b), and their associated elementary reactions modeled in figure 2.
First, the regulations denoted as 'syn DARK ' and 'syn RHG ' induced by external signal in figure 1 The inhibition of DARK by DRONC, denoted as regulation ' IC inh DARK ' in figure 1(b), accounts for the inhibitory proteolytic cleavage of DARK by active DRONC as shown in recent experiments [35]. This inhibition process is modeled by the reactions 14-16 in  The regulation ' RHG inh DIAP ' in figure 1(b) represents the inhibition of DIAP by RHG, through the reversible binding reaction 23 in figure 2. Once DIAP is sequestered into the DIAP·RHG complex as shown in formula ( RHG inh DIAP ) below, it gets degraded, following the same mechanism of SMAC-XIAP interaction [18] DIAP RGH DIAP RHG.
The regulation ' EC inh DIAP ' in figure 1(b) corresponds to the catalytic removal of the DIAP-binding Finally, the regulation ' IC act EC ' in figure 1(b) refers to the catalytic activation of DRICE by DRONC. This process is modeled by reactions 7-10 in figure 2 and formula ( IC act EC ). In particular, reactions 8-10 represent the activation of DRICE by active DRONC dimers. We also model the catalysis of DRICE by zymogen DRONC by reaction 7, assuming its catalytic rate is hundreds of times smaller than those of reactions 8-10 [18,37]. The complete list of reactions and the governing differential equations are described in Supporting Material. It is noteworthy that our model includes a large number of oligomerization processes, a feature that also appears in previous models of mammalian apoptosis pathways [34,38].
Parameter calibration for apoptotic switch The dynamics of the apoptosis pathway has been characterized as a 'slow induction plus fast switching' response of the EC, as observed at both average cellpopulation level and single-cell level [8,10,13,14,39]. Specifically, the activity of EC stays low for a while after the stimulus signal is applied, and this slow induction period is quantified by an apoptotic delay time, denoted by T d . After a slow induction phase, EC is activated, quickly switching from the survival state to the apoptotic state. Finally, EC will maintain at the high apoptotic level. Experimental quantification indicates that the execution of apoptosis generally follows the same fashion of 'slow induction plus fast switching' behavior in both mammalian and fly cells [13,40,41]. Nevertheless, the apoptotic response in Drosophila remains poorly characterized.
Due to the lack of detailed data and kinetic rate constants for the apoptosis process in fly cells, we choose to calibrate the dynamics of our model using data obtained from mammalian study, based on the recognition that the apoptotic pathways and regulatory modules are highly conserved between organisms, and the general trend of caspase response are similar between flies and mammals. Previous experimental and theoretical studies have suggested that the apoptotic delay time T d is an important system feature for mammalian extrinsic apoptosis pathways at the single-cell level [13,14,39,42]. It characterizes the transduction duration of the cell signaling information from apoptotic stimuli to apoptotic response, and is inversely dependent on the strength of stimuli. Therefore, the value of T d and its relationship to apoptotic stimuli offers an excellent quantitative measure for evaluating the system behavior of a single-cell model of Drosophila apoptosis pathway, because of the functional conservation between flies and mammals. To this end, we first adopt the ranges of the parameters used in the Legewie model for mammalian apoptosis pathway [18]. Same as [18], our model contains only elementary reactions formulated based on the mass action kinetics. The model parameters can be categorized into five functional types according to the associated reactions: protein synthesis rate, protein degradation rate, binding rate, unbinding rate, and enzyme catalysis rate. We obtain ranges for the five categories of parameters from the Legewie model and apply them to constrain the parameters of our model. Second, we calibrate the constrained parameters to achieve experimentally quantified switching behavior [14]. The experimental data of T d in relationship to varying apoptotic input measured in mammalian cells has shown that higher input level leads to faster response and hence smaller T d . In particular, singlecell data of T d measured across a population, with quantified mean and standard deviation values (table 1 of [14]) , is used here to calibrate the parameters. The delay time T d is defined as the time it takes to reach the half-maximum EC activity, analogous to the quantity previously used in experiments to characterize the response time of extrinsic apoptosis pathway in single cells ( figure 3(a)) [14,42,43]. The calibration of parameter set is performed via a hybrid searching method combining a stochastic parameter search in the global space, implemented by GA, with a deterministic parameter optimization in the local space, implemented by nonlinear programming (Methods). The resulting parameters are referred as the nominal parameter set of the model hereafter, and their values are listed in the supporting material.
At the nominal parameter set, the difference between the simulated and the experimentally measured values of T d at varying input strength is minimized. The comparison between the results from simulation and experiment shows that the nominal parameter set can generate delay times that fall within one standard deviation of the mean value measured experimentally ( figure 3(b)). In addition, the simulated time trajectories demonstrate that the activation of DRICE has a steep sigmoidal shape, whose activation time from the OFF state to the ON state does not change as the level of input varies ( figure 3(a)). This constant activation time is also in good agreement with previous experimental observation [14,42]. Thus far, we have developed a mathematical model for the average single-cell behavior of apoptotic switching in Drosophila.

Robust bistability of Drosophila apoptosis pathways
The sharp transition of the time course of DRICE from a near-zero state to a sustained high state in the simulation suggests that the apoptosis system in Drosophila is likely bistable. Bistability is a desired systems property for an apoptosis model, as it enforces two discrete stable states, with one representing the OFF state, the other representing the ON state, of the molecular switch. Consequently, the cell fate decision of apoptosis will be executed in an all-or-none bimodal fashion, avoiding potential indefinite failure of cellular suicide that may be harmful to the organisms.
To theoretically confirm that the model of Drosophila apoptosis pathway has two stable states, we perform bifurcation analysis, where the steady states and their stability are calculated in an extended region of a parameter (Methods). Figure 4 shows the bifurcation diagrams of the steady-state active EC versus four single parameters, namely the basal synthesis rate constants of DRONC, DRICE, DIAP and RHG. These diagrams show that the model of single-cell apoptosis pathway in Drosophila is indeed bistable at the nominal parameter set (indicated by the gray dashed line in figure 4), with two stable states and one unstable state in between. The low steady state of DRICE* is near zero (<10 −4 nM) and the high steady state of DRICE has a value greater than 100 nM. Moreover, the bistable behavior persists robustly in extended singleparameter intervals between the two saddle-node points (figure 4), allowing the apoptotic switch to perform reliably in the presence of parametric perturbations.
The different configurations of the bistability curves reflect the function of the parameters associated with the four major signaling proteins. For instance, the synthesis rate constants of DRICE and DRONC both serve to enhance the activity of the caspases. Therefore, the apoptotic switch is always OFF (ON) below (above) the bistable interval (figures 4(a) and (b)), indicating a positive correlation between the abundance of caspase proteins and their function to evoke apoptosis response. On the contrary, the apoptotic switch is always ON (OFF) below (above) the bistable interval of the synthesis rate constant of DIAP  [14]). The black dots represent the measured mean value of T d at specific input signals and the trend of input-output relationship is interpolated by a black curve. The mean ± standard deviation for T d is bounded within a pair of '+' symbols.
( figure 4(c)), which suggests that the amount of DIAP negates the caspase response. The synthesis of RHG family protein also promotes the caspase response, as the switch maintains at an ON state when RHG is above the threshold concentration at the right saddlenode point ( figure 4(d)). The bistability, however, persists when the basal synthesis rate of RHG is zero, meaning that apoptosis can be induced by input signal mediated by the DARK apoptosome even if RHG is not present. Qualitatively speaking, the signaling structure of Drosophila apoptosis activation using only the assembly of apoptosome is in analogy to that of the receptor-clustering model of the mammalian extrinsic apoptosis pathway, where the activation occurs independent of ligand [44].
To further explore the robustness of the bistable behavior, we plot the region of bistability in two-parameter planes. Figures 5(a)-(h) demonstrates the results against eight pairs of parameters. The parameter pairs are of three natures: pairs of bindingunbinding rates ( figure 5(a)), pairs of irreversible catalytic rates (figures 5(b)-(d)), and pairs of synthesis and degradation rates of primary proteins (figures 5(e)-(h)). In all eight cases, the white dot represents the location of the nominal parameter of the system. The bistable behavior persists in reasonably large regions in all of the two-dimensional parameter planes. Therefore, the bifurcation analysis of our model predicts that the Drosophila apoptosis pathway exhibits robust behavior of bistable response. Such systems property of robust bistability is similar to that of the mammalian apoptosis pathways as reported previously [18,43].

Reversible bistability and feedback topologies
Analysis of the Legewie model of intrinsic apoptosis pathways in mammals shows that the bistable caspase response induced by the input signal of APAF1 is irreversible [18]. It is known that irreversible bistability yields unflappable or 'long-term' memory of activation, which is a desired feature for the mechanism underlying a critical cellular decision process [12]. Long-term memory guarantees that once a fate is chosen, the cell will maintain the historic fate despite of future retraction of the inducing agent, thus preventing decision error in the event of transient noise. Since the function of APAF1 as an input to induce the activation of caspase is essential for the operation of mammalian apoptosis pathways, analogously we analyze the response of our model to the input signal of DARK, the homolog of APAF1 in Drosophila. This comparison allows for the theoretical interrogation of whether such an error-proof property, beneficial for the performance of cellular decision making, may be present in flies.
The bifurcation diagram versus the abundance of DARK shows that the two saddle-node points of the bifurcation curve both are located in the first quadrant, implicating that the system is reversibly bistable (figure 6(d)), instead of irreversible as in mammals. That is, in theory, after the level of DRICE has reached the high steady state when passing the right saddlenode point, the removal of the DARK input can reverse the response to the low steady state. According to our model prediction, for Drosophila cells, not only the activation of apoptotis response but also the maintenance of apoptosis phenotype require sufficient amount of DARK. If the DARK protein is retrieved or lost, a Drosophila cell that already enters the apoptosis stage could turn back to be alive again due to deactivated EC. Therefore, the supply of DARK needs to be continuously stable during the whole process of programmed cell death to guarantee its proper completion. This predicted reversible bistability can be experimentally tested by observing the DRICE* activity at increasing level of DARK followed by decreasing level of DARK, where DRICE* should be induced above a critical amount (corresponding to the right saddle-node point) during the DARK-addition phase and then diminish below another critical amount (corresponding to the left saddle-node point) during the DARK-reduction phase. Importantly, the lack of irreversibility in the bistability with respect to DARK input for the Drosophila model suggests a potential distinction from the behavior of the apoptotic switch in mammals.
We continue to computationally probe what mechanism may lead to the disparate bistable behaviorof apoptotic switch in Drosophila. The aforementioned qualitative analysis of molecular interactions has pointed to a distinction in the network organization between flies and mammals in that a direct positive feedback from EC to IC exists in mammals while a double-negative feedback topology exists in Drosophila ( figure 6(a)). Since a number of theoretical analysis has indicated that the signaling interconnections of a positive feedback loop dictate its capacity to generate irreversible bistability, we hypothesize that the different topologies of the positive feedback from EC to IC yields the distinctive bistable behavior in flies and mammals. To test this hypothesis, we study the impact of two topologies of positive feedback from EC to IC on the bistable behavior of Drosophila apoptotic switch: (i) feedback with the direct activation of DRONC by DRICE* (supporting material, reaction 43), which operates in mammals [6,18], and (ii) feedback where the cleavage of DRONC by DRICE* removes the DIAP binding site (figure 2, reactions [11][12][13], which operates in flies [29,30]. To that end, we analyze the Drosophila models with different combinations of topologies (i) and (ii) of the feedback from EC to IC, while keeping all other reactions intact, and plot the steady-state response of DRICE* at varying level of DARK input. In the presence of both (i) and (ii) (i.e. with the positive feedback in figure 6(b)), the left saddle-node point of the bifurcation curve moves to the second quadrant, and thus the bistable activation of DRICE becomes irreversible (figure 6(e)). In the presence of (i) and the absence of (ii) (i.e. with the positive feedback in figure 6(c)), the bistable curve is again irreversible, where the left saddle-node point moves toward further left of the second quadrant (figure 6(f)). The results in figures 6(e) and (f), differing from the reversibility in figure 6(d), suggests that the positive feedback in the format of direct activation of IC by EC is the plausible component responsible for the irreversible bistability. In the irreversible scheme, the direct mutual upregulation enables the two caspase proteins (IC and EC) to sustain each other without any upstream signal of DARK once the circuit is switched on. In contrast, in Drosophila the repressive action of DRICE on DIAP, by inhibiting its sequestering binding to DRONC, forms a double-negative interaction between DRICE and DRONC. Although the net effect of such feedback  figure S1 for details). This type of bifurcation is called subcritical Hopf bifurcation. Such bifurcation, nonetheless, does not affect the region of bistability or the reversibility of the bistability.
loop is still a positive nature with system bistability, it lacks the catalytic conversion of zymogen DRONC to active DRONC. And thus when the DRONC-catalyzing protein DARK is absent the caspase response will switch off, whereby the reversibility arises. These results of single-parameter bifurcation analysis implicate that the reversible versus irreversible bistability of the apoptotic switch is due to the disparate feedback organizations between flies and mammals.
We further evaluate this finding by checking the bistable behavior of the system with the three aforementioned feedback topologies when varying the synthesis rate of two caspase proteins, DRONC and DRICE (i.e. k 25 , k 26 ). Specifically, at varying values of k 25 and k 26 , the bistability curve with respect to the concentration of DARK is checked for reversibility and irreversibility. First, for the original model of Drosophila pathway with the double-negative feedback (ii), the entire gray area plotted in figure 5(e) presents reversible bistability. Second, when the direct activation from EC to IC (i) is added on top of (ii), partial area of the bistable region becomes irreversible ( figure 7(a)). In the third case where we remove the double-negative feedback (ii) while keeping only the direct positive feedback (i), the bistable region is completely irreversible ( figure 7(b)). Hence we infer that the topology (i) is the mechanism underlying the irreversible bistability, whereas the topology (ii) is the mechanism underlying the reversible bistability.
Finally, we check the bistable reversibility/irreversibility presented by the models incorporated with the different feedback topologies listed in figures 6(a)-(c) at varying value of parameters, as it has been suggested that such system property may be parameter value dependent [18]. First, we find that the reversible bistability of DRICE* versus DARK for the model of Drosophila apoptosis pathway is maintained across all the two-dimensional bistable planes in figure 5. Furthermore, we randomly vary the model parameters and gauge their values using comparison of the performance of the Drosophila apoptosis model with the experimental data shown in figure 3(b). Specifically, the apoptotic delay time T d of the model with the randomly chosen parameters is constrained to be distributed within zero to five standard deviations from the experimental mean, where three standard deviations 'nearly certainly' cover all the statistically possible cases [45] based on the experimental observation, and going beyond three standard deviations gives rise to cases significantly far from the nominal behavior. The resulting bifurcation diagrams of DRICE* versus DARK show that, when the model performance is within three standard deviations, the bistability behavior in figure 6 remain essentially the same in that, at all the varied parameter sets, the model with the nominal feedback configuration in figure 6(a) is reversibly bistable, while the models with the feedback configurations in figures 6(b) and (c) are irreversibly bistable, with the latter showing wider regions of irreversibility (see figures S2-S4 in supporting material). When its performance is beyond three standard deviations, the model with the nominal feedback is again reversibly bistable at all the varied parameter sets (first row in figures S5-S6), and the models with direct positive feedback show irreversible bistability, albeit with biologically unattainable (negative) lower stable branch (second and third rows in figures S5-S6). Therefore, the Drosophila apoptosis model exhibits reversible bistability with significant degree of parametric robustness. The irreversibility of the models with the direct positive feedback is at least robust with respect to parameters yielding similar performance to the nominal set.
Taken together, our results show that the known differences in the network connectivity between Drosophila and mammalian network have the potential to generate a functional difference in reversibility of the apoptotic cellular decision. Note that the double negative feedback structure belongs to the Core Activation Module ( figure 1(b)). Therefore, the results also show that the Core Activation Module supplies the central mechanism underlying the system behavior of the Drosophila apoptosis pathway. In addition, since irreversible bistability endows the superior property of long-term memory of a transient stimulus [12], this study predicts that the apoptotic switch in mammals possesses a functionally more advanced systems property for a cell-fate decision process than its counterpart in Drosophila. Our theoretical finding supports that the caspase activation in Drosophila plausibly follows a more primitive mechanism that may be the precursor to the pathways in mammals, a proposition previously made based on qualitative examination of apoptosis pathways [46].

Conclusions
In this study, the molecular mechanisms of the conserved cellular process of apoptosis in Drosophila are probed by the method of theoretical modeling. We find that although the pathway of apoptosis is mostly conserved between flies and mammals, the underlying elementary reactions differ considerably. The massaction model of Drosophila apoptosis pathway gives rise to all-or-none dynamic of caspase response, which is desired for the system to act as a molecular switch. Same as the mammalian case, the model of apoptosis pathway presents robust bistability for Drosophila, with two discrete states designating the life-or-death cell fate. However, the model predicts that, in contrast to the irreversible bistability of the caspase response to the APAF1 induction in mammals, the caspase activation elicited by DARK, the homolog of APAF1, is reversibly bistable. Further analysis shows that the reversible bistability arises from the absence of the direct feedback activation of the initiator caspase by the EC. Due to the essential role played by the irreversibility in the robust apoptosis decision in mammals, our finding highlights a potential mechanistic distinction between the apoptotic switch in flies and that in mammals. In addition, the theoretical results propose a conjecture that the molecular mechanism of a conserved cellular function may acquire properties that are more beneficial for systems performance of regulatory networks through evolution.
Author Contributions L M and R Z designed the research. R Z and L M developed the model. R Z performed the computation. R Z and L M analyzed the data. L M and R Z wrote the paper.