# Design of Fixed Points in Boolean Networks Using Feedback Vertex Sets and Model Reduction

Authors: Koichi Kobayashi
Journal: Complexity
 (2019)

---

## 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.


---

## Body

## 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 2n×2n 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 Figure1 for an overview of the proposed methods.Figure 1
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 Figure1). Based on this result, a class of BNs such that the number of fixed points becomes 1 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 Section2, 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 x1,x2,…,xn and the index set I={i1,i2,…,im}⊆{1,2,…,n}, define [xi]i∈I≔xi1xi2⋯xim⊤. Let In and 0m×n denote the n×n identity matrix and the m×n zero matrix, respectively. For simplicity of notation, we sometimes use the symbol 0 instead of 0m×n, and the symbol I instead of In.


## 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:(1)x1k+1=f1xjkj∈P1,¬xjkj∈N1,x2k+1=f2xjkj∈P2,¬xjkj∈N2,⋮xnk+1=fnxjkj∈Pn,¬xjkj∈Nn,where xi∈{0,1}, i∈{1,2,…,n} is the state (e.g., the expression of genes) and k∈{0,1,2,…} is the discrete time. The set Pi⊆{1,2,…,n} is a given index set of nodes that affecting the state xi positively. The set Ni⊆{1,2,…,n}∖Pi is a given index set of nodes affecting the state xi negatively. Since Pi∩Ni=∅ holds, we consider only unate functions. The function fi:{0,1}|Pi|×{0,1}|Ni|→{0,1}1 is a given Boolean function satisfying the following assumptions:(i)
fi is minimal (i.e., redundant terms such as xi∨¬xi(=1) are not included),(ii)
Logical operators infi are composed of AND (∧) and OR (∨),(iii)
fi is identical 0 or 1 if Pi=Ni=∅ holds.We present a simple example.Example 1.
Consider the following BN of an apoptosis network:(2)x1k+1=x1k,x2k+1=x1k∧¬x3k,x3k+1=¬x2k∧x4k,x4k+1=x1k∨x3k,where we have the following:x1: the concentration level (high or low) of the tumor necrosis factor (TNF, a stimulus)x2: the concentration level of the inhibitor of apoptosis proteins (IAP)x3: the concentration level of the active caspase 3 (C3a)x4: 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, x3(k)=0 implies cell survival, and x2(k)=0,x3(k)=1 imply cell death [28]. In this BN, the following relations hold:P1=1,N1=∅,P2=1,N2=3,P3=4,N3=2,P4={1,3},N3=∅.


### 2.2. Definitions

