Attractor Transformation by Impulsive Control in Boolean Control Network

Abstract

Boolean control networks have recently been attracting considerable interests as computational models for genetic regulatory networks. In this paper, we present an approach of impulsive control for attractor transitions in Boolean control networks based on the recent developed matrix semitensor product theory. The reachability of attractors is estimated, and the controller is also obtained. The general derivation proposed here is exemplified with a kind of gene model, which is the protein-nucleic acid interactions network, on numerical simulations.

1. Introduction

Gene regulatory networks (GRNs) have been offering plenty of holistic approaches to biological processes. They can explicitly represent the causality of developmental processes and exactly describe the state set of biological systems [1]. Waddington and Kacser proposed a metaphor that the developmental process of GRNs can be represented by a ball rolling down along a landscape with peaks and valleys, and the steady states, which were called attractors, were found at the bottom of the basins [2]. In cell model, there is a one-to-one relationship between the attractors and the observed phenotypes. This means that different cell types can be characterized by different attractors [3]. The states of a GRN will stay in an attractor, unless it is perturbed by an outside impact [4].

In several studies on GRNs such as genetic organogenesis and diseases, researchers have considered to make the states of GRN transit from one attractor to another one by using control methods [5–7]. It was found that repression of a single RNA binding polypyrimidine tract-binding protein was sufficient to induce transdifferentiation of fibroblasts into functional neurons in [8]. An approach was presented to evaluate drug targets of GRN inference to ovarian cancer in [9]. The previous studies mainly focused on intervening the system to help it transit to the desirable attractors by controlling a (or some) valid genetic locus. Since most of the existing achievements in related fields were obtained based on the experiments, the actual impact of control on the same GRN is uncertain [10]. In brief, to estimate the effectiveness of the controller for the transformation of GRNs from one attractor to another one, still remains an open crucial theoretical problem [11–15].

A Boolean network (BN) is often used as a model for gene regulation which treats genes as binary nodes that are either expressed or unexpressed [4]. In order to manipulate networks, the control of BNs is an important topic. A Boolean control network (BCN) can be considered as a BN with additional binary inputs. BCNs are attracting considerable interests as computational models for GRNs which use the exogenous inputs. BCN has been widely used in yeast cell-cycle [16], Drosophila melanogaster [17], and other kinds of cells.

In this paper, we propose a theoretical method to estimate the effects of a certain impulsive controller in a BCN and solve the appropriate controller by using semitensor product. Compared with the existing methods, based on the results of the experiments, our mathematics-based approach is more accurate and simpler.

The rest of this paper is organized as follows. Section 2 reviews STP and the model of BN. In Section 3, the reachability of an attractor is realized and the controller is also obtained. Section 4 gives our main results of an example. Section 5 is the conclusion.

2. Preliminaries

2.1. Semitensor Product

Semitensor product (STP) of matrices was firstly proposed by Cheng and Dong. It is the algebraic form and the coordinate transformation of BN and BCN. Based on STP, BNs and BCNs can be converted into equivalent algebraic form of some standard discrete-time system [18]. In this paper, STP is denoted by “⋉”.

Definition 1. Assuming there are two matrices 𝐴∈ℝ^𝑚×𝑛 and 𝐵∈ℝ^𝑝×𝑞, the STP of 𝐴 and 𝐵 is 𝐴⋉𝐵=(𝐴⊗𝐼_𝛼/𝑛)(𝐵⊗𝐼_𝛼/𝑝), where 𝛼 is the least common multiple of 𝑛 and 𝑝, “⊗” is the Kronecker product, and 𝐼_𝑘 is the identity matrix.
The STP of matrices makes all the fundamental properties of the conventional matrix product remain true [19]. With STP, Boolean operation can be converted into matrix product. These two logical values, “true” and “false,” are expressed in vector forms as 𝛿¹₂ and 𝛿²₂, where 𝛿^𝑟_𝑛 denotes the 𝑟th column of the identity matrix 𝐼_𝑘. Some fundamental logical functions are identified as 𝑀=[𝛿^𝑖1_𝑛,𝛿^𝑖2_𝑛,…,𝛿^𝑖𝑠_𝑛], which is also briefly expressed as 𝑀=𝛿_𝑛[𝑖₁,𝑖₂,…,𝑖_𝑠]. And the logic relationships are(1)negation: 𝑀_𝑛=𝛿₂[2,1];(2)disjunction: 𝑀_𝑑=𝛿₂[1,1,1,2];(3)conjunction: 𝑀_𝑐=𝛿₂[1,2,2,2];(4)XOR: 𝑀_𝑝=𝛿₂[2,1,1,2].
The above matrices are called the structure matrices.

