Skip to main content

Advertisement

Identification of control targets in Boolean molecular network models via computational algebra

BMC Systems Biology volume 10, Article number: 94 (2016) | Download Citation

Article metrics

Abstract

Background

Many problems in biomedicine and other areas of the life sciences can be characterized as control problems, with the goal of finding strategies to change a disease or otherwise undesirable state of a biological system into another, more desirable, state through an intervention, such as a drug or other therapeutic treatment. The identification of such strategies is typically based on a mathematical model of the process to be altered through targeted control inputs. This paper focuses on processes at the molecular level that determine the state of an individual cell, involving signaling or gene regulation. The mathematical model type considered is that of Boolean networks. The potential control targets can be represented by a set of nodes and edges that can be manipulated to produce a desired effect on the system.

Results

This paper presents a method for the identification of potential intervention targets in Boolean molecular network models using algebraic techniques. The approach exploits an algebraic representation of Boolean networks to encode the control candidates in the network wiring diagram as the solutions of a system of polynomials equations, and then uses computational algebra techniques to find such controllers. The control methods in this paper are validated through the identification of combinatorial interventions in the signaling pathways of previously reported control targets in two well studied systems, a p53-mdm2 network and a blood T cell lymphocyte granular leukemia survival signaling network. Supplementary data is available online and our code in Macaulay2 and Matlab are available via http://www.ms.uky.edu/~dmu228/ControlAlg.

Conclusions

This paper presents a novel method for the identification of intervention targets in Boolean network models. The results in this paper show that the proposed methods are useful and efficient for moderately large networks.

Background

The dynamics of gene regulatory networks (GRNs) have been studied using different modeling frameworks, with the goal of building computational models of GRNs to get new insights into important cellular processes, e.g., the cell cycle, [1, 2], oscillations in the p53-mdm2 system, [3–5], the phage-lambda system, [6–8], or the T cell large granular lymphocyte (T-LGL) leukemia network, [9, 10]. A generally difficult problem is to design control policies to achieve desired dynamics in GRNs. This is particularly important in the control of cancer cells, [5, 11–14] and cell fate reprogramming, [15, 16]. Thus, developing tools to control mathematical models of GRNs are key to the design of experimental control policies.

There is a rich theory for the control of continuous models such as systems of differential equations, [17–20]. Discrete models such as Boolean networks (BN) have been proposed to study GRNs. In a BN, the genes of a GRN are represented by a set of nodes that can take on only two possible states (ON or OFF). Time is discrete, and the state of a gene at the next time step is determined by a Boolean function over a subset of nodes of the BN. The dependence of a gene on the state of another gene is graphically represented by a directed edge. BN models are widely used in molecular and systems biology to capture coarse-grained dynamics of a variety of regulatory networks and have been shown to provide a good approximation of the dynamics of continuous processes [8, 10, 21–28]. However, control methods for discrete models are still in their infancy, compared to the theory for continuous models.

In this work, we propose a framework to study the control of BNs. Therapeutic interventions are modeled by manipulating the wiring diagram of a BN accordingly. For example, a gene knockout is modeled by fixing the value of one node of the BN to OFF. Similarly, deleting directed edges of a BN models the action of a drug that inactivates the interaction between two gene products. The control problem consists of finding the sequence of actions (node and edge deletions) to get the BN out of an undesirable state, and drive it towards a desirable state. Undesirable states may represent a disease condition of a cell such as, for example, a highly proliferative state of a cancer cell, and a desirable state may represent cell death. Thus, a therapeutic intervention could aim at inducing the death of aberrant tumor cells.

Several approaches to address this problem have been used in recent years. For example, the optimal control techniques developed in [29–32] assume a set of control nodes to derive a control policy that minimizes the likelihood of transitioning into an undesirable state in a stochastic context. Other approaches for the identification of intervention targets include the use of stable motifs in the network to induce the system into a desirable state [33]. In [34], integer programming was used to find a set of nodes to control the states of BNs representing tumor and normal cells. Optimization techniques were used in [35] to determine possible node manipulations for BNs. There are also search algorithms based on genetic and greedy algorithms described in [36–38]. For continuous models, a related control approach is given in [19].

The idea behind our approach is to encode the problem of finding control candidates as a problem of solving a system of polynomial equations. Then we can use computational algebra techniques to solve the system. This approach takes advantage of the rich algorithmic theory of computer algebra (e.g. Gröbner basis techniques) and the theoretical foundations of algebraic geometry with software available for computations (e.g., Macaulay2 [39]). The output of our method is a set of candidate actions with the capacity to induce the GRN towards desirable states. The method also has the ability of identifying combinatorial interventions in the network which are sometimes more effective as will be shown in the “Results” section. The algebraic view of discrete models has been previously used for the development of computational tools to determine the attractors of BNs [40–43], and also for inferring network structure from dynamics [44–47].

Methods

Boolean networks

A Boolean network is a dynamical system that is discrete in time as well as in variable states. More formally, consider a collection x 1,…,x n of variables that take values in the binary set {0,1}. Then a Boolean network in the variables x 1,…,x n is a function

$$ \mathbf{F} = (f_{1},\dots,f_{n}):\{0,1\}^{n}\rightarrow \{0,1\}^{n} $$

where each coordinate function f i :{0,1}n→{0,1} is a Boolean function on a subset of {x 1,…,x n } which represents how the future value of the i-th variable depends on the present values of the other variables.

Example 2.1

For concreteness, we illustrate the definitions using the following toy network

$$ \begin{array}{lll} f_{1} = \neg x_{3} \wedge \neg x_{5}, \qquad \qquad & f_{2} = \neg x_{1} \vee x_{4}, \qquad \qquad & f_{3} = \neg x_{2}, \\ f_{4} = x_{3}, & f_{5} = \neg x_{4}, \end{array} $$

where ∧,∨, and ¬ are the AND, OR, and NOT logical operators, respectively. In the context of modeling biological systems, ∧ corresponds to activation by the combination of regulators (all regulators are necessary), ∨ corresponds to independent activation (one regulator is sufficient), and ¬ corresponds to negative regulation.

Given a Boolean network F=(f 1,…,f n ), a directed graph \(\mathcal {W}\) with n nodes x 1,…,x n is associated with F. There is a directed edge in \(\mathcal {W}\) from x j to x i if x j appears in f i . Notice that the presence of the interaction x j →x i implies that the regulatory function f i depends on x j , say \(f_{i}(x_{k_{1}},\dots,x_{j},\dots,x_{k_{m}})\) with \(x_{j}\in \{x_{k_{1}},\dots,x_{k_{m}}\}\). In the context of a molecular network model, \(\mathcal {W}\) represents the wiring diagram of the network.

Example 2.1 (cont.) The wiring diagram of the toy network is shown in Fig. 1 a.

Fig. 1
figure1

Description of the algebraic approach of identification of control targets. a. A BN model of a molecular network. The control variables are the entries of u. b. In the absence of control policies (u=0), the attractor landscape (L F ) can have undesirable attractors. c. The goal is to choose control values that give a desired attractor landscape, L ∗. d. To use the algebraic approach we first find the polynomial representation of the BN (see Section Control Actions: edge and node manipulations). e. The next step is to set up the desired attractor landscape as a system of equations that the BN has to satisfy, \(L_{\mathcal {F}(x,u)}=L^{*}\) (see Section Control targets in Boolean networksControl targets in Booleannetworks). f. Solving the equation \(L_{\mathcal {F}(x,u)}=L^{*}\) for u will provide the control values to achieve the desired landscape (see Section Identifying control targets). This approach not only finds individual control policies (u=A, single node; u=B, single edge), but also combinatorial control policies (u=C, two nodes; u=D, two edges). In a combinatorial control policy, the desired attractor landscape is achieved by the combination of two or more entries of u

