 Research Article
 Open Access
Identification of control targets in Boolean molecular network models via computational algebra
 David Murrugarra†^{1}Email author,
 Alan VelizCuba†^{2},
 Boris Aguilar^{3} and
 Reinhard Laubenbacher^{4, 5}
https://doi.org/10.1186/s129180160332x
© The Author(s) 2016
Received: 8 September 2015
Accepted: 23 August 2016
Published: 23 September 2016
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 p53mdm2 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.
Keywords
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 p53mdm2 system, [3–5], the phagelambda system, [6–8], or the T cell large granular lymphocyte (TLGL) 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 coarsegrained 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
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
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.
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.
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].
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

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 }.
Definition 2.3

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.
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.
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.
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
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.
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
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.
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
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 \).
In order to avoid regions of the form x _{3}=0, we need to find parameters for which the system has no solutions. This is guaranteed if any of the polynomials is equal to 1. Namely, we select the equations that only have the control parameters \(u_{3}^{+}=0, u_{2}^{}=0\). If we use \(u_{3}^{+}=1\) or \(u_{2}^{}=1\), the system is guaranteed to have no solutions. Thus, \(u_{3}^{+}=1\) or \(u_{2}^{}=1\) independently are sufficient to achieve the control of the system (second control policy illustrated in Fig. 1 f, u=A).
As we will show in Examples 2.4–2.5, the computation of a Gröbner basis allows us to read out all controls for which the system of equations has a solution. Furthermore, our algebraic approach can detect combinatorial control actions (control by the synergistic combination of more than one action).
Results
We test our control methods using two published models where potential intervention targets were identified along with experimental validations. In the first example we focus on control with edge deletions while for the second we use control with node deletions and constant expressions.
Example 2.4
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 selfloops 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 selfloops were omitted in Fig. 2 to be consistent with the original model.
Difference of impact in the combinatorial action of edge deletions
Controllers applied  Ref.  Basin size of y _{0} 