2.2. Attractor

A BN, which is typically formulated as a directed graph, composed of 𝑛 nodes, whose state indicates whether the gene is switched 1 (on) or 0 (off). The state of each node at time 𝑡+1 is determined by the state of its spatial neighbors at time 𝑡. The system can be described by 𝑥₁(𝑡+1)=𝑓₁𝑥₁(𝑡),𝑥₂(𝑡),…,𝑥_𝑛(𝑡),𝑥₂(𝑡+1)=𝑓₂𝑥₁(𝑡),𝑥₂(𝑡),…,𝑥_𝑛(𝑡),⋮𝑥_𝑛(𝑡+1)=𝑓_𝑛𝑥₁(𝑡),𝑥₂(𝑡),…,𝑥_𝑛(𝑡), (1) where 𝑓_𝑖(𝑖=1,2,…,𝑛) is an 𝑛-ary logical function.

The BN is a globally convergence system. An attractor, called the stable state of system, is in the form of either a single state (fixed point) or a repeating set of states (cycle) [20]. Here, we consider how to find the attractors of (1). According to STP, we define 𝐴(𝑡)=⋉^𝑛_𝑖=1𝑥_𝑛(𝑡). (2) Then 𝐴(𝑡+1)=⋉^𝑛_𝑖=1𝑀_𝑖𝐴(𝑡), (3) where 𝑀_𝑖(𝑖=1⋯𝑛) is the structure matrix. Using the properties of STP, (3) can be converted into an algebraic form as 𝐴(𝑡+1)=𝐿𝐴(𝑡), (4) where 𝐿∈Δ₂^𝑛_×2^𝑛 is called the transition matrix. The state of (1) is uniquely determined by the transition matrix. Each column of 𝐿, which is called state number, represents a state of the BN. The attractor has the following definitions:(1)a state 𝑥(𝑡)∈Δ₂^𝑛 is called a fixed point if 𝐿𝑥(𝑡)=𝑥(𝑡); (2){𝑥(𝑡),𝐿𝑥(𝑡),…,𝐿^𝑘𝑥(𝑡)} is called a cycle with length 𝑘 if 𝐿^𝑘𝑥(𝑡)=𝑥(𝑡) and the elements in the set {𝑥(𝑡),𝐿𝑥(𝑡),…,𝐾^𝑘−1𝑥(𝑡)} are distinct.

Theorem 2. In system (1), the number of length 𝑠 cycles, 𝑁_𝑠, is inductively determined by 𝑁₁(𝐿),𝑁=𝑇𝑟𝑎𝑐𝑒_𝑘=𝑇𝑟𝑎𝑐𝑒𝐿^𝑘−∑_{(𝑘)𝑠∈𝒫}𝑠𝑁_𝑠𝑘,2≤𝑘≤2^𝑛, (5) where 𝒫(𝑘) is the set of proper factors of 𝑘. According to (5), one can find all the attractors in the state space of BN (1) [21].

3. Main Results

In GRNs, much attention focuses on making the whole system transit from one attractor to another by control methods [22, 23]. Since the impulsive controller could destroy the cycle structure of the biological system, it is widely used in GRNs [24]. First, we will judge the reachability of an attractor.

Assume 𝒜₁={𝑥_ℎ∣ℎ=1⋯𝑙₁} and 𝒜₂={𝑦_𝑟∣𝑟=1⋯𝑙₂} are two attractors in system (1). 𝑙₁ and 𝑙₂ are the lengths of 𝒜₁ and 𝒜₂, respectively. We will determine the reachability from 𝒜₁ to 𝒜₂.

