The inference of Boolean networks from biological data is at present coming of age. On the one hand, a methodology that has been employed to infer a number of Boolean gene regulatory networks has established itself. This methodology follows the next three steps. The first one is determining the set of genes regulating each gene, in addition to the “sign” of each regulation based on experimental results. Second, a number of constraints are considered. These constraints can be, for example, (a) a set of input-output pairs representing time-series data (Liang et al., 1998; Akutsu et al., 1999; Lähdesmäki et al., 2003; Chueh and Lu, 2012; Haider and Pal, 2012; Berestovsky and Nakhleh, 2013; Han et al., 2014; Barman and Kwon, 2017), (b) a desired set of fixed-point attractors (Albert and Othmer, 2003; Li et al., 2004; Pal et al., 2005; Mendoza, 2006; Davidich and Bornholdt, 2008; Tarissan et al., 2008; La Rota et al., 2011; Azpeitia et al., 2013, 2014), or (c) a set of temporal-logic constraints (Bernot et al., 2004; Chabrier-Rivier et al., 2004; Fages et al., 2004; Calzone et al., 2006a,b; Mateus et al., 2007; Khalis et al., 2009; Richard et al., 2012). Third, search is performed to find Boolean networks consistent with the available information.

On the other hand, there are now many methods and algorithms for the analysis of Boolean networks that have demonstrated their effectiveness. We noticed that some of these algorithms can also be used as auxiliary techniques for the inference of network dynamics. In particular, we have selected Dubrova and Teslenko's (2011) algorithm for finding cyclic attractors in synchronous Boolean networks, Garg et al.'s (2008) method for finding basins of attraction, and Naldi et al.'s (2007) formula for finding fixed-point attractors. (We use the terms “steady state dynamics” and “fixed-point attractors” interchangeably.) Griffin has combined these algorithms with novel techniques, offering a tool automating many steps of the inference of molecular network dynamics.

Boolean networks have variables representing the genes, mRNA, proteins, or any other type of molecule included in the network. Variables have only two values, and time is discrete. A state is a tuple (i.e., ordered list) of values, one for each variable. A Boolean network consists of a (finite) set of states together with a “transition” relation (sometimes also called “update” or “successor” relation).

In asynchronous Boolean networks, states may have more than one successor. In synchronous Boolean networks, by contrast, states have exactly one successor. Previous studies suggest that having multiple successors can provide a closer description of the biological phenomena and can eliminate dynamical artifacts (Garg et al., 2008). There may be advantages, nevertheless, in considering synchronous networks.

The rest of this article only deals with synchronous Boolean networks. For simplicity, we will sometimes drop the “synchronous” adjective and refer to such networks simply as “Boolean.”

In a nutshell, Griffin's methodology is the following. Griffin has as input a (generalization of a) “regulation” graph (sometimes also called “interaction” or “influence” graph), in addition to biological constraints, such as an expected set of either fixed-point or cyclic attractors of both the wild-type and of mutants of an organism. The output is typically a set of Boolean networks satisfying such constraints.

Because in Griffin the regulation graph plays a prominent role, the definition of a regulation is fundamental. Griffin uses the definition of Naldi et al. (2007), Richard et al. (2012), and Mori and Mochizuki (2017).

A regulation graph is normally an underdetermined specification of a Boolean network because even with the formal definition of regulation and sometimes also due to lack of information the regulation graph may be satisfied by many Boolean networks.

Griffin then combines the regulation graph with the biological constraints, to reduce the number of possible networks. This is essentially done by constructing a typically large Boolean formula representing: (1) a formal definition of the regulation graph and (2) the biological constraints. This formula is then given to a “SAT solver.” These solvers employ algorithms finding value assignments to the variables occurring in the formula that make such a formula true. Each such assignment represents a possible network.

Notice that both the regulation graph and many of the biological constraints rely on biological information. Thus, Griffin mainly produces biological meaningful solutions (i.e., networks that are coherent with the available experimental information).