Full size image

If the set {0,1} is given the structure of a finite field with standard addition and multiplication, \(\mathbb F_{2} = \{0,1\}\), then the functions \(f_{i}:\mathbb F_{2}^{n}\rightarrow \mathbb F_{2}\) can be represented as polynomials over \(\mathbb F_{2}\), and the dynamical system \(\mathbf {F} = (f_{1},\dots,f_{n}):\mathbb F_{2}^{n}\rightarrow \mathbb F_{2}^{n}\) becomes a polynomial dynamical system, see [41], which gives access to a range of mathematical tools for the analysis of F.

Example 2.1 (cont.) The rules to transform Boolean functions into polynomials in \(\mathbb {F}_{2}[x_{1},\ldots,x_{n}]\) are given by a∧b=a b,a∨b=a+b+a b, and ¬a=1+a, where the operations are computed modulo 2. For the toy network we obtain

$$ \begin{array}{l} f_{1} = (1+ x_{3})(1+x_{5})=1+x_{3}+x_{5}+x_{3}x_{5},\\ f_{2} = (1+ x_{1})+ x_{4} +(1+ x_{1})x_{4} = 1+x_{1}+x_{1}x_{4},\\ f_{3} = 1+ x_{2}, \\ f_{4} = x_{3},\\ f_{5} = 1+ x_{4}. \end{array} $$

The dynamical properties of a Boolean network are given by the difference equation x(t+1)=F(x(t)); that is, the dynamics is generated by iteration of F. More precisely, the dynamics of F is given by the state space graph S, defined as the graph with vertices in \(\mathbb F_{2}^{n}=\{0,1\}^{n}\), which has an edge from \(x\in \mathbb F_{2}^{n}\) to \(y\in \mathbb F_{2}^{n}\) if and only if y=F(x). The states \(x\in \mathbb F_{2}^{n}\) where the system will get stabilized are of particular importance. These special points of the state space are called attractors of a Boolean network and these may include steady states (fixed points), where F(x)=x, and cycles, where F k(x)=x for some integer k>1. Attractors in Boolean network modeling might represent cell types [48] or cellular states such as apoptosis, proliferation, or cell senescence [49, 50]. Identifying the attractors of a system is an important step towards the control of that system and can be done using tools from computational algebra [40, 41].

Example 2.1 (cont.) The steady states of the toy network are found by solving the system of equations f i =x i ,i=1,…,5. This means we want to find the roots of g i =0, where g i =f i −x i . This gives the system of equations

$$ \begin{array}{l} g_{1} = 1+x_{3}+x_{5}+x_{3}x_{5} + x_{1}=0,\\ g_{2} = 1+x_{1}+x_{1}x_{4} + x_{2}=0,\\ g_{3} = 1+ x_{2} + x_{3}=0, \\ g_{4} = x_{3} + x_{4}=0,\\ g_{5} = 1+ x_{4} + x_{5}=0. \end{array} $$
(1)

Note that, since the system is not linear, we cannot make use of tools such as Gaussian elimination. Computational algebra allows us to solve such systems by encoding the solutions as an algebraic object called an ideal of polynomials, I=〈g 1,g 2,g 3,g 4,g 5〉, and then finding an equivalent, but simpler representation. More precisely, for the system in Eq. 1 we can find its Gröbner basis [51] and obtain

$$ I = \langle \, g_{1},g_{2},g_{3},g_{4},g_{5}\, \rangle = \langle x_{5}+1, x_{4}, x_{3}, x_{2}+1, x_{1} \rangle. $$

This means that the original system has the same solutions as the simpler system

$$ \begin{array}{l} x_{1} = 0, x_{2}+1 =0, x_{3} = 0, x_{4} =0, x_{5}+1 =0, \end{array} $$

which is easily solved to obtain x=01001 as a steady state of the toy network (parentheses omitted for brevity).

Control Actions: edge and node manipulations

This paper considers two types of control actions: 1. Deletion of edges and 2. Deletion (or constant expression) of nodes. The motivation for considering these actions is that they represent the common interventions that prevent a regulation from happening. For instance, an edge deletion can be achieved by the use of therapeutic drugs that target specific gene interactions and node deletions represent the blocking of effects of products of genes associated to these nodes; see [5].

A schematic of our approach is given in Fig. 1 where the dynamics of a Boolean network can be manipulated by a set of controls consisting of deletion or constant expressions of edges and nodes. Below we define and explain all the steps in more detail.

A Boolean network with control is given by a function \(\mathcal {F}:\mathbb {F}_{2}^{n}\times U\rightarrow \mathbb {F}_{2}^{n}\), where U is a set that denotes all possible control inputs (Fig. 1 a). Given a control u in U, the dynamics are given by \(x(t+1)=\mathcal {F}(x(t),u)\).

We consider a BN \(\mathbf {F}=(f_{1},\ldots,f_{n}):\mathbb {F}_{2}^{n}\rightarrow \mathbb {F}_{2}^{n}\) and show how to encode edge and node controls by \(\mathcal {F}:\mathbb {F}_{2}^{n}\times U \rightarrow \mathbb {F}_{2}^{n}\), such that \(\mathcal {F}(x,0)=\mathbf {F}(x)\). That is, the case of no control coincides with the original BN. We remark that these control types can be combined, but for clarity we present them separately.

Definition 2.2

Consider the edge x i →x j in the wiring diagram \(\mathcal {W}\). The function

$$ \mathcal{F}_{j}(x,u_{i,j}) := f_{j}(x_{1},\dots,(u_{i,j}+1)x_{i},\dots,x_{n}) $$
(2)

encodes the control of the edge x i →x j , since for each possible value of \(u_{i,j}\in \mathbb {F}_{2}\) we have the following control settings:

  • If \(u_{i,j}=0, \mathcal {F}_{j}(x,0) = f_{j}(x_{1},\dots,x_{i},\dots,x_{n})\). That is, the control is not active.

  • If \(u_{i,j}=1, \mathcal {F}_{j}(x,1) = f_{j}(x_{1},\dots,x_{i}=0,\dots,x_{n})\). In this case, the control is active, and the action represents the removal of the edge x i →x j .

Definition 2.2 can easily be extended for the control of many edges, so that we obtain \(\mathcal {F}:\mathbb {F}_{2}^{n}\times \mathbb {F}_{2}^{e} \rightarrow \mathbb {F}_{2}^{n}\), where e is the number of edges in the wiring diagram. Each coordinate, u i,j , of u in \(\mathcal {F}(x,u)\) encodes the control of an edge x i →x j .

Example 2.1 (cont.) Incorporating edge control in the toy network results in

$$ \begin{array}{l} \mathcal{F}_{1} = 1+(u_{3,1}+1)x_{3}+(u_{5,1}+1)x_{5}+\\ \quad \quad \quad (u_{3,1}+1)x_{3}(u_{5,1}+1)x_{5},\\ \mathcal{F}_{2} = 1+(u_{1,2}+1)x_{1}+(u_{1,2}+1)x_{1}(u_{4,2}+1)x_{4},\\ \mathcal{F}_{3} = 1+ (u_{2,3}+1)x_{2}, \\ \mathcal{F}_{4} = (u_{3,4}+1)x_{3},\\ \mathcal{F}_{5} = 1+ (u_{4,5}+1)x_{4}. \end{array} $$

