# Boolean network simulations for life scientists

- István Albert
^{1}Email author, - Juilee Thakar
^{2}, - Song Li
^{3}, - Ranran Zhang
^{4}and - Réka Albert
^{1, 2, 3}

**3**:16

**DOI: **10.1186/1751-0473-3-16

© Albert et al; licensee BioMed Central Ltd. 2008

**Received: **24 September 2008

**Accepted: **14 November 2008

**Published: **14 November 2008

## Abstract

Modern life sciences research increasingly relies on computational solutions, from large scale data analyses to theoretical modeling. Within the theoretical models Boolean networks occupy an increasing role as they are eminently suited at mapping biological observations and hypotheses into a mathematical formalism. The conceptual underpinnings of Boolean modeling are very accessible even without a background in quantitative sciences, yet it allows life scientists to describe and explore a wide range of surprisingly complex phenomena. In this paper we provide a clear overview of the concepts used in Boolean simulations, present a software library that can perform these simulations based on simple text inputs and give three case studies. The large scale simulations in these case studies demonstrate the Boolean paradigms and their applicability as well as the advanced features and complex use cases that our software package allows. Our software is distributed via a liberal Open Source license and is freely accessible from http://booleannet.googlecode.com

## Introduction

At the most general level systems biology approaches consist of two steps. The first is building a *model* of the biological system of interest, a representation that incorporates existing knowledge and experimental observations. This model then can be subjected to various conditions and may be allowed to evolve in time, a step typically referred to as *simulation*. These simulations then can be used to generate qualitative or quantitative predictions on the overall behavior of the system.

