Skip to main content
Download PDF

Simulating Quantitative Cellular Responses Using Asynchronous Threshold Boolean Network Ensembles

BMC Systems Biology volume 5, Article number: 109 (2011) Cite this article

Abstract

Background

With increasing knowledge about the potential mechanisms underlying cellular functions, it is becoming feasible to predict the response of biological systems to genetic and environmental perturbations. Due to the lack of homogeneity in living tissues it is difficult to estimate the physiological effect of chemicals, including potential toxicity. Here we investigate a biologically motivated model for estimating tissue level responses by aggregating the behavior of a cell population. We assume that the molecular state of individual cells is independently governed by discrete non-deterministic signaling mechanisms. This results in noisy but highly reproducible aggregate level responses that are consistent with experimental data.

Results

We developed an asynchronous threshold Boolean network simulation algorithm to model signal transduction in a single cell, and then used an ensemble of these models to estimate the aggregate response across a cell population. Using published data, we derived a putative crosstalk network involving growth factors and cytokines - i.e., Epidermal Growth Factor, Insulin, Insulin like Growth Factor Type 1, and Tumor Necrosis Factor α - to describe early signaling events in cell proliferation signal transduction. Reproducibility of the modeling technique across ensembles of Boolean networks representing cell populations is investigated. Furthermore, we compare our simulation results to experimental observations of hepatocytes reported in the literature.

Conclusion

A systematic analysis of the results following differential stimulation of this model by growth factors and cytokines suggests that: (a) using Boolean network ensembles with asynchronous updating provides biologically plausible noisy individual cellular responses with reproducible mean behavior for large cell populations, and (b) with sufficient data our model can estimate the response to different concentrations of extracellular ligands. Our results suggest that this approach is both quantitative, allowing statistical verification and calibration, and extensible, allowing modification and revision as guided by experimental evidence. The simulation methodology is part of the US EPA Virtual Liver, which is investigating the effects of everyday contaminants on living tissues. Future models will incorporate additional crosstalk surrounding proliferation as well as the putative effects of xenobiotics on these signaling cascades within hepatocytes.

Background

Motivation

Thousands of chemicals are used in commerce and evaluating their public health risk remains a challenging problem [1, 2]. Much of our knowledge about mechanisms of toxicity is based on evidence from in vivo animal studies and in vitro experiments, where we can begin to unravel some of the molecular signaling and transcriptional changes induced via chemical perturbation; however, there are three main issues in translating these findings to humans. First, it is often impractical to design experiments with sufficient power to detect the subtle effects of very low environmentally relevant exposures [3]. Hence most observations about chemical effects are at relatively high concentrations that cannot be easily used to quantify the risk of long-term and low-dose exposure to complex mixtures of chemicals [4]. Second, since the molecular response to chemicals is not always conserved between species [5], the effects observed in rodents cannot be directly extrapolated from rodents to humans without additional mechanistic insight [6]. Third, toxicity is a tissue level phenomenon arising from the behaviors of heterogeneous cell populations. Understanding the complex signaling processes between these different cell types is crucial in determining toxicity potential. We are building a cell-based tissue model to estimate the quantitative population-level effects of chemical exposures [7, 8]. Here we describe an asynchronous threshold Boolean network (BN) approach to model signal transduction in individual cells and to estimate tissue level responses using an ensemble of BNs.

Boolean Networks

A BN describes a signaling network as a digital circuit in which logical elements (proteins or genes) are either 'ON' or 'OFF'. The temporal evolution of the signaling network is calculated using a set of Boolean functions (AND, OR, NOT) to model regulatory interactions. Since they offer a biologically relevant and computationally efficient formalism for analyzing the relationship between molecular network topology and function, BNs have been used extensively to simulate the behavior of cells based on their network activity. Genetic regulatory networks have been particularly amenable to this formalism due to the binary nature of gene activation [9]. The availability of large-scale transcriptional profiles spurred more recent applications of deterministic [10] and probabilistic [11] BNs for reconstructing and simulating genetic regulatory networks. Additionally, BNs have been used for modeling the cell cycle [12–15]; cell proliferation [16–18]; apoptosis [19]; and cancer [20]. BNs can be used to represent the binary activity of molecular species across cell populations (in vitro and in vivo). One of the limitations of BNs is that they cannot readily estimate continuous functional responses, i.e., quantitative dose-response, which are fundamentally important in pharmacology and toxicology.

Cancer Pathways

Liver cancer is a frequent outcome in testing the long-term effects of chemicals in rodents [21] and often considered in regulatory decisions [22]. Since the mechanisms of carcinogenesis are poorly understood, it is difficult to translate chemical effects from rodents to humans. Cancer is believed to arise due to a breakdown of the homeostatic processes that maintain balance between cell death and division [23]. Some chemicals (called mutagens) can alter cell phenotypes by damaging DNA resulting in harmful mutations that can spur the activation of oncogenes. Nongenotoxic carcinogens, on the other hand, can act via insidious mechanisms that suppress apoptosis or to stimulate cell proliferation.

It has been suggested that nongenotoxic carcinogens may increase hepatocyte proliferation by perturbing the crosstalk network regulated by growth factors and cytokines [24]. Crosstalk refers to the sequence of protein regulation activated by any one growth factor or cytokine ligand overlapping with the sequences of other ligands, which allows for complex behavior. The presence of crosstalk allows a cell to behave as a multiplexer, integrating multiple signals to select from many possible outcomes, such as cell cycle initiation and progression.

A number of computational models have been proposed for simulating cell proliferation [12–18, 25], however, BNs have not been extensively used in modeling chemical induced toxicity or in hepatic biology. In order for a chemical to produce a chronic or acute tissue level effect, there must be some level of perturbed protein activity in the signal transduction of one or more cells. We are evaluating BNs for modeling early molecular signaling events in hepatocytes that lead to proliferative changes, which are key events in carcinogenesis. Hence, our initial objective is to model some of the normal cues, i.e. homeostatic processes, that can stimulate healthy, quiescent hepatocytes (G0) into entering the cell cycle (G1).

Technological advancements such as flow cytometry and high content screening have made it possible to measure protein levels with single cell resolution. Experimental observations suggest that protein levels within cells may exhibit a switch-like 'all' or 'nothing' ('ON' or 'OFF') response - for example, p53 response to DNA damage [26, 27], TNFα stimulation [28, 29], MAPK signaling events [30], and drug treatment [31]. These types of observations serve as a foundation for the hypothesis that a Boolean representation is sufficient for describing the molecular multiplicities of individual cells in a simulation framework. Next, we assume that the aggregate activity of molecules across a population of hundreds, thousands or millions of cells can be used to estimate tissue level responses.

Key Contributions

