The constructed yeast apoptosis network contains 73 species and 115 reactions (Table 1). Six species do not have successor (Sink/Output) or predecessor (Source/Input; Table 2). In the schematic diagram (Figure 1) species are represented as nodes and reactions as edges. Species include: processes, proteins, and small molecules (metabolites or signals) and Reactions include activation (green arrow) and inhibition (red arrow). Nine nodes depict inputs to the model and are colored blue. Seventeen elements are active in nucleus and are shown in gray, 12 mitochondria species are shown in yellow, and 34 orange boxes represent species residing in cytoplasm. System has only one output node which is called “Apoptosis” and it is in dark blue. Filled blue circles are used as an “AND” gates between two or more reactions to indicate the necessity of presence of two species for activation or inhibition of a reaction (Figure 1).

Network is activated via input signals corresponding to the following species: Acetic Acid, Heat, H₂O₂, Adozelesin, Mg²⁺, Cu²⁺, Salt (NaCl), mating, and osmotic stress.

Like many molecular mechanisms in living cells, apoptosis can be approximated as an outcome of sequential regulation steps which do not occur all at the same time. As a result, upon induction, the state of the cell changes with passing time. In order to capture this feature, each reaction was assigned activation time scale implicating occurrences of different scenarios in sequential order. By introducing time delays to the logical function, it is possible to describe dynamic behavior of the given process using logical networks (Thomas and D’Ari, 1990).

For this study we used Cell Net Analyser (CNA; Klamt et al., 2006) which provides the function to capture signal propagation in a time series to get a snapshot of the network and discrimination of signaling events. This approach has been successfully used to model human apoptotic network revealing new feedback loops (Schlatter et al., 2009). To describe dynamic behavior of the yeast apoptosis five timescales t = [0, 2, 4, 5, 6] are assigned. These timescales are constants and indicate at which timescale each node gets activated or inhibited. Simulation of the network at t = x comprises of all interactions with timescale t ≤ x but not interactions with t ≥ x. This gives the possibility of separation of different functional groups such as different signaling pathways or feedback loops. It should be noted that specified timescales do not refer to real time but are only indicators of sequential regulatory steps (difference between timescale t = 0 and t = 2 is equal to difference between timescales t = 2 and t = 4 or t = 4 and t = 5 or between t = 5 and t = 6). Timescales are not assigned based on speed of reactions (how fast or slow reactions are), they are assigned based on the sequence of events. To allow flexibility in our model and facilitate changes and insertions of new events timescales 1 and 3 are reserved and not used in the current network structure. First timescale t = 0 is assigned to genes which are already present and constantly active in cell (27 reactions). Values for all input arcs to stimuli species are set to 1 at second timescale t = 2 (53 reactions). Further interactions are activated at time t = 4 (17 reactions). Although yeast apoptosis does not contain any feedback loops time point t = 5 is reserved for feedback loops which get activated in response to osmotic and heats shock and will lead to apoptosis. Considering impact of feedback loops on system behavior it is reasonable to assign them a separate timescale. Finally interaction occurring at the very end are assigned timescale t = 6 (18 reactions). Details on each reaction, their time points and species involved in reactions are in Table 3.

The model simulated the induction of apoptosis both in the independent mode (assessing each individual stimulus separately) as well as in the additive mode, where the activation of all inputs were set at the same time (Table 4).

Based on qualitative knowledge we have constructed a discrete Boolean model of yeast apoptosis. While it can capture its essential behavior, the question remained how this model can be used to predict the qualitative behavior of the system. In order to address this question we expanded the model by transforming the discrete Boolean model into a continuous model using Odefy (Wittmann et al., 2009; Krumsiek et al., 2010; see Materials and Methods). In order to test the continuous model and its predictive capacity, we performed three independent case-studies: (i) induction of apoptosis by activation of Hog1p by heat stress, (ii) inhibition by Bir1p of acetic acid-induced apoptosis, and (iii) induction of apoptosis with H₂O₂ by activation of Stm1p. We show here that the behavior that emerges from specific interactions in the model is in agreement with published experimental data.