m d m2→p53  [5]  35581 (54.29 %) 
p53→W i p1  
p53→W i p1  A control set that  39856 (60.82 %) 
M d m2→p21  forces y _{0} to be  
M d m2→p53  a fixed point, from  
p21→C a s p a s e  Eq. 11.  
m d m2→p53  A control set to make  65536 (100 %) 
p53→W i p1  y _{0} a fixed point  
m d m2→p21  and for blocking  
p21→C a s p a s e  the transition in red  
A T M→R b  at Fig. 3.  
m d m2→R b  
m d m x→p53  
R b→E2F1  
B c l2→B a x 
Furthermore, we can aim to destroy the limit cycle in Fig. 3. We target one of the transitions within this limit cycle, let us take the transition x _{5}→x _{6}, see dashed transition in Fig. 3. Blocking this transition gives additional controllers that help to stabilize the system at the desired fixed point; see the last row of Table 1.
The deletion of the edges mdm2→p53 and p53→Wip1 were identified in [5] by deleting one edge at a time and checking if the modified system had the desired attractor. Choi et al. [5] reported that the combinatorial action of these controllers increased the basin of attraction of the desired fixed point to more than 50 %, which was validated experimentally. However, doing this type of search for finding all possible combinations of controls is infeasible. Since for each edge we have 2 possible actions (control or no control), and there are 50 edges, then there are 2^{50} 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
This system has three steady states, one that represents the normal state, x _{0}=0000000000000001, where Apoptosis is ON, and 0001101000101110, 0011101000101110, the disease states in which Caspase and Apoptosis are OFF.
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^{16} networks to be analyzed in total. Thus, even for this small network, exhaustive search is computationally challenging.
Control nodes for the reduced TLGL 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 
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 crossedout 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 p53mdm2 system and the TLGL 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. W911NF1410486.
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
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
 Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cellcycle network is robustly designed. Proc Natl Acad Sci U S A. 2004; 101(14):4781–6. doi:10.1073/pnas.0305937101.PubMedPubMed CentralView ArticleGoogle Scholar
 GevaZatorsky 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
 Batchelor E, Loewer A, Lahav G. The ups and downs of p53: understanding protein dynamics in single cells. Nat Rev Cancer. 2009; 9(5):371–7. doi:10.1038/nrc2604.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: Transientmediated stochastic fate determination in cells infected by bacteriophages. PLoS Comput Biol. 2011; 7(3). doi:10.1371/journal.pcbi.1002006.
 Zeng L, Skinner SO, Zong C, Sippy J, Feiss M, Golding I. Decision making at a subcellular level determines the outcome of bacteriophage infection. Cell. 2010; 141(4):682–91. doi:10.1016/j.cell.2010.03.034.PubMedPubMed CentralView ArticleGoogle Scholar
 Murrugarra D, VelizCuba 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
 Wang W. Therapeutic hints from analyzing the attractor landscape of the p53 regulatory circuit. Sci Signal. 2013; 6(261):5. doi:10.1126/scisignal.2003820.View ArticleGoogle Scholar
 Lee MJ, Ye AS, Gardino AK, Heijink AM, Sorger PK, MacBeath G, Yaffe MB. Sequential application of anticancer drugs enhances cell death by rewiring apoptotic signaling networks. Cell. 2012; 149(4):780–94. doi:10.1016/j.cell.2012.03.031.PubMedPubMed CentralView ArticleGoogle Scholar
 Wang Z, Deisboeck TS. Mathematical modeling in cancer drug discovery. Drug Discov Today. 2014; 19(2):145–50. doi:10.1016/j.drudis.2013.06.015.PubMedView ArticleGoogle Scholar
 Erler JT, Linding R. Network medicine strikes a blow against breast cancer. Cell. 2012; 149(4):731–3. doi:10.1016/j.cell.2012.04.014.PubMedView ArticleGoogle Scholar
 Takahashi K, Yamanaka S. Induction of pluripotent stem cells from mouse embryonic and adult fibroblast cultures by defined factors. Cell. 2006; 126(4):663–76. doi:10.1016/j.cell.2006.07.024.PubMedView ArticleGoogle Scholar
 Young RA. Control of embryonic stem cell state. Cell. 2011; 144(6):940–54. doi:10.1016/j.cell.2011.01.032.PubMedPubMed CentralView ArticleGoogle Scholar
 Iglesias PA, Ingalls BP. Control Theory and Systems Biology. Cambridge, MA [etc.]: MIT. (cop. 2010).Google Scholar
 Shin YJ, Bleris L. Linear control theory for gene network modeling. PLoS One. 2010; 5(9). doi:10.1371/journal.pone.0012785.
 Cornelius SP, Kath WL, Motter AE. Realistic control of network dynamics. Nat Commun. 2013; 4:1942. doi:10.1038/ncomms2939.PubMedPubMed CentralView ArticleGoogle Scholar
 Motter AE. Networkcontrology. Chaos. 2015; 25(9):097621. doi:10.1063/1.4931570.PubMedPubMed CentralView 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
 SaezRodriguez 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, AlvarezBuylla 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
 Davidich MI, Bornholdt S. Boolean network model predicts cell cycle sequence of fission yeast. PLoS One. 2008; 3(2):1672. doi:10.1371/journal.pone.0001672.View ArticleGoogle Scholar
 VelizCuba A, Stigler B. Boolean models can explain bistability in the lac operon. J Comput Biol. 2011; 18(6):783–94. doi:10.1089/cmb.2011.0031.PubMedView ArticleGoogle Scholar
 VelizCuba 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/s1153801297781.PubMedView ArticleGoogle Scholar
 VelizCuba A, Kumar A, Josić K. Piecewise linear and boolean models of chemical reaction networks. Bull Math Biol. 2014; 76(12):2945–84. doi:10.1007/s115380140040x.PubMedPubMed CentralView ArticleGoogle Scholar
 Yousefi MR, Datta A, Dougherty ER. Optimal intervention strategies for therapeutic methods with fixedlength duration of drug effectiveness. Signal Process IEEE Trans. 2012; 60(9):4930–44.View ArticleGoogle Scholar
 Yousefi MR, Dougherty ER. Intervention in gene regulatory networks with maximal phenotype alteration. Bioinformatics. 2013; 29(14):1758–67. doi:10.1093/bioinformatics/btt242. http://bioinformatics.oxfordjournals.org/content/29/14/1758.full.pdf+html.PubMedView ArticleGoogle Scholar
 Yousefi MR, Datta A, Dougherty ER. Optimal intervention in markovian gene regulatory networks with randomlength 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/1687415320146.View ArticleGoogle Scholar
 Zañudo JGT, Albert R. Cell fate reprogramming by control of intracellular network dynamics. PLoS Comput Biol. 2015; 11(4):1004193. doi:10.1371/journal.pcbi.1004193.View ArticleGoogle Scholar
 Qiu Y, Tamura T, Ching WK, Akutsu T. On control of singleton attractors in multiple boolean networks: integer programmingbased method. BMC Syst Biol. 2014; 8 Suppl 1:7. doi:10.1186/175205098S1S7.View ArticleGoogle Scholar
 Shmulevich I, Dougherty ER, Zhang W. Gene perturbation and intervention in probabilistic boolean networks. Bioinformatics. 2002; 18(10):1319–31. doi:10.1093/bioinformatics/18.10.1319. http://bioinformatics.oxfordjournals.org/content/18/10/1319.full.pdf+html.PubMedView ArticleGoogle Scholar
 Xiao Y, Dougherty ER. The impact of function perturbations in boolean networks. Bioinformatics. 2007; 23(10):1265–73. doi:10.1093/bioinformatics/btm093.PubMedView ArticleGoogle Scholar
 VeraLicona P, Bonnet E, Barillot E, Zinovyev A. Ocsana: optimal combinations of interventions from network analysis. Bioinformatics. 2013; 29(12):1571–3. doi:10.1093/bioinformatics/btt195.PubMedView ArticleGoogle Scholar
 Poret A, Boissel JP. An in silico target identification using boolean network attractors: Avoiding pathological phenotypes. Comptes rendus biologies. 2014; 337(12):661–78. doi:10.1016/j.crvi.2014.10.002.PubMedView ArticleGoogle Scholar
 Grayson DR, Stillman ME. Macaulay2, a Software System for Research in Algebraic Geometry. http://www.math.uiuc.edu/Macaulay2/.
 VelizCuba 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/1471210515221.View ArticleGoogle Scholar
 VelizCuba A, Jarrah AS, Laubenbacher R. Polynomial algebra of discrete models in systems biology. Bioinformatics. 2010; 26(13):1637–43. doi:10.1093/bioinformatics/btq240.PubMedView ArticleGoogle Scholar
 Hinkelmann F, Brandon M, Guang B, McNeill R, Blekherman G, VelizCuba A, Laubenbacher R. Adam: analysis of discrete models of biological systems using computer algebra. BMC Bioinforma. 2011; 12:295. doi:10.1186/1471210512295.View ArticleGoogle Scholar
 VelizCuba A. Reduction of Boolean network models. J Theor Biol. 2011; 289:167–72.PubMedView ArticleGoogle Scholar
 Curto C, Itskov V, VelizCuba 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
 VelizCuba 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
 Laubenbacher R, Stigler B. A computational algebra approach to the reverse engineering of gene regulatory networks. J Theor Biol. 2004; 229(4):523–37. doi:10.1016/j.jtbi.2004.04.037.PubMedView 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: SpringerVerlag: 1998.Google Scholar
 Shmulevich I, Dougherty ER, Kim S, Zhang W. Probabilistic boolean networks: a rulebased 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
 Thomas R, D’Ari R. Biological Feedback. Boca Raton: CRC Press; 1990. https://lccn.loc.gov/89023874.Google Scholar
 Murrugarra D, Laubenbacher R. The number of multistate nested canalyzing functions. Physica D: Nonlinear Phenomena. 2012; 241(10):929–38. doi:10.1016/j.physd.2012.02.011.View ArticleGoogle Scholar
 Murrugarra D, Laubenbacher R. Regulatory patterns in molecular interaction networks. J Theor Biol. 2011; 288(0):66–72. doi:10.1016/j.jtbi.2011.08.015.PubMedView ArticleGoogle Scholar
 Ahmad J, Niazi U, Mansoor S, Siddique U, Bibby J. Formal modeling and analysis of the malassociated biological regulatory network: insight into cerebral malaria. PLoS One. 2012; 7(3):33532. doi:10.1371/journal.pone.0033532.View ArticleGoogle Scholar
 Didier G, Remy E, Chaouiya C. Mapping multivalued onto boolean dynamics. J Theor Biol. 2011; 270(1):177–84. doi:10.1016/j.jtbi.2010.09.017.PubMedView ArticleGoogle Scholar
 Murrugarra D, Dimitrova ES. Molecular network control through boolean canalization. EURASIP J Bioinform Syst Biol. 2015; 2015(1):9. doi:10.1186/s1363701500292.PubMedPubMed CentralView ArticleGoogle Scholar
 Li R, Yang M, Chu T. Controllability and observability of boolean networks arising from biology. Chaos. 2015; 25(2):023104. doi:10.1063/1.4907708.PubMedView 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