Griffin's methodology required the development of multiple innovations, including the generalization of the concept of a regulation graph that we call R-graph, the gradual application of biological constraints at runtime using clause learning, and the handling of both mutations and partially known fixed-point attractors. Some of Griffin's innovations will be described below.

Griffin relies heavily on “symbolic” methods. Such methods do not represent sets explicitly, but rather codified as Boolean expressions. This allows for the representation of large sets of states, as well as the pruning of large fragments of the search space. It is important to observe that numerous off-the-shelf symbolic computational tools exist.

The original prototype of Griffin appeared in Rosenblueth et al. (2014), where the authors presented a formal framework for the inference of Boolean networks from a standard regulation graph by direct application of a SAT solver to Boolean formulas codifying regulations, fixed-point attractors, and single-point mutations.

Examples illustrating the use of Griffin already appeared in: Rosenblueth et al. (2014), employing sets of desired fixed-points for mutations; Weinstein et al. (2015), looking for interactions that are necessary for the existence of a cyclic attractor; García-Gómez et al. (2017), verifying if a set of regulations was enough to obtain a desired set of expected attractors; and Azpeitia et al. (2017), making extensive use of hypotheses and partially defined sets of attractors for multi-point mutations.

Features illustrated in previous articles mentioning Griffin have been left out from this work. We demonstrate in this article Griffin's new functionalities through three cases studies in section 2.2.

The rest of this article is organized as follows. The section “sect:results” has two parts: first we give a description of Griffin, and then we turn our attention to three case studies. Next follows the section “sect:discussion,” relating our work with other approaches and summarizing our results. Section 4 gives formal definitions appearing in the pseudo-code of the algorithms. Finally, the Supplementary Material includes all query files for the case studies, and a detailed explanation of the syntax used to formulate partially defined fixed-point constraints.

We demonstrate Griffin's functionalities with three different Boolean network models taken from the literature. The first example allows us to illustrate queries for (a) finding Boolean networks given an R-graph, (b) finding the sets of fixed-point attractors of the Boolean networks of a given R-graph, and (c) finding Boolean networks given an R-graph and a desired set of fixed-point attractors.

The second example is devoted to hypothetical regulations. Such regulations will enable us to exemplify Griffin's whole-query vs. query-splitting approaches.

The third example shows how uncertainty in the steady state behavior of a system can be effectively expressed by combining partially known fixed-point attractors and explicit exclusions of fixed-point attractors.

Some features, such as the use of constraints on mutants and specification of desired cyclic attractors, have been left out from this article. However, we have illustrated the use of these capabilities in previous articles (Rosenblueth et al., 2014; Weinstein et al., 2015; Azpeitia et al., 2017; García-Gómez et al., 2017).

In the inference of Boolean networks form regulation graphs, the fact that multiple networks might be consistent with the same (given) regulation graph is sometimes neglected. Often, publications reporting the inference of Boolean networks consider a single network. Even if several such networks are analyzed, only a fragment of the possibly vast number of consistent networks is usually considered. Griffin, by contrast, employs a formal definition of regulation (Naldi et al., 2007; Richard et al., 2012; Mori and Mochizuki, 2017) to build an implicit representation of all possible networks satisfying the regulation graph. Griffin then augments such a representation with biological constraints, resulting in a uniform representation of all potential networks and at the same time reducing such potential number of solution networks. This representation is then given to a powerful, symbolic search mechanism.

A key contribution to Model Checking (Clarke et al., 1999) was the incorporation of symbolic representations of data (Burch et al., 1992). Such representations allowed Model Checking to verify systems, first with 10²⁰ states, and subsequently with many more. As with model checkers, Griffin is able to handle large systems because of also using symbolic data structures. There exist two main symbolic techniques: binary-decision diagrams (BDDs) and SAT solvers. Griffin employs both, for solving different subproblems. Griffin finds cyclic attractors with a SAT solver (Dubrova and Teslenko, 2011), whereas it computes basins of attraction (Garg et al., 2008) with BDDs. The top level in Griffin, which builds a Boolean formula representing all biological information and constraints, also employs a SAT solver. Griffin thus capitalizes on the phenomenal recent developments of such solvers. This dramatic progress has “resulted in speed-ups of many orders of magnitude that have turned many problems that were considered intractable in the 1980s into trivially solved problems now.” (Franco and Martin, 2009).