We simulated the activation of the mitogen-activated protein kinase that is a key component of the HOG pathway (Albertyn and Hohmann, 1994; Van Wuytswinkel et al., 2000; Hohmann, 2002; de Nadal et al., 2004), which is an osmoregulatory signal transduction cascade (Hohmann, 2009). Upon stress, Hog1p (encoded by HOG1/YLR113W) regulates the expression of almost 600 genes by phosphorylating several different transcription factors (Hohmann, 2002; Westfall et al., 2004). The activity and nuclear localization of Hog1p is regulated by its phosphorylation state, and that in turn is regulated by the kinase MAPKK Pbs2p and the phosphatases Ptc1p, Ptc2p, Ptc3p, Ptp2p, and Ptp3p (Brewster et al., 1993; Ferrigno et al., 1998; Warmka et al., 2001; Young et al., 2002). Besides the induction by osmotic stress, the HOG pathway can be induced by heat stress, via a Sho1p-dependent sensory mechanism (Winkler et al., 2002), thus we used heat stress as an input to activate the HOG pathway and simulate the cell death response. Heat stress activates Slt2 which then activates Rlm1. Consequently Rlm1 triggers expression of Slt2 and PTP2 forming a feedback loop. On the other hand heat shock inhibits activity of PKA leading to release the inhibition effect of PKA on MSN2-4. MSN2-4 activates Sdp1 which then along with PTP2 inhibit Slt2. We started by setting all other stimuli to zero and then by transforming the Boolean apoptosis model into Hill Cube continuous model we performed the simulations. Steady states predicted from continuous Hill Cube (Figure 2) and synchronous Boolean (Figure 3) model are in perfect agreement with each other for all nodes apart from the nodes corresponding to the feedback loops in response to heat and osmotic stress. Feedback loop in heat activated pathway includes activation of Rlm1 by Slt2 and expression of Slt2 by Rlm1. Osmotic shock induced Hog1 activates Rlm1 to regulate Slt2 which was inhibited by osmotic shock and PTP2. Unknown Hog1 dependent transcription factor triggers transcription of PTP3 and PTP2 resulting in dephosphorylation of Hog1 phosphotyrosine which inhibits Hog1 activities. All nodes in our model were connected to activation of apoptosis via the heat shock pathway and were indeed activated in our simulation and are colored in red; other pathways are colored in blue indicating that they are not activated during the simulation (Figure 2). Model predicts that upon heat induction, concentration of Hog1p changes trough time never reaching its maximum level of concentration and as an intensity of stimulus decreases level of Hog1p also decreases. Our simulations also proposed that the maximal level of activation of Hog1p during heat stress is 70% of the total activation and as the heat stimulus continues over time its activity decreases and reaches a plateau at 40% of total activation (Figure 4).

This model prediction is in agreement with an experimental measurement performed by Winkler et al. (2002) which shows the induced activity of Hog1p upon heat induction (Figure 4). The difference between the maximal activity values and the lag phase of the activation is due to the fact that we did not use any quantitative data as input to the model.