Definition 2.3

Consider the node x i in the wiring diagram \(\mathcal {W}\). The function

$$ \mathcal{F}_{j}(x,u^{-}_{i},u^{+}_{i}) := (u^{-}_{i}+u^{+}_{i}+1)f_{j}(x) + u^{+}_{i} $$
(3)

encodes the control (knock-out or constant expression) of the node x i , since for each possible value of \(\left (u^{-}_{i},u^{+}_{i}\right)\in \mathbb {F}_{2}^{2}\) we have the following control settings:

  • For \(u^{-}_{i}=0, u^{+}_{i}=0, \mathcal {F}_{j}(x,0,0) = f_{j}(x)\). That is, the control is not active.

  • For \(u^{-}_{i}=1, u^{+}_{i}=0, \mathcal {F}_{j}(x,1,0) = 0\). This action represents the knock out of the node x j .

  • For \(u^{-}_{i}=0, u^{+}_{i}=1, \mathcal {F}_{j}(x,0,1) = 1\). This action represents the constant expression of the node x j .

  • For \(u^{-}_{i}=1, u^{+}_{i}=1, \mathcal {F}_{j}(x,1,1) = f_{j}\left (x_{t_{1}},\dots,x_{t_{m}}\right)+1\). This action changes the Boolean function to its negative value and might not be a relevant case of control.

Example 2.1 (cont.) Incorporating node control in the toy network results in

$$ \begin{array}{l} \mathcal{F}_{1} = (u_{1}^{-}+u_{1}^{+} + 1)(1+x_{3}+x_{5}+x_{3}x_{5}) + u_{1}^{+},\\ \mathcal{F}_{2} = (u_{2}^{-}+u_{2}^{+} + 1)(1+x_{1}+x_{1}x_{4}) + u_{2}^{+},\\ \mathcal{F}_{3} = (u_{3}^{-}+u_{3}^{+} + 1)(1+ x_{2}) + u_{3}^{+}, \\ \mathcal{F}_{4} = (u_{4}^{-}+u_{4}^{+} + 1)x_{3} + u_{4}^{+},\\ \mathcal{F}_{5} = (u_{5}^{-}+u_{5}^{+} + 1)(1+ x_{4}) + u_{5}^{+}. \end{array} $$

Control targets in Boolean networks

We consider a BN with control \(\mathcal {F}:\mathbb F_{2}^{n}\times U \rightarrow \mathbb F_{2}^{n}\), and denote by F the BN with no control (\(\mathcal {F}(x,0) = \mathbf {F}(x)\)). We remark that in each case of interest, both edge and node control could be analyzed simultaneously.

Generating new steady states

Suppose that \(\textbf {y}_{0} = (y_{01},\dots,y_{0n})\in \mathbb F_{2}^{n}\) is a desirable cell state (for instance, it could represent the state of cell senescence) but is not a fixed point, i.e., F(y 0)≠y 0. The problem is then to choose a control u such that \(\mathcal {F}(\mathbf {y}_{0},u)= \mathbf {y}_{0}\). We now show how this can be achieved in our framework.

After encoding our BN with control as a polynomial system \(\mathcal {F}_{j}(x,u)\in \mathbb {F}_{2}[x,u]\) (see Section “Control Actions: edge and node manipulations”), we consider the system of polynomial equations in the u parameters:

$$ \mathcal{F}_{j}(\textbf{y}_{0},u)- y_{0j}=0, j = 1,\dots,m. $$
(4)

Example 2.1 (cont.) Here we consider a toy network and assume we are interested in controlling edges to make y 0=01111 a steady state. For simplicity we only consider control using the edges x 1→x 2,x 4→x 2,x 2→x 3,x 4→x 5. The network is

$$ \begin{array}{l} \mathcal{F}_{1} = 1+x_{3}+x_{5}+x_{3}x_{5},\\ \mathcal{F}_{2} = 1+(u_{1,2}+1)x_{1}+(u_{1,2}+1)x_{1}(u_{4,2}+1)x_{4},\\ \mathcal{F}_{3} = 1+ (u_{2,3}+1)x_{2}, \\ \mathcal{F}_{4} = x_{3},\\ \mathcal{F}_{5} = 1+ (u_{4,5}+1)x_{4}. \end{array} $$

Evaluating at y 0=01111 we obtain

$$ \begin{array}{l} \mathcal{F}_{1} = 0, \ \ \mathcal{F}_{2} = 1, \ \ \mathcal{F}_{3} = u_{2,3}, \ \ \mathcal{F}_{4} = 1, \ \ \mathcal{F}_{5} = u_{4,5}. \end{array} $$

Then, 01111 will be a steady state if and only if \(\mathcal {F}_{i}=\textbf {y}_{0i}\) for i=1,…,5. This gives the following solution u 2,3=1,u 4,5=1. Thus, in this case there is a unique control, u 2,3=u 4,5=1, that guarantees y 0=01111 is a steady state (this control policy is illustrated in Fig. 1 f, u=D). In general, the equations can be more complex and nonlinear, so computational algebra is needed. See the “Results” section for applications to more complex models.

Example 2.1 (cont.) We now show how the problem of creating new steady states can also be solved using node controls. We again use the toy network, and assume that we want to make y 0=11110 a steady state. For simplicity we only consider control of nodes x 1 (constant expression), x 3 (knock-out and constant expression), and node x 4 (constant expression). The network is

$$ \begin{array}{l} \mathcal{F}_{1} = (u_{1}^{+} + 1)(1+x_{3}+x_{5}+x_{3}x_{5}) + u_{1}^{+},\\ \mathcal{F}_{2} = 1+x_{1}+x_{1}x_{4},\\ \mathcal{F}_{3} = (u_{3}^{-}+u_{3}^{+} + 1)(1+ x_{2}) + u_{3}^{+}, \\ \mathcal{F}_{4} = (u_{4}^{+} + 1)x_{3} + u_{4}^{+},\\ \mathcal{F}_{5} = 1+ x_{4}. \end{array} $$

Evaluating at y 0=11110 we obtain

$$ \begin{array}{l} \mathcal{F}_{1} = u_{1}^{+},\ \ \mathcal{F}_{2} = 1,\ \ \mathcal{F}_{3} = u_{3}^{+}, \ \ \mathcal{F}_{4} = 1,\ \ \mathcal{F}_{5} = 0. \end{array} $$

Then, 11110 will be a steady state if and only if \(u_{1}^{+}=1\) and \(u_{3}^{+}=1\). That is, neither control by itself is sufficient, but together they create the steady state (this control policy is illustrated in Fig. 1 f, u=C).

Destroying existing steady states, or, in general, blocking transitions

Suppose that \(\textbf {x}_{0}\in \mathbb F_{2}^{n}\) is an undesirable steady state of F(x), that is, F(x 0)=x 0 (for instance, it could represent a tumor proliferative cell state that needs to be avoided). The goal here is to find a set of controls such that \(\mathcal {F}(\textbf {x}_{0},u) \neq \textbf {x}_{0}\). More generally, we may want to avoid a transition between two states x 0 and z 0. That is, we want to find controls such that \(\mathcal {F}(\textbf {x}_{0},u) \neq \textbf {z}_{0}\). To solve this problem consider the following equation,

