A Minimal Model of Signaling Network Elucidates Cell-to-Cell Stochastic Variability in Apoptosis

Background Signaling networks are designed to sense an environmental stimulus and adapt to it. We propose and study a minimal model of signaling network that can sense and respond to external stimuli of varying strength in an adaptive manner. The structure of this minimal network is derived based on some simple assumptions on its differential response to external stimuli. Methodology We employ stochastic differential equations and probability distributions obtained from stochastic simulations to characterize differential signaling response in our minimal network model. Gillespie's stochastic simulation algorithm (SSA) is used in this study. Conclusions/Significance We show that the proposed minimal signaling network displays two distinct types of response as the strength of the stimulus is decreased. The signaling network has a deterministic part that undergoes rapid activation by a strong stimulus in which case cell-to-cell fluctuations can be ignored. As the strength of the stimulus decreases, the stochastic part of the network begins dominating the signaling response where slow activation is observed with characteristic large cell-to-cell stochastic variability. Interestingly, this proposed stochastic signaling network can capture some of the essential signaling behaviors of a complex apoptotic cell death signaling network that has been studied through experiments and large-scale computer simulations. Thus we claim that the proposed signaling network is an appropriate minimal model of apoptosis signaling. Elucidating the fundamental design principles of complex cellular signaling pathways such as apoptosis signaling remains a challenging task. We demonstrate how our proposed minimal model can help elucidate the effect of a specific apoptotic inhibitor Bcl-2 on apoptotic signaling in a cell-type independent manner. We also discuss the implications of our study in elucidating the adaptive strategy of cell death signaling pathways.


Introduction
Cellular signaling networks are designed to sense an environmental stimulus and respond in a strength dependent manner. In this manuscript, we develop and study a minimal model of a signaling network that can respond to an external stimulus in a manner such that the activation is fast under a strong stimulus but slow if the stimulus is weak. We derive the minimal network by assuming that the cell-to-cell variability dominates the slow signaling activation, under weak stimuli, in order to adapt to a fluctuating environment. In such a scenario, population average over many cells cannot capture cell-to-cell variability in signaling. We employ stochastic differential equations and stochastic simulations to study the signaling response in our proposed minimal signaling network. We carry out a sensitivity analysis of the minimal model with respect to parameter variations that also provides simple quantitative relations connecting different parameter values. We further use probability distributions of signaling molecules to characterize differential signaling response in the minimal network and such distributions show very distinct types of behavior depending on the strength of the stimulus. Specifically, for the case of a weak stimulus, a characteristic bimodal distribution is obtained for the activation of a downstream signaling molecule indicating large cell-to-cell fluctuations.
Interestingly, the results from our minimal stochastic signaling model capture the essential stochastic signaling behavior observed in simulations of complex apoptotic cell death signaling pathways [1]. Details of the apoptotic signaling response vary depending on the cell type under consideration and also on the type of apoptotic stimulus applied [2][3][4][5][6][7]. Our developed minimal signaling network demonstrates that large cell-to-cell stochastic variability increases as the strength of the stimulus is decreased, a feature also observed in large scale simulations and experiments of apoptosis, irrespective of cell types and the stimulus types used in those studies [1,2,7]. Hence, the minimal signaling network developed here captures cell-type independent features of apoptosis signaling and thus can serve as a general signaling model of apoptosis. We also discuss how the study of such a minimal network can provide crucial insights into the process of biological evolution of apoptosis signaling pathways.

