# CABeRNET: a Cytoscape app for augmented Boolean models of gene regulatory NETworks

- Andrea Paroni†
^{1}, - Alex Graudenzi†
^{1, 2}Email author, - Giulio Caravagna†
^{1, 3}, - Chiara Damiani
^{1, 4}, - Giancarlo Mauri
^{1, 2, 4}and - Marco Antoniotti
^{1, 5}

**17**:64

https://doi.org/10.1186/s12859-016-0914-z

© Paroni et al. 2016

**Received: **20 May 2015

**Accepted: **27 January 2016

**Published: **4 February 2016

## Abstract

### Background

Dynamical models of gene regulatory networks (GRNs) are highly effective in describing complex biological phenomena and processes, such as cell differentiation and cancer development. Yet, the topological and functional characterization of real GRNs is often still partial and an exhaustive picture of their functioning is missing.

### Results

We here introduce CABeRNET, a Cytoscape app for the generation, simulation and analysis of Boolean models of GRNs, specifically focused on their augmentation when a only partial topological and functional characterization of the network is available. By generating large ensembles of networks in which user-defined entities and relations are added to the original core, CABeRNET allows to formulate hypotheses on the missing portions of real networks, as well to investigate their generic properties, in the spirit of complexity science.

### Conclusions

CABeRNET offers a series of innovative simulation and modeling functions and tools, including (but not being limited to) the dynamical characterization of the gene activation patterns ruling cell types and differentiation fates, and sophisticated robustness assessments, as in the case of gene knockouts. The integration within the widely used Cytoscape framework for the visualization and analysis of biological networks, makes CABeRNET a new essential instrument for both the bioinformatician and the computational biologist, as well as a computational support for the experimentalist. An example application concerning the analysis of an augmented T-helper cell GRN is provided.

### Keywords

Gene regulatory networks Cell differentiation Attractors Cancer development Network augmentation Robustness analysis## Background

Consistently with the increasing availability of *big data* regarding biological systems, is the need of *mathematical* and *computational models* aimed at their effective analysis and interpretation [1, 2]. Many methodologies aim at *inferring such models from the data*, with the final goal of selecting a unique descriptive model for a phenomenon of interest. Other techiques, instead, explore the space of all possible models via a *generative approach*, with the aim of *identifying the common characteristics*, properties and regularities of those models that are “consistent” with the data. The rationale underlying the latter approach is that, although lack of reliable data might prevent us to infer the exact true model, the statistical analysis of *ensembles* of *plausible* models can provide fundamental insights about the reference phenomenon, its origin and, in specific cases, its evolutionary history. This methodology is typical of *complex systems science*, which borrows ideas and techniques from *statistical physics*, in order to focus on the *emergent* dynamical behaviours and the so-called *generic* or *universal* properties of real-life phenomena, often by means of simplified mathematical and computational models [3, 4].

Within this context, one of the best examples, involving genomic data, is provided by Boolean models of *gene regulatory networks* (GRNs), which have repeatedly proved fruitful in describing key properties of real systems, as well as in providing cues and hints for wet-lab experiments (see, e.g., [3, 5–12]). The simulation of partially characterized regulatory architectures with a Boolean approach, in particular, has recently gained attention (see, e.g., [13–21]). A first motivation lies in the inherently “dynamical” nature of gene (de)regulation processes, and in the clear limitations of a “static” analysis capturing only a partial picture of such complex processes. For example, a structural analysis of the genomic interactions might preclude to determine the influence of a *target-selective therapy* on the overall GRN interplay ruling *tumorigenesis*. Moreover, as compared to other quantitative approaches, such as ODEs models (see [22] for a review on GRN modeling), the Boolean abstraction allows for a clear and effective characterization of the gene activation patterns (or *attractors* in the terminology of dynamical systems) characterizing the different phenotypic functions, such as cell types, modes and fates, under the metaphor of “emergent collective behaviours” [3, 23–31]. In this context, phenomena such as tumorigenesis can be explained as rare emergent pattern, triggered by *signals*, stochastic fluctuations and *biological noise* (see, e.g., [32, 33]).

In line with this approach, we here continue on a recent research strand that has involved the development of simulation and analysis tools for dynamical Boolean GRNs (see, e.g., [34–36]). The foundations of our theoretical framework lay in the seminal work by Stuart Kauffman on *Random Boolean Networks* (RBNs) [37] and, more recently, on the dynamical model of cell differentiation introduced in [38, 39] and based on *Noisy Random Boolean Networks* (NRBNs). In this theoretical framework, cell types are associated to dynamical gene activation patterns and differentiation to cell-specific noise-resistance mechanisms, the underlying hypothesis being that such mechanisms refine along with differentiation stages (see [40] and references therein). In the Background Section the main features of this approach are outlined.