Mathematical models of biological systems range from continuous to discrete (based on the representation of the status of the system's components) and from deterministic to stochastic (based on their incorporation of randomness and noise) [1–4]. The simplest models assume that each element of the system has a binary (Boolean) state, and are therefore discrete, deterministic and parameter-free. However, several extensions of the Boolean modeling formalism also allow for iterative parameterization and the incorporation of continuous and stochastic elements [5–7]. The most important barrier precluding a more general use of Boolean models consists of the difficulties of the computational implementation of a given model. This implementation needs to perform in a consistent, correct and extensible manner to allow for the integration of all empirical knowledge as well as the unfettered exploration of the dynamical behaviors allowed by the model. There are very few tools that focus on qualitative modeling, for example the Genetic Network Analyzer [8] supports qualitative predictions via a piece-wise linear model and provides advanced visualization capabilities. There are several software packages that focus on quantitative modeling via differential equations [9–11], and discrete modeling [7, 12] but these tools are less suited for exploratory analyses of biological systems in which the majority of kinetic parameters are unknown.

In this paper we present a software toolbox (from now on referred to as BooleanNet) that greatly facilitates the implementation and study of Boolean dynamic models of biological systems. It is a tool that can simulate a Boolean model based on a very simple text based input describing the interactions and regulatory relationships in the system. The main distinguishing feature of our software compared to previous efforts is that we aim to provide support for modeling the dynamic behavior of well defined biological sub-systems, rather than focusing on a larger scale network inference, analysis or modeling based on high throughput data. Once the rules are expressed the software can employ several simulation strategies: synchronous iterations, stochastic updates or hybrid modeling via a system of piecewise linear differential equations. More importantly our system allows the integration of non-boolean mechanisms into the simulation thus expanding its applicability to a wider domain. Every aspect of the simulation process may be customized: node states may be overridden at different stages of the operation, updating rules may be altered, and differential equations may be augmented or replaced. We place particular focus on the documentation and demonstrate use cases from simple usage samples to advanced and realistic examples. A series of 6 tutorials guide users through various aspects of the simulations. Our code, written in the Python programming language, is available as a software library and is distributed via a liberal Open Source license, freely accessible through the Google project hosting service at http://booleannet.googlecode.com

### Boolean network models

A Boolean network model is a directed graph (network) whose nodes represent the elements of a system, edges represent regulatory relationships between elements, and every node is characterized by a True (1) or False (0) state [12–14]. Thus a network with N nodes will have 2^{N} possible states. As time passes the state of each node is determined by the states of its neighbors (nodes that point to it), through a rule called a **transfer function**. In the simplest case the transfer function can be represented as a statement acting on the inputs via a logical function using the logical operators **NOT**, **AND**, **OR**; this statement also returns a True/False state. Depending on the output of the transfer function, the state of the node either stays the same or changes. A change in the state of a given node generally triggers changes in the state of nodes regulated by it (the nodes it points to). This way the state of the network goes through a dynamic trajectory.

This type of Boolean representation is very common in biology, although hardly if ever is named as so. For example in the case of a gene regulatory network scientists routinely describe genes and pathways as being activated or repressed (On/Off), the analysis of gene expression levels leads to a classification of genes as being differentially expressed or not, during an infection a pathogen is either cleared or it persists. This classification is a very natural desire to order and categorize observed phenomena, and intuitively estimate the outcome of regulatory mechanisms. However, the complexity of biological networks makes even the enumeration of network components and interactions a daunting task, and human intuition must be complemented by formal modeling. Boolean models and their extensions (collectively called qualitative models) have the common feature of very closely reflecting the topology of the regulatory network (i.e. the upstream/downstream regulatory relationships invoked in verbal pathway descriptions). In addition they readily incorporate inhibitory relationships (through the NOT operator) and combinatorial regulation of a target node by several regulators, whether these regulators act independently (represented by the OR operator) or they are conditionally dependent (represented by the AND operator). Different types of qualitative models offer a range of choices for representing the passing of time, from regular time steps to sampling over relative timescales and to explicit incorporation of known decay times.

In its simplest formalism a Boolean model assumes that the processes represented as edges in the network have similar durations, and correspondingly node states are updated at multiples of a fixed time step [15]. In our software we call this mode of operation as **synchronous update**. During each iteration nodes are updated to their new values only after all rules have been applied, thus the order of updates does not affect the outcome. Synchronous Boolean models starting with a given initial condition will always reach the same state after the same number of steps, i.e. the system is entirely deterministic. Moreover, due to the finite number of states attainable by the system, the system's trajectory in state space will ultimately converge into an attractor, either a single state called a **steady state** or a repeating sequence of states called a **cycle** and characterized by a cycle length. Any given Boolean model has one or several attractors, and each attractor is associated with a set of states (called its basin of attraction) that if used as an initial condition, converge into that attractor. In synchronous Boolean models the basins of attraction of each attractor are non-overlapping. Even in this simple case much can be learned from varying the initial conditions and analyzing the model's steady states and cycles (or lack thereof). We may be able to indentify nodes (variables) that are more important and drive the early appearance of steady states, or we may find nodes whose effect matters very little, either because of the redundancy of the pathway topology or because their effect is being cancelled out. Our software implements a steady state and cycle detection routine that can be used to quickly determine their length and first occurrence.

The synchronous model cannot properly account for the different time scales over which various events take place in the biological system. Most often these time scales are not known, nonetheless imposing the equality of all time scales, as the synchronous model does, introduces an artificial constraint. We can extend the base model to account for timescales by allowing the node updates to affect the state instantaneously while performing the updates in a random order within each iteration [5, 2]. This extension introduces stochasticity into the evolution of the system and allows it to sample all timescales. In this mode of operation (named **asynchronous update** in BooleanNet) the study of the network most often consists of performing a large number of replicate simulations that start from the same initial condition, then tallying the attractors reachable from this initial condition (there can be several due to the stochasticity of the asynchronous update) and the states of the nodes at different time steps. This model can be further expanded by introducing rank ordered updates, essentially grouping nodes into separate ranks based on information regarding the relative duration of their synthesis or decay. The update orders in each rank group are randomized but are executed before the rules contained in a higher rank group. This mode is also supported by our software, where the ranks are numbers affixed to each rule, and rules sharing a certain number are considered a rank group.

*AND*C will be replaced by a differential equation of the form:

Above an inequality will numerically evaluate to 1 if the statement is True and 0 otherwise. Here the Boolean (first) term corresponds to the *regulated synthesis* of A while the decay (second) term corresponds to its *free (unregulated) dissociation*. When the Boolean term is *True* the equation is of the form $\frac{dcon{c}_{A}}{dt}=1-deca{y}_{A}\times con{c}_{A}$, leading to an asymptotic increase to *1*/*delay*_{
A
}, or maintenance of an initial concentration equal to this value. When the Boolean term is *False* the equation is of the form $\frac{dcon{c}_{A}}{dt}=-deca{y}_{A}\times con{c}_{A}$, leading to an exponential decrease to zero, or maintenance of an initial concentration equal to zero. The limiting values 0 and *1*/*delay*_{
A
}of the continuous variable represent, respectively, the absence of species *A* and maximal concentration of species *A*, and correspond to the discrete values False (0) and True (1). Note that the steady states of a Boolean model are the same, regardless of the mechanism of update, and coincide with the steady state of the Boolean variables of the corresponding piecewise linear model.

## Case studies

Below we present a number of realistic scenarios that apply the BooleanNet library to several different problem domains. The rule sets, data and the code that produce the plots are part of the software package and can be found in the 'examples' directory. The case studies are based on published models; our main goal for presenting these is to demonstrate the modeling aspect of the approaches, to illustrate and explain customizations and adaptations that are often necessary to capture the biological phenomena of interest. Biological modeling is rarely "pure" in the sense of being able to fully and rigorously adhere to the mathematical formalisms. There are always elements, actors and processes do not completely fit the model, and need to be altered to accommodate some observed behaviors or limits. We believe that our software allows for an unprecedented customization, where all aspects of the simulation process can be interacted with, nodes, rules and states can be observed and modified in every single timestep. The rules, simulation and visualization code for each of the examples below can be found in the main source code distribution. Each of the three projects below have been published as separate article [18–20] that should be consulted for a more in depth description on the details of the biological systems and of the dynamic models. On our website we also distribute a series of accessible tutorials that gradually introduce the concepts used in the simulations.

### Case study 1: Modeling abscisic acid (ABA) – induced stomatal closure in plants

Plants take up carbon dioxide for photosynthesis through microscopic pores called stomata. The size of the pores is determined by the two cells (called guard cells) that flank the stomata, and the shape of the guard cells is in turn determined by their water content and turgor pressure. Plants also lose water by transpiration through the pores, thus they need to regulate the size of the stomata to balance carbon gain with water loss. During drought conditions plants synthesize a hormone called abscisic acid (ABA). ABA acts as a signal to a complex signal transduction pathway leading to the closure of the stomata. Following an extensive curation of the experimental literature, in [17] we have synthesized a network for ABA-induced closure that includes 43 nodes and 69 regulatory edges (see website or source distribution for the rules). The nodes of the network include proteins, ion channels and secondary messengers, as well as more abstract concepts such as membrane depolarization and stomatal closure. The state of 38 of the 43 nodes is regulated by other nodes in the network and in [17] we expressed this regulation as Boolean rules. Here we use these rules to compare the asynchronous and piece-wise linear frameworks of BooleanNet in exploring the dynamics of the system.

*θ*equal to 0.5, and we performed 300 simulations with random initial states. For both modes, we tested the wild type as well as a number of important knockouts and plotted the average state of the node "Closure" as a function of time steps. Specifically, in the asynchronous Boolean simulation we plot the percentage of simulations that have Closure = 1, while in the piecewise linear simulation we plot the mean and mean ± standard deviation of the continuous variable corresponding to closure. In the wild type (WT) simulation, all nodes are updated as specified by the Boolean rules. We also performed knockout simulations of sphingosine one-phosphate (S1P), phosphatidic acid (PA), cytosolic pH (pH

_{c}) and abscisic acid insensitive (ABI1), setting the corresponding node's state to False (0). We found that for WT simulations, the two model variants give the same closure response (ABA Figure 1a/b, blue curves). In addition, the two simulations produced qualitatively similar results for knockouts. We also observed that the variation of the closure response was smaller in WT than in the mutants, whereas the mutants showed similar, higher than WT variations. These results suggest that knockouts not only affect the mean closure responses but also their variation [17].

### Case study 2: T-cell large granular lymphocyte leukemia simulation

### Case study 3: Modeling the mammalian immune response to B. bronchiseptica infection

The dynamic interplay between a pathogen (e.g. virus or bacterium) and its host's defenses decides whether the pathogen will be cleared or will establish a niche in the host. The gram-negative bacterium *Bordetella bronchiseptica* persists within its mammalian hosts by interfering with the hosts' immune responses [25]. In [26] we synthesized a network of interactions among the immune effectors of a mammalian host and the virulence factors of *B. bronchiseptica* from the experimental literature. The network contains 34 nodes (see website for rules) that represent immune cell types, signaling molecules (cytokines), antibodies and bacterial factors such as antigens. The regulatory edges of the network represent immune processes such as antigen presentation, activation of immune cells as well as the modulation of immune components by bacteria. Based on our best interpretation of the experimental literature we formulated a piece-wise linear model and we used BooleanNet to simulate the interplay between host immune components and bacterial factors.

where *DC* represents the concentration of dendritic cells and subscripts 1 and 2 indicate the concentrations in different compartments. Since cytokines flow to various places through lymph and blood their inter-compartmental dynamics were given by the random variables *rc* and *r* so that a fraction *rc* - *r* <*f* <*rc* + *r* of the total concentration of the cytokine is present in the compartment in which the cytokine is produced and a fraction *1*- *f* is transported to the other compartment. We changed the fraction *f* after a variable expiration time several times during the simulation.

We also simulated previously known network perturbations e.g. B cell, T cell deletion and the infections by Type III secretion system (TTSS) defective Bordetella strain. The model could reproduce the persistence of bacteria in the absence of B cells (see Figure 3, DEL) and T cells [27, 28]. It can also simulate the earlier clearance of bacteria in Δ*bscN* mutant (TTSS defective strain) infection [29].

## Conclusion

The success of qualitative models in describing specific cellular systems and processes such as flower development [30, 31], the yeast cell cycle [32], and Drosophila embryonic development [33–35] indicates that the topology of regulatory networks has a significant role in restricting their dynamical behavior. Our software allows straightforward implementation of qualitative models for systems where the network topology or pathway is at least partially known.

## Declarations

### Acknowledgements

This work was partially supported by NSF grants MCB-0618402 and CCF-0643529 (CAREER), NIH grants 1R55AI065507 – 01A2 and 1 R01 GM083113-01 and USDA grant 2006-35100-17254.

## Authors’ Affiliations

## References

- Fall CP, Marland ES, Wagner JM, Tyson JJ: Computational Cell Biology. 2002, New York: SpringerGoogle Scholar
- Voit EO: Computational Analysis of Biochemical Systems. 2000, Cambridge: Cambridge University PressGoogle Scholar
- de Jong H: Modeling and simulation of genetic regulatory systems: a literature review. J Comput Biol. 2002, 9: 67-103. 10.1089/10665270252833208.View ArticlePubMedGoogle Scholar
- Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits. 2006Google Scholar
- Chaves M, Albert R, Sontag ED: Robustness and fragility of Boolean models for genetic regulatory networks. J Theor Biol. 2005, 235: 431-449. 10.1016/j.jtbi.2005.01.023.View ArticlePubMedGoogle Scholar
- Chaves M, Sontag ED, Albert R: Methods of robustness analysis for Boolean models of gene control networks. Syst Biol (Stevenage). 2006, 153: 154-167.View ArticleGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- de Jong H, Geiselmann J, Hernandez C, Page M: Genetic Network Analyzer: qualitative simulation of genetic regulatory networks. Bioinformatics. 2003, 19: 336-344. 10.1093/bioinformatics/btf851.View ArticlePubMedGoogle Scholar
- Mendes P: GEPASI: a software package for modelling the dynamics, steady states and control of biochemical and other systems. Comput Appl Biosci. 1993, 9: 563-571.PubMedGoogle Scholar
- Loew LM, Schaff JC: The Virtual Cell: a software environment for computational cell biology. Trends Biotechnol. 2001, 19: 401-406. 10.1016/S0167-7799(01)01740-1.View ArticlePubMedGoogle Scholar
- Olivier BG, Rohwer JM, Hofmeyr JH: Modelling cellular systems with PySCeS. Bioinformatics. 2005, 21: 560-561. 10.1093/bioinformatics/bti046.View ArticlePubMedGoogle Scholar
- Kauffman SA: Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969, 22: 437-467. 10.1016/0022-5193(69)90015-0.View ArticlePubMedGoogle Scholar
- Thomas R: Boolean formalization of genetic control circuits. J Theor Biol. 1973, 42: 563-585. 10.1016/0022-5193(73)90247-6.View ArticlePubMedGoogle Scholar
- Glass L, Kauffman SA: The logical analysis of continuous, non-linear biochemical control networks. J Theor Biol. 1973, 39: 103-129. 10.1016/0022-5193(73)90208-7.View ArticlePubMedGoogle Scholar
- Kauffman SA: The origins of order: self organization and selection in evolution. 1993, New York: Oxford University PressGoogle Scholar
- Glass L: Classification of biological networks by their qualitative dynamics. J Theor Biol. 1975, 54: 85-107. 10.1016/S0022-5193(75)80056-7.View ArticlePubMedGoogle Scholar
- Li S, Assmann SM, Albert R: Predicting essential components of signal transduction networks: a dynamic model of guard cell abscisic acid signaling. PLoS Biol. 2006, 4: e312-10.1371/journal.pbio.0040312.PubMed CentralView ArticlePubMedGoogle Scholar
- Sokol L, Loughran TP: Large Granular Lymphocyte Leukemia. Oncologist. 2006, 11: 263-273. 10.1634/theoncologist.11-3-263.View ArticlePubMedGoogle Scholar
- Klebanoff CA, Gattinoni L, Restifo NP: CD8+ T-cell memory in tumor immunology and immunotherapy. Immunological Reviews. 2006, 211: 214-224. 10.1111/j.0105-2896.2006.00391.x.PubMed CentralView ArticlePubMedGoogle Scholar
- Lamy T, Liu JH, Landowski TH, Dalton WS, Loughran TP: Dysregulation of CD95/CD95 Ligand-Apoptotic Pathway in CD3+ Large Granular Lymphocyte Leukemia. Blood. 1998, 92: 4771-4777.PubMedGoogle Scholar
- Epling-Burnette PK, Bai F, Wei S, Chaurasia P, Painter JS, Olashaw N, Hamilton A, Sebti S, Djeu JY, Loughran TP: ERK couples chronic survival of NK cells to constitutively activated Ras in lymphoproliferative disease of granular lymphocytes (LDGL). Oncogene. 2004, 23: 9220-9229.PubMedGoogle Scholar
- Schade AE, Powers JJ, Wlodarski MW, Maciejewski JP: Phosphatidylinositol-3-phosphate kinase pathway activation protects leukemic large granular lymphocytes from undergoing homeostatic apoptosis. Blood. 2006, 107: 4834-4840. 10.1182/blood-2005-08-3076.PubMed CentralView ArticlePubMedGoogle Scholar
- Epling-Burnette PK, Liu JH, Catlett-Falcone R, Turkson J, Oshiro M, Kothapalli R, Li Y, Wang J-M, Yang-Yen H-F, Karras J, Jove R, Loughran TP: Inhibition of STAT3 signaling leads to apoptosis of leukemic large granular lymphocytes and decreased Mcl-1 expression. J Clin Invest. 2001, 107: 351-362. 10.1172/JCI9940.PubMed CentralView ArticlePubMedGoogle Scholar
- 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. Proc Natl Acad Sci USA. 2008, 105: 16308-16313. 10.1073/pnas.0806447105.PubMed CentralView ArticlePubMedGoogle Scholar
- Mills KH: Regulatory T cells: friend or foe in immunity to infection?. Nat Rev Immunol. 2004, 4: 841-855. 10.1038/nri1485.View ArticlePubMedGoogle Scholar
- Thakar J, Pilione M, Kirimanjeswara G, Harvill ET, Albert R: Modeling systems-level regulation of host immune responses. PLoS Comput Biol. 2007, 3: e109-10.1371/journal.pcbi.0030109.PubMed CentralView ArticlePubMedGoogle Scholar
- Kirimanjeswara GS, Mann PB, Harvill ET: Role of antibodies in immunity to Bordetella infections. Infect Immun. 2003, 71: 1719-1724. 10.1128/IAI.71.4.1719-1724.2003.PubMed CentralView ArticlePubMedGoogle Scholar
- Harvill ET, Cotter PA, Miller JF: Pregenomic comparative analysis between bordetella bronchiseptica RB50 and Bordetella pertussis tohama I in murine models of respiratory tract infection. Infect Immun. 1999, 67: 6109-6118.PubMed CentralPubMedGoogle Scholar
- Pilione MR, Harvill ET: The Bordetella bronchiseptica type III secretion system inhibits gamma interferon production that is required for efficient antibody-mediated bacterial clearance. Infect Immun. 2006, 74: 1043-1049. 10.1128/IAI.74.2.1043-1049.2006.PubMed CentralView ArticlePubMedGoogle Scholar
- Mendoza L, Thieffry D, Alvarez-Buylla ER: Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis. Bioinformatics. 1999, 15: 593-606. 10.1093/bioinformatics/15.7.593.View ArticlePubMedGoogle Scholar
- Espinosa-Soto C, Padilla-Longoria P, Alvarez-Buylla ER: A gene regulatory network model for cell-fate determination during Arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles. Plant Cell. 2004, 16: 2923-2939. 10.1105/tpc.104.021725.PubMed CentralView ArticlePubMedGoogle Scholar
- Li F, Long T, Lu Y, Ouyang Q, Tang C: The yeast cell-cycle network is robustly designed. Proc Natl Acad Sci USA. 2004, 101: 4781-4786. 10.1073/pnas.0305937101.PubMed CentralView ArticlePubMedGoogle Scholar
- Sanchez L, Thieffry D: A logical analysis of the Drosophila gap-gene system. J Theor Biol. 2001, 211: 115-141. 10.1006/jtbi.2001.2335.View ArticlePubMedGoogle Scholar
- Sanchez L, Thieffry D: Segmenting the fly embryo: a logical analysis of the pair-rule cross-regulatory module. J Theor Biol. 2003, 224: 517-537. 10.1016/S0022-5193(03)00201-7.View ArticlePubMedGoogle Scholar
- Ghysen A, Thomas R: The formation of sense organs in Drosophila: a logical approach. Bioessays. 2003, 25: 802-807. 10.1002/bies.10311.View ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.