$$ (\mathcal{F}_{1}(\mathbf{x}_{0},u)-\mathbf{z}_{01} +1)\ldots (\mathcal{F}_{n}(\mathbf{x}_{0},u)-\mathbf{z}_{0n} +1)=0. $$
(5)

Equation 5 defines a polynomial equation in the u parameters. It can be shown that \(\mathcal {F}(\textbf {x}_{0},u) \neq \textbf {z}_{0}\) if and only if Eq. 5 is satisfied.

Example 2.1 (cont.) Here we focus on finding edges to block the transition from x 0=00101 to F(x 0)=01111. For simplicity we only consider control using the edges x 3→x 1,x 5→x 1,x 2→x 3,x 3→x 4. The network is

$$ \begin{array}{l} \mathcal{F}_{1} = 1+(u_{3,1}+1)x_{3}+(u_{5,1}+1)x_{5}+\\ \quad \quad \quad (u_{3,1}+1)x_{3}(u_{5,1}+1)x_{5},\\ \mathcal{F}_{2} = 1+x_{1}+x_{1}x_{4},\\ \mathcal{F}_{3} = 1+ (u_{2,3}+1)x_{2}, \\ \mathcal{F}_{4} = (u_{3,4}+1)x_{3},\\ \mathcal{F}_{5} = 1+ x_{4}. \end{array} $$

Evaluating at x 0=00101 we obtain

$$ \begin{array}{lllll} \mathcal{F}_{1} = u_{3,1}u_{5,1}, \mathcal{F}_{2} = 1, \mathcal{F}_{3} = 1, \mathcal{F}_{4} = u_{3,4}, \mathcal{F}_{5} = 1. \end{array} $$

Then, Eq. 5 becomes (u 3,1 u 5,1+1)(u 3,4+1)=0 which has two solutions: u 3,4=1 and u 3,1=u 5,1=1 (first control policy illustrated in Fig. 1 f, u=B). The second solution is what we refer to as a combinatorial control policy. Neither u 3,1=1 or u 5,1=1 is sufficient, but combined they block the transition. In general, the equations can be more complex and nonlinear, so computational algebra is needed.

Blocking regions in the state space

We now consider the case where we want the dynamics to avoid certain regions. For example, if a particular value of a variable, \(x_{k}=a\in \mathbb {F}_{2}\) triggers an undesirable pathway, or is the signature of an abnormal cell, then we want all steady states of the system to satisfy x k ≠a. In this case, we consider the systems of equations

$$ \begin{aligned} \mathcal{F}_{j}(x,u)- x_{j}=0, j = 1,\dots,m,\\ x_{k}-a=0. \end{aligned} $$
(6)

Note that, in contrast to Sections “Generating new steady states and Destroying existing steady states, or, in general, blocking transitions”, we are now using variables for x instead of specific values. Since the steady states with x k =a are to be avoided, we want to find controls u for which Eq. 6 has no solution.

Example 2.1 (cont.) Here we consider the toy network and focus on controlling nodes to avoid regions of the form x 3=0. For simplicity we only consider control using the nodes x 2 (knock-out), x 3 (constant expression), x 4 (constant expression). The network is

$$ \begin{array}{l} \mathcal{F}_{1} = 1+x_{3}+x_{5}+x_{3}x_{5},\\ \mathcal{F}_{2} = (u_{2}^{-} + 1)(1+x_{1}+x_{1}x_{4}),\\ \mathcal{F}_{3} = (u_{3}^{+} + 1)(1+ x_{2}) + u_{3}^{+}, \\ \mathcal{F}_{4} = (u_{4}^{+} + 1)x_{3} + u_{4}^{+},\\ \mathcal{F}_{5} = 1+ x_{4}. \end{array} $$

Then, Eq. 6 becomes

$$ \begin{array}{l} 1+x_{3}+x_{5}+x_{3}x_{5} +x_{1}=0,\\ (u_{2}^{-} + 1)(1+x_{1}+x_{1}x_{4}) +x_{2}=0,\\ (u_{3}^{+} + 1)(1+ x_{2}) + u_{3}^{+} + x_{3}=0, \\ (u_{4}^{+} + 1)x_{3} + u_{4}^{+} + x_{4}=0,\\ 1+ x_{4}+x_{5}=0,\\ x_{3}=0. \end{array} $$

In contrast with the previous examples, this system of equations cannot be analyzed by hand. In Section “Identifying control targets” we will show how computational algebra gives an equivalent, but simpler, system of equations.

Identifying control targets

In each case of Section “Control targets in Boolean networks” we obtained a system of equations (or a single equation) that we need to solve to find the appropriate controls. This can be done using computational algebra tools. For instance, we can compute the Gröbner basis of the ideal associated to Eq. 4, see [51],

$$ I=\left\langle \mathcal{F}_{1}(\textbf{y}_{0},u)- y_{01},\ldots,\mathcal{F}_{n}(\textbf{y}_{0},u)- y_{0n}\right\rangle. $$
(7)

Example 2.1 (cont.) Now we continue the previous example where the goal was to avoid regions of the form x 3=0. We encode the system of equations as \( I =\langle 1+x_{3}+x_{5}+x_{3}x_{5} +x_{1}, (u_{2}^{-} + 1)(1+x_{1}+x_{1}x_{4}) +x_{2}, (u_{3}^{+} + 1)(1+ x_{2}) + u_{3}^{+} + x_{3}, (u_{4}^{+} + 1)x_{3} + u_{4}^{+} + x_{4}, 1+ x_{4}+x_{5}, x_{3}\rangle \). Using computational algebra we find a Gröbner basis of this ideal: \(I=\langle u_{3}^{+}, u_{2}^{-}, x_{5}+u_{4}^{+} + 1, x_{4}+u_{4}^{+}, x_{3}, x_{2}+1, x_{1}+u_{4}^{+} \rangle \).

Thus, the original system of equations has the same solutions as the system

$$ \begin{array}{lllll} {}u_{3}^{+}\!=0,& \quad u_{2}^{-}\!=0,&\quad x_{5}+u_{4}^{+} \,+\, 1\,=\,0,&\quad x_{4}\,+\,u_{4}^{+}\,=\,0, \\ x_{3}=0,& \quad x_{2}+1=0, &\quad x_{1}+u_{4}^{+}\,=\,0. \end{array} $$

In order to avoid regions of the form x 3=0, we need to find parameters for which the system has no solutions. This is guaranteed if any of the polynomials is equal to 1. Namely, we select the equations that only have the control parameters \(u_{3}^{+}=0, u_{2}^{-}=0\). If we use \(u_{3}^{+}=1\) or \(u_{2}^{-}=1\), the system is guaranteed to have no solutions. Thus, \(u_{3}^{+}=1\) or \(u_{2}^{-}=1\) independently are sufficient to achieve the control of the system (second control policy illustrated in Fig. 1 f, u=A).

As we will show in Examples 2.4–2.5, the computation of a Gröbner basis allows us to read out all controls for which the system of equations has a solution. Furthermore, our algebraic approach can detect combinatorial control actions (control by the synergistic combination of more than one action).

Results

We test our control methods using two published models where potential intervention targets were identified along with experimental validations. In the first example we focus on control with edge deletions while for the second we use control with node deletions and constant expressions.

