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.
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 x1,…,xn of variables that take values in the binary set {0,1}. Then a Boolean network in the variables x1,…,xn is a function
where each coordinate function fi:{0,1}n→{0,1} is a Boolean function on a subset of {x1,…,xn} 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
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=(f1,…,fn), a directed graph \(\mathcal {W}\) with n nodes x1,…,xn is associated with F. There is a directed edge in \(\mathcal {W}\) from xj to xi if xj appears in fi. Notice that the presence of the interaction xj→xi implies that the regulatory function fi depends on xj, 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. 1a.
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 (LF) 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=ab,a∨b=a+b+ab, and ¬a=1+a, where the operations are computed modulo 2. For the toy network we obtain
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 Fk(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 fi=xi,i=1,…,5. This means we want to find the roots of gi=0, where gi=fi−xi. This gives the system of equations
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=〈g1,g2,g3,g4,g5〉, 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
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. 1a). 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 xi→xj in the wiring diagram \(\mathcal {W}\). The function
encodes the control of the edge xi→xj, 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 xi→xj.
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, ui,j, of u in \(\mathcal {F}(x,u)\) encodes the control of an edge xi→xj.
Example 2.1 (cont.) Incorporating edge control in the toy network results in
encodes the control (knock-out or constant expression) of the node xi, 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 xj.
For \(u^{-}_{i}=0, u^{+}_{i}=1, \mathcal {F}_{j}(x,0,1) = 1\). This action represents the constant expression of the node xj.
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
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(y0)≠y0. 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:
Example 2.1 (cont.) Here we consider a toy network and assume we are interested in controlling edges to make y0=01111 a steady state. For simplicity we only consider control using the edges x1→x2,x4→x2,x2→x3,x4→x5. The network is
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 u2,3=1,u4,5=1. Thus, in this case there is a unique control, u2,3=u4,5=1, that guarantees y0=01111 is a steady state (this control policy is illustrated in Fig. 1f, 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 y0=11110 a steady state. For simplicity we only consider control of nodes x1 (constant expression), x3 (knock-out and constant expression), and node x4 (constant expression). The network is
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. 1f, 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(x0)=x0 (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 x0 and z0. 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,
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 x0=00101 to F(x0)=01111. For simplicity we only consider control using the edges x3→x1,x5→x1,x2→x3,x3→x4. The network is
Then, Eq. 5 becomes (u3,1u5,1+1)(u3,4+1)=0 which has two solutions: u3,4=1 and u3,1=u5,1=1 (first control policy illustrated in Fig. 1f, u=B). The second solution is what we refer to as a combinatorial control policy. Neither u3,1=1 or u5,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 xk≠a. In this case, we consider the systems of equations
Example 2.1 (cont.) Here we consider the toy network and focus on controlling nodes to avoid regions of the form x3=0. For simplicity we only consider control using the nodes x2 (knock-out), x3 (constant expression), x4 (constant expression). The network is
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],
Example 2.1 (cont.) Now we continue the previous example where the goal was to avoid regions of the form x3=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
In order to avoid regions of the form x3=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. 1f, 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)
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 xi into a polynomial fi, the function fi might depend on xi. 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,
where x16=caspase is active. In order to make y0 a steady state with the approach described in the “Methods” section, we form the following system of polynomial equations,
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 y0 as well as the ten sets that give the smallest basin for y0; 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
Control edges that increase the basin of attraction of cell death represented by y0 in Eq. 9. There are 216=65536 possible states. The number in parentheses is the ratio between the basin size and the total number of states
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 x5→x6, 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
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
This system has three steady states, one that represents the normal state, x0=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,
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
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 FLIPcombined 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
Solution
Control targets
Attractor
Basin size
\(u^{+}_{8}=1\)
Ceramide=ON
0000000100000001
100 %
\(u^{+}_{9}=1\)
DISC=ON
0000000010000001
100 %
\(u^{+}_{10}=1\)
Caspase=ON
0000000001000001
100 %
\(u^{+}_{12}=1\)
BID=ON
0000000000010001
100 %
\(u^{-}_{14}=1\)
MCL1=OFF
0000000000000001
100 %
\(u^{-}_{15}=1\)
S1P=OFF
0000000000000001
100 %
\(u^{+}_{6}=1\)
Fas=ON
0000010000000001
100 %
\(u^{-}_{11}=1\)
FLIP=OFF
\(u^{-}_{7}=1\)
sFas=OFF
0000000000000001
100 %
\(u^{-}_{11}=1\)
FLIP=OFF
The last two rows represent combinatorial actions of two nodes. All attractors are steady states, and the basin sizes include the steady states themselves. Notice that node x16=Apoptosis is a conceptual node in this model, thus it is not a relevant solution for network control
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.
Declarations
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.
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.
Authors’ Affiliations
(1)
Department of Mathematics, University of Kentucky
(2)
Department of Mathematics, University of Dayton
(3)
Institute for Systems Biology
(4)
Center for Quantitative Medicine, University of Connecticut Health Center
(5)
Jackson Laboratory for Genomic Medicine
References
Tyson JJ, Chen K, Novak B. Network dynamics and cell physiology. Nat Rev Mol Cell Biol. 2001; 2(12):908–16.PubMedView ArticleGoogle Scholar
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.PubMedPubMed CentralView ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.
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.View ArticleGoogle Scholar
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.PubMedPubMed CentralView ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.PubMedView ArticleGoogle Scholar
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.PubMedPubMed CentralView ArticleGoogle Scholar
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.
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.View ArticleGoogle Scholar
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.PubMedView ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.PubMedView ArticleGoogle Scholar
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.View ArticleGoogle Scholar
Jarrah AS, Laubenbacher R, Stigler B, Stillman M. Reverse–engineering of polynomial dynamical systems. Adv Appl Math. 2007; 39(4):477–89.View ArticleGoogle Scholar
Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969; 22(3):437–67.PubMedView ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.
Cox D, Little J, O’shea D. Using algebraic geometry, volume 185 of Graduate Texts in Mathematics. New York: Springer-Verlag: 1998.Google Scholar
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.PubMedView ArticleGoogle Scholar
Thomas R. Regulatory networks seen as asynchronous automata: a logical description. J Theor Biol. 1991; 153(1):1–23.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.PubMedView ArticleGoogle Scholar