Another strong point of our tool is the concept of an R-graph, a generalization of ordinary regulation graphs. R-graphs represent advantages from both the user and the computational viewpoints. The user has available the whole repertoire of combinations involving the sign of a regulation and whether a regulation is compulsory or simply hypothetical. For example, the user can specify a hypothetical regulation that, if present, should be negative and cannot be positive. From a computational point of view, an R-graph encompasses all hypothetical regulations in one formula, thus avoiding having to perform many separate tests. Computer tools employing ordinary regulation graphs, by contrast, usually test all possible combinations of presence and absence of all hypothetical regulations, resulting in numerous analyses.

Griffin is hence not only an innovative but also powerful computer tool for the inference of Boolean networks. It is possible, for example, to detect erroneous inferences due to the experimental, error-prone nature of biological data. Regulations that were believed to be positive might in fact turn out to be negative or vice versa. Griffin has detected these situations, as reported in Azpeitia et al. (2017).

As a first comparison between our work and related articles, it is important to point out a difference in the type of input data. In this work, Griffin input is composed of (partial) information about the network topology (R-graphs) along with other data representing biological constraints. R-graphs contain information on genetic connectivity that was inferred from data obtained by direct measurement of gene expression or protein interaction. By contrast, inference methods proposed by other authors (e.g., Laubenbacher and Stigler, 2004) have an input composed of time series together with other data that capture information on the network dynamics. Time series, unlike the R-graphs, represent experimental data obtained by direct measurement of gene expression or protein interaction.