Here, we consider a BCN with 𝑚 impulsive inputs at time 𝑡, and it is defined as 𝑥₁(𝑡+1)=𝑓₁𝑥₁(𝑡),…,𝑥_𝑛(𝑡),𝑢₁(𝑡),…,𝑢_𝑚(𝑡),⋮𝑥_𝑛(𝑡+1)=𝑓_𝑛𝑥₁(𝑡),…,𝑥_𝑛(𝑡),𝑢₁(𝑡),…,𝑢_𝑚(𝑡), (6) where 𝑓_𝑖(𝑖=1,2,…,𝑛) is an 𝑛-ary logical function and 𝑢_𝑗(𝑗=1⋯𝑚) is the impulsive input. 𝑢_𝑗 is described as 𝑢_𝑗(𝑡)=input,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘, (7) where input is 0 or 1. When system (6) is at time 𝑡=𝑡_𝑘, according to STP, we define 𝐴(𝑡)=⋉^𝑛_𝑖=1𝑥_𝑖(𝑡)⋉⋉^𝑚_𝑗=1𝑢_𝑗(𝑡), (8) then 𝐴(𝑡+1)=⋉^𝑛_𝑖=1𝑀_𝑖𝐴(𝑡), (9) where 𝑀_𝑖(𝑖=1⋯𝑛) is the structure matrix. So (9) can be converted into an algebraic form as follows: 𝐴(𝑡+1)=𝐿𝐴(𝑡), (10) where 𝐿∈Δ₂^𝑛_×2^𝑛+𝑚 is called the state transition matrix of system (6).

Theorem 3. Consider (6). The transformation from attractor 𝒜₁ to attractor 𝒜₂ is reachable with controllers 𝑢₁(𝑡)⋯𝑢_𝑚(𝑡) if and only if 𝐿𝒜₁∩𝒜₂≠𝜙, (11) where 𝜙 is the null set.

Proof. Since 𝐿 is the linear representation of matrix 𝐿 with the inputs 𝑢₁(𝑡)⋯𝑢_𝑚(𝑡), 𝐿𝒜₁ is the reachable state set of 𝒜₁. The intersection set of 𝐿𝒜₁ and 𝒜₂ is the destination states from 𝒜₁ to 𝒜₂ with controllers 𝑢₁(𝑡)⋯𝑢_𝑚(𝑡).

Definition 4. Assume 𝑦_𝑟 is the destination state for the transition from 𝒜₁ to 𝒜₂ if and only if 𝑦_𝑟∈𝐿𝒜₁∩𝒜₂, (12) then 𝑥_ℎ is the source state if and only if 𝑦_𝑟=𝐿𝑥_ℎ, (13) where 𝑥_ℎ∈𝒜₁ and 𝑦_𝑟∈𝒜₂.
Next, we will find the existence of the controller for attractor transition.

Theorem 5. Consider system (6). Assume 𝐿𝑥_ℎ={𝑒_𝑘∣𝑘=1⋯2^𝑚}; if 𝑦_𝑟=𝑒_𝑘, we have 𝛿^𝑘₂^𝑚. The impulsive controllers are obtained by 𝛿^𝑘₂^𝑚=𝑢₁(𝑡)⋉𝑢₂(𝑡)⋉⋯𝑢_𝑚(𝑡), (14) where 𝑢_𝑗(𝑡)(𝑗=1⋯𝑚) is the impulsive input at time 𝑡.

Proof. Since 𝐿𝑥_ℎ is the destination state from 𝑥_ℎ with impulsive inputs 𝑢₁(𝑡)⋯𝑢_𝑚(𝑡), each 𝑒_𝑘 represents each state which is from 𝑥_ℎ by inputs 𝑢₁(𝑡)⋯𝑢_𝑚(𝑡). Equation (14) is based on the properties of STP.
We can solve the input values for a deterministic target by using (14).

4. Example

In order to illustrate our approach, an example is given in this section. It is an idealized protein-nucleic acid interaction involved in gene regulation model in cells [25]. The type of unit component which we will study is shown as follows: 𝑥₁(𝑡+1)=1+𝑥₃(𝑡)+𝑥₆(𝑡)+𝑥₃(𝑡)𝑥₆(𝑡),𝑥₂(𝑡+1)=𝑥₁(𝑡),𝑥₃(𝑡+1)=𝑥₂(𝑡),𝑥₄(𝑡+1)=1+𝑥₃(𝑡)+𝑥₆(𝑡)+𝑥₃(𝑡)𝑥₆(𝑡),𝑥₅(𝑡+1)=𝑥₄(𝑡),𝑥₆(𝑡+1)=𝑥₅(𝑡), (15) where 𝐴⋅𝐵 represents 𝐴 Conjunction 𝐵 and 𝐴+𝐵 represents the XOR operation between 𝐴 and 𝐵. Based on (2)–(4), the 𝐿 matrix of system (15) is 𝐿 = 𝛿₆₄[37 37 38 38 39 39 40 40 37 1 38 2 39 3 40 4 45 45 46 46 47 47 48 48 45 9 46 10 47 11 48 12 53 53 54 54 55 55 56 56 53 17 54 18 55 19 56 20 61 61 62 62 63 63 64 64 61 25 62 26 63 27 64 28].