Example 2.4

P53-mdm2 network. A Boolean network model for the signaling response of DNA damage in a p53 network was developed in [5]; the wiring diagram of this model is reproduced in Fig. 2. The tumor suppressor protein p53 can trigger either cell cycle arrest or apoptosis in response to DNA damage in normal cells. The authors of [5] extended their analysis to a breast cancer cell line, MCF7, with the purpose of identifying potential therapeutic interventions in the network for the cancer cell. Thus, in this example we will focus on the cancer cell model where PTEN and p14ARf are always inactive (fixed to zero) and cyclinG is always active (fixed to 1), see Table 5 of [5]. This system can be represented as a discrete dynamical system \(\mathbf {F}=(f_{1},\ldots,f_{16}):\mathbb F_{2}^{16}\rightarrow \mathbb F_{2}^{16}\) with 16 nodes and 50 edges. We represent the nodes by

Fig. 2
figure2

The p53-mdm2 network adapted from [5]. Arrows in green represent activation while hammerhead arrows (in red) represent inhibition. Self loops were omitted, see text in Example 2.4 for an explanation. For the cancer cell model, PTEN and p14ARf are always inactive (fixed to zero) and cyclinG is always active (fixed to 1)

Full size image
$$ \begin{aligned} &x_{1} = {ATM}, x_{2} = {p53}, \\ &x_{3} = {Mdm2}, x_{4} = {MdmX}, \\ &x_{5} = {Wip1}, x_{6} = {cyclinG},\\ &x_{7} = {PTEN}, x_{8} = p21, \\ &x_{9} = {AKT}, x_{10} = {cyclinE}, \\ &x_{11} ={Rb}, x_{12} = {E2F1}, \\ &x_{13} = {p14ARf}, x_{14} = {Bcl2}, \\ &x_{15} = {Bax}, x_{16} = {caspase}. \end{aligned} $$
(8)

The polynomial functions for this network are listed in the Additional file 1. We remark that the original model in [5] considers threshold functions for all the regulatory rules and this type of functions might introduce self-loops in the wiring diagram. That is, when we translated the threshold rule for x i into a polynomial f i , the function f i might depend on x i . Notice that the self-loops were omitted in Fig. 2 to be consistent with the original model.

For this model, in the presence of DNA damage, the system has a single limit cycle representing the state of cell cycle arrest, where p53 and p21 are oscillating; see Fig. 3.

Fig. 3
figure3

States of limit cycle representing cell cycle arrest in the p53 model. The order of the vector entries follows the indexing in Eq. 8. The dashed edge represents the transition target to destroy the limit cycle

Full size image

Suppose that we want to find out which set of edges one can manipulate in order to destroy this limit cycle and to lead the system into a different attractor, one that represents a desirable cell state. Let us start by considering a desirable state that represents cell death,

$$ \textbf{y}_{0} = (1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1), $$
(9)

where x 16=caspase is active. In order to make y 0 a steady state with the approach described in the “Methods” section, we form the following system of polynomial equations,

$$ \mathcal{F}_{i}(\textbf{y}_{0},\mu) = y_{0i}, j = 1,\dots,16. $$
(10)

The solutions for the system of Eqs. 10 are given by the nonzero generators of the ideal associated to the system 10, given in Eq. 11,

$$ \begin{aligned} &\left\{u_{2,5}+1, u_{16,11}(u_{1,11}\,+\,1),u_{8,8}(u_{3,8}\,+\,1),u_{2,2}(u_{3,2}\,+\,1),\right.\\ &\left.u_{1,2}(u_{3,2}+1), u_{1,1}u_{12,1},u_{12,16}u_{16,16}(u_{8,16}+1)\right\} \end{aligned} $$
(11)

There are 945 solutions for Eq. 11. Each solution gives a controller that guarantees the desired steady state, but each controller may have a different impact on the dynamics of the network. Table 1 (middle row) shows one of the controllers from Eq. 11. In Additional file 1: Table S1, we list the top ten control combinations from Eq. 11, that is, those that give the largest basin of y 0 as well as the ten sets that give the smallest basin for y 0; see Additional file 1: Table S2. For each solution in Additional file 1: Tables S1-S2, we crossed-out the edges that become nonessential after the other controllers in the set have been applied; see the “Discussion” for more about nonessential edges.

Table 1 Difference of impact in the combinatorial action of edge deletions
Full size table

Furthermore, we can aim to destroy the limit cycle in Fig. 3. We target one of the transitions within this limit cycle, let us take the transition x 5→x 6, see dashed transition in Fig. 3. Blocking this transition gives additional controllers that help to stabilize the system at the desired fixed point; see the last row of Table 1.

The deletion of the edges mdm2→p53 and p53→Wip1 were identified in [5] by deleting one edge at a time and checking if the modified system had the desired attractor. Choi et al. [5] reported that the combinatorial action of these controllers increased the basin of attraction of the desired fixed point to more than 50 %, which was validated experimentally. However, doing this type of search for finding all possible combinations of controls is infeasible. Since for each edge we have 2 possible actions (control or no control), and there are 50 edges, then there are 250 networks to be analyzed in total. In contrast, the computational algebra approach of this paper allows to obtain all combinations of edges that guarantee the desired steady state of the system in one process.

Example 2.5

T-LGL network. A Boolean network model for the blood cancer large granular lymphocyte (T-LGL) leukemia was developed in [9] and further analyzed in [10, 33]. T-LGL leukemia is characterized by escaping cell death through abnormal mechanisms, which are insensitive to Fas-induced apoptosis, [9]. This network has 60 nodes but was reduced to a subnetwork of 16 nodes for steady state analysis, see Fig. 4. Here we use the 16-nodes network to identify potential control targets.

Fig. 4
figure4

Reduced T-LGL network adapted from [10]. Arrows in green represent activation while hammerhead arrows (in red) represent inhibition. All the negative edges from Apoptosis were omitted for simplicity

Full size image

We represent the nodes by

$$\begin{aligned} x_{1}\,=\,{CREB}, & x_{2}\,=\,{IFNG}, x_{3}\,=\,{P2},\\ x_{4}\,=\,{GPCR}, & x_{5}\,=\,{SMAD}, x_{6}\,=\,{Fas},\\ x_{7}\,=\,{sFas}, & x_{8}\,=\,{Ceramide}, x_{9}\,=\,{DISC},\\ x_{10}\,=\,{Caspase}, & x_{11}\,=\,{FLIP}, x_{12}\,=\,{BID},\\ x_{13}\,=\,{IAP}, & x_{14}\,=\,{MCL1}, x_{15}\,=\,{S1P},\\ x_{16}\,=\,{Apoptosis}. & \end{aligned} $$

The polynomial rules for this system are listed in Eq. 12,