Our work is based on two extensions of asynchronous BNs, which employ a non-deterministic updating scheme. First, we use threshold functions to calculate the activation of each protein in our model. This technique has been applied to other systems [32, 33], and it provides a simple representation and adjustable parameter for investigating the interactions between signaling molecules. Second, we model an ensemble of BNs to simulate the quantitative responses of thousands of cells. As such, we can estimate dose dependent responses of cell populations. We defined the topology of the BN semi-automatically using structured information about canonical signaling network from a public pathway repository. Here we describe our initial results on the reproducibility of asynchronous threshold Boolean network ensembles and their potential utility for estimating quantitative time- and concentration-dependent biological responses.

Results

Cell Signaling Network Construction

We used the Science Signaling database (or STKE) [34–37] to construct the protein signaling network. The canonical pathways in the network include: Epidermal Growth Factor (EGF) signaling, Insulin (INS) signaling, Insulin like Growth Factor type 1 (IGF-1) signaling, and Tumor Necrosis Factor alpha (TNFα) signaling. The number of proteins and molecular interactions in each of these pathways are summarized in Table 1. We performed several steps to systematically build a crosstalk network from these canonical pathways. First, we combined all of the proteins and interactions from the four pathways into one integrated molecular interaction network. After filtering for uniqueness among proteins and interactions, we produced a non-redundant crosstalk network with 102 proteins and 150 interactions. Second, we excluded the proteins and interactions that did not lead to c-Jun and c-Fos activity, which are important components in the formation of the activator protein 1 (AP-1) transcription factor complex.

Table 1 Pathway Information for Building Crosstalk Model
Full size table

In our initial model, we focused on early signaling events in cell proliferation and did not consider gene expression changes which lead to mitosis. Hence, we assumed that AP-1 formation, encoded as a c-Jun/c-Fos dimer, is an early marker of cell cycle progression. This allowed us to further simplify the network by removing all proteins and interactions that are not on a pathway from one of the four receptors to either c-Jun or c-Fos. Furthermore, we manually removed an additional six proteins with in degree less than 2. We did, however, leave some proteins with in degree less than two: the extracellular ligands and their receptors, as well as Rat Sarcoma (RAS), ribosomal s6 kinase (RSK), v-erb-b2 erythroblastic leukemia viral oncogene homolog 2 (ErbB2) and homolog 3 (ErbB3), Rho GTPase (RHO), p55 gamma (p55g), Vav proto-oncogene (VAV), c-Jun, mitogen- and stress-activated protein kinase 2 (MSK2), mitogen activated protein kinase kinase (MAPKK), phosphoinositide-dependent kinase 1 (PDK1), and 1,2-Diacylglycerol (DAG). These molecules are implicated in the EGF signaling pathway, which was simulated and compared to in vitro data, except p55g which is involved INS pathway. Finally, we included in the model a molecular species representing AP-1 transcription factor complex formation, by adding two additional interactions involving the c-Jun and c-Fos dimerization. The final biochemical interaction network contained 46 proteins and 77 interactions. The protein signaling network in Figure 1 was drawn with Cytoscape [38], an open source tool convenient for visualizing large scale networks.

Figure 1
figure1

Putative Crosstalk Network for Simulation. The network shows signaling interactions due to extracellular ligands including: EGF, IGF-1, INS, and TNFα. The pathways share common adapter proteins (e.g., SHC and GRB2) as well as downstream signals (ERK1/2 and c-Fos/c-Jun activation).

Full size image

Simulating individual cellular responses

We used the biochemical interaction network in Figure 1 to describe the response of an individual hepatocyte to the growth factors (EGF, IGF-1 and INS) and the inflammatory cytokine (TNFα). In order to simulate the dynamics of signal transduction, we translated the biochemical interaction network into a threshold BN. As in a traditional BN approach, we assumed that: (a) proteins in the network are described by one of two states, active (ON) or inactive (OFF) and, (b) interactions result in either the activation or inhibition of output proteins by input proteins. Our approach deviates from traditional BNs in three important ways. First, we replace the logical operators with an integer summation function that incorporates an activation threshold. This allows us to adjust the activating potential of each protein in the network. Second, we employ a nondeterministic, asynchronous updating scheme (see Methods), which can simulate biological 'noise' observed in protein signaling cascades. Third, we provide a method for using Hill functions for calibrating the probability of activation for proteins in the network, which can be calibrated with concentration-response data.

In our methodology, the signaling network in a single cell is represented as one asynchronous threshold BN. Figure 2(A) illustrates the model of an individual cell as a BN and its discrete dynamic response following treatment with INS. The BN is constructed from the model shown in Figure 1 (see Methods). The temporal evolution of protein activity in one BN is visualized as a heatmap in Figure 2(A) (right panel). The abscissa of the heatmap shows the simulation time steps (denoted as τ). At τ = 0 the cell was 'treated' with INS by switching the ligand from OFF to ON. Each column in the heatmap shows the dynamic changes in the state of proteins (given in the ordinate) at time steps following INS treatment. The simulation continues until (τ = 369) when it reaches a steady state, which involves the activation of the AP-1 transcription factor complex. The discrete profile of each protein shows the asynchronous dynamics of signal transduction through the insulin receptor (IR) in our crosstalk model (shown in Figure 1). In other words, Figure 2a depicts a putative sequence of signaling events that could occur in a single hepatocyte after insulin stimulation. This result is qualitatively concordant with in vitro observations on AP-1 formation following insulin treatment [39]. Since the output of a single BN is binary, however, it is difficult to evaluate the activation of AP-1 to different concentrations of INS or other ligands for a single cell.

Figure 2
figure2

Overview of Boolean Network Ensembles. (A) A single synthetic cell as a Boolean network (left) along with discrete dynamic response of proteins in one simulated cell to INS treatment (right). (B) The aggregate activity profiles of four proteins (right) across 1000 simulated cells (a replicate) in response to INS treatment. The abscissa shows time steps and the simulation terminates when all cells have reached steady state.

Full size image

Simulating cell population responses

