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 ith 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.
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 (knockout 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} (knockout 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} (knockout), 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).