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

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 ₁,…,x _n of variables that take values in the binary set {0,1}. Then a Boolean network in the variables x ₁,…,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}ⁿ→{0,1} is a Boolean function on a subset of {x ₁,…,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 ₁,…,f _n ), a directed graph $\mathcal {W}$ with n nodes x ₁,…,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

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

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 ₁,g ₂,g ₃,g ₄,g ₅〉, 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 ₀)≠y ₀. 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 ₀=01111 a steady state. For simplicity we only consider control using the edges x ₁→x ₂,x ₄→x ₂,x ₂→x ₃,x ₄→x ₅. 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 ₀=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 ₀=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 ₀=11110 a steady state. For simplicity we only consider control of nodes x ₁ (constant expression), x ₃ (knock-out and constant expression), and node x ₄ (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 ₀=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 ₀)=x ₀ (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 ₀ and z ₀. 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 ₀=00101 to F(x ₀)=01111. For simplicity we only consider control using the edges x ₃→x ₁,x ₅→x ₁,x ₂→x ₃,x ₃→x ₄. 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 ₀=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 ₃=0. For simplicity we only consider control using the nodes x ₂ (knock-out), x ₃ (constant expression), x ₄ (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 ₃=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 ₃=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

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)

$$ \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

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

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 ₁₆=caspase is active. In order to make y ₀ 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 ₀ as well as the ten sets that give the smallest basin for y ₀; 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

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 ₅→x ₆, 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 2⁵⁰ 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

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

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 ₀=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 3¹⁶ 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

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.