Computing preimages of Boolean networks

Abstract

In this paper we present an algorithm based on the sum-product algorithm that finds elements in the preimage of a feed-forward Boolean networks given an output of the network. Our probabilistic method runs in linear time with respect to the number of nodes in the network. We evaluate our algorithm for randomly constructed Boolean networks and a regulatory network of Escherichia coli and found that it gives a valid solution in most cases.

Introduction

In systems and computational biology Boolean networks (BN) are widely used to model regulative dependencies of organisms [1, 2]. We consider networks, which map a set of environmental conditions to the presence of proteins and finally to actual chemical reactions, which are often modeled as fluxes of a flux-balance analysis [3]. Hence, these networks are used to make in silico predictions of behavior of organisms in a certain environment [4].

In this paper we address the inverse problem, i.e., we want to predict environmental conditions that allow certain reactions to take place, and others not. Hence, in general, we need to find a set of possible inputs that lead to a given output. This so called predecessor problem or preimage problem has been addressed by Wuensche in [5] and has been shown to NP-hard in general [6], which makes it infeasible to solve it for large networks. In [7] an algorithm with reduced complexity for BNs with canalizing Boolean functions has been introduced. However, the problem is infeasible under certain conditions. Both algorithms are designed to find the whole set of preimages, i.e., all inputs to the BN with lead to a certain, desired, output.

In some applications, knowledge of the whole preimage set is not important, merely it can be sufficient to know a subset of the preimage set. Here, we propose a probabilistic algorithm, which solves this problem in linear time with respect to the number of nodes in the network, based on a variation of the well known Sum-Product Algorithm (SPA) [8], which is used for a variety of tasks, including decoding error correction codes in communication engineering [9].

Methods

Boolean networks and main idea

We consider networks like shown in Figure 1, mapping the values of the N in-nodes $\mathbb{I}=\left\{1,2,3\right\}$ to the M out-nodes $\mathbb{O}=\left\{12,13,14,15,16\right\}$, i.e., we can represent this BN as a function mapping the N input values uniquely to the M output values:

Figure 1

Example of a feed-forward network

$$f:{\left\{0,1\right\}}^N\to {\left\{0,1\right\}}^M.$$

The network itself consists of n nodes, and a set of directed edges connecting these nodes. Each node i has a certain state, which can be either zero or one, represented by a variable x _i . Its value is determined by evaluating a Boolean function (BF) f _i . Further, lets define the set ñ(f _j ) as the incoming nodes of node j. For example in Figure 1, ñ(f₅) = {1, 3}. The BF f _j is a function mapping k _j = |ñ(f _j )| values of {0, 1}^k to {0, 1}, where k is also called the in-degree of node j. The number of edges emerging from a node is called out-degree.

Given a vector of input values x ∈ {0, 1}^N , x = (x₁, x₂, . . . , x _N ) the corresponding output of f is y = f (x), y ∈ {0, 1}^M. In general there does not exist a unique inverse function f⁻¹. Instead the cardinality of the set Ω _y := {x : f(x) = y} will be larger one. We call Ω _y the set of preimages of y. In this paper we are interested to find at least parts of Ω _y . Suppose there is a probability distribution P _y on {0, 1}^N such that