The minimal signaling network
We designed our proposed minimal signaling network to respond to an external stimulus in a manner dependent on the strength of the stimulus.
(i) Under a strong stimulus, activation of the signaling network is fast and there is no need for cell-to-cell stochastic fluctuations; (ii) whereas under a weak stimulus, activation of the signaling network is slow with large cell-to-cell variability so that it can adapt to a fluctuating environment.
We attempted to find a minimal set of signaling reactions that can achieve the above response, with the nature of signaling reactions being that of irreversible chemical reactions of the type X ?X Ã . Rapid deterministic activation under a strong stimulus can be achieved through a one-step direct activation of downstream signaling molecules (Figure 1a (i)). However, this one-step pathway would generate a slow deterministic activation of the signaling network under weak stimuli. Hence, for the case of a weak stimulus, we need a minimum of two separate rate constants to generate slow activation with large cell-to-cell variability. First, a low-rate constant step is required to generate slow stochastic activation of signaling molecules as the low kinetic rate constant will lead to a low probability of signaling activation with enhanced stochastic effects. However, one slow step cannot generate large cell-to-cell stochastic fluctuations and therefore it is essential to add a second fast step. This fast reaction will activate the signaling network rapidly once a few molecules have been activated by the low-probability event of a slow step (Figure 1a (ii)). This slow-fast combination will presumably lead to the slow stochastic response in a manner such that large cell-to-cell variability with all-or-none type activation of the signaling network is achieved. The same external stimulus activates both pathways of the minimal network. Hence, we need to combine the two pathways in such a manner that, under a strong stimulus the one-step deterministic (denoted by type A) pathway wins, and,s under a weak stimulus the stochastic (denoted by type B) pathway dominates. This is achieved by adding a fast step to the second stochastic pathway that connects the second pathway to the external stimulus before the slow-fast steps. The rate constant for that fast step is chosen to be greater than that of the deterministic reaction of the first pathway to ensure that the stochastic pathway dominates under weak stimuli.The final downstream signaling molecule in the two pathways are considered to be the same molecule so that the two pathways are joined in a loop structure (Figure 1b). We used a combination of stochastic differential equations and stochastic simulations to study the behavior of the signaling network as the strength of an external stimulus is varied.
Signaling in the minimal network: stochastic differential equations We derive a stochastic differential equation for the one-step irreversible reaction X ?X Ã by equating the corresponding Fokker-Planck equation with the continuum limit of the Master equation for that activation reaction [8]. The stochastic equation is given by x denotes the number of molecules of the signaling species X that are activated to X Ã . k is the reaction rate constant and dW denotes the differential white noise (W (t) is a Brownian process). The equation that governs the number of activated species is given by dx Ã~{ dx~kxdtz ffiffiffiffiffi ffi kx p dW , which shows that the total number xzx Ã of a signaling species remains conserved. Eq. (1) is a nonlinear stochastic differential equation (SDE) with interesting structure of the noise term: for a low rate constant (k%1), the noise term may not vanish even in the presence of a large number of molecules, which is in contrast to earlier considerations where stochastic effects were thought to be important only when a small number of signaling molecules were involved [9,10]. We are now able to write the set of stochastic differential equations that describes all four reactions in type A and type B pathways: The number of molecules of the activated species corresponding to x i is denoted by x Ã i and obeys the equation dx Ã i~{ dx i . The reaction rates in Eq. (2) are given by k 1~k Thus the reaction rates depend on the numbers of activated signaling molecules from the previous reaction step. x 0 is taken as the external stimulus that activates the signaling network under consideration. We denote the initial concentrations of all the inactivated signaling species by x 1 (0)~x 0 1 , and x 2 (0)~x 0 2 , x 3 (0)~x 0 3 , and the initial number of activated molecules for all three species are set to zero. Hence, the total number of molecules for the i th signaling species is taken to be x 0 i . Same variables are used in our Monte Carlo study of the signaling network. We first consider the solutions of the above set of equations for the type A and type B pathways separately and then study their combined behavior.
(a) Deterministic signaling through type A pathway. This is the one step signaling reaction described by the dx 3 reaction in Eq. (2) with k 32 for the type B pathway being set to zero. This equation has the generic form of the Eq. (1): dx~{kxdt{ ffiffiffiffiffi ffi kx p dW and is amenable to analytical solution. We make a change of variable z~ffi ffiffiffiffi ffi kx p . The solution to the equation in terms of our new variable z is given as s is a variable that is integrated over and dW is a stochastic differential (W (s) is a Brownian process). The number of signaling molecules of x 3 is calculated from the above solution using x 3~z 2 =k. The number of activated molecules is then estimated simply by using x Ã (b) Stochastic signaling through type B pathway. Here we use the three steps signaling reactions described by the equation set (2) with k 31 for the type A pathway set to zero in the last equation. The signaling reactions in each step again have the generic form of Eq. (1) which is amenable to an analytical solution. However, such a scheme with step-by-step solutions where the solution in one step is fed into the reaction constant of the next step quickly leads to complicated expressions that are not easy to analyze. Hence we numerically solved the set of stochastic differential equations given for the type B pathway. In Figure 2b, we show the activation of the end point of the signaling network, namely the activated x 3 molecules, as obtained from the numerical solution of the type B reactions. Clearly, in the type B signaling, the average behavior does not capture cell-to-cell stochastic fluctuations. We also probe the robustness of such a behavior in terms of parameter variations. The upper bound of the rate constant for the intermediate step k 0 2 is constrained by two factors: (a) it has to be small (k 2 x 2 v1) to generate stochastic fluctuations in the type B pathway, and (b) it also has to be significantly smaller than the rate constant for the final activation step (k 0 32 &k 0 2 ) to generate all-or-none type activation with large cell-to-cell variability. Even though there is no lower bound needed to be defined for the constant k 0 2 , it sets the time-scale of the signaling activation. In addition, lower the value of k 0 2 , lower the amount of stimulus needed for the type B pathway to get activated (or equivalently longer the time to switch to the type B pathway). Hence the slowness of the stochastic pathway under weak stimuli also determines the relative dominance between type A and type B pathways.
(c) Transition from deterministic to stochastic signaling in a combined network of type A and type B pathways. By comparing the two coefficients in the equation for the x 3 activation, we can derive a condition for which the deterministic type A pathway will dominate over the stochastic type B pathway (over the entire time-course of activation). The relevant coefficient for the type A pathway is *k 31~k 0 31 x 0 and the same for the type B pathway is *k 32~k 0 32 x Ã 2 . We approximate the intermediate reaction step of the type B pathway by the stochastic activation term such that dx 2 &{ ffiffiffiffiffiffiffiffiffi ffi Activation of the first step of the type B pathway occurs in a rapid manner and we could approximate Thus comparing the coefficients k 0 , we derive a condition for which the deterministic type A pathway will dominate over the stochastic type B pathway: Figure S1). Hence, below a certain threshold value of the stimulus x 0 , the stochastic type B pathway plays a role and above that value the deterministic type A signaling leads to rapid activation of x 3 . Note, the observed transition from deterministic to stochastic behavior in signaling, depending on the strength of an external stimulus, confirms that the signaling network under consideration has a suitable minimal structure. However, several simplifying assumptions were made in the calculations based on stochastic differential equations. We use Monte Carlo simulations for a detailed study and also to corroborate the results obtained from our theoretical analysis.

Signaling in the minimal network: Monte Carlo simulations
We use Gillespies Monte Carlo sampling technique [11] to simulate the Master equations that describe the set of chemical reactions given in the equation Eq. (2). The initial condition is set as x 0 1~1 00, x 0 2~5 0, x 0 3~5 0. We consider activation of the downstream signaling molecule x Ã 3 (for k 0 31 %k 0 1 ) as we vary the external stimulus x 0 .
Slow stochastic activation under a weak stimulus. For a weak external stimulus (x 0 v10), large cell-to-cell fluctuations dominate the signaling behavior, where x 3 activation at the single cell level is rapid compared to the time-scale range for x 3 activation for all cells (Figure 3: left panel). The plot was generated for x 0~1 and with reaction constants set to: k 0 31~1 0 {5 (for the type A pathway), and k 0 1~1 :0, k 0 2~1 0 {7 , and k 0 32~1 :0 (for the type B pathway) in suitable units. We derive a condition for which the signaling will be entirely dominated by the type B pathway by comparing k 0 31 x 0 x 0 3 *k 0 2 x 0 1 x 0 2 . For the given parameter values, we obtain x 0 *1 (or less) for pure type B signaling, which agrees well with our numerical simulation results (Figure 3). Decreasing k 0 2 would demand increasingly lower x 0 for the type B pathway to dominate at early times of signaling ( Figure S1). The small rate constant in the intermediate step of the type B pathway leads to slow activation of the pathway where different cells initiate activation at very different time-points and then a subsequent fast rate constant leads to a spike in x 3 activation. We call this phenomenon all-ornone type activation with large cell-to-cell variability. Decreasing stimulus slows down the downstream signaling and increases cell-tocell stochastic variability. We can estimate the variance of x 3 activation at the single cell level using a sharp-rising approximation scheme for x 3 activation (ignoring the spread along the time axis for a single cell). Suppose, m(t) is the number of cells, out of a population N cells, in which x 3 got activated within a given time period t. Assuming rapid activation of x 3 , the average of x 3 is given by vx 3 w~x 0 3 (1{m(t)=N) and the average of x 3 square is vx 2 3 w~(x 0 3 ) 2 (1{m(t)=N). The variance can then be estimated as . Hence a simple scaling relation is obeyed by the Fano factor, i.e. variance(t)/ average(t),~x 0 3 (m(t)=N), which remains &1 and increases with time implying presence of large cell-to-cell stochastic fluctuations. Calculation of the Fano factor from actual simulation data (such as in Figure 3: left panel) also shows that Fano factor at a given time increases with increasing x 0 3 and remains larger than one as predicted by the scaling relation. However, the scaling relation is obtained only within a sharp-rising approximation and thus slightly deviates from the variance/average ratio obtained from simulations.
The stochastic signaling observed in response to a weak external stimulus is essentially dominated by the type B pathway, as a similar signaling response is generated if the kinetic constant for the type A pathway is set equal to zero. This result is also consistent with our theoretical analyses of the minimal network that showed the presence of large stochastic fluctuations under a weak stimulus. Stochastic signaling behavior through the type B pathway in response to a weak stimulus is maintained as long as the fast-slow-fast kinetics of signaling reactions is also maintained in the type B pathway ( Figure S1). For a stimulus of intermediate strength (x 0^1 0{1000), a mixed signaling behavior is observed with a sudden change from type A to type B at the level of single cells (Figure 3: middle panel). Also note, for pure type B signaling (when the kinetic constant for type A is set to zero), presence of a large number of molecules of an initial stimulus (x 0 &1) does not change the observed stochastic signaling behavior. Rapid deterministic activation under a strong stimulus: Under a strong external stimulus, i.e. for large values of x 0 (w1000), the type A pathway dominates the signaling behavior ( Figure 3: right panel), as also observed in the theoretical study using stochastic equations. The plot was generated for x 0~1 000 and with reaction constants set to: k 0 31~1 0 {5 (for the type A pathway), and k 0 1~1 :0, (order of magnitude approximation), which matches well with our stochastic simulation results. Similar deterministic signaling is observed when the kinetic constant for the type B pathway is set to zero to assure a pure type A response. As mentioned earlier, characteristics of such type A signaling is small cell-to-cell fluctuations where the average (over a population of cells) behavior dominates the signaling. As we decrease the strength of an external stimulus increased cell-to-cell stochastic fluctuations are observed with a gradual switch to the type B activation.
Probability distribution approach to characterize stochastic signaling. Large (w1) variance/average ratio for weak stimuli indicates large stochastic variations in signaling. To characterize such stochastic fluctuations in the minimal network we determine the probability distribution of activated x Ã 3 . A histogram of activated x Ã 3 is generated from many runs (i.e. many cells) of our simulations that defines the probability distribution of x Ã 3 . In the case of a strong stimulus x 0~1 000, i.e. when the type A pathway dominates the signaling, the probability distribution shows a gradual shift (along the x axis) over time. Thus the average behavior is representative of the entire cell population and cell-tocell stochastic fluctuations are not significant when compared to the average (Figure 4, top panel). The plot was generated for x 0~1 000 and with reaction constants set to: k 0 31~1 0 {5 (for the type A pathway), and k 0 1~1 :0, k 0 2~1 0 {7 , and k 0 32~1 :0 (for the type B pathway) in suitable units. The initial condition is set as x 0 1~1 00, x 0 2~5 0, x 0 3~1 00. On the contrary, if the stimulus is weak (x 0~1 ), large cell-to-cell stochastic fluctuations lead to a bimodal distribution (Figure 4, bottom panel) for activated x Ã 3 molecules where one peak gets taller while the other peak gets shorter. Such bimodal distribution arises because activation of x Ã 3 in an individual cell is completed in a period of time that is orders of magnitude smaller than the time over which activation occurs for a large population of cells. We define such signaling response as all-or-none type activation with large cell-to-cell variability. Hence, such a bimodal probability distribution can be considered as a characteristic of signaling through the type B pathway with large stochastic variations. For intermediate stimuli, the probability distribution shows characteristics of both type A and type B activation (Figure 4, middle panel).

The minimal network elucidates signaling mechanisms in apoptosis
Apoptotic cell death signaling pathway is one of the most complex intracellular signaling networks and consists of a large number of components [12,13]. Large number of signaling components, molecules, and complicated network structure all contribute to the enormous complexity of the apoptotic cell death signaling pathways, and thus mask the essential design principles of apoptosis signaling. However, one may want to know if there exists a minimal signaling network structure that can capture the essential qualitative features of the cell death signaling through the apoptotic pathways.
Apoptotic cell death signaling is generally mediated by two distinct pathways that are connected in a unique loop structure through the final downstream activation of caspase-3 molecules [13]. Activation of the apoptotic cell death signaling pathway has been shown to have two disparate time scales: (a) fast (minutes) activation in some cases such as ( [14]), and (b) orders of magnitude slower (w10 hrs) activation under certain conditions [5,14]. However, the presence of a specific signaling molecule and its concentrations are often cell-type specific [15], calling into question the robustness of the results obtained from experiments studying apoptosis signaling. Our proposed minimal model of stochastic signaling network can capture the experimentally observed slow activation of the apoptotic signaling pathway in a cell-type independent and robust manner. The slow intermediate step in the type B pathway of the minimal network captures the slow reaction of Bax activation or apoptosome formation (or both of those) in apoptosis signaling. The time-to-death and its cell-to-cell variability in apoptosis signaling is limited by stochastic Bax activation induced cytochrome c release and apoptosome formation induced caspase-9 activation [1,6,7]. The relative contribution of those two slow stochastic events in apoptosis can vary significantly depending on the cell type and apoptoptic stimuli used. We could effectively simulate the effect of both slow steps in apoptosis, and their stochastic variations, by varying the rate constant of the intermediate slow step in the type B pathway of the minimal model ( Figure 5). Thus the proposed minimal network with stochastic signaling behavior seems to be an appropriate minimal model for the apoptosis signaling where large cell-to-cell fluctuations cause slow cell death. This minimal network could also be derived from the full apoptotic signaling network by (i) grouping functionally redundant proteins and (ii) replacing a set of reactions by a single effective reaction with modified rate constants. However, the approach taken in this work has the advantage that prior knowledge of the apoptotic pathway is not needed, instead, a simple set of assumptions on the response of the signaling network is sufficient to generate the minimal network structure.
Results of our minimal model indicate that a slow rate constant followed by a faster kinetics in our minimal model is enough to generate large cell-to-cell stochastic variations as observed in the large scale simulations, as well as in single cell experiments of apoptosis signaling [1][2][3]6,7]. Variations in the initial number of signaling molecules in the type B pathway will enhance the cell-tocell stochastic variability, a feature also observed in apoptosis signaling [7]. We have further connected our minimal model of stochastic signaling with a specific model of apoptosis inhibition mediated by anti-apoptotic Bcl-2 protein. The loop structure of tBid-Bcl-2-Bax along with the Bcl-2 level determine the time course of Bax2 activation and thus apoptotic activation through the mitochondrial pathway of apoptosis [7,13]. The rate constants for the tBid-Bcl-2, tBid-Bax and Bax-Bcl-2 reactions are not small and the number of molecules initially present are also large. However, inhibition by Bcl-2 molecules reduces the effective rate for Bax conversion to activated Bax and subsequent Bax-2 complex formation. In cancer cells, Bcl-2 inhibits over-expressed BH3 molecule such as Bid, which can directly activate Bax with a very small rate constant [7,16], and thus through the Bid-Bcl-2-Bax loop stochastic Bax activation is generated. Such Bcl-2 inhibition of Bax activation dynamically generates stochastic variability in apoptosis even when the initial number of molecules present are not small [9,10]. In our minimal model, we have an effective X ?X Ã reaction with a slow reaction rate constant which can be decreased to simulate an increase in Bcl-2 levels. Recent experimental and theoretical studies elucidated that increase in Bcl-2 levels increases time-to-death and its cell-to-cell variability in apoptosis signaling [7]. For very high Bcl-2 levels, as observed in some cancer cells, apoptotic activation is strongly inhibited with rare stochastic activations of a few cells [7]. In our minimal model, reduction in the rate constant of the slow intermediate step of stochastic pathway generates slower activation with increased cellto-cell variability and thus captures the signaling behavior in apoptosis with over-expressed Bcl-2 levels ( Figure 5).

Discussion
In this article, we proposed a minimal signaling network that is capable of sensing an environmental stimulus and generating an appropriate signaling response in accordance with the strength of the stimulus. We derived this network structure starting with a simple set of assumptions related to the differential response of this signaling network to varying external stimuli. Specifically, under a weak stimulus, the minimal network is designed to use slow stochastic signaling to adapt to a fluctuating environment. We think the proposed minimal network can serve as a minimal signaling model for the complex cell death signaling pathway. Having large cell-to-cell stochastic fluctuations could be a strategy which cells use to respond to a weak environmental stimulus and thus diversify their options for adapting to a fluctuating environment. We hope our results can be tested using synthetic networks and will guide the engineering of biosynthetic mimics designed to sense the environment in the proposed manner. We also used a probability distribution based approach for the characterization of stochastic signaling that shows non-trivial bimodal distribution under weak stimuli.
Large-scale complex signaling pathways, such as cell growth, cell death or immune response signaling pathways, were thought to behave in a deterministic fashion due to the presence of a large numbers of molecules [9,10,17,18]. In this work, we show that stochastic fluctuations can persist even in the presence of a large number of molecules. This result becomes even more important in light of our recent simulations of the apoptotic cell death pathway, where stochastic signaling behavior is also observed. The study of minimal models prove the robustness (cell-type independence) of the results obtained from our Monte Carlo simulation of the apoptosis signaling network and thus can elucidate the essential design principles of apoptotic cell death signaling in normal and diseased cells.
In addition, the proposed minimal signaling network can provide crucial insights into the adaptive evolutionary strategy of complex cell death signaling network. The CED3-CED4 death pathway found in C. Elegans resembles the type B pathway of the proposed minimal signaling network [19]. The large-scale type 2 apoptotic network could have evolved such complexity in response to the needs of increasingly complex and higher level species. Also note, the type 1 apoptotic pathway, which induces fast deterministic activation under strong apoptotic stimuli, first appeared in vertebrates concurrent with the first appearance of adaptive immunity. This type 1 apoptotic pathway resembles the deterministic type A pathway of our minimal network and could have first appeared in order to control adaptive immune response through rapid activation of the apoptotic programmed cell death. Thus the proposed minimal stochastic signaling network could be an evolutionarily conserved network motif and thus could serve as a general model for signaling during apoptosis. Such a minimal model based approach to elucidate design principles of cell signaling pathways is very general and can also be applied to study other large-scale signaling systems. Figure S1 Differential signaling response of the minimal network as k 2 0 and x 0 are varied. This phase diagram is generated using two relations obtained from a parameter sensitivity analysis of the minimal model: (i) type A signaling dominates for all times when