In this paper we introduce CABERNET, a Cytoscape [41] app to generate, import, simulate and analyse Boolean models of GRNs. In addition to the simulation of specific regulatory architectures, a key feature of this new tool lies in the possibility of generating and analysing ensembles of GRN models sharing user-defined structural and dynamical features (such as, for instance, a specific degree distribution, a set of plausible regulatory functions, or a particular differentiation scheme), in order to investigate the generic properties of classes of GRNs, in the spirit of the previously mentioned ensemble-based generative approach. Most importantly, CABERNET allows to *augment* partially-characterized GRNs by adding user-defined entities and relations. The underlying motivation is that, despite the increasing knowledge on gene regulation in real organisms, the topological and functional characterization of real networks is stil far from being comprehensive, thus preventing to capture the complexity of the overall interplay. CABERNET is conceived to randomly generate the missing portions of these partially-characterized GRNs, in order to accommodate the advantages of the ensemble when we want to study GRNs that share a common core and other structural and functional parameters. Clearly, this shall allow to test hypotheses on the yet unknown missing portions of real networks.

In addition, CABERNET guides the user through various robustness analyses of the GRNs, which can be matched against genomic experimental data. Notice that CABERNET can also generate completely random GRNs with defined structural and functional constraints, as well as simulate the dynamics of completely characterized GRNs.

The paper is organized as follows. In the Background Section a quick overview of the GRN dynamical model is presented. In the Implementation Section the main features and functionalities of CABERNET are described. As a proof of principle, in the Results Section we present the augmentation of the T-helper cell signaling network and describe the analysis of its dynamics, with particular focus on the emergent differentiation scheme and its robustness against the knockout of specific genes. Conclusions are presented in the last Section.

## Background: the GRN model

In CABERNET, GRNs are represented as Noisy Random Boolean Networks (NRBNs) [39] (for an extended description of the model and a recent review on Boolean modeling please refer to [42]). Canonical RBNs are oriented graphs in which Boolean nodes represent (either active or inactive) genes, and edges stand for regulatory interactions, modeled as *boolean functions*. The dynamics of RBNs is synchronous and deterministic, and the activation value of each node is updated at every discrete time step, according to a node-specific Boolean function depending on the value of the input nodes. Accordingly, the overall dynamics eventually ends up in (at least unitary) state cycles, termed *attractors*, which represent the emerging gene activation patterns displayed by the GRN^{1}

In the more recent NRBN model [39], transitions among attractors are made possible as a consequence of random modifications of the node’s activation value, lasting a defined time span (i.e., *flips*). Thus, it is possible to determine the so-called *Attractors Transition Network* (ATN, or *matrix*, ATM), i.e., a *stability matrix* displaying the noise-induced transition probability among attractors.

*Threshold Ergodic Sets*, TESs) represent the gene patterns of specific cell types. Besides, given that it is hypothesized that less differentiated cell types are characterized by less refined noise-control mechanisms (see, e.g., [43–46]), different thresholds are introduced to exclude those transitions that are unlikely to occur for a cell in a specific differentiation stage. Accordingly, each threshold determines a set of (disconnected) TESs representing the distinct cell types of a specific differentiation stage, starting from toti-/multi-potent stem cells (generally characterize by a single TES including many attractors) to intermediate states, to fully differentiated cells (generally characterized by one attractor-TESs). Thus, such an emergent hierarchical structure resembles and can be compared to that of differentiation trees in real cells. Notice that the approach is general as it does not refer to any specific organism (see Fig. 1 for a simplified representation of the model).

## Implementation

Main functionalities and parameters of CABERNET. A schematic representation of the various functions and parameters of CABERNET is provided, as explicitly described in the Implementation Section in the main text. For a thorough explanation please refer to the user manual (see Availability)

### Input GRNs (generation, import and augmentation)

**[Random network generation]**
CABERNET can randomly generate and simulate ensembles of GRNs with certain structural parameters such as: (*i*) number of genes; (*i*
*i*) ingoing/outgoing GRN’s topology, e.g., *fixed*, *Erd*
\(\ddot {o}\)
*s-Rényi’s random* [47], *Barabasi-Alberts’ scale-free* (either via *preferential attachment* or *power-law* generation) [48] or *Watts-Strogatz’s small-world* [49]; (*i*
*i*
*i*) type of regulation functions (node-specific Boolean functions). Concerning the latter, these can be set by the user or randomly generated to accomodate: (*i*)*bias-based random* functions, (*i*
*i*)*canalyzing* [9] or (*i*
*i*
*i*)*logical* functions - expressed in the canonical AND/OR notation. Network and simulation parameters (e.g., samples size) can be defined either via an input form or a textual file^{2}.

**[Import of characterized GRNs]** GRNs that are characterized with respect to both the topology and the regulatory functions, can be loaded in CABERNET via input textual file. Limited to the network topology, CABERNET can also directly import GRNs from any public database that Cytoscape interfaces with or from any database that allow export in SBML format.

**[Augmentation of GRNs]** Any GRN loaded in the tool (see above) can be augmented by CABERNET, by randomly generating a chosen number of augmented networks, in which entities and relations are added to the input network according to user-defined structural and functional parameters. Parameters for augmentation are the same that must be defined for the random generation (see above). The resulting ensemble of augmented networks will share the input topological and functional core and will differ for the randomly generated portion.

For instance, in the Results Section, the *T-helper* cell signaling network curated from [20] (40 genes and 51 regulatory interactions), is augmented with 160 further nodes and 349 edges according to a Erdos-Rényi random topology. Assessment of the functional effect of this augmented network is also discussed.