$$ \begin{aligned} f_{1}&=x_{2}(1+x_{16}), f_{2}=(x_{5}+1)(x_{3}+1)(1+x_{16}),\\ f_{3}&= (x_{2}x_{3}+x_{2}+x_{3})(1+x_{16}),\\ f_{4}&= x_{15}(1+x_{16}), f_{5}= x_{4}(1+x_{16}),\\ f_{6}&= (x_{7}+1)(1+x_{16}), f_{7}= x_{15}(1+x_{16}),\\ f_{8}&= (x_{15}+1)x_{6}(1+x_{16}),\\ f_{9}&= (x_{6}x_{8}x_{11}+x_{6}x_{8}+x_{6}x_{11}+x_{6}+x_{8})(1+x_{16}),\\ f_{10}&= (x_{9}x_{12}x_{13}\!+x_{9}x_{12}+x_{12}x_{13}+x_{9}+x_{12})(1\!+x_{16}),\\ f_{11}&= (x_{9}+1)(1+x_{16}), f_{12}= (x_{14}+1)(1+x_{16}),\\ f_{13}&= (x_{12}+1)(1+x_{16}), f_{14}= (x_{9}+1)(1+x_{16}),\\ f_{15}&= (x_{8}+1)(1+x_{16}),\\ f_{16}&= x_{10}x_{16}+x_{16}+x_{10}. \end{aligned} $$
(12)

This system has three steady states, one that represents the normal state, x 0=0000000000000001, where Apoptosis is ON, and 0001101000101110, 0011101000101110, the disease states in which Caspase and Apoptosis are OFF.

To find the potential node deletions (or constant expressions) that will block the disease states we use Eq. 6,

$$ \begin{aligned} (u^{-}_{1}+u^{+}_{1}+1)f_{1}(\mathbf{x})+u^{+}_{1} - x_{1}=0\\ \vdots\\ (u^{-}_{16}+u^{+}_{16}+1)f_{16}(\mathbf{x})+u^{+}_{16} - x_{16}=0\\ x_{16}=0, x_{10}=0 \end{aligned} $$
(13)

Equation 13 encodes all parameters for which there is a steady state where Caspase and Apoptosis are OFF. Thus, we are interested in the parameters for which this system of equations has no solution. Since for each node we have 3 possible actions (no control, node deletion, constant expression), there are 316 networks to be analyzed in total. Thus, even for this small network, exhaustive search is computationally challenging.

On the other hand, computational algebra allows us to obtain the parameter combinations that guarantee that the disease states are not fixed points of the system (see the Additional file 1 for details). The parameter combinations (enclosed in brackets, where entries not shown are equal to zero) are

$$ \begin{aligned} \{u^{+}_{8}=1\}, \{u^{+}_{9}=1\}, \{u^{+}_{10}=1\}, \{u^{+}_{12}=1\}, \{u^{-}_{14}=1\},\\ \{u^{-}_{15}=1\}, \{u^{+}_{6}=1,u^{-}_{11}=1\}, \{u^{-}_{7}=1,u^{-}_{11}=1\}. \end{aligned} $$
(14)

Thus we obtain that the constant expression of Ceramide, DISC, Caspase, or BID, or the deletion of MCL1 or S1P, will guarantee that the disease states are not steady states of the system (Table 2). These controls could also be found by trying one control at a time as was done in [10]. Importantly, our computational algebra approach shows that there are two additional control policies that consist of a combination of different controls. It can be shown that neither Fas, sFas, FLIP individually can eliminate the disease states, but the deletion of FLIP combined with the constant expression of Fas or the deletion of sFas will work (Table 2). Furthermore, those are the only combinations that guarantee the disease states will not be attractors.

Table 2 Control nodes for the reduced T-LGL network
Full size table

We note that, in the worst case, computing the Gröbner basis for a system of polynomial equations has doubly exponential complexity in the number of solutions. However, for the type of networks discussed in this paper, namely, biological networks where most of the nodes are regulated by only a small subset of the other nodes, Gröbner bases can be computed in a reasonable time, see [42].

Discussion

The design of control policies for gene regulatory networks is an important challenge in systems biology. The method described here exploits the interplay between the structure and the dynamics of the network to identify potential control interventions that will drive the system towards desired dynamics. The formulation of the problem as that of finding all solutions to a system of polynomial equations provides an alternative to exhaustive search of all possible combinations of interventions, which often is not feasible.

One shortcoming of the method is that, in its current form, it requires the Boolean network model to be updated synchronously, with deterministic dynamics. While steady states of the network do not depend on the update order used, general limit cycles and attractor basins do, however. Thus, some of the methods described here might not be applicable in a stochastic setting. For instance, if the dynamics is generated using the asynchronous update or more generally using the stochastic settings in [8, 52–54], then encoding the controllers for blocking transitions will need to have a different setup than the one described here. However, the method for producing a new steady state is still valid for all variants of stochastic updates because the system will maintain the steady state. Another shortcoming is that the control methods described in this paper were developed for Boolean networks only. However, many published discrete models of GRNs are multistate. These methods could potentially be extended to a more general setting, where the network variables might attain more than two states [41, 54–57]. We remark that it is possible to map a multistate model into a Boolean model (see [58]) where our methods can be applied and then it would be possible to recover a control set for the multistate model.

An important challenge in the process of identification of control targets in a network is to develop methods that can identify controllers that can guarantee global reachability of a desired steady state. This problem is not addressed in this paper, and remains to be studied in the future. For instance, the control strategies do not give any information about the basin size of a fixed point generated by the methods of this paper. However, we remark that some algebraic methods allow to estimate the change in the basin size after an edge deletion, see [59]. Nonetheless, the control targets identified by the algebraic techniques described here could be used for further analysis of the system, such as for studying reachability [60], or for designing optimal control policies in a stochastic setting [29–32].

Finally, it is worth pointing out that the methods of this paper might produce a large number of control strategies, which all give the desired result, but many of these might be biologically meaningless or infeasible as actual interventions. Eliminating those might be challenging. For instance, some solution sets might contain nonessential [61] or nonfunctional [62] edges. That is, an edge could become nonessential after the other controllers in the set have been applied. In Additional file 1: Table S1, we list the top ten control combinations for Example 2.4 where we crossed-out the edges that become nonessential after the other controllers have been applied. Moreover, one can group the solution sets by considering only minimal sets where all the control edges are functional. For instance, we can group the first four solutions of Additional file 1: Table S1 into one group with the minimal representative set given in the middle row of Table 1.

Conclusions

This paper presents a novel approach to the identification of potential interventions in Boolean molecular networks. The methods use the theoretical foundations of algebraic geometry to encode the structure of a network by a set of polynomials and then, with the use of computer algebra techniques, find a set of nodes and edges to perform interventions in silico. The methods were validated using two published models where dynamic network interventions were identified, the p53-mdm2 system and the T-LGL leukemia model. It was shown that the methods in this paper can identify the controllers that were already reported and also find new potential targets. Some of these new control targets are combinatorial in nature and might result in more efficient strategies as was shown in the “Results” section using the change in the basin size of the system as a measure of efficiency.