$$P_{\mathsf{\text{y}}}\left\{\mathsf{\text{x}}\right\}=\left\{\begin{array}{cc}\hfill \frac{1}{\left|{\text{Ω}}_y\right|}\hfill & \hfill \mathsf{\text{if x}}\in {\text{Ω}}_y\hfill \\ 0\hfill & \hfill \mathsf{\text{else}}\hfill \end{array}\right..$$

If we knew the probability distribution P _y , we would have solved the problem. But as explained, this is too difficult in general. Our main idea now is to approximate P _y by the product of the marginal distributions P _i on the individual x _i , i.e.,

$$P_{\mathsf{\text{y}}}\approx {\displaystyle \prod _{i=1}^N}{P}_i,$$

as the well-known SPA can be used to compute the marginals efficiently. If the approximation is good enough sampling out the product of the marginals will yield an element in Ω _y with reasonable probability.

Proposed algorithm

In this section we will first discuss the basic principles of factor graphs and the SPA. Then we will describe the BN as factor graph and will formulate the actual algorithm to find the marginals. Finally, the sampling is described.

Factor graphs and sum-product algorithm

Assume some function g(x₁, . . . , x _n ) defined on some domain ${\mathbb{A}}^n$, which can be factorized in m local functions h _j , j ∈ [m] := {1, 2, . . . , m}, i.e.,

$$g\left(x_1,\dots ,x_n\right)={\displaystyle \prod _j}{h}_j\left(X_j\right),$$

where X _j is the subset of [n] containing the argument of h _j . We can then define a factor graph [8] as a bipartite graph consisting of n nodes representing variables {x₁, . . . , x _n } (variable nodes) and of m nodes representing functions {h₁, . . . h _m } (function node). Edges exist between a function node and a variable node if and only if x _i is an input to function h _j .

The marginal function g _i (x _i ) is defined as [8]

$$g_i\left(x_i\right)={\displaystyle \sum _{~\left\{x_j\right\}}}g\left(x_1,\dots ,x_n\right),$$

where ${\sum }_{~\left\{x_j\right\}}g\left(x_1,\dots ,x_n\right)$ is defined as

$$\begin{array}{l}{\displaystyle \sum _{~\left\{x_i\right\}}}g\left(x_1,\dots ,x_n\right)\\ ={\displaystyle \sum _{x_1\in A}}\dots {\displaystyle \sum _{x_{i-1}\in A}}{\displaystyle \sum _{x_{i+1}\in A}}\dots {\displaystyle \sum _{x_n\in A}}g\left(x_1,\dots ,x_n\right),\end{array}$$

In general the computation of the g _i is difficult, but due to the factorization of g the task can be efficiently solved using the the so called Sum-Product Algorithm (SPA) [8]. The algorithm iteratively passes messages between the nodes of the graph. At each iteration the messages µ are sent from the function nodes to the variable nodes, containing the corresponding marginal function of the local function. These messages are computed as follows [8]:

Function to variable node

$${\mu }_{h\to x}\left(x\right)={\displaystyle \sum _{~\left\{x\right\}}}\left(h\left(n\left(h\right)\right){\displaystyle \prod _{y\in n\left(h\right)\backslash \left\{x\right\}}}{\lambda }_{y\to h}\left(y\right)\right),$$

where n(i) give the set of neighboring nodes of node i.

At the variable nodes, these messages are then combined to a marginal function λ and sent back to the function nodes [8]:

Variable to function node

$${\lambda }_{x\to h}\left(x\right)={\displaystyle \prod _{q\in n\left(x\right)\backslash \left\{h\right\}}}{\mu }_{q\to x}\left(x\right).$$

The Boolean network as factor graph

We apply the concept of factor graphs to BNs. Each node in the network represents one variable x _i ∈ {0, 1}, i ∈ [n] of the factor graph, hence we have n variable nodes. Each BF f _j of the BN $\left(j\in \left[n\right]\backslash \mathbb{I}\right)$ is a function node and is connected to the node j and the incoming nodes ñ(f _j ). Lets to define ${\tilde{X}}_j$ as the variables of the incoming nodes of node j, i.e. the argument of the BN f _j . Further, we define ${\tilde{X}}_j^{\left(i\right)}$ as ${\tilde{X}}_j$ without the node i.

Finally, if we consider the variables as each node as random variables, we have a common distribution of all variables nodes described by the density function,

$$g_{x_1,\dots ,x_n}\left(x_1,\dots ,x_n\right)\equiv g\left(x_1,\dots ,x_n\right),$$

For sake of readability we will omit the subscripts of the density function, if they are obvious from context. We are interested in finding the marginal distributions of the in-nodes, which can be described by the density functions

$$g_{x_i}\left(x_i\right)={\displaystyle \sum _{~x_i}}{g}_{x_1,...,x_n}\left(x_1,\dots ,x_n\right)\forall i\in I.$$

This problem is an instance of the problem described in Section Factor Graphs and Sum-Product Algorithm, hence we apply the same methods here.

Update rule: function to variable node

If we focus on one function node $j\in \left[n\right]\backslash I$ there exists a common distribution of all variables relevant for this node. Namely, these relevant variables are the ones located in ${\tilde{X}}_j$ of the BF f _j , and the value of node j. We can write the density of this distribution as:

$$p\left(x_j,{\tilde{X}}_j\right).$$

Lets define ñ(f _j ) as the set of indices of the input nodes of the BF f _j .

We need to send the local marginal distribution of each variable i ∈ {j} ∪ ñ(f _j ) back to the variable node, or more formally:

$${\mu }_{j\to i}\left(x_i\right)={\displaystyle \sum _{~\left\{x_i\right\}}}p\left(x_j,{\tilde{X}}_j\right)={\displaystyle \sum _{~\left\{x_i\right\}}}p\left(x_j,x_i,{\tilde{X}}_j^{\left(i\right)}\right)$$

(1)

If i = j, i.e. if the message is designated for the node containing the output of the BF, the density of the marginal distribution becomes:

$$\begin{array}{c}{\mu }_{j\to j}\left(x_j\right)={\displaystyle \sum _{~\left\{x_j\right\}}}p\left(x_j|{\tilde{X}}_j\right)\cdot \left({\tilde{X}}_j\right)\\ ={\displaystyle \sum _{~\left\{x_j\right\}}}{f}_j\left({\tilde{X}}_j\right)\cdot \left({\tilde{X}}_j\right)\end{array}$$

which is the probability distribution of the functions output. We can assume that the elements of ${\tilde{X}}_j$ are pairwise independent, hence we can write:

$$p\left({\tilde{X}}_j\right)={\displaystyle \prod _{l\in ñ\left(f_j\right)}}{\lambda }_l\left(x_l\right),$$

where λ _l is the probability distribution of variable node l and is defined in Eq. 3.

In the other cases, i.e., i ≠ j, Eq. (1) becomes:

$${\mu }_j\to i\left(x_i\right)={\displaystyle \sum _{~\left\{x_i\right\}}}p\left(x_i|x_j,{\tilde{X}}_j^{\left(i\right)}\right)\cdot \left(x_j,{\tilde{X}}_j^{\left(i\right)}\right).$$

We still can assume that the elements of ${\tilde{X}}_j^{\left(i\right)}$ are pairwise independent, hence we can write:

$$\begin{array}{ll}\hfill p\left(x_j,n\left(f_j\right)\backslash x_i\right)& =p\left(x_j|{\tilde{X}}_j^{\left(i\right)}\right)\cdot \left({\tilde{X}}_j^{\left(i\right)}\right)\\ =p\left(x_j|{\tilde{X}}_j^{\left(i\right)}\right){\displaystyle \prod _{l\in ñ\left(f_j\right)\backslash \left\{i\right\}}}{\lambda }_l\left(x_l\right).\end{array}$$

If the Boolean functions output $x_j=f_j\left({\tilde{X}}_j\right)$ is already completely determined by ${\tilde{X}}_j^{\left(i\right)}$, i.e., if the variable x _i has no influence on the output for this particular choice of the other variables, we assume x _i to be uniformly distributed:

$$p\left(x_i|x_j,{\tilde{X}}_j^{\left(i\right)}\right)=\frac{1}{2}{p}_{x_j}\left(f\left({\tilde{X}}_j^{\left(i\right)},x_i\right)=x_j\right)$$

and since x _j is completely determined by ${\tilde{X}}_j^{\left(i\right)}$

$$p(x_j,{\tilde{X}}_j^{(i)})={\displaystyle \prod _{l\in \tilde{n}(f_j)/\left\{i\right\}}{\lambda }_l(x_l)}.$$

Otherwise, x _i is totally determined by x _j and the other variables, i.e., x _i is 0 or 1 depending on BF. Hence, we can write

$$p\left(x_i|x_j,n\left(f_j\right)\backslash x_i\right)=p_{x_j}\left(f\left({\tilde{X}}_j^{\left(i\right)},x_i\right)=x_j\right),$$

where $p_{x_j}\left(f\left({\tilde{X}}_j^{\left(i\right)},x_i\right)=x_j\right)$ is either 0 or 1. Further we can assume x _j independent of ${\tilde{X}}_j^{\left(i\right)}$, hence

$$p\left(x_j,{\tilde{X}}_j^{\left(i\right)}\right)={\lambda }_j\left(x_j\right){\displaystyle \prod _{l\in ñ\left(f_j\right)\backslash \left\{i\right\}}}{\lambda }_l\left(x_l\right).$$

Finally, we can summarize for i ≠ j:

$${\mu }_{j\to i}\left(x_i\right)={\displaystyle \sum _{~\left\{x_i\right\}}}{\xi }_{i,j}{p}_{x_j}\left(f\left({\tilde{X}}_j^{\left(i\right)},x_i\right)=x_j\right){\displaystyle \prod _{l\in ñ\left(f_j\right)\backslash \left\{i\right\}}}{\lambda }_l\left(x_l\right),$$

(2)

with

$${\xi }_{i,j}=\left\{\begin{array}{cc}\frac{1}{2}\hfill & \hfill ,\mathsf{\text{if }}{f}_j\left({\tilde{X}}_j^{\left(i\right)},x_i=0\right)=f_j\left({\tilde{X}}_j^{\left(i\right)},x_i=1\right)\hfill \\ \hfill {\lambda }_j\left(x_j\right)\hfill & \hfill ,\mathsf{\text{else}}\hfill \end{array}\right..$$

Update rule: variable to function node

The update rule is the same for all variable nodes j ∈ [n] and is independent of the function node to which they are directed.

$${\lambda }_j\left(x_j\right)={\displaystyle \prod _{l\in {\mathbb{S}}_j}}{\mu }_{l\to j}\left(x_j\right),$$

where $S_j$ is the set of all function nodes, which have node j as input.

Finding the input distributions

In our algorithm, we use the well known log-likelihood ratio (LLR) to represent the probability distribution of binary variables [10]. It is defined as:

$$L_X=\text{ln}\frac{p\left(x=0\right)}{p\left(x=1\right)}.$$

(4)

A scheme of the algorithm is given in Algorithm 1.

The probability distribution of each node j ∈ [n] at iteration t is given as $L_j^{\left(i\right)}$ and are initialized with $L_j^{\left(0\right)}=0$, which is equivalent to the uniform distribution. Then we set the LLRs for the out-nodes to either −∞ or +∞ depending on the desired output y of the BN. At each iteration the algorithm can be split in two steps. The first step iterates over all function nodes $\left(j\in \left[n\right]\backslash \mathbb{I}\right)$ and all input variables i ∈ ñ(f _j ) calculating the LLR $L_{j\to i}^{\left(t\right)}$ using Eq. (2) and Eq. (4).

In the second step we update all variables-nodes, where the LLRs L _j represents the distributions λ _j and, hence the product of Eq. 3 becomes a summation. Please note, that the LLR of the previous iteration is also added to the sum, in order to prevent rapid changes of the distributions.

After performing a certain number of iterations t _max , the desired marginal distributions of the input variables are found.

Algorithm 1

    Initialize $L_j^{\left(0\right)}=0$ for all nodes

    Set the desired LLRs of the out-nodes, i.e., $L_j^{\left(0\right)}$ is either −∞ or +∞, for all out-nodes $j\in O$.

    t = 0

    repeat

      t=t+1

      for each non-in-node $\left(j\in \left[n\right]\backslash I\right)$do

        for each input variable i ∈ ñ(f _j ) do

          calculate L _j→i using Eq. (2) and Eq. (4)

        end for

      end for

      for each non-out-node v do

$L_j^{\left(t\right)}=L_j^{\left(t-1\right)}+{\sum }_{l\in S_j}{L}_{l\to j}^{\left(t\right)}$

      end for

    until maximum number of iterations reached

Sampling

The sampling part of our approach is straightforward. Using the marginal distributions $L_j^{\left(t_{max}\right)}$, $j\in \mathbb{I}$ we randomly draw vectors x and check if they fulfill y = f(x). If so, they are added to the set ${\widetilde{\text{Ω}}}_y$. This procedure is repeated for a certain number of samples.

Simulation results and discussion

We tested our algorithm with randomly generated networks and the regulatory network of Escherichia coli (E-coli) [2]. The random networks consist of 2400 nodes with N = 200 and M = 1200. We have chosen the BFs from:

· all functions with k ≤ 15 (Type A)

unate, i.e. locally monotone, functions with k ≤ 15 (Type B)

After generating a network we draw a certain number T of uniformly distributed input vectors x and obtain y = f(x). For each y we applied then Algorithm 1 to obtain the marginal distributions$L_j^{\left(t_{max}\right)}$,$j\in \mathbb{I}$ . To investigate the convergence behavior with respect to t _max we first apply hard-decision to evaluate a good choice for t _max , i.e., we generate an estimate $\tilde{x}$ by setting

$${\tilde{x}}_j=\left\{\begin{array}{cc}\hfill 0\hfill & \hfill \mathsf{\text{if}}{L}_j^{\left(t_{max}\right)}>0\hfill \\ \hfill 1\hfill & \hfill \mathsf{\text{if}}{L}_j^{\left(t_{max}\right)}<0\hfill \end{array}\right.$$

Then we evaluate the network $ỹ=\text{f}\left(\tilde{x}\right)$, and measure the similarity between y and $ỹ$ by counting the equal entries and divide them by the length of y. We did so for 100 networks of Type A and B, and set T = 100. The averaged results can be seen in Figure 2.

Figure 2

Similarity of y and $ỹ$ versus t _max

One can see, that for t _max ≥ 14 there is almost no improvement in the similarity. This number is equal to two times the number of nodes between input and output, i.e., it seems to be sufficient that the messages travel once through the network and back. Thus, the following simulations have been perform setting t _max = 14.

Next, we apply sampling as described in Section Sampling. We did so for 100 different networks of Type A and B, and the E-coli network. For each random network we did T = 100 runs, for E-coli T = 1000. The results can be viewed in Table 1. We depict the percentage of solved networks, i.e. the portion of networks we found at least one valid x ∈ Ω _y . Further, we give the average number of valid x and the average number of unique x.

Table 1 Simulation results for different networks

One can see from the results, that in general for most networks and y s at least one preimage can be found. It is worth mentioning, that for the E-coli network every sampled solution was unique. This is due to the fact, that there exist a few inputs, who completely determine the output. The other input variables have then no influence and hence a marginal distribution of 0.5. Further, the results for the network of type B are much better than for type A. It seems that the marginal distributions for unate functions give better estimation of the actual distribution than the marginal distributions for non-unate functions.

Conclusions

In this work, we proposed a probabilistic algorithm to address the preimage problem of Boolean networks. This is of interest when designing experiments, in which certain regulators are supposed to be in a specific state. Performing a series of simulations with Random networks we showed, that the algorithm works not only for unate functions, of which most biologically motivated networks consist, but for any kind of Boolean functions. By replacing the fixed output values of the network by probabilities one can simply apply the algorithm to networks, whose designated output is described by probability distributions. Further, the algorithm may be easily adjusted to work on stochastic, e.g. Bayesian, networks, where the nodes contain only transition probabilities instead of Boolean function. Therefore, it is needed to adapt the update rules accordingly. It remains an open question, which influence topographic properties, such as number of layers and number of nodes in these layers, have to the performance of the proposed algorithms, since we only investigated networks which are similar to the regulatory network of E-coli.