Using (5), we obtain that there are the following two attractors in the state space: 𝒜₁=()⟶()⟶(),𝒜194619₂=(1)⟶()⟶()⟶()⟶()⟶()⟶(1).3755642810 (16) The attractors of this system represent different quantities of the generation of a metabolic species.

Next, assume the system is already in attractor 𝒜₁; we want to transit the whole system from 𝒜₁ to 𝒜₂ with some impulsive controllers.

The BCN is expressed as 𝑥₁(𝑡+1)=1+𝑥₃(𝑡)+𝑥₆(𝑡)+𝑥₃(𝑡)𝑥₆(𝑡),𝑥₂(𝑡+1)=𝑥₁(𝑡)+𝑢₁(𝑡),𝑥₃(𝑡+1)=𝑥₂(𝑡),𝑥₄(𝑡+1)=1+𝑥₃(𝑡)+𝑥₆(𝑡)+𝑥₃(𝑡)𝑥₆(𝑡),𝑥₅(𝑡+1)=𝑥₄(𝑡)+𝑢₂(𝑡),𝑥₆(𝑡+1)=𝑥₅(𝑡), (17) where 𝑢₁, 𝑢₂ are controllers.

Step 1. Using (8)–(10), we obtain the matrix 𝐿. Based on the computing of 𝐿, we obtain 𝐿𝒜₁=𝛿₆₄[64624846131719]. According to (11), the intersection of two sets is 𝐿𝒜₁∩𝒜₂=𝛿₆₄[641]≠𝜙. So, the transformation from 𝒜₁ to 𝒜₂ is reachable, and 𝑦_𝑟1=𝛿₆₄[64] and 𝑦_𝑟2=𝛿₆₄[1] are the destination states.

Step 2. Since 𝒜₁ is made up of two states, which are 𝒜₁₁=𝛿₆₄[19] and 𝒜₁₂=𝛿₆₄[46], we have 𝐿𝒜₁₁=𝛿₆₄[64624846],𝐿𝒜₁₂=𝛿₆₄[131719], (18) then 𝐿𝒜₁₁∩𝑦_𝑟1≠𝜙 and 𝐿𝒜₁₂∩𝑦_𝑟2≠𝜙.
So 𝑥_ℎ1=𝒜₁₁=𝛿₆₄[19] and 𝑥_ℎ2=𝒜₁₂=𝛿₆₄[46] are source states.

Step 3. Letting 𝐿𝒜₁₁={𝑒_𝑘∣𝑘=1,…,4}, 𝑦_𝑟1=𝑒₁, we have 𝛿¹₄. Based on Theorem 5, the impulsive controllers are described as 𝑢₁(𝑡)=1,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘,𝑢₂(𝑡)=1,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘. (19) Similarly, we can obtain other impulsive controllers which are described as 𝑢₁(𝑡)=1,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘,𝑢₂(𝑡)=1,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘, (20) whose source state is 𝒜₁₂=𝛿₆₄[46].

The conclusion of this example is that there are two kinds of controllers which can transform the state of system (17) from attractor 𝒜₁ to attractor 𝒜₂. They can be described as follows.

The First. The source state number is 19, the destination state number is 64, and the impulsive controllers are 𝑢₁(𝑡)=1,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘,𝑢₂(𝑡)=1,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘. (21)

The Second. The source state number is 46, destination state number is 1, and the impulsive controllers are 𝑢₁(𝑡)=1,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘,𝑢₂(𝑡)=1,𝑡=𝑡_𝑘,noinput,𝑡≠𝑡_𝑘. (22)

5. Conclusion

This paper explores the problem of attractor transformation by impulsive control in BCN. We propose an effective algorithm which allows us to realize the transformation among different attractors of the BCN. Although the protein-nucleic acid gene network has relatively simple structure compared with those exhibited by metazoans, the attractors transformation by impulsive control is impressively significant. Our findings open a new perspective in the attractor transformation by impulsive control which is of utmost importance in areas as diverse as drug target and gene regulation and so forth. Developing more effective algorithms or approximate techniques for the present approaches will be the future work.

Acknowledgments

This paper is supported by the National Natural Science Foundation of China (Grant nos. 61100204, 61070209, 61121061, and 61272402), the Asia Foresight Program under NSFC Grant (Grant no. 61161140320), and Inner Mongolia Colleges and Universities Scientific and Technological Research Projects (Grant no. NJZY12176).