References

  1. 1

    Tyson JJ, Chen K, Novak B. Network dynamics and cell physiology. Nat Rev Mol Cell Biol. 2001; 2(12):908–16.

  2. 2

    Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cell-cycle network is robustly designed. Proc Natl Acad Sci U S A. 2004; 101(14):4781–6. doi:10.1073/pnas.0305937101.

  3. 3

    Geva-Zatorsky N, Rosenfeld N, Itzkovitz S, Milo R, Sigal A, Dekel E, Yarnitzky T, Liron Y, Polak P, Lahav G, Alon U. Oscillations and variability in the p53 system. Mol Syst Biol. 2006; 2:2006–0033. doi:10.1038/msb4100068.

  4. 4

    Batchelor E, Loewer A, Lahav G. The ups and downs of p53: understanding protein dynamics in single cells. Nat Rev Cancer. 2009; 9(5):371–7. doi:10.1038/nrc2604.

  5. 5

    Choi M, Shi J, Jung SH, Chen X, Cho KH. Attractor landscape analysis reveals feedback loops in the p53 network that control the cellular response to dna damage. Sci Signal. 2012; 5(251):83.

  6. 6

    Joh RI, Weitz JS. To lyse or not to lyse: Transient-mediated stochastic fate determination in cells infected by bacteriophages. PLoS Comput Biol. 2011; 7(3). doi:10.1371/journal.pcbi.1002006.

  7. 7

    Zeng L, Skinner SO, Zong C, Sippy J, Feiss M, Golding I. Decision making at a subcellular level determines the outcome of bacteriophage infection. Cell. 2010; 141(4):682–91. doi:10.1016/j.cell.2010.03.034.

  8. 8

    Murrugarra D, Veliz-Cuba A, Aguilar B, Arat S, Laubenbacher R. Modeling stochasticity and variability in gene regulatory networks. EURASIP J Bioinforma Syst Biol. 2012; 2012(1):5.

  9. 9

    Zhang R, Shah MV, Yang J, Nyland SB, Liu X, Yun JK, Albert R, Jr Loughran TP. Network model of survival signaling in large granular lymphocyte leukemia. Proc Natl Acad Sci U S A. 2008; 105(42):16308–13. doi:10.1073/pnas.0806447105.

  10. 10

    Saadatpour A, Wang RS, Liao A, Liu X, Loughran TP, Albert I, Albert R. Dynamical and structural analysis of a t cell survival network identifies novel candidate therapeutic targets for large granular lymphocyte leukemia. PLoS Comput Biol. 2011; 7(11):1002267. doi:10.1371/journal.pcbi.1002267.

  11. 11

    Wang W. Therapeutic hints from analyzing the attractor landscape of the p53 regulatory circuit. Sci Signal. 2013; 6(261):5. doi:10.1126/scisignal.2003820.

  12. 12

    Lee MJ, Ye AS, Gardino AK, Heijink AM, Sorger PK, MacBeath G, Yaffe MB. Sequential application of anticancer drugs enhances cell death by rewiring apoptotic signaling networks. Cell. 2012; 149(4):780–94. doi:10.1016/j.cell.2012.03.031.

  13. 13

    Wang Z, Deisboeck TS. Mathematical modeling in cancer drug discovery. Drug Discov Today. 2014; 19(2):145–50. doi:10.1016/j.drudis.2013.06.015.

  14. 14

    Erler JT, Linding R. Network medicine strikes a blow against breast cancer. Cell. 2012; 149(4):731–3. doi:10.1016/j.cell.2012.04.014.

  15. 15

    Takahashi K, Yamanaka S. Induction of pluripotent stem cells from mouse embryonic and adult fibroblast cultures by defined factors. Cell. 2006; 126(4):663–76. doi:10.1016/j.cell.2006.07.024.

  16. 16

    Young RA. Control of embryonic stem cell state. Cell. 2011; 144(6):940–54. doi:10.1016/j.cell.2011.01.032.

  17. 17

    Iglesias PA, Ingalls BP. Control Theory and Systems Biology. Cambridge, MA [etc.]: MIT. (cop. 2010).

  18. 18

    Shin YJ, Bleris L. Linear control theory for gene network modeling. PLoS One. 2010; 5(9). doi:10.1371/journal.pone.0012785.

  19. 19

    Cornelius SP, Kath WL, Motter AE. Realistic control of network dynamics. Nat Commun. 2013; 4:1942. doi:10.1038/ncomms2939.

  20. 20

    Motter AE. Networkcontrology. Chaos. 2015; 25(9):097621. doi:10.1063/1.4931570.

  21. 21

    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):1–18.

  22. 22

    Kauffman S, Peterson C, Samuelsson B, Troein C. Random boolean network models and the yeast transcriptional network. Proc Natl Acad Sci. 2003; 100(25):14796–9.

  23. 23

    Saez-Rodriguez J, Simeoni L, Lindquist JA, Hemenway R, Bommhardt U, Arndt B, Haus U, Weismantel R, Gilles ED, Klamt S, Schraven B. A logical model provides insights into T cell receptor signaling. PLoS Comput Biology. 2007; 3(8). doi:10.1371/journal.pcbi.0030163.

  24. 24

    Balleza E, Alvarez-Buylla ER, Chaos A, Kauffman S, Shmulevich I, Aldana M. Critical dynamics in genetic regulatory networks: examples from four kingdoms. PLoS One. 2008; 3(6):2456. doi:10.1371/journal.pone.0002456.

  25. 25

    Davidich MI, Bornholdt S. Boolean network model predicts cell cycle sequence of fission yeast. PLoS One. 2008; 3(2):1672. doi:10.1371/journal.pone.0001672.

  26. 26

    Veliz-Cuba A, Stigler B. Boolean models can explain bistability in the lac operon. J Comput Biol. 2011; 18(6):783–94. doi:10.1089/cmb.2011.0031.

  27. 27

    Veliz-Cuba A, Arthur J, Hochstetler L, Klomps V, Korpi E. On the relationship of steady states of continuous and discrete models arising from biology. Bull Math Biol. 2012; 74(12):2779–92. doi:10.1007/s11538-012-9778-1.

  28. 28

    Veliz-Cuba A, Kumar A, Josić K. Piecewise linear and boolean models of chemical reaction networks. Bull Math Biol. 2014; 76(12):2945–84. doi:10.1007/s11538-014-0040-x.

  29. 29

    Yousefi MR, Datta A, Dougherty ER. Optimal intervention strategies for therapeutic methods with fixed-length duration of drug effectiveness. Signal Process IEEE Trans. 2012; 60(9):4930–44.

  30. 30

    Yousefi MR, Dougherty ER. Intervention in gene regulatory networks with maximal phenotype alteration. Bioinformatics. 2013; 29(14):1758–67. doi:10.1093/bioinformatics/btt242. http://bioinformatics.oxfordjournals.org/content/29/14/1758.full.pdf+html.

  31. 31

    Yousefi MR, Datta A, Dougherty ER. Optimal intervention in markovian gene regulatory networks with random-length therapeutic response to antitumor drug. Biomed Eng IEEE Transac. 2013; 60(12):3542–52. doi:10.1109/TBME.2013.2272891.

  32. 32

    Yousefi MR, Dougherty ER. A comparison study of optimal and suboptimal intervention policies for gene regulatory networks in the presence of uncertainty. EURASIP J Bioinforma Syst Biol. 2014; 2014(1):6–6. doi:10.1186/1687-4153-2014-6.

  33. 33

    Zañudo JGT, Albert R. Cell fate reprogramming by control of intracellular network dynamics. PLoS Comput Biol. 2015; 11(4):1004193. doi:10.1371/journal.pcbi.1004193.

  34. 34

    Qiu Y, Tamura T, Ching WK, Akutsu T. On control of singleton attractors in multiple boolean networks: integer programming-based method. BMC Syst Biol. 2014; 8 Suppl 1:7. doi:10.1186/1752-0509-8-S1-S7.

  35. 35

    Shmulevich I, Dougherty ER, Zhang W. Gene perturbation and intervention in probabilistic boolean networks. Bioinformatics. 2002; 18(10):1319–31. doi:10.1093/bioinformatics/18.10.1319. http://bioinformatics.oxfordjournals.org/content/18/10/1319.full.pdf+html.

  36. 36

    Xiao Y, Dougherty ER. The impact of function perturbations in boolean networks. Bioinformatics. 2007; 23(10):1265–73. doi:10.1093/bioinformatics/btm093.

  37. 37

    Vera-Licona P, Bonnet E, Barillot E, Zinovyev A. Ocsana: optimal combinations of interventions from network analysis. Bioinformatics. 2013; 29(12):1571–3. doi:10.1093/bioinformatics/btt195.

  38. 38

    Poret A, Boissel JP. An in silico target identification using boolean network attractors: Avoiding pathological phenotypes. Comptes rendus biologies. 2014; 337(12):661–78. doi:10.1016/j.crvi.2014.10.002.

  39. 39

    Grayson DR, Stillman ME. Macaulay2, a Software System for Research in Algebraic Geometry. http://www.math.uiuc.edu/Macaulay2/.

  40. 40

    Veliz-Cuba A, Aguilar B, Hinkelmann F, Laubenbacher R. Steady state analysis of boolean molecular network models via model reduction and computational algebra. BMC Bioinforma. 2014; 15:221. doi:10.1186/1471-2105-15-221.

  41. 41

    Veliz-Cuba A, Jarrah AS, Laubenbacher R. Polynomial algebra of discrete models in systems biology. Bioinformatics. 2010; 26(13):1637–43. doi:10.1093/bioinformatics/btq240.

  42. 42

    Hinkelmann F, Brandon M, Guang B, McNeill R, Blekherman G, Veliz-Cuba A, Laubenbacher R. Adam: analysis of discrete models of biological systems using computer algebra. BMC Bioinforma. 2011; 12:295. doi:10.1186/1471-2105-12-295.

  43. 43

    Veliz-Cuba A. Reduction of Boolean network models. J Theor Biol. 2011; 289:167–72.

  44. 44

    Curto C, Itskov V, Veliz-Cuba A, Youngs N. The neural ring: An algebraic tool for analyzing the intrinsic structure of neural codes. Bull Math Biol. 2013; 75(9):1571–611.

  45. 45

    Veliz-Cuba A. An algebraic approach to reverse engineering finite dynamical systems arising from biology. SIAM J Appl Dyn Syst. 2012; 11(1):31–48.

  46. 46

    Jarrah AS, Laubenbacher R, Stigler B, Stillman M. Reverse–engineering of polynomial dynamical systems. Adv Appl Math. 2007; 39(4):477–89.

  47. 47

    Laubenbacher R, Stigler B. A computational algebra approach to the reverse engineering of gene regulatory networks. J Theor Biol. 2004; 229(4):523–37. doi:10.1016/j.jtbi.2004.04.037.

  48. 48

    Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969; 22(3):437–67.

  49. 49

    Huang S. Gene expression profiling, genetic networks, and cellular states: an integrating concept for tumorigenesis and drug discovery. J Mol Med (Berl). 1999; 77(6):469–80.

  50. 50

    Shmulevich I, Dougherty ER. Probabilistic Boolean Networks - The Modeling and Control of Gene Regulatory Networks: SIAM; 2010. http://dx.doi.org/10.1137/1.9780898717631.

  51. 51

    Cox D, Little J, O’shea D. Using algebraic geometry, volume 185 of Graduate Texts in Mathematics. New York: Springer-Verlag: 1998.

  52. 52

    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.

  53. 53

    Thomas R. Regulatory networks seen as asynchronous automata: a logical description. J Theor Biol. 1991; 153(1):1–23.

  54. 54

    Thomas R, D’Ari R. Biological Feedback. Boca Raton: CRC Press; 1990. https://lccn.loc.gov/89023874.

  55. 55

    Murrugarra D, Laubenbacher R. The number of multistate nested canalyzing functions. Physica D: Nonlinear Phenomena. 2012; 241(10):929–38. doi:10.1016/j.physd.2012.02.011.

  56. 56

    Murrugarra D, Laubenbacher R. Regulatory patterns in molecular interaction networks. J Theor Biol. 2011; 288(0):66–72. doi:10.1016/j.jtbi.2011.08.015.

  57. 57

    Ahmad J, Niazi U, Mansoor S, Siddique U, Bibby J. Formal modeling and analysis of the mal-associated biological regulatory network: insight into cerebral malaria. PLoS One. 2012; 7(3):33532. doi:10.1371/journal.pone.0033532.

  58. 58

    Didier G, Remy E, Chaouiya C. Mapping multivalued onto boolean dynamics. J Theor Biol. 2011; 270(1):177–84. doi:10.1016/j.jtbi.2010.09.017.

  59. 59

    Murrugarra D, Dimitrova ES. Molecular network control through boolean canalization. EURASIP J Bioinform Syst Biol. 2015; 2015(1):9. doi:10.1186/s13637-015-0029-2.

  60. 60

    Li R, Yang M, Chu T. Controllability and observability of boolean networks arising from biology. Chaos. 2015; 25(2):023104. doi:10.1063/1.4907708.

  61. 61

    Li Y, Adeyeye JO, Murrugarra D, Aguilar B, Laubenbacher R. Boolean nested canalizing functions: A comprehensive analysis. Theor Comput Sci. 2013; 481(0):24–36. doi:10.1016/j.tcs.2013.02.020.

  62. 62

    Comet JP, Noual M, Richard A, Aracena J, Calzone L, Demongeot J, Kaufman M, Naldi A, Snoussi EH, Thieffry D. On circuit functionality in boolean networks. Bull Math Biol. 2013; 75(6):906–19.