There are a number of tutorials on Boolean-network inference (D'haeseleer et al., 2000; Markowetz and Spang, 2007; Karlebach and Shamir, 2008; Hecker et al., 2009; Hickman and Hodgman, 2009; Berestovsky and Nakhleh, 2013). From these tutorials we can classify algorithms according to (1) the expected input, (2) the kind of model inferred, and (3) the search strategy.

Much of the effort in Boolean-network inference has been aimed at having a binarized time-series data as input. Hence, multiple methods have been proposed and some of them have been compared with each other (Berestovsky and Nakhleh, 2013). An influential method in this category is reveal (Liang et al., 1998) (employing Shannon's mutual information between all pairs of molecules to extract an influence graph; the truth tables of the update functions for each molecule are simply taken from the time series). The use of mutual information and time-series for the inference of Boolean networks continues to develop (Barman and Kwon, 2017), as well as alternative methods based on time series. Example are Han et al. (2014) (using a Bayesian approximation), Shmulevich et al. (2003), Lähdesmäki et al. (2003), Akutsu et al. (1999), Laubenbacher and Stigler (2004), and Layek et al. (2011) (using a generate-and-test method, generating all possible update functions for one gene and testing with the input data). Extra information can be included in addition to time-series data. For example, an expected set of stable states (Layek et al., 2011), previously known regulations (Haider and Pal, 2012) and gene expression data (Chueh and Lu, 2012) are used as an aid to curtail the number of possible solutions.

Griffin belongs to a second family of methods taking as input a possibly partial regulation graph (perhaps obtained form the literature). There are also approaches employing both time-series data and a regulation graph, such as Ostrowski et al. (2016).

A third important area of research is the development of algorithms taking as input temporal-logic specifications, based on Model Checking (Clarke et al., 1999). Works following this approach are: Calzone et al. (2006b), Mateus et al. (2007), and Streck and Siebert (2015).

As for the kind of model inferred, here we would be concerned with Boolean networks and similar formalisms. Among the nonprobabilistic approaches, we find synchronous Boolean networks, asynchronous networks based on Thomas's formalism (Bernot et al., 2004; Khalis et al., 2009; Corblin et al., 2012; Richard et al., 2012), and polynomial dynamical systems (Laubenbacher and Stigler, 2004). Typically, methods based on temporal logic infer Kripke structures, which are closely related to Boolean networks. Probabilistic models, one the other hand, include Bayesian networks, and have the advantage of being able to deal with noise and uncertainty.

From the search-strategy point of view, Boolean-network inference methods may employ a simple random-value assignment (Pal et al., 2005), exhaustive search (Akutsu et al., 1999) or more elaborate algorithms. The work of Chueh and Lu (2012) is based on p-scores and that of Liang et al. (1998) guides search with Shannon's mutual information. Genetic algorithms are used by Saez-Rodriguez et al. (2009) and Ghaffarizadeh et al. (2017). Linear programming is the basis of Tarissan et al. (2008). The methods of Ostrowski et al. (2016) and Corblin et al. (2012), are based on Answer Set Programming (Brewka et al., 2011). Algebraic methods use reductions (of polynomials over a finite field) modulo an ideal of vanishing polynomials.

Approaches based on temporal logic (sometimes augmented with constraints), such as Calzone et al. (2006b) and Mateus et al. (2007), normally employ Model Checking (Clarke et al., 1999). Model Checking, in turn, is often based on symbolic approaches: BDDs and SAT solvers. Biocham (Calzone et al., 2006b) and SMBioNet (Bernot et al., 2004; Khalis et al., 2009; Richard et al., 2012) use a model checker as part of a generate-and-test method.

Having classified various approaches, we now mention a work similar to ours in spirit: Pal et al. (2005). These authors also propagate fixed-point constraints onto the truth tables of the update functions of each variable (molecule species). There is, however, no search technique, other than randomly giving values to the remaining entries of such truth tables. There is a random assignment of values to such entries, and a check for unwanted attractors. Neither is there a formal definition of regulation.

In contrast with the logical approach of our work, we devote our attention now to an algebraic approach to the problem of inference (reverse engineering) of Boolean networks. Instead of using a Boolean network to model the dynamics of a gene network, the algebraic approach (Laubenbacher and Stigler, 2004; Jarrah et al., 2007; Veliz-Cuba, 2012) uses a polynomial dynamical system F: Sⁿ → Sⁿ where S has the structure of a finite field (ℤ/p for an appropriate prime number p). The first benefit of this algebraic approach is that each component of F can be expressed by a polynomial (in n variables with coefficients in S) such that the degree of each variable is at most equal to the cardinality of the field S.

Following a computational algebra approach, in a framework of modeling with polynomial dynamical systems, Laubenbacher and Stigler (2004) propose a reverse-engineering algorithm. This algorithm takes as input a time series s1,…,sm∈Sn of network states (where S is a finite field), and produces as output a polynomial dynamical system F=(F1,…,Fn):Sn→Sn such that ∀i ∈ {1, …, n}: F_i ∈ S[x₁, …, x_n] and F_i(s_j) = s_j+1,i for j ∈ {1, …, m}. An advantage of the algebraic approach of this algorithm is that there exists a well-developed theory of algorithmic polynomial algebra, with a variety of procedures already implemented, supporting the implementation task. The time complexity of this algorithm is quadratic in the number of variables (n) and exponential in the number of time points (m).

Comparing Griffin with the reverse-engineering algebraic algorithms proposed by Laubenbacher and Stigler (2004) and Veliz-Cuba (2012), we found three basic differences. (1) Algebraic algorithms can handle discrete multi-valued variables, while Griffin only handles Boolean (two-valued) variables. Multi-valued variables give more flexibility and detail to the modeling process, but Boolean variables (of Boolean networks) lead to simpler models (see section 1). (2) The input of the algebraic algorithms, typically time series, provides simple information coming directly from experimental measurements, while input of Griffin, R-graphs and Griffin queries, provides structured information allowing a more precise specification of the required Boolean network. (3) The algebraic algorithm of Veliz-Cuba (2012) uses a formal definition of regulation, but this definition does not match the definition of regulation used by Griffin. While Griffin allows for R-regulations based on Boolean combinations of positive and negative regulations, Veliz-Cuba (2012) uses regulations restricted to unate functions h such that, for all variable x: h does not depend on x, or h depends positively on x, or h depends negatively on x.

Finally, we observe that sometimes results might have been reported overoptimistically. There have been some doubts cast upon the effectiveness of a number of methods of inference of network dynamics (Wimburly et al., 2003), especially those based on more general-purpose learning methods. It is therefore important to establish tests such as the dream challenges (Stolovitzky et al., 2007) emphasizing reproducibility.

Both Griffin and model checkers are based on symbolic search-algorithms, either SAT solvers or BDDs. Hence, a natural possibility to consider is the use of an off-the-shelf model checker as a search mechanism for Griffin. Model checkers take as input a temporal-logic formula. We did not find any advantage, however, in using a temporal logic, such as computation-tree logic (CTL) or linear-time logic (LTL), for expressing regulation over ordinary Boolean logic. Once a problem is expressed in Boolean logic, a SAT solver is precisely a method for finding a solution to such a problem. Hence, we disposed of temporal logic and directly expressed the regulation graph in Boolean logic.

Lacking a temporal logic makes the regulation graph the main means of communication of the user with Griffin. This suggested enriching the concept of a regulation graph as much as possible, which resulted in our generalization: the R-graph. In addition, the value of a number of parameters can be established with Griffin's query language. Nevertheless, as shown in our case studies, from a Biology point of view, it is easy to fall into situations where Griffin's query language is not expressive enough. In such situations, the user can employ Griffin's API, thus bypassing the query language. Moreover, Griffin's query language (apart from the R-graph, which can be viewed as a shorthand of a Boolean-logic formula) does not have a logical basis. Therefore, compared with a model checker, it would be desirable that Griffin have a temporal logic.

Temporal logic has been shown to be effective for the analysis of Boolean networks (Arellano et al., 2011; Carrillo et al., 2012; Klarner and Siebert, 2015). At the same time, our case studies have shown that clause learning could be effective in pruning the search space. This suggests combining Model Checking with clause learning. Much more complex queries than are currently possible could be formulated to Griffin in this manner. This is, therefore, an important avenue of research.

Another possibility for future research would explore the combination of Griffin with different approaches of Boolean network inference. If the search space is too large for Griffin to give useful solutions because of exhausting available resources, Griffin could incorporate algorithms based on other approaches, such as genetic algorithms, through its API.

Yet another direction of future work would extend Griffin to include either (a) biological information relating species with each other, such as protein complexes, thus reducing the size of the search space, or (b) partially defined Boolean functions for a particular gene.

Finally, we mention an important improvement that would allow the computation and the output of sets of Boolean networks instead of a single network at a time. Currently, the SAT solver employed by Griffin returns one assignment of all the variables appearing in the Boolean formula representing the R-graph together with the constraints (i.e., a “minterm”) at a time. By contrast, it would be more useful to have a SAT solver returning a set of assignments represented by a partial assignment (i.e., a term that is not necessarily a minterm). This would allow Griffin, in turn, to output sets of Boolean networks.

ℕ⁺ is the set of positive natural numbers. Unless differently stated, we assume that n ∈ ℕ⁺. ℕn+ is an initial segment of ℕ⁺, ℕn+={x∈ℕ+∣x≤n}. 𝔹 = {0, 1} is a set of Boolean values. If b ∈ 𝔹, then b′ is the complement of b. If x ∈ 𝔹ⁿ, we say that x is a state and x_i denotes the i-th component of x. We sometimes omit parentheses and commas when writing a state. If f:X→X, and k ∈ ℕ⁺, we write f^k to denote the iterated composition of f with itself, i.e., f¹ = f and fⁿ⁺¹ = f ◦ fⁿ. If v ∈ 𝔹^k, p∈(ℕn+)k, and p_i ≠ p_j for i ≠ j, then x[v/p] is the state resulting from replacing, for all i∈ℕk+, the p_i-th component of x by v_i. We denote as x ~ i the vector resulting from replacing the i-th component of x with the complement of x_i, that is, x~i=x[xi′/i]. Given k < n, v ∈ 𝔹^k, and p∈(ℕn+)k, 𝔹ⁿ[v/p] is a subspace of 𝔹ⁿ defined by 𝔹ⁿ[v/p] = {x[v/p]∣x ∈ 𝔹ⁿ}. Alternatively, we describe subspaces by means of strings of zeros, ones, and asterisks used as wildcards indicating the free (unknown value) components. For example, 𝔹⁶[(0, 1, 0)/(1, 4, 6)] can be written simply as 0^**1^*0.

Definition 1. We define a synchronous Boolean network with n components as a function f:𝔹ⁿ→𝔹ⁿ. The i-th component of f is a function fi:𝔹n→𝔹 such that f_i(x) = f(x)_i. We use BN to denote the class of Boolean networks with n components.

To relate a synchronous Boolean network with a molecular network dynamics, we interpret each component of a state x as representing the activation state of a particular variable denoting a molecule included in the network (which can be a gene, a protein, a hormone, for example). Given a set of variables V, we use a bijective function v:ℕn+→V to relate the components of a state x with their respective variables. The molecule represented by the variable v(i), also denoted v_i, is said to be active if x_i = 1 and inactive otherwise.

Definition 2. If p∈(ℕn+)k, and v ∈ 𝔹^k with k∈ℕn+, the mutation of f by p = v is defined as f^{p = v}: 𝔹ⁿ → 𝔹ⁿ, where, for all x ∈ 𝔹ⁿ and j∈ℕk+, fpjp=v(x)=vj, and fip=v(x)=fi(x) if i∈ℕn+-{p1,…,pk}. We say that f^{p = v} is a single-point mutation of f if k = 1, and a multi-point mutation if k > 1.

Definition 3. (Naldi et al., 2007; Richard et al., 2012; Mori and Mochizuki, 2017). A positive regulation between variables v_i and v_j is a function R⁺: V × BN × V → 𝔹 such that R+(vi,f,vj)={1 if ∃x∈𝔹n:fj(x)≠fj(x~i) and xi=fj(x)0 otherwise

Similarly, a negative regulation between v_i and v_j is a function R⁻: V × BN × V → 𝔹 such that R−(vi,f,vj)={1  if ∃x∈𝔹n:fj(x)≠fj(x~i) and xi≠fj(x)0  otherwise

Definition 4. We say that G is a regulation graph if G = (V, I⁺, I⁻), where: V is the set of vertices, I⁺ ⊆ V × V is the set of positive regulations, and I⁻ ⊆ V × V is a set of negative regulations. If i ∈ I⁺ ∩ I⁻, we say that i is an ambiguous regulation.

Definition 5. The regulation graph of f is Gf=(V,I+,I-) where:

Intuitively, the regulation graph of f, G_f, describes the structure of f and its interpretation as a molecular network where edges represent molecular regulations. Note that G_f may have both a positive and a negative regulation from j to i. Observe that whereas a Boolean network f has a unique regulation graph, a regulation graph may have more than one Boolean network.

Definition 6. Given a regulation graph G, we say that f is a satisfying Boolean network of G, denoted by f ⊨ G, if G = G_f. The set of all satisfying Boolean networks of G is denoted by F_G = {f∣f ⊨ G}.

Definition 7. The state graph of a Boolean network f is the graph Ĝf=(𝔹n,{(x,f(x))∣x∈𝔹n}).

Definition 8. We say that ω ⊆ 𝔹ⁿ is an attractor of f if ω is a terminal strongly connected component of Ĝ_f (Ruet, 2017). That is, ω is an attractor if ω is strongly connected (each state in ω is reachable from any other state in ω) and no state or edge of Ĝ_f can be added to ω without causing ω to be no longer strongly connected. In particular, ω is a fixed point (or stationary state) of f if |ω| = 1, it is a cyclic attractor otherwise. The size k of an attractor is defined as its cardinality. The set of attractors of f is denoted as A(f). The set of fixed points of f is denoted as FP(f). The set of cyclic attractors of f is denoted as CA(f). If ω ∈ A(f), the basin of attraction of ω is BA(ω)={x∣x∈𝔹nand∃k∈ℕ2n-1+:fk(x)∈ω}.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.