### GRN’s dynamics simulation

Given that the space of the possible configurations of a NRBN can be dramatically large (there are 2^{
N
} possible configurations for a network with *N* nodes), CABERNET can simulate the dynamics of a network by either: (*i*) uniformly sampling the initial conditions to test (via a required user-defined parameter) or (*i*
*i*) performing an exhaustive search (for small networks only). CABERNET allows to investigate key statistics of the emerging attractors such as, e.g., number, length, robustness and reachability. The stability of any pattern to perturbations can be assessed either via *temporary flips* (with duration of 1 step) or via *permanent* gene knock-in/knock-out.

Following the simulation, CABERNET can compute and display the threshold-dependent ATN (here called the *TES network*) for specific threshold values. Different views on such a network are available in the tool and, for instance, allow to display the genes’ configurations in a pattern, and the variation of the number of TESs alongside thresholds, as proposed in [51].

### GRN-selection constrained by differentiation scheme

One might look for a GRN giving rise to a specific differentiation tree, as done e.g. in [35]. Trees can be inputed to CABERNET in textual format or from the Cytoscape active window; CABERNET can select those NRBNs that display a differentiation tree structurally similar to the loaded one, where the measure of similarity is defined by the user. This feature is implemented as a batch process scanning among generated NRBNs. Notice that, usually, a single NRBN exhibits various emerging trees, according to the various possible combinations of thresholds thus this feature is the computationally most demanding in CABERNET.

Besides, the statistically *representative* differentiation tree(s) of each specific network, defined as the most frequent emerging tree for different (sampled) threshold combinations, provided a specific tree depth, can be computed.

### Visualization

Each computational task is tracked by a progress bar. Once the simulation of the dynamics is completed, the powerful visualization capabilities of Cytoscape can be used to analyze the topological and dynamical properties of the networks.

In particular, with CABERNET it is possible to visualize: (*i*) the NRBNs, (*i*
*i*) the *attractor graph network*, in which all the states of the attractors and the transitions among them are displayed, (*i*
*i*
*i*) the threshold dependent ATNs and (*i*
*v*) the representative differentiation tree. By clicking on a specific network, it gets visualized within Cytoscape, so that it can be further analyzed.

Different network styles have been defined and can be selected: (*i*)CABERNET
*network*, aimed at visualizing the properties of the NRBN: the color of each node being related to the Boolean function bias and the size of each node proportional to its degree, (*i*
*i*)CABERNET
*attractors*, for the visualization of the attractor graph network, (*i*
*i*
*i*)CABERNET
*TES*, for the visualization of the ATN and (*i*
*v*)CABERNET
*collapsed TES*, for the collapsed visualization of the ATN: in these last two styles, the edge size is proportional to the transition probability.

### Robustness analysis

Different kind of perturbations can be simulated to assess the robustness of the attractors of a network. In particular, it is possible to perform a user-defined number of (*i*) temporary (i.e., flips) or (*i*
*i*) permanent (i.e., knock-in/knock-out) perturbations on (*i*) a chosen number of randomly selected nodes or (*i*
*i*) specific nodes. The robustness analysis can be performed on single networks or on the whole ensemble of simulated GRNs.

Network’s stability is assessed via robustness analyses, by means of standard measures such as *avalanches* (i.e., the number of nodes whose activation pattern is different in a perturbation experiment with respect to the wild type scenario) and *sensitivity* (i.e., the number of perturbation experiments in which a certain node’s pattern is affected) [11]. The results of the analyses can be exported in csv files.

### Network analysis

CABERNET offers a wide range of network-specific statistics. These include (*i*) the distribution of the attractors’ lengths, (*i*
*i*) the basins of attraction, (*i*
*i*
*i*) the proportion of frozen and oscillating nodes, plus other classical network measures such as *clustering coefficient*, *network diameter* and *average path length*. All the statistics can be visualized and exported; further network measures are accessible via the network analysis tools included in Cytoscape.

### Outputs

All the networks and the relative topological, functional and dynamical properties can be exported as textual files, from both the Wizard and the Function menu. For instance, the complete topological and functional description of the networks can be exported so that it can be used in simulation environments external to Cytoscape such as, e.g., CHASTE [52], CompuCell3D [53] or the simulator described in [35] (see [54] for a recent review on multiscale models of multicellular systems). Also, it is possible to export the complete description of all the attractor states, as well as information of their basins on attraction.

### Comparison with other tools

Software comparison. Comparison of the main features implemented in the principal tools for the qualitative simulation of the dynamics of GNRs [34, 55–58, 61–65]. Green cells and the symbol ‘V’ indicate feature that are implemented, as opposed to red cells and the symbol ‘X’. When the assignation is not neat a footnote provides further remarks

We conclude by noting that CABERNET is a follow-up of our previous software GESTODIFFERENT [34], an earlier GRN tool developed by our group. Novelties include, but are not limited to, the following key functionalities and simulation options: the network augmentation module, the simulation of specific regulatory architectures, the possibility of querying publicly available databases, new network topologies, functions and constraints, the robustness analysis and the network structural analysis modules, new visualization styles and new I/O functions.