Download references

Acknowledgements

The authors thank the referees for their very insightful comments that have improved the manuscript. RL was supported by a grant from the U.S. Army Research Office, Grant Nr. W911NF-14-1-0486.

Authors’ contributions

DM, AVC, BA, and RL conceived the study. DM and AVC performed the numerical experiments and theoretical analysis. BA designed code. All authors helped in the writing of the manuscript. All authors approved the final version of the manuscript.

Competing interests

The authors declare that they have no competing interests.

Author information

Author notes

    Affiliations

    1. Department of Mathematics, University of Kentucky, Lexington, 40506-0027, KY, USA

      • David Murrugarra

    2. Department of Mathematics, University of Dayton, Dayton, 45469, OH, USA

      • Alan Veliz-Cuba

    3. Institute for Systems Biology, Seattle, 98109-5263, WA, USA

      • Boris Aguilar

    4. Center for Quantitative Medicine, University of Connecticut Health Center, Farmington, 06030-6033, CT, USA

      • Reinhard Laubenbacher

    5. Jackson Laboratory for Genomic Medicine, Farmington, 06030, CT, USA

      • Reinhard Laubenbacher

    Authors

    1. Search for David Murrugarra in:

    2. Search for Alan Veliz-Cuba in:

    3. Search for Boris Aguilar in:

    4. Search for Reinhard Laubenbacher in:

    Corresponding author

    Correspondence to David Murrugarra.

    Additional file

    Additional file 1

    Supporting material for: Identification of control targets in Boolean molecular network models via computational algebra. (PDF 214 KB)

    Rights and permissions

    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.

    Reprints and Permissions

    About this article

    Verify currency and authenticity via CrossMark

    Keywords

    • Boolean networks
    • Algebraic control
    • Network control
    • Edge deletions
    • Blocking transitions

    BMC Systems Biology

    ISSN: 1752-0509

    Contact us