Next, some definitions are given. A fixed point and a cyclic attractor are defined as follows.Definition 2.
The statex(k) is called afixed point if x(k+1)=x(k) holds.A fixed point is also called asingleton attractor.Definition 3.
A set of states withp elements, i.e., {x(k),x(k+1),…,x(k+p-1)} (p≥2), is called acyclic attractor if x(k+p)=x(k) holds for all k.We present a simple example.Example 4.
Consider the BN in Example1 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 (0000⊤, 0011⊤, 1011⊤, and 1101⊤) and two cyclic attractors ({0001⊤,0010⊤} and {1001⊤,1111⊤}).Figure 2
State transition diagram, where four fixed points and two cyclic attractors exist.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 graphG=(V,E,L), where V={1,2,…,n} is the set of vertices corresponding to xi, i∈{1,2,…,n}, E={(j,i)∈{1,2,…,n}×{1,2,…,n}∣j∈Pi∪Ni} is the set of arcs, and L:E→{+,-} is the labeling function such that L(j,i)=+ if j∈Pi while L(j,i)=- if j∈Ni.In the next example, an interaction graph is explained.Example 6.
Consider the BN in Example1 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.Figure 3
Interaction graph.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 afeedback vertex set if removal of vertices results in an acyclic graph. In particular, a feedback vertex set is called aminimum 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 aninput vertex if it corresponds to the Boolean function fi satisfying either xi(k+1)=xi(k) (i.e., Pi={i}) or xi(k+1)=0(1) (i.e., Pi=∅). In both cases, Ni=∅ holds.In the case ofxi(k+1)=xi(k), if xi(0) can be arbitrarily determined, then xi(k) can be regarded as a constant input. In the case of xi(k+1)=0(1), even if xi(0) can be arbitrarily determined, xi(k), k≥1 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 anoutput vertex if it corresponds to the Boolean function fi satisfying Pj∩{i}=∅, Nj∩{i}=∅, j≠i.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 Figure3. The minimum vertex set is given by {1,3}. The input vertex is given by node 1. 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 byxi(k+1)=0(1). 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 Example1 again. The input vertex is given by 1 (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 x1(0). In the case of x1(0)=1 (i.e., x1(k)=0 for k≥1), fixed points are given by 1011⊤ and 1101⊤. In the case of x1(0)=0, fixed points are given by 0000⊤ and 0011⊤.Hereafter, under Assumption11, 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.


## 2.1. Boolean Networks

Consider the following BN:(1)x1k+1=f1xjkj∈P1,¬xjkj∈N1,x2k+1=f2xjkj∈P2,¬xjkj∈N2,⋮xnk+1=fnxjkj∈Pn,¬xjkj∈Nn,where xi∈{0,1}, i∈{1,2,…,n} is the state (e.g., the expression of genes) and k∈{0,1,2,…} is the discrete time. The set Pi⊆{1,2,…,n} is a given index set of nodes that affecting the state xi positively. The set Ni⊆{1,2,…,n}∖Pi is a given index set of nodes affecting the state xi negatively. Since Pi∩Ni=∅ holds, we consider only unate functions. The function fi:{0,1}|Pi|×{0,1}|Ni|→{0,1}1 is a given Boolean function satisfying the following assumptions:(i)
fi is minimal (i.e., redundant terms such as xi∨¬xi(=1) are not included),(ii)
Logical operators infi are composed of AND (∧) and OR (∨),(iii)
fi is identical 0 or 1 if Pi=Ni=∅ holds.We present a simple example.Example 1.
Consider the following BN of an apoptosis network:(2)x1k+1=x1k,x2k+1=x1k∧¬x3k,x3k+1=¬x2k∧x4k,x4k+1=x1k∨x3k,where we have the following:x1: the concentration level (high or low) of the tumor necrosis factor (TNF, a stimulus)x2: the concentration level of the inhibitor of apoptosis proteins (IAP)x3: the concentration level of the active caspase 3 (C3a)x4: 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, x3(k)=0 implies cell survival, and x2(k)=0,x3(k)=1 imply cell death [28]. In this BN, the following relations hold:P1=1,N1=∅,P2=1,N2=3,P3=4,N3=2,P4={1,3},N3=∅.


## 2.2. Definitions

Next, some definitions are given. A fixed point and a cyclic attractor are defined as follows.Definition 2.
The statex(k) is called afixed point if x(k+1)=x(k) holds.A fixed point is also called asingleton attractor.Definition 3.
A set of states withp elements, i.e., {x(k),x(k+1),…,x(k+p-1)} (p≥2), is called acyclic attractor if x(k+p)=x(k) holds for all k.We present a simple example.Example 4.
Consider the BN in Example1 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 (0000⊤, 0011⊤, 1011⊤, and 1101⊤) and two cyclic attractors ({0001⊤,0010⊤} and {1001⊤,1111⊤}).Figure 2
State transition diagram, where four fixed points and two cyclic attractors exist.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 graphG=(V,E,L), where V={1,2,…,n} is the set of vertices corresponding to xi, i∈{1,2,…,n}, E={(j,i)∈{1,2,…,n}×{1,2,…,n}∣j∈Pi∪Ni} is the set of arcs, and L:E→{+,-} is the labeling function such that L(j,i)=+ if j∈Pi while L(j,i)=- if j∈Ni.In the next example, an interaction graph is explained.Example 6.
Consider the BN in Example1 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.Figure 3
Interaction graph.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 afeedback vertex set if removal of vertices results in an acyclic graph. In particular, a feedback vertex set is called aminimum 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 aninput vertex if it corresponds to the Boolean function fi satisfying either xi(k+1)=xi(k) (i.e., Pi={i}) or xi(k+1)=0(1) (i.e., Pi=∅). In both cases, Ni=∅ holds.In the case ofxi(k+1)=xi(k), if xi(0) can be arbitrarily determined, then xi(k) can be regarded as a constant input. In the case of xi(k+1)=0(1), even if xi(0) can be arbitrarily determined, xi(k), k≥1 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 anoutput vertex if it corresponds to the Boolean function fi satisfying Pj∩{i}=∅, Nj∩{i}=∅, j≠i.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 Figure3. The minimum vertex set is given by {1,3}. The input vertex is given by node 1. 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 byxi(k+1)=0(1). 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 Example1 again. The input vertex is given by 1 (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 x1(0). In the case of x1(0)=1 (i.e., x1(k)=0 for k≥1), fixed points are given by 1011⊤ and 1101⊤. In the case of x1(0)=0, fixed points are given by 0000⊤ and 0011⊤.Hereafter, under Assumption11, 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

LetW 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 BNsStep 1.
For a given interaction graph, find a minimum feedback vertex set and a set of output vertices. Letx^i and f^i(·),i∈W denote the state and the Boolean function corresponding to vertices obtained, respectively.Step 2.
Replacexj(k+1)=fj(·), j∈{1,2,…,n}∖W with xj(k)=fj(·).Step 3.
Replacexj(k) appearing in f^i(·) with fj(·). Repeat until there exists no xj(k) appearing in f^i(·).Step 4.
Simplify the Boolean functionsf^i,i∈W.We remark here that, under Assumption11, 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 Example1 is considered. In Step 1, we can obtain W={1,3}. In Step 2, we can obtain x2(k)=x^1(k)∧¬x^3(k) and x4(k)=x^1(k)∨x^3(k). In Step 3, we can obtain(3)x^1k+1=x^1k,x^3k+1=¬x^1k∧¬x^3k∧x^1k∨x^3k.Finally, in Step 4, we can obtain(4)x^1k+1=x^1k,x^3k+1=x^3k.Hereafter, the reduced model obtained by the above procedure is denoted by(5)x^1k+1=f^1x^jkj∈P^1,¬x^jkj∈N^1,x^2k+1=f^2x^jkj∈P^2,¬x^jkj∈N^2,⋮x^n^k+1=f^n^x^jkj∈P^n^,¬x^jkj∈N^n^,where n^≔W.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 Example1 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 (00⊤, 01⊤, 10⊤, and 11⊤). 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 α1,α2,…,αnd∈{0,1}n^-m and undesired fixed points β1,β2,…,βnu∈{0,1}n^-m are given, where m 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δ1 and δ2. Then, the following relations hold:
(i)¬δ1 is equivalent to 1-δ1.
(ii)δ1∧δ2 is equivalent to δ1δ2.
(iii)δ1∨δ2 is equivalent to δ1+δ2-δ1δ2.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δj∈{0,1}, j∈J are given, where J is some index set. Then z=∏j∈J‍δj is equivalent to the following linear inequalities:(6)∑j∈Jδj-z≤J-1,-∑j∈Jδj+Jz≤0,where J is the cardinality of J.Using Lemmas18 and 19, the reduced BN (5) can be equivalently rewritten as the following pair of a linear state equation and a linear inequality:(7)x^k+1uk+1=A1A20Ix^kuk+B0zk,Cx^kuk+Dzk≤E,where x^∈{0,1}n^-m is the vector consisting of states corresponding to vertices except for input vertices, and u∈{0,1}m consists of states corresponding to input vertices. The vector z∈{0,1}p is the auxiliary binary variable, and the dimension p 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.
Problem17 is equivalent to the following problem.Problem 21.
Findu∈{0,1}m, z1d,z2d,…,zndd∈{0,1}p, and z1u,z2u,…,znuu∈{0,1}p subject to(8)αi=A1αi+A2u+Bzid,i∈1,2,…,nd,(9)Cαiu+Dzid≤E,i∈1,2,…,ndand(10)βi≠A1βi+A2u+Bziu,i∈1,2,…,nu,(11)Cβiu+Dziu≤E,i∈1,2,…,nu.Finally, consider (10). Noting that the left-hand side value and the right-hand side value are binary, (10) is equivalent to(12)A1-Iβi+A2u+Bziu=δ1-δ2,(13)δ1+δ2=1,where δ1,δ2∈{0,1}. Thus, we can obtain the following theorem.Theorem 22.
Problem17 is equivalent to the following ILP problem.Problem 23.
Findu∈{0,1}m, z1d,z2d,…,zndd∈{0,1}p, z1u,z2u,…,znuu∈{0,1}p, and δ1,δ2∈{0,1} subject to (8), (9), (11), (12), and (13).An ILP problem can be solved by using a free/commercial solver. In this paper, Problem17 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].


## 3.1. Model Reduction of Boolean Networks

LetW 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 BNsStep 1.
For a given interaction graph, find a minimum feedback vertex set and a set of output vertices. Letx^i and f^i(·),i∈W denote the state and the Boolean function corresponding to vertices obtained, respectively.Step 2.
Replacexj(k+1)=fj(·), j∈{1,2,…,n}∖W with xj(k)=fj(·).Step 3.
Replacexj(k) appearing in f^i(·) with fj(·). Repeat until there exists no xj(k) appearing in f^i(·).Step 4.
Simplify the Boolean functionsf^i,i∈W.We remark here that, under Assumption11, 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 Example1 is considered. In Step 1, we can obtain W={1,3}. In Step 2, we can obtain x2(k)=x^1(k)∧¬x^3(k) and x4(k)=x^1(k)∨x^3(k). In Step 3, we can obtain(3)x^1k+1=x^1k,x^3k+1=¬x^1k∧¬x^3k∧x^1k∨x^3k.Finally, in Step 4, we can obtain(4)x^1k+1=x^1k,x^3k+1=x^3k.Hereafter, the reduced model obtained by the above procedure is denoted by(5)x^1k+1=f^1x^jkj∈P^1,¬x^jkj∈N^1,x^2k+1=f^2x^jkj∈P^2,¬x^jkj∈N^2,⋮x^n^k+1=f^n^x^jkj∈P^n^,¬x^jkj∈N^n^,where n^≔W.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 Example1 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 (00⊤, 01⊤, 10⊤, and 11⊤). 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 α1,α2,…,αnd∈{0,1}n^-m and undesired fixed points β1,β2,…,βnu∈{0,1}n^-m are given, where m 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δ1 and δ2. Then, the following relations hold:
(i)¬δ1 is equivalent to 1-δ1.
(ii)δ1∧δ2 is equivalent to δ1δ2.
(iii)δ1∨δ2 is equivalent to δ1+δ2-δ1δ2.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δj∈{0,1}, j∈J are given, where J is some index set. Then z=∏j∈J‍δj is equivalent to the following linear inequalities:(6)∑j∈Jδj-z≤J-1,-∑j∈Jδj+Jz≤0,where J is the cardinality of J.Using Lemmas18 and 19, the reduced BN (5) can be equivalently rewritten as the following pair of a linear state equation and a linear inequality:(7)x^k+1uk+1=A1A20Ix^kuk+B0zk,Cx^kuk+Dzk≤E,where x^∈{0,1}n^-m is the vector consisting of states corresponding to vertices except for input vertices, and u∈{0,1}m consists of states corresponding to input vertices. The vector z∈{0,1}p is the auxiliary binary variable, and the dimension p 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.
Problem17 is equivalent to the following problem.Problem 21.
Findu∈{0,1}m, z1d,z2d,…,zndd∈{0,1}p, and z1u,z2u,…,znuu∈{0,1}p subject to(8)αi=A1αi+A2u+Bzid,i∈1,2,…,nd,(9)Cαiu+Dzid≤E,i∈1,2,…,ndand(10)βi≠A1βi+A2u+Bziu,i∈1,2,…,nu,(11)Cβiu+Dziu≤E,i∈1,2,…,nu.Finally, consider (10). Noting that the left-hand side value and the right-hand side value are binary, (10) is equivalent to(12)A1-Iβi+A2u+Bziu=δ1-δ2,(13)δ1+δ2=1,where δ1,δ2∈{0,1}. Thus, we can obtain the following theorem.Theorem 22.
Problem17 is equivalent to the following ILP problem.Problem 23.
Findu∈{0,1}m, z1d,z2d,…,zndd∈{0,1}p, z1u,z2u,…,znuu∈{0,1}p, and δ1,δ2∈{0,1} subject to (8), (9), (11), (12), and (13).An ILP problem can be solved by using a free/commercial solver. In this paper, Problem17 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 x^i(k+1)=0 or x^i(k+1)=1; 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 x^i(k+1)=x^i(k). 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:(14)x1k+1=x1k,x2k+1=¬x1k,x3k+1=x2k∧x4k,x4k+1=x1k.In the interaction graph of this BN, the input vertex is only 1. The fixed points are derived as 0100⊤ and 1001⊤. Hence, the fixed points can be uniquely specified by using x1(0).


### 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 α∈{0,1}n-m is given, where m 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 x1,x2,…,xn-m, and states corresponding to input vertices are given by xn-m+1,xn-m+2,…,xn. Let αi denote the i-th element of α.The design procedure is given as follows.Procedure of Designing Initial Values of States Corresponding to Input Vertices from an Interaction GraphStep 1.
Enumerate all combinations of initial values ofxn-m+1(k)xn-m+2(k)⋯xn(k)⊤. Let γi, i∈{1,2, …,2m} denote combinations of initial values obtained.Step 2.
Seti=1.Step 3.
Substituteα and γi into the BN (1).Step 4.
If all arguments offi are equal to αi, then a solution of Problem 27 is obtained as γi, and the procedure is terminated. If this condition does not hold, and i<2m, then set i≔i+1 and go to Step 3. If this condition does not hold, and i=2m, 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:(15)x1k+1=¬x2k∧x4k,x2k+1=¬x1k∧x2k,x3k+1=x2k∨x5k,x4k+1=x4k,x5k+1=x5k.Since Boolean functions are not explicitly used in this method, this BN is denoted by(16)x1k+1=f1¬x2k,x4k,x2k+1=f2¬x1k,x2k,x3k+1=f3x2k,x5k,where x4(k+1)=x4(k) and x5(k+1)=x5(k) are omitted. In the interaction graph of this BN, the input vertices are 4 and 5. Suppose that α is given by α=011⊤. First, consider the case of x4(0)=x5(0)=0. In Step 4 of the above procedure, we can obtain(17)0=f10,0,1=f21,1,1=f31,0.Hence, the condition of Step 4 is not satisfied. Next, consider the case of x4(0)=0 and x5(0)=1. In Step 4, we can obtain(18)0=f10,0,1=f21,1,1=f31,1.Hence, the condition of Step 4 is satisfied. Thus, we see that α becomes one of the fixed points by setting x4(0)=0 and x5(0)=1.
The number of combinations ofx4 and x5 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 x4(0)=0 and x5(0)=1. 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 α=100⊤. In the case of x4(0)=1 and x5(0)=0, we can obtain(19)1=f11,1,0=f20,0,0=f30,0.We see that α becomes one of the fixed points by setting x4(0)=1 and x5(0)=0. The conditions x4(0)=1 and x5(0)=0 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.


## 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 x^i(k+1)=0 or x^i(k+1)=1; 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 x^i(k+1)=x^i(k). 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:(14)x1k+1=x1k,x2k+1=¬x1k,x3k+1=x2k∧x4k,x4k+1=x1k.In the interaction graph of this BN, the input vertex is only 1. The fixed points are derived as 0100⊤ and 1001⊤. Hence, the fixed points can be uniquely specified by using x1(0).


## 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 α∈{0,1}n-m is given, where m 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 x1,x2,…,xn-m, and states corresponding to input vertices are given by xn-m+1,xn-m+2,…,xn. Let αi denote the i-th element of α.The design procedure is given as follows.Procedure of Designing Initial Values of States Corresponding to Input Vertices from an Interaction GraphStep 1.
Enumerate all combinations of initial values ofxn-m+1(k)xn-m+2(k)⋯xn(k)⊤. Let γi, i∈{1,2, …,2m} denote combinations of initial values obtained.Step 2.
Seti=1.Step 3.
Substituteα and γi into the BN (1).Step 4.
If all arguments offi are equal to αi, then a solution of Problem 27 is obtained as γi, and the procedure is terminated. If this condition does not hold, and i<2m, then set i≔i+1 and go to Step 3. If this condition does not hold, and i=2m, 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:(15)x1k+1=¬x2k∧x4k,x2k+1=¬x1k∧x2k,x3k+1=x2k∨x5k,x4k+1=x4k,x5k+1=x5k.Since Boolean functions are not explicitly used in this method, this BN is denoted by(16)x1k+1=f1¬x2k,x4k,x2k+1=f2¬x1k,x2k,x3k+1=f3x2k,x5k,where x4(k+1)=x4(k) and x5(k+1)=x5(k) are omitted. In the interaction graph of this BN, the input vertices are 4 and 5. Suppose that α is given by α=011⊤. First, consider the case of x4(0)=x5(0)=0. In Step 4 of the above procedure, we can obtain(17)0=f10,0,1=f21,1,1=f31,0.Hence, the condition of Step 4 is not satisfied. Next, consider the case of x4(0)=0 and x5(0)=1. In Step 4, we can obtain(18)0=f10,0,1=f21,1,1=f31,1.Hence, the condition of Step 4 is satisfied. Thus, we see that α becomes one of the fixed points by setting x4(0)=0 and x5(0)=1.
The number of combinations ofx4 and x5 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 x4(0)=0 and x5(0)=1. 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 α=100⊤. In the case of x4(0)=1 and x5(0)=0, we can obtain(19)1=f11,1,0=f20,0,0=f30,0.We see that α becomes one of the fixed points by setting x4(0)=1 and x5(0)=0. The conditions x4(0)=1 and x5(0)=0 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 Example1, 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 Figure4):(20)x1k+1=¬x8k∧x12k,x2k+1=¬x6k∧¬x9k∧x12k,x3k+1=¬x5k,x4k+1=x3k∧¬x5k,x5k+1=¬x2k∧x4k∧x12k∨¬x2k∨x4k∧¬x12k,x6k+1=x4k∧x12k,x7k+1=x4k∧¬x9k∧x12k∨x4k∨¬x9k∧¬x12k,x8k+1=x4k,x9k+1=¬x7k∧x10k,x10k+1=x1k∨x9k∧¬x11kx11k+1=x4k∧¬x9k∧x12k,∨x4k∨¬x9k∧¬x12k,x12k+1=x12k,Figure 4
Interaction graph of an apoptosis network addressed in [23]. Some edges that do not have a sign imply that both signs are assigned.wherex1,x2,…,x12 are as follows:x1: a second complex (T2),x2: inhibitor of IκB kinases (IKKa),x3: nuclear factor κB in the cytoplam (NFκB),x4: nuclear factor κB in the nucleus (NFκBnuc),x5: inhibitor of NFκB (IκB),x6: a protein regulating IKK activity (A20a),x7: inhibitor of apoptosis proteins (IAP),x8: a protein associated with inhibitation of complex T2 (FLIP),x9: active caspase 3 (C3a),x10: active caspase 8 (C8a),x11: caspase-8 and -10-associated RING proteins (CARP),x12: tumor necrosis factor (TNF).The state corresponding to the input vertex is given byx1. Furthermore, in this model, x7=0 and x9=1 imply cell death, and x7=1 and x9=0 imply cell survival. This BN has two fixed points, i.e., 000010001100⊤ and 000010100010⊤, and has no fixed point for which x12=1 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 x5 by x5(k+1)=(¬x2(k)∨x4(k))∧¬x12(k). The resulting BN has the four fixed points 000010001100⊤, 000010100010⊤, 001101011101⊤, and 001101110011⊤.Consider applying the method in Section3 to this apoptosis model. From Figure 4, we see that one of the minimum vertex sets is given by {4,9,12} (i.e., NFκBnuc, C3a, and TNF). Then, we can obtain the following reduced BN:(21)x^4k+1=x^12k,x^9k+1=¬x^4k∧x^9k∨x^12k∨x^4k∧x^9k∧x^12k,x^12k+1=x^12k.This reduced BN has the four fixed points 000⊤, 010⊤, 101⊤, and 111⊤, and these are in a one-to-one correspondence with the fixed points of the modified BN above.According to [23], x^12(k) (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 settingx^12(0)=0, fixed points can be placed to x^4x^9⊤=00⊤ and x^4x^9⊤=01⊤.(ii)
By settingx^12(0)=1, fixed points can be placed to x^4x^9⊤=10⊤ and x^4x^9⊤=11⊤.Also from this result, we see that it is impossible to choose either cell death or cell survival by controlling onlyx^12(0). These properties hold for the original BN. We suggest that additional control inputs are required.Next, consider solving Problem17 (i.e., the ILP problem). Here, we derive (7) under u(k)=x^12(k). By a simple calculation, we can obtain(22)x^4k+1x^9k+1uk+1=001011001x^4kx^9kuk+001011001z1kz2kz3k,where z1(k)=x^4(k)x^9(k), z2(k)=x^4(k)u(k), and z3(k)=x^9(k)u(k), which can be transformed into a set of linear inequalities with respect to x^4(k), x^9(k), u(k) by using Lemma 19. In Problem 17, we set α1=00⊤ (the desired fixed point) and β1=11⊤ (the undesired fixed point). This setting is artificially given. By solving the ILP problem (Problem 23), we can obtain x^12(0)=0.


## 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 Section4.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.



---

*Source: 9261793-2019-03-04.xml*