### Performance evaluation and computational complexity

Performance evaluation

Ordered networks ( | Critical networks ( | Chaotic networks ( | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Nodes | 100 | 500 | 1000 | 5000 | 100 | 500 | 1000 | 5000 | 100 | 500 | 1000 | 5000 |

Avg. (sec) | 10.20 | 53.66 | 101.82 | 603.56 | 9.61 | 52.53 | 109.65 | 746.80 | 13.62 | 73.22 | 150 | 1341.94 |

St. Dev. (sec) | 1.07 | 1.76 | 2.45 | 29.40 | 0.40 | 1.74 | 3.82 | 203.43 | 0.48 | 2.18 | 6.91 | 39.58 |

We also remark that it is not possible to provide a precise upper bound for the size of the networks to simulate CABERNET, because, as it is known from RBN/NRBN literature, the dynamical behaviour of even small networks can be dramatically heterogeneous, strongly depending on the *dynamical regime*, which is defined by a series of key structural parameters (for a discussion on this topic pleas refer to, e.g., [3]). This aspect, which is particularly relevant when dealing with critical and chaotic networks, is reflected on the computation of the attractors and their stability, as well as on the optional tree-matching procedure. Accordingly, computation time cut-offs can be set to avoid excessively long simulation. Clearly, the generation of large ensemble of GRNs, as opposed to the simulation of single architectures, increases the risk of encountering a dynamically complex network. As a general comment, CABERNET can generate and simulate networks with a few thousands nodes, yet networks with higher average connectivity are more likely to display very heterogeneous behaviours.

In regard to the network augmentation function, its computational complexity corresponds with that of simple incremental approaches. Indeed, our general task is that of augmenting a network with *n* nodes and *m* edges, regardless its current topology, to one with *N*>*n* nodes and *M*>*m* edges. Thus its cost is that of adding *N*−*n* new nodes and *M*−*m* new edges; the way we select edges determines the Cabernet’s network type. No matter what, we can exploit standard generative algorithms with no computational overload^{3}.

## Results

Our group has recently been focusing on the investigation of the dynamical properties of multicellular systems via multiscale simulations, with particular attention to the conditions that would favor the emergence and development of tumors. To this end, CABERNET was recently used to generate, simulate and visualize the GRNs ruling the behaviour of an intestinal crypt in CHASTE’s multiscale simulation engine, allowing to identify conditions for cancer’s emergence and crypt’s colonization [36]. In the following, we propose a further example to show some of the potential applications of CABERNET.

**[Augmentation of T-helper signaling network]** The signaling network of human T-helper cells was recently characterized with respect to both the topology and the regulatory functions. In [59] the dynamics of such network was simulated with a Boolean approach and it was shown that the attractors actually reproduce real gene activation patterns of distinctly differentiated T-helper cell types.

*i*) generate a large ensemble of random networks with the T-helper functional core, (

*i*

*i*) select only those networks in which the emergent dynamical behaviour is in accordance with the hematopoietic differentiation tree, in which the T-helper cell type is supposed to be one of the leaves, and (

*i*

*i*

*i*) eventually detect and analyze the generic properties of the selected networks, which result as

*plausible*GRNs ruling the hematopoietic differentiation fate, specifically investigating their robustness (see Fig. 2).

The distinct networks differ for the augmented portion, which is randomly generated with structural parameters (i.e., Erdos-Renyi random topology, average connectivity =2, random Boolean functions with bias =0.5) that are classically used in similar studies (see, e.g., [11]). Accordingly, only certain networks will eventually display the desired emergent dynamical behavior. The underlying idea is that the *matching* networks might allow for the formulation of hypotheses on the missing portions of the relevant GRN ruling the overall hematopoiesis process. Besides, the characterization of the attractors could be matched with the real gene activation patterns driving the functioning of the various cell types^{4}.

In the experiment, a large number of distinct augmented networks was generated and simulated with CABERNET: only 1 on a total of 600 augmented NRBNs actually displayed the expected dynamical behavior, i.e., the correct differentiation tree, hinting at the complexity of the tuning process driven by the evolutionary pressure that led to the topology of current GRNs^{5}.

In Fig. 2 one can see the original T-helper signaling network, originally mapped in [20], and the NRBN that was selected as correct, in which the augmented portion is highlighted. Notice that, in the augmented network, both the topology and the functions of the original core are slightly different from those of the T-helper GRN, as a consequence of adding new relations linking the new and the original portions of the net. The visualization of the network is provided via CABERNET, by applying the suitable styles (see above).

In Fig. 2 one can notice that this specific network exhibits 8 distinct attractors (each one characterized by a length equal to 8 NRBN time steps). By pruning the ATN with increasingly larger thresholds, the TES at the higher level, including all the 8 attractors connected by noise-induced transitions and representing multi-/toti-potent cells, progressively splits in TESs enclosing an increasingly lower number of attractors, up to the 7 TESs at the lower level, which correspond to single attractors, when the threshold is equal to 1. This network was specifically selected as plausible because the resulting emergent differentiation tree matches that of hematopoietic cells (taken from [60], see Fig. 2)^{6}.

Note that, in case augmented networks were more likely to display a matching emergent tree, one may exploit CABERNET to perform ensemble-level analyses on the matching set, aimed at the formulation of hypotheses on the generic properties of real networks.

*knockout experiments*(KO), by forcing the specific Boolean function of each gene to inactivation (i.e., 0 output for any regulatory input), and we tried to match the resulting differentiation trees with that of hematopoietic cells. Remarkably, in 35 cases (88 %), the KO experiment resulted in a mismatching tree, hinting at the role of those specific genes in the interplay leading to the emergence of the hematopoietic tree. The similarity among the hematopoietic tree,

*h*, and a tree

*T*resulting from a KO experiment was measured as follows:

where *l*
^{∗} and *k*
^{∗} are the maximum depth and the maximum number of a node’s children in both *h* and *T*. Function *n*
_{
x
}(*k,l*) returns the number of nodes at level *l* with *k* children in tree *x*; thus, this quantity measures the structural *level-by-level similarity* of two trees by assessing the number of parents with *k* children, per level. Since we focus on differentiation trees, this can be interpreted as a measure of the ability of a certain cell type, a progenitor, to differentiate in a set of distinct subtypes.

In Fig. 2 we show values of \(\widehat {d}\) in our experiments; the lower the value the closer is *T* to *h*. Values of \(\widehat {d}\) range around 9, with a maximum of 17 and minimum 0; in 8 cases, the value lower than 5 suggests a close similarity between the emergent and the hematopoietic tree. Besides, the dynamics turned out to be completely insensitive (i.e., \(\widehat {d}= 0\)) to the KO of 5 specific genes, i.e., *CRE*, *Ca*, *PLCg-a*, *Fos* and *Gads*, which, accordingly, might be not relevant in the differentiation process. Clearly, further investigations are needed to corroborate this hypothesis.

We finally remark that the generated networks could be used within any multiscale simulation frameworks, in order to investigate, e.g., the processes of homeostasis and clonal expansion, as proposed in [35, 36].

## Conclusions

In this work we introduced CABERNET - a new Cytoscape app for the generation, simulation and analysis of augmented Boolean models of gene regulatory networks - and described some of its key functionalities, as well as an example application to real GRN data.

CABERNET is the result of a long-time effort aimed at bridging different fields and disciplines, such as computer science, statistics and complex systems science, for the effective study of complex biological systems. The numerous modeling and simulation functionalities, the various effective analysis tools and the fine integration within the widely used Cytoscape framework, might settle the ground for CABERNET becoming a powerful instrument for bioinformaticians and computational biologists, especially in providing a computational support for experimentalists.

In particular, CABERNET can provide an essential tool to effectively investigate key and still partially undeciphered biological phenomena, such as, e.g., gene regulation, cell differentiation and tumorigenesis, with particular focus on the properties of dynamical gene activation patterns and their relation with biological noise.

## Availability and requirements

**Project name:**
CABERNET: a Cytoscape app for the generation and the Analysis of Boolean models of gene Regulatory NETworks **Version:** 1.1 **Plugin website:**
http://bimib.disco.unimib.it/index.php/CABERNET
**Operating systems:** platform independent **Software requirement:** Cytoscape 3.x (http://www.cytoscape.org/) **Programming language:** Java **License:** BSD-like license (see website)

## Endnotes

^{1} We remark that one of the key advantages of employing a Boolean modeling approach lies in the possibility of including entities involved in the regulatory interplay other than genes and proteins, such as, e.g., non coding RNAs. In fact, any gene product can be considered as a Boolean entity (i.e., present/absent) interacting with other genes/products. For instance, miRNAs - non coding RNAs characterized by inhibitory functions that are able to modulate gene expression and are supposed to confer robustness against biological noise - might be represented by associating a canalyzing inhibitory function [9] to their target genes.

^{2} An extension of the software to import/export GRN models in the SBML Qual file format [50] is currently underway.

^{3} Erdos–Renyi random network are generated regardless the initial structures. For the Barabasi-Albert preferential attachment model the initial network determines the attachment probability, via the degree. For the power-law network, we sample the number of edges to assign to each node from a power law, and augment the original network accordingly. For the small word, we operate similarly. Constraint such as fixing the number of edges for a node can be set with no computational overload of this generative procedure – see the Manual for details.

^{4} Notice that the quest for an effective characterization of the relation between the functional modules of biochemical networks (i.e., sub-graphs in wider GRN models) and the phenotypic functions, in normal and aberrant cells, is nowadays central in biomedical research. This is not in contrast with the complexity of the regulation/differentiation complexity, whereas, instead, highlights the importance of the evolutionary pressure-based mechanisms that led to the specialization, modularization, hierarchical organization and plasticity of current GRNs. To this end, CABERNET can provide an effective instrument to investigate the properties of such modules, as in the presented example.

^{5} The analysis of the possible relation between the topology of the augmented portion, e.g., random, scale-free or small-world, and the likelihood of displaying the correct differentiation tree could be easily performed with the tool. We leave this interesting further application of CABERNET to future works.

^{6} The hematopoietic differentiation scheme is characterized by a multi-potent progenitor (*MPP*; antecedent hematopoietic stem cells, HSCs, are not shown in the scheme), with the potential to differentiate into two lineages, i.e., common myeloid progenitor (*CMP*) and common lymphoid progenitor (*CLP*). CMP further divide into megakaryocyte-erythroid progenitor (*MEP*) and granulocyte/monocyte progenitor (*GMP*), finally committing to mature blood cells including erythrocytes (*EC*), megakaryocyte (*MK*), monocyte (*M*) and granulocytes, i.e. neutrophils, eosinophils, basophils (*N/E/B*). Conversely, CLP further differentiate into B-cell progenitors (*B PROG*) and T-cell and natural killer cell progenitors (*T/NK PROG*), with a final commitment to mature B cells (*B*), T cells (*T*) and NK cells (*NK*) [60].

## Declarations

### Acknowledgements

This project was partially supported by the ASTIL program, project “RetroNet”, grant n. 12-4-5148000-40; U.A 053, and by NEDD Project [ID14546A Rif SAL-7] Fondo Accordi Istituzionali 2009.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## Authors’ Affiliations

## References

- Kitano H. Foundations of Systems Biology. Cambridge: MIT Press; 2001, pp. 1–36.Google Scholar
- Kitano H. Computational systems biology. Nature 420.6912. 2002:206–210.Google Scholar
- Kauffman SA. At Home in the Universe: Oxford University Press; 1995.Google Scholar
- Kaneko K. Life: An Introduction to Complex Systems Biology. Berlin Heidelberg: Springer-Verlag; 2006.Google Scholar
- Shmulevich I, Dougherty ER, Kim S, Zhang W. Probabilistic boolean networks: a rule-based uncertainty model for gene regulatory networks. Bioinformatics. 2002; 18(2):261–74.View ArticlePubMedGoogle Scholar
- Shmulevich I, Dougherty ER, Zhang W. From boolean to probabilistic boolean networks as models of genetic regulatory networks. Proc IEEE. 2002; 90(11):1778–92.View ArticleGoogle Scholar
- Kauffman SA, Peterson C, Samuelsson B, Troein C. Random boolean network models and the yeast transcriptional network. Proc Nat Acad Sci USA. 2003; 100:14796–9.PubMed CentralView ArticlePubMedGoogle Scholar
- Serra R, Villani M, Semeria A. Genetic network models and statistical properties of gene expression data in knock-out experiments. J Theor Biol. 2004; 227:149–57.View ArticlePubMedGoogle Scholar
- Kauffman SA, Peterson C, Samuelsson B, Troein C. Genetic networks with canalyzing boolean rules are always stable. Proc Nat Acad Sci USA. 2004; 101:17102–7.PubMed CentralView ArticlePubMedGoogle Scholar
- Ramo P, Kesseli Y, Yli-Harja O. Perturbation avalanches and criticality in gene regulatory networks. J Theor Biol. 2006; 242:164–70.View ArticlePubMedGoogle Scholar
- Serra R, Villani M, Graudenzi A, Kauffman SA. Why a simple model of genetic regulatory networks describes the distribution of avalanches in gene expression data. J Theor Biol. 2007; 249:449–60.View ArticleGoogle Scholar
- Cheng D, Qi H, Li Z. Analysis and Control of Boolean Networks: a Semi-tensor Product Approach. London: Springer-Verlag; 2010.Google Scholar
- Sanchez L, van Helden J, Thieffry D. Establishement of the dorso-ventral pattern during embryonic development of drosophila melanogasater: a logical analysis. J Theor Biol. 1997; 189:377–89.View ArticlePubMedGoogle Scholar
- Akutsu T, Miyano S, Kuhara S. Identification of genetic networks from a small number of gene expression patterns under the Boolean network model. In: Pacific symposium on biocomputing. Vol. 4. 1999., p. 17–28. http://psb.stanford.edu/psb-online/proceedings/psb99/.
- Sanchez L, Thieffry D. A logical analysis of the drosophila gap-gene system. J Theor Biol. 2001; 211:115–41.View ArticlePubMedGoogle Scholar
- Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in drosophila melanogaster. J Theor Biol. 2003; 223:1–18.View ArticlePubMedGoogle Scholar
- Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cell-cycle network is robustly designed. Proc Nat Acad Sci USA. 2004; 101(14):4781–786.PubMed CentralView ArticlePubMedGoogle Scholar
- Chaos A, Aldana M, Espinosa-Soto C, de León BGP, Arroyo AG, Alvarez-Buylla ER. From genes to flower patterns and evolution: dynamic models of gene regulatory networks. J Plant Growth Regul. 2006; 25(4):278–89.View ArticleGoogle Scholar
- Fauré 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(14):124–31.View ArticleGoogle Scholar
- Klamt S, Saez-Rodriguez J, Lindquist JA, Simeoni L, Gilles ED. A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinformatics. 2006; 7(1):56.PubMed CentralView ArticlePubMedGoogle Scholar
- Davidich MI, Bornholdt S. Boolean network model predicts cell cycle sequence of fission yeast. PloS one. 2008; 3(2):1672.View ArticleGoogle Scholar
- Karlebach G, Shamir R. Modelling and analysis of gene regulatory networks. Nat Rev Mol Cell Biol. 2008; 9(10):770–80.View ArticlePubMedGoogle Scholar
- Huang S, Ingber DE. Shape-dependent control of cell growth, differentiation, and apoptosis: switching between attractors in cell regulatory networks. Exp Cell Res. 2000; 261(1):91–103.View ArticlePubMedGoogle Scholar
- Huang S, Eichler G, Bar-Yam Y, Ingber DE. Cell fates as high-dimensional attractor states of a complex gene regulatory network. Phys Rev Lett. 2005; 94(12):128701.View ArticlePubMedGoogle Scholar
- Forgacs G, Newman SA. Biological Physics of the Developing Embryo: Cambridge University Press; 2006.Google Scholar
- Huang S, Ernberg I, Kauffman SA. Cancer attractors: a systems view of tumors from a gene network dynamics and developmental perspective.Semin Cell Dev Biol. 2009; 20(7):869–76.PubMed CentralView ArticlePubMedGoogle Scholar
- Canham MA, Sharov AA, Ko MSH, JM B. Functional heterogeneity of embryonic stem cells revealed through translational amplification of an early endodermal transcript. PLoS Biol. 2010; 8(5):1000379.View ArticleGoogle Scholar
- Felli N, Cianetti L, Pelosi E, Carè A, Liu CG, Calin GA, et al. Hematopoietic differentiation: a coordinated dynamical process towards attractor stable states. BMC Syst Biol. 2010; 4(1):85.PubMed CentralView ArticlePubMedGoogle Scholar
- Furusawa C, Kaneko K. A dynamical-systems view of stem cell biology. Science. 2012; 338(6104):215–7.View ArticlePubMedGoogle Scholar
- Cheng WY, Yang T-HO, Anastassiou D. Biomolecular events in cancer revealed by attractor metagenes. PLoS Comput Biol. 2013; 9(2):1002920.View ArticleGoogle Scholar
- Nikolov S, Wolkenhauer O, Vera J. Tumors as chaotic attractors. Mol BioSystems. 2014; 10(2):172–9.View ArticleGoogle Scholar
- Eldar A, Elowitz MB. Functional roles for noise in genetic circuits. Nature. 2010; 467:167–73.PubMed CentralView ArticlePubMedGoogle Scholar
- Tsimring LS. Noise in biology. Reports Prog Phys. 2014; 77(2):026601.View ArticleGoogle Scholar
- Antoniotti M, Bader GD, Caravagna G, Crippa S, Graudenzi A, Mauri G. Gestodifferent: a cytoscape plugin for the generation and the identification of gene regulatory networks describing a stochastic cell differentiation process. Bioinformatics. 2013; 29(4):513–14.PubMed CentralView ArticlePubMedGoogle Scholar
- Graudenzi A, Caravagna G, De Matteis G, Antoniotti M. Investigating the relation between stochastic differentiation, homeostasis and clonal expansion in intestinal crypts via multiscale modeling. PLoS ONE. 2014; 9(5):e97272. http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0097272.PubMed CentralView ArticlePubMedGoogle Scholar
- Rubinacci S, Graudenzi A, Caravagna G, Mauri G, Osborne J, Pitt-Francis J, et al. Cognac: a chaste plugin for the multiscale simulation of gene regulatory networks driving the spatial dynamics of tissues and cancer. Cancer Informatics. 2015; 14(Suppl 4):53.PubMed CentralPubMedGoogle Scholar
- Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969; 22:437–67.View ArticlePubMedGoogle Scholar
- Serra R, Villani M, Barbieri A, Kauffman SA, Colacci A. On the dynamics of random boolean networks subject to noise: attractors, ergodic sets and cell types. J Theor Biol. 2010; 265:185–93.View ArticlePubMedGoogle Scholar
- Villani M, Barbieri A, Serra R. A dynamical model of genetic networks for cell differentiation. PLoS ONE. 2011; 6(3):17703–1013710017703.View ArticleGoogle Scholar
- Hoffman M, Chang HH, Huang S, Ingber DE, Loeffler M, Galle J. Noise driven stem cell and progenitor population dynamics. PLoS ONE. 2008; 3:2922.View ArticleGoogle Scholar
- Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramadge D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13(11):2498–504.PubMed CentralView ArticlePubMedGoogle Scholar
- Albert R, Thakar J. Boolean modeling: a logic-based dynamic approach for understanding signaling and regulatory networks and for making useful predictions. Wiley Interdiscip Rev Syst Biol Med. 2014; 6(5):353–69.View ArticlePubMedGoogle Scholar
- Hu M, Krause D, Greaves M, Sharkis S, Dexter M, Heyworth C, et al. Multilineage gene expression precedes commitment in the hemopoietic system. Genes Dev. 1997; 11:774–85.View ArticlePubMedGoogle Scholar
- Hayashi K, Lopes SM, Surani MA. Dynamic equilibrium and heterogeneity of mouse pluripotent stem cells with distinct functional and epigenetic states. Cell Stem Cell. 2008; 3:391–440.View ArticlePubMedGoogle Scholar
- Furusawa C, Kaneko K. Chaotic expression dynamics implies pluripotency: when theory and experiment meet. Biol Direct. 2009; 4:17.PubMed CentralView ArticlePubMedGoogle Scholar
- Yamanaka H. Elite and stochastic models for induced pluripotent stem cell generation. Nature. 2009; 460:49–52.View ArticlePubMedGoogle Scholar
- Erdős P, Rényi A. On random graphs. Publicationes Mathematicae. 1959; 6:290–7.Google Scholar
- Barabasi A, Albert R. Emergence of scaling in random networks. Science. 1999; 286:509–12.View ArticlePubMedGoogle Scholar
- Watts DJ, Strogatz SH. Collective dynamics of ‘small-world’networks. Nature. 1998; 393(6684):440–2.View ArticlePubMedGoogle Scholar
- Chaouiya C, Bérenguier D, Keating SM, Naldi A, Van Iersel MP, Rodriguez N, et al. Sbml qualitative models: a model representation format and infrastructure to foster interactions between qualitative modelling formalisms and tools. BMC Syst Biol. 2013; 7(1):135.PubMed CentralView ArticlePubMedGoogle Scholar
- Graudenzi A, Damiani C, Paroni A, Filisetti A, Villani M, Serra R, et al. Investigating the role of network topology and dynamical regimes on the dynamics of a cell differentiation model. In: Advances in Artificial Life and Evolutionary Computation. Switzerland: Springer International Publishing: 2014. p. 151–68.Google Scholar
- Pitt-Francis J, Pathmanathan P, Bernabeu MO, Bordas R, Cooper J, Fletcher AG. Chaste: A test-driven approach to software development for biological modelling. Comput Phys Commun. 2009; 180(12):2452–71.View ArticleGoogle Scholar
- Swat M, Thomas GL, Belmonte JM, Shirinifard A, Hmeljak D, Glazier JA. 2012. Multi-Scale Modeling of Tissues Using CompuCell3D, Vol. 110.Google Scholar
- De Matteis G, Graudenzi A, Antoniotti M. A review of spatial computational models for multi-cellular systems, with regard to intestinal crypts and colorectal cancer development. J Math Biol. 2013; 66(7):1409–62.View ArticlePubMedGoogle Scholar
- Bock M, Scharp T, Talnikar C, Klipp E. Boolesim: an interactive boolean network simulator. Bioinformatics. 2014; 30(1):131–2.View ArticlePubMedGoogle Scholar
- Helikar T, Kowal B, McClenathan S, Bruckner M, Rowley T, Madrahimov A, et al. The cell collective: toward an open and collaborative approach to systems biology. BMC Syst Biol. 2012; 6(1):96.PubMed CentralView ArticlePubMedGoogle Scholar
- Naldi A, Berenguier D, Fauré A, Lopez F, Thieffry D, Chaouiya C. Logical modelling of regulatory networks with ginsim 2.3. Biosystems. 2009; 97(2):134–9.View ArticlePubMedGoogle Scholar
- Di Cara A, Garg A, De Micheli G, Xenarios I, Mendoza L. Dynamic simulation of regulatory networks using squad. BMC Bioinformatics. 2007; 8(1):462.PubMed CentralView ArticlePubMedGoogle Scholar
- Mendoza L, Xenarios I. A method for the generation of standardized qualitative dynamical systems of regulatory networks. Theor Biol Med Modelling. 2006; 3:13.View ArticleGoogle Scholar
- Lim WF, Inoue-Yokoo T, Tan KS, Lai MI, Sugiyama D. Hematopoietic cell differentiation from embryonic and induced pluripotent stem cells. Stem Cell Res Ther. 2013; 4(3):71.PubMed CentralPubMedGoogle Scholar
- Terfve C, Cokelaer T, Henriques D, MacNamara A, Goncalves E, Morris MK, et al. Cellnoptr: a flexible toolkit to train protein signaling networks to data using multiple logic formalisms. BMC Syst Biol. 2012; 6(1):133.PubMed CentralView ArticlePubMedGoogle Scholar
- Zheng J, Zhang D, Przytycki PF, Zielinski R, Capala J, Przytycka TM. Simboolnet - cytoscape plugin for dynamic simulation of signaling networks. Bioinformatics. 2010; 26(1):141–2.PubMed CentralView ArticlePubMedGoogle Scholar
- Albert I, Thakar J, Li S, Zhang R, Albert R. Boolean network simulations for life scientists. Source Code Biol Med. 2008; 3(1):1–8.View ArticleGoogle Scholar
- Hinkelmann F, Brandon M, Guang B, McNeill R, Blekherman G, Veliz-Cuba A, et al. Adam: analysis of discrete models of biological systems using computer algebra. BMC Bioinformatics. 2011; 12(1):295.PubMed CentralView ArticlePubMedGoogle Scholar
- Le DH, Kwon YK. Netds: a cytoscape plugin to analyze the robustness of dynamics and feedforward/feedback loop structures of biological networks. Bioinformatics. 2011; 27(19):2767–8.View ArticlePubMedGoogle Scholar