Abstract
Fixed points in Boolean networks (BNs) represent cell types or states of cells and are important to decide characteristics of cells. As the control problem on fixed points, it is important to consider the problem of changing fixed points by using external stimuli (i.e., control inputs). In this paper, we propose two methods for designing fixed points. First, a design method using model reduction is proposed. Using the reduced model, the problem of placing fixed points can be rewritten as an integer linear programming problem. Next, we consider the design problem using only the graph structure of a given BN and derive some results. In both methods, a feedback vertex set of a directed graph plays an important role. Finally, a biological example is presented.
1. Introduction
One of the central problems in systems biology is to develop control theory of gene regulatory networks. Control of gene regulatory networks is closely related to therapeutic interventions, which are realized by radiation, chemotherapy, and so on. Towards gene therapy technologies in the future (see, e.g., [1]), establishment of theoretical foundations on modeling, analysis, and control of gene regulatory networks is important. Furthermore, experimental studies on control of gene regulatory networks have been done (see, e.g., [2–5]). Thus, much attention in control of gene regulatory networks has been attracted from both theoretical and practical viewpoints.
In control methods, a mathematical model should be chosen depending on characteristics of dynamics. A gene regulatory network is modeled by discrete dynamical systems (e.g., Boolean networks (BNs) [6]), continuous dynamical systems (e.g., differential/difference equations), or hybrid dynamical systems (e.g., piecewise affine models) (see, e.g., [7] for further details). In the last decade, BNs are widely used as a mathematical model for control methods of gene regulatory networks (see, e.g., [8–13]). In a BN, gene expression is modeled by a binary value (0 or 1), and interactions between genes are modeled by a set of Boolean functions.
While this model is quite simple, it is still useful in developing a control method for gene regulatory networks. A BN has been extended to a probabilistic BN (PBN) and a context-sensitive PBN (see, e.g., [14–17]).
In this paper, we focus on the control problem related to fixed points. Fixed points in BNs represent cell types or states of cells [18] and are important to decide characteristics of cells. Many results on analysis of fixed points have been obtained (see, e.g., [9, 19–22]). Change of fixed points by external stimuli has been analyzed in, e.g., [23]. However, to the best of our knowledge, the problem of designing fixed points by control has not been formally formulated and studied.
Furthermore, the semitensor product method has been widely used in analysis and control of BNs (see, e.g., [10, 11]). This method is powerful, but there is a serious technical issue. That is, matrices with at least size must be manipulated. Hence, it is important to develop the methods for simplifying a given BN.
In this paper, the problem of designing fixed points by external stimuli is studied, and two design methods are proposed. See also Figure 1 for an overview of the proposed methods.
First, the design method using model reduction is proposed. The model reduction method used in this paper is based on the methods in [24, 25]. Based on these methods, we propose a more sophisticated procedure using a minimum feedback vertex set of an interaction graph [20] obtained from a given BN. A feedback vertex set of a graph is a set of vertices whose removal results an acyclic graph. If the number of its elements is minimal, then it is called a minimum feedback vertex set. A feedback vertex set is well known as an important element to characterize the behavior of biological systems modeled by differential equations [26, 27]. In this paper, it is shown that a feedback vertex set is significant in model reduction of BNs. Using the reduced BN, the problem of designing fixed points by external stimuli can be rewritten as an integer linear programming (ILP) problem (see also Figure 1).
Next, the design method using an interaction graph obtained from a given BN is proposed. In this method, Boolean functions are not used, and only graph structure of a given BN is used. First, the number of fixed points is characterized by using a minimum feedback vertex set (see also Figure 1). Based on this result, a class of BNs such that the number of fixed points becomes by external stimuli is clarified. Next, the design procedure using an interaction graph is proposed (see also Figure 1).
This paper is organized as follows. In Section 2, properties of BNs are explained. In Section 3, the design method using model reduction is explained. In Section 4, the design method using only an interaction graph is explained. In Section 5, a biological example on an apoptosis network is presented. In Section 6, we conclude this paper.
Notation. For the numbers and the index set , define . Let and denote the identity matrix and the zero matrix, respectively. For simplicity of notation, we sometimes use the symbol instead of , and the symbol instead of .
2. Preliminaries
In this section, first, the outline of BNs is explained. Next, some definitions are given.
2.1. Boolean Networks
Consider the following BN:where , is the state (e.g., the expression of genes) and is the discrete time. The set is a given index set of nodes that affecting the state positively. The set is a given index set of nodes affecting the state negatively. Since holds, we consider only unate functions. The function is a given Boolean function satisfying the following assumptions:(i) is minimal (i.e., redundant terms such as are not included),(ii)Logical operators in are composed of AND () and OR (),(iii) is identical or if holds.
We present a simple example.
Example 1. Consider the following BN of an apoptosis network:where we have the following: : the concentration level (high or low) of the tumor necrosis factor (TNF, a stimulus) : the concentration level of the inhibitor of apoptosis proteins (IAP) : the concentration level of the active caspase 3 (C3a) : the concentration level of the active caspase 8 (C8a) This model is described in [28] and is a simplified version of an apoptosis network model in [23]. In this model, implies cell survival, and imply cell death [28]. In this BN, the following relations hold: , , ,
2.2. Definitions
Next, some definitions are given. A fixed point and a cyclic attractor are defined as follows.
Definition 2. The state is called a fixed point if holds.
A fixed point is also called a singleton attractor.
Definition 3. A set of states with elements, i.e., (), is called a cyclic attractor if holds for all .
We present a simple example.
Example 4. Consider the BN in Example 1 again. Then, the state transition diagram can be obtained as Figure 2, where each node corresponds to the value of the state. From this figure, we see that there exist four fixed points (, , , and ) and two cyclic attractors ( and ).
An interaction graph obtained from a given BN is defined as follows (see, e.g., [20]).
Definition 5. An interaction graph of BNs is defined by a signed directed graph , where is the set of vertices corresponding to , , is the set of arcs, and is the labeling function such that if while if .
In the next example, an interaction graph is explained.
Example 6. Consider the BN in Example 1 again. Then, the interaction graph can be obtained as Figure 3. From this graph, we see that an interaction graph of a given BN can model interactions between genes.
For a given interaction graph, a feedback vertex set is defined as follows (see, e.g., [29, 30]).
Definition 7. A set of vertices of an interaction graph is called a feedback vertex set if removal of vertices results in an acyclic graph. In particular, a feedback vertex set is called a minimum feedback vertex set if the number of its elements is minimum.
We remark here that, in the above definition, the sign (+ or −) in an interaction graph is ignored. The computational complexity of finding a minimum feedback vertex set is NP-complete [30], but an approximate algorithm of finding it has been developed (see, e.g., [29]).
Finally, in this paper, input vertices and output vertices for a given interaction graph are newly defined as follows.
Definition 8. A vertex of a given interaction graph is called an input vertex if it corresponds to the Boolean function satisfying either (i.e., ) or (i.e., ). In both cases, holds.
In the case of , if can be arbitrarily determined, then can be regarded as a constant input. In the case of , even if can be arbitrarily determined, , cannot be controlled. In gene regulatory networks, an input vertex corresponds to an external stimulus from out of the cell. In other words, the state corresponding to the input vertex is not influenced by other states. For simplicity, we may consider that the state corresponding to the input vertex is regarded as a constant.
Definition 9. A vertex of a given interaction graph is called an output vertex if it corresponds to the Boolean function satisfying , , .
In other words, the state corresponding to the output vertex does not influence other states.
We present a simple example.
Example 10. Consider the interaction graph in Figure 3. The minimum vertex set is given by . The input vertex is given by node . In this interaction graph, there is no output vertex.
Hereafter, the following assumption is made for input vertices.
Assumption 11. There exists no input vertex such that the Boolean function is given by . In addition, for the state corresponding to the input vertex, its initial state can be arbitrarily controlled.
In other words, initial states corresponding to input vertices are regarded as a control input.
Example 12. Consider the BN in Example 1 again. The input vertex is given by (see also Figure 3). The state transition diagram of this BN is shown in Figure 2. From Figure 2, we see that fixed points are changed depending on the value of . In the case of (i.e., for ), fixed points are given by and . In the case of , fixed points are given by and .
Hereafter, under Assumption 11, we consider the problem of designing fixed points using initial values of states corresponding to input vertices. Input vertices correspond to external stimuli from out of cells (see, e.g., Example 1 and Section 5). Furthermore, it is difficult to change the values of external stimuli dynamically. Hence, it is appropriate to consider the design problem under Assumption 11.
3. Design of Fixed Points Using Model Reduction
In this section, we consider designing fixed points using reduced BNs. First, a model reduction method for BNs is explained. Model reduction of BNs has been studied so far (see, e.g., [24, 25]). In this paper, using these results and the notion of minimum feedback vertex sets, we introduce a sophisticated procedure of model reduction. In the procedure below, the number of states in the reduced model can be determined at the first step. Hence, this procedure is visible and simple. Next, we consider the problem of designing fixed points using reduced BNs.
3.1. Model Reduction of Boolean Networks
Let denote a set of vertices consisting of a minimum feedback vertex set and a set of output vertices. The procedure of model reduction introduced in this paper is given as follows.
Procedure of Model Reduction of BNs
Step 1. For a given interaction graph, find a minimum feedback vertex set and a set of output vertices. Let and denote the state and the Boolean function corresponding to vertices obtained, respectively.
Step 2. Replace , with .
Step 3. Replace appearing in with . Repeat until there exists no appearing in .
Step 4. Simplify the Boolean functions .
We remark here that, under Assumption 11, all input vertices are included in a minimum feedback vertex set. In addition, from the definition of minimum feedback vertex sets, it is guaranteed that Step 3 is terminated. We explain this procedure using a simple example.
Example 13. The BN in Example 1 is considered. In Step 1, we can obtain . In Step 2, we can obtain and . In Step 3, we can obtainFinally, in Step 4, we can obtain
Hereafter, the reduced model obtained by the above procedure is denoted bywhere .
For the reduced model obtained, the following theorem has been obtained [25].
Theorem 14. The set of fixed points for the BN (1) and the set of fixed points for the reduced BN (5) are one-to-one correspondence.
We present a simple example.
Example 15. Consider the BN in Example 1 again. From Figure 2, we see that there exist four fixed points. On the other hand, the reduced model for (2) is given by (4). From (4), we see that there exist four fixed points (, , , and ). Thus, we see that fixed points for (2) and fixed points for (4) are one-to-one correspondence.
Furthermore, from the above procedure and the definitions of feedback vertex sets and output vertices, we can obtain the following theorem immediately.
Theorem 16. The minimal number of states for the reduced BN is given by a sum of the number of elements in a minimum feedback vertex set and the number of output vertices.
If it is difficult to find a minimum feedback vertex set; then instead of it, we may use a feedback vertex set. In such a case, the minimality of the reduced model is lost. However, in large-scale networks, there are situations that the reduced model with no minimality is useful.
3.2. Design of Fixed Points Using Reduced Boolean Networks
Using the reduced BN, consider designing fixed points. It is difficult to uniquely determine a fixed point. We consider the following problem.
Problem 17. Consider the reduced BN (5). Suppose that for states corresponding to vertices except for input vertices; desired fixed points and undesired fixed points are given, where is the number of input vertices. Find initial values of the states corresponding to the input vertices such that the reduced BN has desired fixed points and has no undesired fixed points.
Fixed points represent cell types or states of cells [18]. Hence, the above problem is important to characterize the property of cells. Since, in this problem, we focus on only fixed points, we can use the reduced BN.
This problem can be equivalently rewritten as an integer linear programming (ILP) problem. We explain this fact below. First, as a preparation, the following lemma [31] is introduced.
Lemma 18. Consider two binary variables and . Then, the following relations hold:
(i) is equivalent to .
(ii) is equivalent to .
(iii) is equivalent to .
Using this lemma, (5) can be equivalently rewritten as a polynomial system with binary variables.
Furthermore, the following lemma [32] is also introduced.
Lemma 19. Suppose that binary variables , are given, where is some index set. Then is equivalent to the following linear inequalities:where is the cardinality of .
Using Lemmas 18 and 19, the reduced BN (5) can be equivalently rewritten as the following pair of a linear state equation and a linear inequality:where is the vector consisting of states corresponding to vertices except for input vertices, and consists of states corresponding to input vertices. The vector is the auxiliary binary variable, and the dimension can be determined from Boolean functions. See Section 5 for a method to derive (7). See also [33]. Then, we can obtain the following lemma.
Lemma 20. Problem 17 is equivalent to the following problem.
Problem 21. Find , , and subject toand
Finally, consider (10). Noting that the left-hand side value and the right-hand side value are binary, (10) is equivalent towhere . Thus, we can obtain the following theorem.
Theorem 22. Problem 17 is equivalent to the following ILP problem.
Problem 23. Find , , , and subject to (8), (9), (11), (12), and (13).
An ILP problem can be solved by using a free/commercial solver. In this paper, Problem 17 is transformed into an ILP problem, but as another approach, we may utilize an SMT (Satisfiability Modulo Theories) solver such as the Yices SMT solver [34].
4. Design of Fixed Points Using Minimum Feedback Vertex Sets and Interaction Graphs
In the previous subsection, we use Boolean functions of the reduced BN (5). However, there is a possibility that, even if a given BN is reduced, then solving the ILP problem (Problem 23) is hard. For such cases, structural analysis using only the interaction graph of a given BN is effective. In this section, first, we consider characterizing the number of fixed points using minimum feedback vertex sets. Next, we consider a general case using interaction graphs.
4.1. Special Cases
Consider two special cases on a minimum feedback vertex set. First, we can obtain the following theorem.
Theorem 24. Consider the BN (1). Assume that a minimum feedback vertex set is empty. Then, the number of fixed points is 1.
Proof. From the assumption, the states in the Boolean functions in (5) consist of the states corresponding to the output vertices. Since a minimum feedback vertex set is empty, the output vertex has no self-loop. Hence, the Boolean functions in the reduced BN (5) are given by either or ; that is, the fixed point is uniquely determined.
The assumption in this theorem is equivalent to that a given interaction graph is acyclic. The same result has been obtained in [9, 20], but it is here reformulated using a minimum feedback vertex set.
Next, we consider the case where a minimum feedback vertex set is not empty.
Theorem 25. Consider the BN (1). Assume that all elements of a minimum feedback vertex set are given by only input vertices. Then, the fixed point is uniquely determined by specifying the initial values of the states corresponding to the input vertices.
Proof. From the assumption, the Boolean functions in (5) are given by a function of the states corresponding to the input vertices. Furthermore, from Assumption 11, the Boolean functions corresponding to input vertices are given by the form of . Hence, Theorem 25 can be obtained.
From this theorem, we see that in the case where a minimum feedback vertex set is given by a set of input vertices; the fixed point can be controlled by the initial values of the states corresponding to the input vertices.
We present a simple example.
Example 26. Consider the following BN:In the interaction graph of this BN, the input vertex is only . The fixed points are derived as and . Hence, the fixed points can be uniquely specified by using .
4.2. General Case
For two special cases, the number of fixed points was characterized. However, the concrete values of fixed points are not discussed. Here, we consider deriving a solution of the following problem.
Problem 27. Consider the BN (1). Suppose that, for states corresponding to vertices except for input vertices, a desired fixed point is given, where is the number of input vertices. Then, using only the interaction graph of the BN, find initial values of the states corresponding to the input vertices such that the BN has a desired fixed point.
In general, it is necessary to determine from the biological viewpoint. Hereafter, for simplicity of notation, states corresponding to vertices except for input vertices are given by , and states corresponding to input vertices are given by . Let denote the -th element of .
The design procedure is given as follows.
Procedure of Designing Initial Values of States Corresponding to Input Vertices from an Interaction Graph
Step 1. Enumerate all combinations of initial values of . Let , , denote combinations of initial values obtained.
Step 2. Set .
Step 3. Substitute and into the BN (1).
Step 4. If all arguments of are equal to , then a solution of Problem 27 is obtained as , and the procedure is terminated. If this condition does not hold, and , then set and go to Step 3. If this condition does not hold, and , then a solution of Problem 27 cannot be found, and the procedure is terminated.
We remark here that, in this procedure, Boolean functions are not computed, and only information obtained from an interaction graph is used.
We present a simple example.
Example 28. Consider the following BN:Since Boolean functions are not explicitly used in this method, this BN is denoted bywhere and are omitted. In the interaction graph of this BN, the input vertices are and . Suppose that is given by . First, consider the case of . In Step 4 of the above procedure, we can obtainHence, the condition of Step 4 is not satisfied. Next, consider the case of and . In Step 4, we can obtainHence, the condition of Step 4 is satisfied. Thus, we see that becomes one of the fixed points by setting and .
The number of combinations of and is four. For all combinations, is one of the fixed points. However, the condition of Step 4 of the above procedure holds in only the case of and . Since Boolean functions are not used, the above procedure provides us a sufficient condition such that is one of the fixed points.
As the other example, suppose that is given by . In the case of and , we can obtainWe see that becomes one of the fixed points by setting and . The conditions and are the unique condition such that is a fixed point.
Thus, the proposed procedure corresponds to a sufficient condition such that is one of the fixed points, but computation of this procedure is easy. Since the proposed procedure is a sufficient condition, there is a possibility that even if is one of the fixed points, the proposed procedure cannot be found this fact. It is one of the future efforts to derive a better sufficient condition based on structural analysis.
5. Biological Example
In Example 1, we used a very simple model of an apoptosis network. In this section, a more detailed model [23] is used.
Consider the following BN expressing an apoptosis network (see also Figure 4):
where are as follows: : a second complex (T2), : inhibitor of IB kinases (IKKa), : nuclear factor B in the cytoplam (NFB), : nuclear factor B in the nucleus (NF), : inhibitor of NFB (IB), : a protein regulating IKK activity (A20a), : inhibitor of apoptosis proteins (IAP), : a protein associated with inhibitation of complex T2 (FLIP), : active caspase 3 (C3a), : active caspase 8 (C8a), : caspase-8 and -10-associated RING proteins (CARP), : tumor necrosis factor (TNF).
The state corresponding to the input vertex is given by . Furthermore, in this model, and imply cell death, and and imply cell survival. This BN has two fixed points, i.e., and , and has no fixed point for which holds. We remark that in [23], the asynchronous BN dynamics are considered, and its fixed points are different to those in the synchronous BN dynamics studied in this paper. Based on observations made for the simplified model in Example 1, we next consider modifying the dynamics of by . The resulting BN has the four fixed points , , , and .
Consider applying the method in Section 3 to this apoptosis model. From Figure 4, we see that one of the minimum vertex sets is given by (i.e., NF, C3a, and TNF). Then, we can obtain the following reduced BN:This reduced BN has the four fixed points , , , and , and these are in a one-to-one correspondence with the fixed points of the modified BN above.
According to [23], (TNF) is regarded as a constant control input of the whole system. The above network is composed of a proapoptotic and an antiapoptotic pathway, which are activated by stimulation of death receptor by a factor such as TNF. Binding of TNF to a death receptor activates the antiapoptotic pathway and, after a certain delay, the proapoptotic pathway is also activated. In the above network, we cannot choose either cell death or cell survival by controlling only TNF. The choice of being either cell death or survival depends on also the initial value of other states. Since this BN is too simple, we can obtain the following facts from the truth table:(i)By setting , fixed points can be placed to and .(ii)By setting , fixed points can be placed to and .
Also from this result, we see that it is impossible to choose either cell death or cell survival by controlling only . These properties hold for the original BN. We suggest that additional control inputs are required.
Next, consider solving Problem 17 (i.e., the ILP problem). Here, we derive (7) under . By a simple calculation, we can obtainwhere , , and , which can be transformed into a set of linear inequalities with respect to , , by using Lemma 19. In Problem 17, we set (the desired fixed point) and (the undesired fixed point). This setting is artificially given. By solving the ILP problem (Problem 23), we can obtain .
6. Conclusion
In this paper, we studied the design problem of fixed points in BNs. First, the design problem using model reduction was considered. This problem is rewritten as an ILP problem. Next, the design problem using only interaction graphs was considered. Finally, a biological example is presented.
There are several open problems.
In the proposed methods, a minimum feedback set of the interaction graph obtained from a given BN plays an important role. We will consider that this fact does not depend on mathematical models and is the design principle in biological systems. However, further discussion from several viewpoints is required.
Since we consider only a simple BN, one of the future efforts is to apply the proposed methods to several biological systems. Reducing conservativeness of the procedure in Section 4.2 using additional information is one of the future efforts. In this paper, we focus on only fixed points, and we do not consider cyclic attractors. Analysis of cyclic attractors is also one of the future efforts.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was partly supported by JSPS KAKENHI Grant no. 17K06486.