In order to estimate the quantitative response to treatments, we assume that cell populations can be modeled as an ensemble of asynchronous BNs. This allows us to estimate the dynamic response across a simulated biological sample as the aggregate activity of each protein across thousands of BNs (see Methods). Hence, an ensemble of BNs can be considered a simulated 'replicate' as illustrated in Figure 2(B). We investigated the response of an ensemble of 1000 asynchronous BNs to treatment with INS (including the BN depicted in Figure 2(A) until all BNs reached a steady state. The resulting aggregate activity profiles of IRS, c-Jun, c-Fos, and ERK1/2 are shown in Figure 2b (right panel). These trends captured by the simulated BN ensembles appear similar to experimental data from molecular assays performed on in vivo and in vitro replicates (which contain a large number of cells). While this requires additional quantitative and mechanistic evaluation, it is important to note that such continuous protein activity profiles cannot be generated using traditional BNs. Before further evaluation with experimental data we analyzed the reproducibility of our approach with respect to the network depicted in Figure 1.

Reproducibility of Protein Activity Profiles

We systematically evaluated the reproducibility of the asynchronous BN ensembles of the model shown in Figure 1 by analyzing their response to different treatment conditions. For each treatment condition, we simulated 100 replicates with 1000 cells per replicate (i.e., 100,000 cells per treatment condition). Each treatment condition is defined by combinatorial stimulation of the four extracellular ligand cues: (i) EGF, TNFα, IGF-1, and INS combined; (ii) TNFα and INS; (iii) TNFα and IGF-1; (iv) TNFα and EGF; (v) EGF only; (vi) IGF-1 only; (vii) INS only; (viii) TNFα only; (ix) and no active extracellular ligands (the control group). We assumed that each cell is exposed to enough ligand in order to activate a sufficient number of receptors for signal propagation. Hence, for each of the simulated treatment conditions 100% of the cells receive stimulation. Moreover, following the logic of Boolean abstraction of protein concentration, we assumed that the ligand is switched 'ON' in every cell upon initialization.

Figure 3(A) shows the dynamic responses of the simulated replicates in the treatment group (i), the combined stimulation of all extracellular ligands. Each of the 12 plots shows the activity profile for one protein from a random sampling of eight replicates. Even though the activity profile of each replicate is noisy, the overall trend across the eight replicates in the (i) treatment group appears to be reproducible. To analyze this further, we calculated the coefficient of variation (CV) for every treatment group (see Methods). These results are summarized in the heatmap in Figure 4(A). The rows of the heatmap correspond to the treatment groups and the columns to proteins in our model. Each cell shows the normalized CV across all proteins and treatment where increasing color saturation signifies decreasing reproducibility. For instance, the simulated treatment with all ligands produces highly reproducible changes in steady state protein activities, whereas there is considerable variation in the absence of any treatment. Overall, the protein activation across replicates is generally reproducible and well within the range of actual experiments [40].

Figure 3
figure3

Simulated Replicate Protein Activity Profiles. The protein activity profiles after treatment with all ligands for 1000 simulated cells. The graphs show the time course response of a subset of proteins in the model including: (A) receptor tyrosine kinases, (B) select intermediate signaling proteins, and (C) select transcription factor proteins including c-Fos and c-Jun. For instance, EGFR shows an initial increase in activity followed by inactivation due to feedback inhibition, which is consistent with our knowledge of receptor internalization and ubiquitination. We have not included similar feedback for the other receptor in the current model. Similarly, IRS is a hub signaling protein with a number of potential inhibitors and some of these are represented in our model. The activation of the MAPKs (ERK1/2) and transcription factors (AKT, c-Fos, c-Jun) in our simulation highlights the putative signaling cascade responsible for activating immediate early genes, which is a key step in cell cycle progression into S phase.

Full size image
Figure 4
figure4

Reproducibility of Protein Activities Across Simulated Replicates. (A) The heatmap shows the coefficient of variation for protein activity (columns) following different treatments (rows). The color intensity increases with decreasing reproducibility. The nine treatments shown in the rows are (i) EGF, TNFα, IGF-1, and INS combined; (ii) TNFα and INS; (iii) TNFα and IGF-1; (iv) TNFα and EGF; (v) EGF only; (vi) IGF-1 only; (vii) INS only; (viii) TNFα only; (ix) and no active extracellular ligands (Control). Each cell in the heatmap represents the coefficient of variation for the activity of one protein across 100 replicates with 1000 cells/replicate (a total of 100,000 Boolean network simulations). (B) The distribution of CV across each treatment condition. (C) The distribution of CV across each protein.

Full size image

In Figure 4(B), we show the distribution of the CV for the steady state protein activities across all treatment groups as a box and whisker plot. Similarly, Figure 4(C) shows the reproducibility across the proteins for each treatment group. Whereas the heatmap of Figure 4(A) shows information on the CV per protein per treatment condition, the plots in Figures 4(B) and 4(C) visualize the overall behavior of the model across each treatment condition and protein, respectively. We found that treatment conditions (iii), (vi), (viii) and (ix) - IGF-1 and TNFα, IGF-1 only, TNFα only, and the untreated control group - were the least reproducible in comparison to all other treatment conditions. For the control group, a possible explanation for the reproducibility result is that the median activity of proteins in the control group is very low. As a result, signaling molecules other than ligands and receptors have a very low probability of being active at initialization, in order to simulate a background level of hepatocyte proliferation. This very low mean value has the effect of inflating the CV. In the case of IGF-1 or TNFα, the low reproducibility could be due to the inclusion of fewer reactions than the two other growth factors (EGF and INS). Hence, stimulating with IGF-1 or TNFα may not sufficiently stimulate the individual BNs for the entire ensemble to synchronize in response to treatment. Similar logic governs the model reproducibility following combined stimulation with IGF-1 and TNFα. These results also help to illuminate the sensitivity of our simulation approach to the topology of the signaling network. Importantly, the key endpoint of the model, AP-1 formation, is very reproducible across all treatment conditions.

Comparison with Experimental Data

We used experimental data on primary hepatocytes in culture [41] for a preliminary evaluation of our simulation approach and putative crosstalk model. In this experiment, rat primary hepatocytes were treated with varying concentrations of EGF and/or TNFα, and then the proportion of cells entering S phase (DNA Synthesis) was measured using Bromodeoxyuridine (5'-bromo-2'-deoxyuridine, BrdU). Although we do not explicitly model S phase in our network, the formation of the AP-1 transcription factor complex is believed to precede S phase in cell cycle progression. Hence, we assumed the formation of AP-1 to be a potential surrogate of S phase and, therefore, correlated with BrdU incorporation. We adjusted the probabilities of activation for proteins in our network in order to closely approximate the levels of BrdU incorporation in the absence of any treatment (control group). Further details on the calibration of the model are described in the Methods section.

We simulated the effects of different treatments on AP-1 formation.

Figure 5 shows the results of simulating 10 replicates for each of the treatment conditions including: combined EGF and TNFα, EGF only and TNFα only. The graphs in Figure 5 show the predicted activity profile of the AP-1 transcription factor complex across simulation time steps with the probability of activation for the treatment molecules set to 100%. Next, we simulated the concentration dependent effects of EGF stimulation. The steady state levels of AP-1 activity are shown in Figure 6 with the experimental data on BrdU incorporation for different treatment conditions. For each treatment condition, we simulated 100 replicates with 1000 cells per replicate. The plot in Figure 6 has a solid black line representing the mean of the fold change of AP-1 activity relative to mean of the control activity (at steady state) across all replicates. Based on the replicate data, we also report the 95% confidence interval for each plot -- the shaded blue region. The experimental data from [41] on BrdU incorporation is shown as points with the standard deviation. The simulation is able to reproduce the experimentally observed trends in DNA synthesis. As EGF is known to activate the AP-1 transcription factor complex (as a c-Jun/c-Fos dimer) in hepatocytes and other cell types [42], this result is consistent with the literature. Finally, the model did not predict the synergistic effect of EGF and TNFα stimulation on S phase. We believe this suggests mechanistic limitations in our crosstalk model that could be improved by incorporating additional mechanistic information about the downstream interactions between TNFα and EGF pathways.

Figure 5
figure5

AP-1 Activity Response Profiles Across Simulated Replicates to Different Treatment Conditions. Each replicate line represents the activity of AP-1 across 1000 simulated cells. There are three different treatment conditions: (A) EGF only, (B) TNFα only or (C) EGF and TNFα.

Full size image
Figure 6
figure6

Quantitative Comparison Between Simulation and Experimental Results. The quantitative comparison of simulated concentration response data with experimental data on EGF. The solid line represents the mean steady state fold change of AP-1 activity (relative to control) across 100 replicates with 1000 cells per replicate, simulated over a range of values from 0 to 25 ng/ml EGF. The 95% confidence interval across replicates is in blue. The experimental results of BrdU incorporation are plotted with standard deviation as points.

Full size image

Discussion

BNs have been used extensively to model the dynamics of molecular signaling and genetic regulatory networks. Because they require the discretization of molecular activity levels, however, a BN cannot be readily used to model the continuous concentration- and time-dependent effects of treatments. To address this issue we extended traditional BNs in three ways. First, we model the molecular response of an individual cell using a BN. Hence, we assumed that the average activity of signaling molecules in individual cells exhibits a switch-like (ON or OFF) response. Although this may not always be the case, we believe it is more biologically plausible than a binary representation of tissue level molecular activities. Second, we assumed that cellular response to stimulation is nondeterministic. It has been suggested synchronous updating schemes for simulating BNs can produce spurious attractors that disappear in the presence of noise [43, 44]. An asynchronous updating scheme, on the other hand, allows for variability (or noise) in cellular responses: two identical cells with the same initial configurations respond differentially to the same stimulus over time. Although asynchronous updating is not new in BNs [44], we define a probability of activation for each molecular species that is amenable to calibration and evaluation using experimental data. Third, we simulate an ensemble of asynchronous BNs to estimate the aggregate activity of each molecular species across cell populations. Our results show that this approach produces continuous responses similar to experimental observations from tissues. We believe this opens up new possibilities for estimating quantitative dose- and time-dependent responses in toxicology and in disease progression using knowledge of molecular mechanisms.

For this work, we used this simulation methodology to analyze the dynamics of a specific biochemical interaction network, which was constructed to investigate early molecular events surrounding hepatocellular proliferation. This is important because sustained cell proliferation is one of the key events in the progression of liver cancer. We find that (a) our extension of BNs yields highly reproducible results that have variability consonant with biological data and (b) our pathway-driven preliminary cytokine and growth factor protein signaling network is concordant with experimental observations on DNA synthesis in hepatocytes.

We investigated the effects of protein deletion from the network. In Additional file 1, we show a heatmap of these results. From this analysis, we believe our network is robust to the deletion of single proteins. Additionally, this information helps illustrate important signaling nodes in the network. For example, ERK1/2 and JNK appear to play important roles, since the removal of these changes the steady state values of other signaling molecules. To our surprise, the removal of IRS activity did not have a large effect on the signaling processes relative to other experimental conditions, even though the signaling molecule has a high connectivity in the graph.

A number of formalisms have been used to model the dynamics of eukaryotic cell cycle initiation/progression. Tyson and coauthors used ordinary differential equations (ODEs) to describe key cell division in frog oocytes [45]. Zielinski and colleagues used fuzzy Boolean logic to simulate receptor mediated crosstalk preceding cell proliferation in SKOV3 human epithelial cell line [46]. Similar fuzzy models have been proposed by others [47]. One advantage common to both techniques - i.e., ODEs and fuzzy logic - is that they can represent continuous or multivalued treatment concentrations. While these methods are powerful, they predict the population behavior of molecular species without emphasizing individual cellular protein activity. Our objective is to model the heterogeneous response of cell populations in order to estimate histological effects due to different treatments, which necessitates a modeling paradigm with individual cell clarity.

In order to help unravel the mechanisms of toxicity, we are compelled to investigate a simulation framework with a strong emphasis on network topology and reduced parameter space. ODE-based methods involve rate constants which are often difficult to quantify. Testing perturbations to systems of ODEs is not always straightforward with limited data. With BN dynamics, we lose some resolution of time but significantly reduce the number of parameters. We believe our technique is amenable to high throughput modeling of perturbations over a diverse chemical landscape, where the calibration of parameters can be limited due to scarce data on large numbers of environmental chemicals. We hope the modeling framework proposed will be useful in testing chemical perturbations from HTS data to generalized networks on signaling events assessed at the tissue level.

A comparable modeling framework was proposed [18] which also considered populations of networks. There are several differences in the modeling approaches. First, we use a binary representation of protein activity in a single cell while [18] used a ternary description to capture the level of protein expression relative to a population average under normal conditions. An ON/OFF representation of protein activity may not always be sufficient (e.g., caspase 8 activity) but it is generally consistent with single cell level observations (e.g. flow cytometry or high-content screening). The ternary representation used by [18] is based on western blot data on individual proteins with a comparison between treatment and control groups. While a mathematical transformation could relate one approach to another, the two techniques use a different abstraction for describing single cell biology.

Second, there is a subtle distinction between the abstraction of population level behavior between the two approaches. The authors of [18] calculated the average behavior of proteins across a set of BN whereas we use summation. This allows us to compare dose-dependent differences in potency and efficacy between treatments. Furthermore, it also enables the quantitative evaluation of population level "up-regulation" or "down-regulation" between treatment and control groups without using a ternary representation.

Third, the authors of [18] evaluate the effects of knockouts by maintaining some of the proteins in at 'control' or 'below control'. On the other hand, we consider the effects of dose-dependent perturbations in protein activities (e.g. extracellular ligands or intracellular signaling molecules) by using Hill functions to define probability of activation for certain protein(s) across the cell population. To our knowledge, no one has used this approach to incorporate concentration-response data in a BN modeling framework. The ability to reproduce and predict concentration-response data is essential for toxicological applications, bridging data from toxicity studies with systems biology to anticipate adverse outcomes.

We did not consider edge weights as an adjustable parameter for the model. All edges are weighted equally (set to 1.0). Modifying the edge weights would change the dynamics of the simulation. For example, in [48], the authors use a sigmoid function of the weighted sum to determine the probability per node in the propagation of the signal. Their technique offers a unique method for additional stochasticity to a threshold modeling framework. Modifying edge weights would require careful consideration of the updating scheme (describe in Methods), which would be affected if the edge weights were allowed to vary throughout simulation.

The thresholds described in the methods section provide a tunable parameter for investigating the signaling interactions. Each protein has a threshold value which defines the biochemical interaction surrounding its activity - that is, the logic underlying the interaction of the activating and inhibiting molecules. To illustrate the effects of modifying the threshold of a molecule in a network, we provide Figure 7: a truth table for variations on threshold values. The truth table displays the differences in the activity of a molecule, P, after one update (time step), as a function of the input value - the sum of the states of the inhibiting proteins subtracted from the sum of the states of the activating proteins - and the threshold value. Setting the threshold to an integer value allows for the molecule to maintain its current state whereas, following the discussion of edge weights, setting the threshold to a noninteger value will force a decisions for a (new) value of 0 or 1.

Figure 7
figure7

Truth Table for Update Scheme. The truth table illustrates how we calculate the activity of a protein (P) after one time step based on the aggregate input (column 1) and differential threshold values (columns 2 - 7). For example, row 2 shows that the activity of P remains the same as the previous time step () when the aggregate input is 0 and θ P = 0, and for all other threshold values. Similarly, the other row show of the consequence of varying the threshold for other values for the aggregate input. Note the distinction between integer and non-integer choices for the threshold.

Full size image

For the crosstalk network model investigated in this paper, randomly increasing the threshold of a single molecule from an integer to a noninteger value has little effect on the steady state levels of AP-1 (results not shown). These changes can effect the activity profile of the proteins over time - allowing sustained (threshold = 0.0) or transient activity (0.0 < threshold < 1.0) after activation, or recapitulating protein lability and potential loss of signal. For example,in the case of EGF stimulation, the receptor has negative feedback encoded in the model. Therefore, increasing the thresholds of proteins downstream of EGFR to noninteger values affects the overall activity profiles of the molecules across time (results not shown). Increasing thresholds by integer values can have the same effect as node deletion when the number of activating (input) molecules is equal to the threshold. In general, increasing the threshold increases the required input value for signal propagation.

There are some important limitations in our approach. First, the model presented in this work did not consider any communication between the cells (autocrine or paracrine signals), or between the cells and the extracellular matrix. Contact inhibition and matrix attachment, and cell communication are important factors in cell cycle progression and we are considering their role in a cell-based model of liver tissues [7]. Second, we also recognize that the biochemical interaction network we used in this work is incomplete. Our primary goal was to develop a flexible modeling approach that can incrementally accommodate additional mechanistic information as it becomes available.

The BN used in this work was constructed semi-automatically from a molecular network topology, which was defined with little manual intervention using curated information on pathways. In future work we will evaluate the effects of additional pathways in the crosstalk network, for example, the interleukin protein family as well as the behavior of mito-inhibitors like Transforming Growth Factor Beta (TGF-β). For this work, we did not investigate alternative hypotheses on the signaling mechanisms of individual hepatocytes. The model development in [16] provides a method for investigating signaling differences between cell types. Understanding the signaling differences between cell types, such as, kinetic differences, is undoubtedly important in extrapolating cell line in vitro data to acute/chronic in vivo responses. Leveraging information from alternate pathway databases [49–51] should increase the descriptive power of our model, and is critical to decipher the role of molecular crosstalk in cellular responses to endogenous ligands and xenobiotics.

AP-1 is among several transcription factors important for cell cycle initiation/progression. In the current model we did not include gene expression regulated by the activation of transcription factors. Therefore, we chose AP-1 formation as the surrogate for downstream events, since it promotes cell cycle progression through increased expression of proteins, such as cyclin D1 [52], and antagonizes the function other molecules, such as p53 and p21 [53]. As we continue to develop this model, we can incorporate the activity of additional important transcription factors, such as Nuclear Factor kappa B (NF-κB) and Forkhead Box (FOX) proteins. Finally, most BN based approaches cannot directly relate simulation time steps to physiologic time intervals. This is an open problem in BN simulation and our approach is not immune to this issue, but we hope to address this in future work.

An advantage of asynchronous BN ensembles is that they can be simulated very rapidly. A single simulated treatment group (1000 cells/replicate) can be executed on a computer in minutes, while some of the more complex simulations, such as the reproducibility investigation (Figure 4) can be simulated on the scale of hours. However, since each cell is initialized/simulated to steady state individually, the approach is amenable to parallelization. We believe this efficiency will allow us to simultaneously investigate the role of molecular network topology using background knowledge on quantitative tissue level responses from experimental data.

Conclusions

Systems biology approaches are vital for efficiently analyzing the effects of environmental contaminants on living tissues in order to evaluate the potential risk to public health. We developed asynchronous threshold Boolean network ensembles for translating mechanisms to in vitro and in vivo observations on molecular activity. The reproducibility of our modeling framework was evaluated by systematically analyzing the variability of our predictions across multiple simulations. We also showed that populations of hepatocytes can be simulated in this manner to predict experimentally observed quantitative responses. We believe that ensembles of Boolean networks can allow us to probe deeper mechanistic questions about the mode of action for chronic liver injury. We are testing this modeling approach as part of a broader computational and experimental effort aimed at estimating the putative effects of xenobiotics on the human liver by integrating chemical concentration, molecular pathways, cellular responses, and the role of cell-cell communication.

Methods

Threshold Boolean Network

We begin the definition of the cellular model with a biochemical interaction network as a signed, directed graph, G(V,E), where V is the set of all vertices (or proteins/molecules) and E is the set of all edges (or reactions). Let v ∈ V. Then, we define the set of all predecessors of v:

(1)

For each edge e uv ∈ E we have Sign(e uv ) ∈ {+,-} where '+' indicates u is involved in the activation of v and '-' indicates u is involved in the inhibition of v. Now we define A v ⊆ P v as the set of all activators of v. More formally,

(2)

Likewise, we define I v ⊆ P v as the set of all inhibitors of v,

(3)

Furthermore, we let and , and let n be the number of proteins in the graph. Additionally, we store the binary vector, , of the state (active or inactive) of every vertex at time τ. The state of a protein is dependent on the states of its predecessors. Therefore, we define a vector, , representing the threshold of activation for each vertex, a biologically inspired variable guiding the interplay predecessor vertices and protein activation. For the model, all thresholds were set to 0.0 with the exception of AP-1 formation, which is set to θAP-1∈ (1,2). This modification to the AP-1 threshold reflects the underlying biochemistry in that both c-Jun and c-Fos must be active for the activation of the AP-1.

Finally, we define the vector , which represents the probability of activation for each protein. For most models/proteins, there is a basal level of activity. We assume that individual BNs can have different protein activity profiles upon initialization (τ = 0), which allows for biological variability across the cell population.

Now, for any direct, signed graph, we formally define our model as follows:

(4)

Temporal Evolution of the Boolean Network

All other components of the model, Cτ , are fixed during the simulation except for one, Δτ, which is the state vector for v∈ V at a given time step. This value is determined by using the following threshold based logic:

(5)

The steady state protein activity in a BN is expressed as the following state vector:

(6)

Model Calibration

In this model, the probabilities of activation for proteins, Φ, are considered a tunable parameter. The probabilities determine the state of each BN at τ = 0, which influences the dynamics of protein activation across the ensemble. We adjusted the values of Φ in three steps using experimental data where available. First, we assumed that in the absence of any treatment (i.e., experimental controls) the ligands, receptors, the adaptor proteins (GRB2 and SHC), the insulin receptor substrates (IRS and IRS-2), and non-ligand receptors (ErbB2 and ErbB3) have φ = 0.0. Second, the values of Φfor the remaining proteins in the network were increased until the predicted activity of AP-1 was close to the experimental level of BrdU incorporation in the control group (~1.5% DNA Synthesis [41]). For these proteins, the probability of activation φ = 0.0025 gave 1.49% ± 0.05 of AP-1 formation. Third, we assumed that the probability of activation of ligands in the model was related to the experimental concentration of the ligand. For EGF, we used the Hill function (Equation 8) to describe the relationship between probability of activation, φ EGF , and treatment concentration (in ng/ml).

(8)

We used the BrdU concentration response data [41] to estimate the parameters for Equation 8. In order to use the Nelder-Mead algorithm to make a maximum likelihood estimation of Vmax, K A and n, we assumed that EGF activity corresponds to AP-1 activation. Although the network was encoded with negative feedback for the EGF receptor, representing the internalization and ubiquitination of this receptor, we make this assumption based on the simulation results with φ EGF = 1.0 (Figure 5a). The maximum likelihood estimates we found are n = 0.7, K A = 5.9ng/ml and Vmax = .28 (probability of activation).

Simulating Populations of Cells

The ensemble of asynchronous threshold Boolean networks at a time step is represented as:

The aggregate activity of each protein across the ensemble at one time step, denoted as vτ , is calculated across Cτ , as follows:

(9)

Similarly the steady state activity of a protein across the ensemble is denoted as vT .

Hence, the coefficient of variation of the steady state protein activity is calculated as follows:

(10)

Implementation and Input Formats

The simulator is implemented in Python and the network model was produced interactively using Cytoscape. The molecular interaction network topology is defined in the Cytoscape exported format and protein information as well as quantitative parameters can be defined in the node attributes file (e.g., protein name, probabilities of activation, and threshold). The Python code as well as the network model are included (see Additional file 2 and Additional file 3).

References

  1. 1.

    Andersen ME, Krewski D: Toxicity Testing in the 21st Century: Bringing the Vision to Life. Toxicological Sciences. 2009, 107: 324-330.

    Article  CAS  PubMed  Google Scholar 

  2. 2.

    Holsapple MP, Afshari CA, Lehman-McKeeman LD: Forum Series: The "Vision" for Toxicity Testing in the 21st Century: Promises and Conundrums. Toxicological Sciences. 2009, 107: 307-308.

    Article  CAS  PubMed  Google Scholar 

  3. 3.

    White RH, Cote I, Zeise L, Fox M, Dominici F, Burke TA, White PD, Hattis DB, Samet JM: State-of-the-Science Workshop Report: Issues and Approaches in Low-Dose-Response Extrapolation for Environmental Health Risk Assessment. Environ Health Perspect. 2008, 117:

    Google Scholar 

  4. 4.

    De Rosa C, El-Masri H, Pohl H, Cibulas W, Mumtaz M: Implications of Chemical Mixtures in Public Health Practice. Journal of Toxicology and Environmental Health Part B: Critical Reviews. 2004, 7: 339-350. 10.1080/10937400490498075.

    Article  CAS  PubMed  Google Scholar 

  5. 5.

    Lawrence JW, Li Y, Chen S, DeLuca JG, Berger JP, Umbenhauer DR, Moller DE, Zhou G: Differential Gene Regulation in Human Versus Rodent Hepatocytes by Peroxisome Proliferator-activated Receptor (PPAR) α. Journal of Biological Chemistry. 2001, 276: 31521-31527. 10.1074/jbc.M103306200

    Article  CAS  PubMed  Google Scholar 

  6. 6.

    Cohen SM, Lawson TA: Rodent bladder tumors do not always predict for humans. Cancer Letters. 1995, 93: 9-16. 10.1016/0304-3835(95)03785-U

    Article  CAS  PubMed  Google Scholar 

  7. 7.

    Wambaugh J, Shah I: Simulating Microdosimetry in a Virtual Hepatic Lobule. PLoS Comput Biol. 2010, 6: e1000756- 10.1371/journal.pcbi.1000756

    PubMed Central  Article  PubMed  Google Scholar 

  8. 8.

    Shah I, Wambaugh J: Virtual Tissues in Toxicology. Journal of Toxicology and Environmental Health Part B: Critical Reviews. 2010, 13: 314-328. 10.1080/10937404.2010.483948.

    Article  CAS  Google Scholar 

  9. 9.

    Kauffman S: Metabolic stability and epigenesis in randomly constructed genetic nets. Journal of Theoretical Biology. 1969, 22: 437-467. 10.1016/0022-5193(69)90015-0

    Article  CAS  PubMed  Google Scholar 

  10. 10.

    Zoltan S, Shoudan L: Modeling the normal and neoplastic cell cycle with 'realistic Boolean genetic networks': Their application for understanding carcinogenesis and assessing therapeutic strategies. Book Modeling the normal and neoplastic cell cycle with 'realistic Boolean genetic networks': Their application for understanding carcinogenesis and assessing therapeutic strategies. 1998, City: World Scientific Publishing,

    Google Scholar 

  11. 11.

    Shmulevich I, Dougherty ER, Kim S, Zhang W: Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks. Bioinformatics. 2002, 18: 261-274. 10.1093/bioinformatics/18.2.261

    Article  CAS  PubMed  Google Scholar 

  12. 12.

    Sánchez L, Thieffry D: A Logical Analysis of the Drosophila Gap-gene System. Journal of Theoretical Biology. 2001, 211: 115-141. 10.1006/jtbi.2001.2335

    Article  PubMed  Google Scholar 

  13. 13.

    Davidich MI, Bornholdt S: Boolean Network Model Predicts Cell Cycle Sequence of Fission Yeast. PLoS ONE. 2008, 3: e1672- 10.1371/journal.pone.0001672

    PubMed Central  Article  PubMed  Google Scholar 

  14. 14.

    Li F, Long T, Lu Y, Ouyang Q, Tang C: The yeast cell-cycle network is robustly designed. Proceedings of the National Academy of Sciences of the United States of America. 2004, 101: 4781-4786. 10.1073/pnas.0305937101

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  15. 15.

    Faure A, Naldi A, Chaouiya C, Thieffry D: Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle. Bioinformatics. 2006, 22: e124-131. 10.1093/bioinformatics/btl210

    Article  CAS  PubMed  Google Scholar 

  16. 16.

    Saez-Rodriguez J, Alexopoulos LG, Epperlein J, Samaga R, Lauffenburger DA, Klamt S, Sorger PK: Discrete logic modelling as a means to link protein signalling networks with functional analysis of mammalian signal transduction. Mol Syst Biol. 2009, 5:

    Google Scholar 

  17. 17.

    Samaga R, Saez-Rodriguez J, Alexopoulos LG, Sorger PK, Klamt S: The Logic of EGFR/ErbB Signaling: Theoretical Properties and Analysis of High-Throughput Data. PLoS Comput Biol. 2009, 5: e1000438- 10.1371/journal.pcbi.1000438

    PubMed Central  Article  PubMed  Google Scholar 

  18. 18.

    Wu M, Yang X, Chan C: A Dynamic Analysis of IRS-PKR Signaling in Liver Cells: A Discrete Modeling Approach. PLoS ONE. 2009, 4: e8040- 10.1371/journal.pone.0008040

    PubMed Central  Article  PubMed  Google Scholar 

  19. 19.

    Philippi N, Walter D, Schlatter R, Ferreira K, Ederer M, Sawodny O, Timmer J, Borner C, Dandekar T: Modeling system states in liver cells: Survival, apoptosis and their modifications in response to viral infection. BMC Systems Biology. 2009, 3: 97- 10.1186/1752-0509-3-97

    PubMed Central  Article  PubMed  Google Scholar 

  20. 20.

    Zhang R, Shah MV, Yang J, Nyland SB, Liu X, Yun JK, Albert R, Loughran TP: Network model of survival signaling in large granular lymphocyte leukemia. Proceedings of the National Academy of Sciences. 2008, 105: 16308-16313. 10.1073/pnas.0806447105.

    Article  CAS  Google Scholar 

  21. 21.

    Martin MT, Judson RS, Reif DM, Kavlock RJ, Dix DJ: Profiling Chemicals Based on Chronic Toxicity Results from the U.S. EPA ToxRef Database. Book Profiling Chemicals Based on Chronic Toxicity Results from the US EPA ToxRef Database. 2009, 117: 392-399. City: National Institute of Environmental Health Sciences,

    Google Scholar 

  22. 22.

    Chhabra RS, Bucher JR, Wolfe M, Portier C: Toxicity characterization of environmental chemicals by the US National Toxicology Program: an overview. International Journal of Hygiene and Environmental Health. 2003, 206: 437-445. 10.1078/1438-4639-00240

    Article  CAS  PubMed  Google Scholar 

  23. 23.

    Hanahan D, Weinberg RA: The Hallmarks of Cancer. Cell. 2000, 100: 57-70. 10.1016/S0092-8674(00)81683-9

    Article  CAS  PubMed  Google Scholar 

  24. 24.

    Butterworth BE: Consideration of both genotoxic and nongenotoxic mechanisms in predicting carcinogenic potential. Mutation Research/Reviews in Genetic Toxicology. 1990, 239: 117-132. 10.1016/0165-1110(90)90033-8.

    Article  CAS  Google Scholar 

  25. 25.

    Csikász-Nagy A, Battogtokh D, Chen KC, Novák B, Tyson JJ: Analysis of a Generic Model of Eukaryotic Cell-Cycle Regulation. 2006, 90: 4361-4379.

    Google Scholar 

  26. 26.

    Lahav G, Rosenfeld N, Sigal A, Geva-Zatorsky N, Levine AJ, Elowitz MB, Alon U: Dynamics of the p53-Mdm2 feedback loop in individual cells. Nat Genet. 2004, 36: 147-150. 10.1038/ng1293

    Article  CAS  PubMed  Google Scholar 

  27. 27.

    Batchelor E, Loewer A, Lahav G: The ups and downs of p53: understanding protein dynamics in single cells. Nat Rev Cancer. 2009, 9: 371-377. 10.1038/nrc2604

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  28. 28.

    Tay S, Hughey JJ, Lee TK, Lipniacki T, Quake SR, Covert MW: Single-cell NF-[kgr]B dynamics reveal digital activation and analogue information processing. Nature. 2010, 466: 267-271. 10.1038/nature09145

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  29. 29.

    Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK: Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature. 2009, 459: 428-432. 10.1038/nature08012

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  30. 30.

    Luke N, Devito M, Portier C, El-Masri H: Employing a Mechanistic Model for the MAPK Pathway to Examine the Impact of Cellular All or None Behavior on Overall Tissue Response. Dose Response. 2010, 8: 347-367. 10.2203/dose-response.09-017.Luke

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  31. 31.

    Cohen AA, Geva-Zatorsky N, Eden E, Frenkel-Morgenstern M, Issaeva I, Sigal A, Milo R, Cohen-Saidon C, Liron Y, Kam Z, et al.: Dynamic Proteomics of Individual Cancer Cells in Response to a Drug. Science. 2008, 322: 1511-1516. 10.1126/science.1160165

    Article  CAS  PubMed  Google Scholar 

  32. 32.

    Derrida B: Dynamical phase transition in nonsymmetric spin glasses. Journal of Physics A: Mathematical and General. 1987, 20: L721-10.1088/0305-4470/20/11/009.

    Article  Google Scholar 

  33. 33.

    Bornholdt S: Boolean network models of cellular regulation: prospects and limitations. Journal of The Royal Society Interface. 2008, 5: S85-S94. 10.1098/rsif.2008.0132.focus.

    PubMed Central  Article  CAS  Google Scholar 

  34. 34.

    Robert J, Lefkowitz SKS: IGF-1 Receptor Signaling through betat-Arrestin. Sci Signal. 2005,

    Google Scholar 

  35. 35.

    Schlessinger J: Epidermal Growth Factor Receptor Pathway. Book Epidermal Growth Factor Receptor Pathway. 2005, City,

    Google Scholar 

  36. 36.

    White MF: Insulin Signaling Pathway. Sci Signal. 2005,

    Google Scholar 

  37. 37.

    David V, Goeddel GC: Tumor Necrosis Factor Pathway. Sci Signal. 2005,

    Google Scholar 

  38. 38.

    Shannon P, Markiel A, Ozier O, Baliga N, Wang J, Ramage D, Amin N, Schwikowski B, Ideker T: Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Research. 2003, 13: 2498-2504. 10.1101/gr.1239303

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  39. 39.

    Ayala J, Boustead J, Chapman S, Svitek C, Oeser J, Robey R, O'Brien R: Insulin-mediated activation of activator protein-1 through the mitogen-activated protein kinase pathway stimulates collagenase-1 gene transcription in the MES 13 mesangial cell line. J Mol Endocrinol. 2004, 33: 263-280. 10.1677/jme.0.0330263

    Article  CAS  PubMed  Google Scholar 

  40. 40.

    Newman JRS, Ghaemmaghami S, Ihmels J, Breslow DK, Noble M, DeRisi JL, Weissman JS: Single-cell proteomic analysis of S. cerevisiae reveals the architecture of biological noise. Nature. 2006, 441: 840-846. 10.1038/nature04785

    Article  CAS  PubMed  Google Scholar 

  41. 41.

    Chevalier S, Macdonald N, Roberts R: Induction of DNA replication by peroxisome proliferators is independent of both tumour necrosis factor (alpha) priming and EGF-receptor tyrosine kinase activity. J Cell Sci. 1999, 112: 4785-4791.

    CAS  PubMed  Google Scholar 

  42. 42.

    Nadori F, Lardeux B, Rahmani M, Bringuier A, Durand-Schneider A, Bernuau D: Presence of distinct AP-1 dimers in normal and transformed rat hepatocytes under basal conditions and after epidermal growth factor stimulation. Hepatology. 1997, 26: 1477-1483. 10.1002/hep.510260614

    Article  CAS  PubMed  Google Scholar 

  43. 43.

    Greil F, Drossel B: Dynamics of Critical Kauffman Networks under Asynchronous Stochastic Update. Physical Review Letters. 2005, 95: 048701-

    Article  PubMed  Google Scholar 

  44. 44.

    Klemm K, Bornholdt S: Stable and unstable attractors in Boolean networks. Physical Review E. 2005, 72: 055101-

    Article  Google Scholar 

  45. 45.

    Zwolak JW, Tyson JJ, Watson LT: Parameter Estimation for a Mathematical Model of the Cell Cycle in Frog Eggs. Journal of Computational Biology. 2005, 12: 48-63. 10.1089/cmb.2005.12.48

    Article  CAS  PubMed  Google Scholar 

  46. 46.

    Zielinski R, Przytycki P, Zheng J, Zhang D, Przytycka T, Capala J: The crosstalk between EGF, IGF, and Insulin cell signaling pathways - computational and experimental analysis. BMC Systems Biology. 2009, 3: 88- 10.1186/1752-0509-3-88

    PubMed Central  Article  PubMed  Google Scholar 

  47. 47.

    Huang Z, Hahn J: Fuzzy modeling of signal transduction networks. Chemical Engineering Science. 2009, 64: 2044-2056. 10.1016/j.ces.2009.01.041.

    Article  CAS  Google Scholar 

  48. 48.

    Rotroff D, Beam A, Dix D, Farmer A, Freeman K, Houck K, Judson R, LeCluyse E, Martin M, Reif D, Ferguson S: Xenobiotic-Metabolizing Enzyme and Transporter Gene Expression in Primary Cultures of Human Hepatocytes Modulated by Toxcast Chemicals. Journal of Toxicology and Environmental Health, Part B: Critical Reviews. 2010, 13: 329-346. 10.1080/10937404.2010.483949.

    Article  CAS  Google Scholar 

  49. 49.

    Kanehisa M, Goto S, Furumichi M, Tanabe M, Hirakawa M: KEGG for representation and analysis of molecular networks involving diseases and drugs. Nucleic Acids Research. 2010, 38: D355-D360. 10.1093/nar/gkp896

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  50. 50.

    Joshi-Tope G, Gillespie M, Vastrik I, D'Eustachio P, Schmidt E, de Bono B, Jassal B, Gopinath GR, Wu GR, Matthews L, et al.: Reactome: a knowledgebase of biological pathways. Nucleic Acids Research. 2005, 33: D428-D432.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  51. 51.

    Schaefer CF, Anthony K, Krupa S, Buchoff J, Day M, Hannay T, Buetow KH: PID: the Pathway Interaction Database. Nucleic Acids Research. 2009, 37: D674-D679. 10.1093/nar/gkn653

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  52. 52.

    Bassingthwaighte J, Liebovitch LS, West BJ: Fractal Physiology. 1994, Oxford University Press, New York,

    Google Scholar 

  53. 53.

    Hunt CA, Yan L, Ropella G, Park S, Sheikh-Bahael S: The multiscale in silico liver. Journal of Critical Care. 2007, 22: 348-349. 10.1016/j.jcrc.2007.10.031.

    Article  Google Scholar 

Download references

Acknowledgements

This work was completed at the National Center for Computational Toxicology at the US Environmental Protection Agency. We would like to thank R. Woodrow Setzer and David Reif of the EPA and Julio Saez-Rodriguez of European Bioinformatics Insititute for useful conversations and critical review.

Disclaimer

The views expressed in this article are those of the authors and do not necessarily reflect the views or policies of the U.S. Environmental Protection Agency. Mention of trade names or commercial products does not constitute endorsement or recommendation for use.

Author information

Affiliations

  1. National Center for Computational Toxicology, Office of Research and Development, U.S. Environmental Protection Agency, Research Triangle Park, North Carolina, USA

    John Jack, John F Wambaugh & Imran Shah

Authors
  1. John Jack
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. John F Wambaugh
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Imran Shah
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Imran Shah.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

JJ designed and implemented the methods, performed the analysis, and wrote the manuscript; JW participated in developing the methods and writing the manuscript; IS supervised the development of the methods, implementation, analysis, and writing the manuscript.

Electronic supplementary material

Additional file 1:Evaluation of Network Behavior for Protein Knockouts. The heatmap shows the simulation results for deleting individual proteins from the network. Each cell in the heatmap represents the mean protein activity at steady state relative to control across 20 replicates with 100 cells per replicate. The color intensity indicates the protein (x-axis) behavior at steady state relative to the baseline simulation (no protein knockout). The y-axis indicates the protein deletion. (DOC 124 KB)

Additional file 2:network_model.zip. The Cytoscape export files for the crosstalk network model on AP-1 formation. These files are to be used in conjunction with the python_simulation_code. (ZIP 2 KB)

Additional file 3:python_simulation_code.zip. The source code to load the network model and produce the simulation results reported in the paper. Refer to the README file for instruction. (ZIP 16 KB)

Authors’ original submitted files for images

Rights and permissions

Reprints and Permissions

About this article

Cite this article

Jack, J., Wambaugh, J.F. & Shah, I. Simulating Quantitative Cellular Responses Using Asynchronous Threshold Boolean Network Ensembles. BMC Syst Biol 5, 109 (2011). https://doi.org/10.1186/1752-0509-5-109

Download citation

Keywords

  • Epidermal Growth Factor
  • Boolean Network
  • Genetic Regulatory Network
  • Aggregate Activity
  • Epidermal Growth Factor Stimulation

BMC Systems Biology

ISSN: 1752-0509

Contact us