Even though there are several genes and pathways that are involved in yeast pro-apoptotic response, there are only few anti-apoptotic genes. For many years Bir1p (encoded by BIR1/YJR089W) was thought to be the only apoptotic inhibitor in yeast, but recently a TSC22-protein family was found and two of the proteins with TSC22-motif have been shown to have anti-apoptotic roles Sno1p and Fvy10p (Khoury et al., 2008). Bir1 belongs to the IAP family and phylogenetic analysis (Uren et al., 1998) showed similarity to Schizosaccharomyces pombe BIR1, human survivin, and Caenorhabditis elegance BIR-1 and BIR-2 proteins. Yeast’s Bir1p has been previously intensively studied in chromosome segregation (as a component of the Aurora kinase complex (chromosome passenger complex, CPC; Ruchaud et al., 2007) but recently studies of its role as a negative regulator of apoptosis have gained momentum. It has been shown that Bir1p is a target for degradation by Nma111p (Owsianowski et al., 2008), an apoptotic serine protease and yeast cells lacking Bir1 are more sensitive to apoptosis, while overexpression of Bir1 reduces apoptosis (Walter et al., 2006).

When building the model we assumed that addition of acetic acid induces cytochrome c release from mitochondria and its translocation to the cytosol. Another assumption is that the execution of apoptosis is Yca1p-dependent (encoded by YCA1/MCA1/YOR197W) in this context and that the activation of the caspase is downstream from the cytochrome c translocation (Ludovico et al., 2001, 2002). Thus, taking these assumptions into account acetic acid was used as an inducer in the simulations with direct induction of apoptosis. The model confirmed experimental evidence that Bir1 is indeed apoptosis inhibitor in yeast (Walter et al., 2006) and as acetic acid is applied as a pulse stimulus (that is then decreasing over time) and cytochrome c present in the cytosol increases, the decrease of Bir1p due to Nma111p degradation promotes apoptosis that then reaches the maximum (Figure 5A). We then performed an artificial and biologically irrelevant simulation: after apoptosis has occurred and no more inducer was present (acetic acid is not added again and the previously added amount has been “used”) the eventual theoretical accumulation of Bir1p did inhibit apoptosis and revert it to zero. Obviously this scenario is biologically implausible, since once the cell has undergone apoptosis it cannot be revived, but has a purpose to show mathematical correctness of the developed model.

Additionally, we have converted discreet apoptosis model to continuous model with constant presence of Bir1 in order to test if this type of model would predict the same outcome as the experimental approach in which Bir1p is overexpressed and provides protection from induction of apoptosis. The continuous model was created taking into account translocation effect of each element from one compartment to another. This effect has been applied to translocation of Nuc1 from mitochondria to nucleus, MCD1 from nucleus to mitochondria, Stm1 from cytosol to nucleus Aif1 (from mitochondria to nucleus, NMA111 from cytosol to nucleus, Ste20 from cytosol to nucleus, and cytochrome c from mitochondria to cytosol. Initially, in this model, translocation of each species is modeled as one node in order to represent the re-localization which introduced self-regulatory loops to the system which is impossible since to our knowledge yeast apoptosis regulatory network does not have such loops. One solution implies intuitive approach in defining single variable for single species to mimic the effect of transferring from one compartment to another. Technically this was solved by introducing two variables for a single species representing both compartments they can belong to. This approach, also, confirmed experimental finding (Walter et al., 2006) that constant presence of Bir1, inhibits apoptosis, validating Bir1 anti-apoptotic role (Figure 5B).

It has been shown that degradation of short-lived pro-apoptotic proteins via proteasomal – ubiquitine pathway plays important role in mammalian apoptosis (Drexler, 1997). One of the yeast proteasomal substrates – Stm1 (YLR150W) was identified in the study by Ligr et al. (2001) as an activator of the cell death process. Conserved orthologs of Stm1 are detected in several species (highly conserved ortholog in Schizosaccharomyces pombe and a putative ortholog in Drosophila melanogaster (Nelson et al., 2000) suggesting that regulation of apoptosis via this protein can be evolutionary conserved process.

We have in silico tested two scenarios confirming the role of Stm1 in apoptosis (Ligr et al., 2001). Accumulation of Stm1 will induce apoptotic behavior followed by DNA fragmentation which is a known marker of cell death, while the Stm1 knock out will promote survival and consequently deterring DNA fragmentation (Figures 6A,B).

S. cerevisiae is a model organism that has conserved genes, proteins and pathways that are similar to the ones in the human cells. This allows for using yeast as a host in which human genes can be expressed, so called “humanized yeast” (Winderickx et al., 2008) and subsequently the physiological roles and molecular mechanisms can be studied. In order to test if our system could be used in a similar way (in silico humanized yeast) and provide reliable simulations and predictions, we inserted three human genes with complete downstream pathways into our initial yeast apoptosis Boolean model. Yeast apoptotic network was “humanized” by in silico insertion of genes belonging to Bcl-2 protein family and Valosin-contain protein (VCP) – distinctive representatives of human apoptosis (Figure 7).

This protein family consists of apoptotic agonist and anti-agonist proteins. Among members of Bcl-2 family Bax as an apoptotic inducer and Bcl-xL as an anti-apoptotic factor were included in our extended model. Bax is located in cytosol in mammalian cells and has vital role in mitochondria morphogenesis. Heterologously expressed Bax causes growth arrest and rapid cell death in S. cerevisiae (Sato et al., 1994; Greenhalf et al., 1996). Also expression of Bax has been linked to the release of cytochrome c from mitochondria (Manon et al., 1997). Yeast cells with Bax expression accumulate ROS and show other apoptotic hallmarks such as DNA fragmentation (Madeo et al., 1999).

Since members of this protein family were previously successfully expressed in yeast, we first validate the “humanized” yeast model by inserting Bcl-xL, Bax, Bad, and Bcl-2 genes (Sato et al., 1994; Greenhalf et al., 1996; Madeo et al., 1999). As an input signal human Akt signaling was used (Akt Signaling), as an extrinsic regulatory switch (since this signaling cascade is not present in yeast). The Akt Serine/Threonine-kinase promotes cell survival by phosphorylating the pro-apoptotic protein BAD (member of the Bcl-2 family), which is the cause of dissociation of BAD from the Bcl-2/Bcl-X complex, and promotion of cell survival (Datta et al., 1997). Besides cell survival, Akt signaling is related to the cell cycle, metabolism, and angiogenesis and therefore a target for anticancer drug development (Falasca, 2010; Hers et al., 2011).

Upon activation of Akt Signaling at the first time step (t = 0, node value Akt Signaling = 1), Bcl-2 and Bcl-xL are activated in the following step (t = 4, node value Bcl-xL = 1 and Bcl-2 = 1) and remain active throughout the simulation and as the inhibitors of apoptosis promote the survival (node value Apoptosis = 0; Table 5). Simulation results suggest the anti-apoptotic role of two Bcl-2 family members: Bcl-2 and Bcl-xL which is consistent with experimental findings (Kharbanda et al., 1997).

The evolutionary conserved Valosin-containing protein (VCP) is a mammalian ortholog of yeast Cdc48, which is the first apoptotic mediator found in S. cereivisae (Braun and Zischka, 2008). VCP is the member AAA-ATPase family which is ubiquitously expressed (Braun and Zischka, 2008). It consists of four domains: N-terminal domain, two ATPase domains D1 and D2 and C-terminal domain (Wang et al., 2004). The major ATPase activity of VCP is carried on by D2 domain (Wang et al., 2004; Song et al., 2003). VCP is involved in different cellular process such as protein degradation, membrane fusion and chaperone activity (Braun and Zischka, 2008). Role of VCP/Cdc48 in fluctuating number of death cell in various types of disease naming cancer and protein deposition diseases is not well understood. Increase in expression of VCP is correlated to the development of cancer and metastasis therefore detecting the level of VCP expression is proposed as cancer progression marker. Also VCP is known as detector of aggregated proteins causing neurodegenerative disease such as Parkinson and Alzheimer (Hirabayashi et al., 2001; Mizuno et al., 2003; Ishigaki et al., 2004). It has been observed that yeast cell carrying mutation of Cdc48 is showing morphological markers of apoptosis (Madeo et al., 1997) and that depletion of this protein has apoptotic effect in other organisms as well (Imamuraab et al., 2003).

It had been shown that VCP can both promote or inhibit apoptosis (Braun and Zischka, 2008). Deletion of VCP triggers ER stress which is followed by unfolded protein response and cell death via ER-associated degradation (ERAD) pathway – a highly conserved pathway between mammalian and yeast cells (Braun and Zischka, 2008). Wild type VCP has pro-apoptotic role in the cells undergoing apoptosis upon ER stress which in turn triggers ER-associated apoptotic pathway. Finally VCP can trigger survival pathway in response to NFkB which is a pro-survival molecule in cells with over expressed level of VCP/Cdc48(Braun and Zischka, 2008).

The dual functionality of VCP in survival and apoptosis mechanism makes it an attractive candidate to be inserted in initial yeast apoptosis Boolean network. Since VCP includes expression of a handful of heterologous proteins, its in vivo insertion to yeast plasmid would be difficult to perform and thus represents a good candidate for further investigation and exploration of model capabilities.

In our model, yeast Cdc48 was replaced with human VCP gene and its downstream pathway effector caspases (caspases 3, 6 and 7), initiator caspases (caspase9 and 12), IAP family, IkB-alpha inhibitory protein which when degraded by proteasome cause the release, and nuclear translocation of active NFkB (represented in the model as NFkB_Cyt and NFkB_Nuc) and gp130/Stat3 pathway (represented as a single node).

Upon activation of Akt Signaling pro-survival role of VCP is observed (node Survival = 1 at t = 6, Table 3). This is achieved when IkBα is dissociated from NFkB (node NFkB/IkBα = 1 at t = 4), which in turn gets activated (node NFkB = 1, t = 5) and is translocated to the nucleus (nodes NFkB_Cyt = 1 and NFkB_Nuc = 1, t = 6) inducing survival of the cell (node Survival = 1, t = 6). Simultaneously, activation of NFkB, promotes activity of IAP family of proteins which inhibits the activity of both effector and initiator caspases (nodes c-3-6-7 = 0 and c-9-12 = 0, t = [0,4,5,6], thus disabling apoptosis-dependent caspase pathway. As an independent, but parallel process, activation of VCP via Akt Signaling activates gp130-Stat3 (node gp130-Stat3 = 1, t = 4) pathway and consequently leads to survival (Table 6).

This result suggests that yeast carrying mammalian VCP gene exhibits the same behavior as it is known from the human apoptotic model (Vandermoere et al., 2006).

Based on the extensive literature study and databases search (Saccharomyces Genome Database-SGD; Cherry et al., 2011) the apoptosis network consisting of 73 nodes and 115 edges is constructed. These nodes are selected based on their interactions and their substrates involved in apoptosis. Cell Designer (Funahashi et al., 2003) was used to visualize genes and their connection. Cell Designer allowed us to represent network using comprehensive graphical notation. Moreover, Cell Designer is able to connect to online data bases such as KEGG (Kanehisa and Goto, 2000), BioModels (Li et al., 2010) and PubMed which expands its connectivity, visualization and model building. The most prominent advantage of Cell Designer is its ability to support SBML (Systems Biology Markup Language; Hucka et al., 2003) format. Creating SBML file avoids creating logic rule file which is more error prone when rules have to be defined by user in simple text and it is not in form of easily drawing boxes.

Identification of Logical Steady States (LSSs) in a Boolean network is an important task as they comprise the states in which a gene-regulatory network resides most of the time. Strong biological implication can be carried out by LSSs of the network. LSSs can even be linked to phenotype (Kauffman, 1969). CNA (Klamt et al., 2006) is used to calculate all possible LLSs based on specified initial value for each gene and signal flow in network. LSSs are used to evaluate the network behavior under perturbation and changes in network structure.

Considering a pair of nodes (a, b) dependency is defined as the influence of node a on node b and vice versa, node a influences node b as follow: a is a total activator of b if there is an activating path between a and b, a is a total inhibitor of b is there is an inhibition path from a to b or a is an ambivalent effecter of b meaning that there is an intermediate node which is involved in a negative feedback loop. Using CNA the interaction matrix and dependency matrix were drawn from the network (Figures 8 and 9).

In order to convert the schematic network into Boolean model we used SQUAD (Mendoza and Xenarios, 2006), a user friendly graphical software which is suitable for modeling signaling network where kinetic reactions are not available. Simulation in SQUAD consists of three steps: (1) the network is first loaded from the SBML file. Components of the network are presented as nodes and value of each node represents the state of that node. SQUAD converts the network to discrete dynamic model. Using Boolean algorithms all steady sates in network are calculated. (2) network is converted to a continues dynamic model generating sets of ordinary differential equations (ODEs) and steady states achieved in pervious step. (3) SQUAD allows perturbation in order to understand the role of each node within the network.

Using Reduced Order Binary Decision Diagram (ROBDD), it is possible to calculate steady states (Di Cara et al., 2007). ROBDD is a memory efficient data structure which is widely used in electronic field and has been proven to work for large binary networks. Moreover, ROBDD computes steady states for large networks (n > 50) in matter of seconds. Another advantage of this algorithm is its ability to identify the cyclic steady states. These oscillating states are reachable when system identifies a cyclic pattern instead of one single state.

In the first step, schematic network is supplied to SQUAD resulting in a discrete network, generating a set of either cyclic or single steady states. These steady states are then used to convert discrete model to continues. Given all calculated steady states or cycles we examined all possible outcomes for each input. Depending on the desire form of output, discrete to continuous conversion can be carried out either as a complete or progressive mode. Our simulations are performed in the complete mode since cell undergoing apoptosis should die after certain time and it is expected to maintain a constant level of apoptosis or survival (in case when apoptosis is not activated) at the end of each run. As an opposite to complete mode, the progressive model allows the user to stop the simulation at any time even before reaching the steady states.

Reconstructing Boolean model of yeast apoptosis from qualitative knowledge never gives details about concentration of molecules in different time points. For this purpose the discrete Boolean model is transformed to continuous model using Odefy. Odefy uses the multivariate polynomial interpolation in order to transform the logical rules into sets of ODEs. Yeast apoptosis Boolean model is converted to continuous model using Hill Cube and normalized Hill Cube where the Hill function is normalized to the unit interval. Behavior of biochemical reactions can be seen as a sigmoid Hill function represented as f(x̄̄)=x̄̄n∕(x̄̄n+kn). Where, n is a Hill coefficient and determines the slope of the curve and is a measurement of cooperativity of the interactions, and parameter k corresponds to values 1 and 0 in the Boolean model in the following manner: threshold value above given k in Boolean model is considered as 1 (on state) and below as 0 (off state).