Controllability and observability of Boolean networks arising from biology Rui Li, Meng Yang, and Tianguang Chu Citation: Chaos 25, 023104 (2015); doi: 10.1063/1.4907708 View online: http://dx.doi.org/10.1063/1.4907708 View Table of Contents: http://scitation.aip.org/content/aip/journal/chaos/25/2?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Controllability of asynchronous Boolean multiplex control networks Chaos 24, 033108 (2014); 10.1063/1.4887278 Synchronization of coupled large-scale Boolean networks Chaos 24, 013115 (2014); 10.1063/1.4863858 Generalized Boolean Descriptors for Biological Macromolecules AIP Conf. Proc. 963, 552 (2007); 10.1063/1.2836138 Complete synchronizability of chaotic systems: A geometric approach Chaos 13, 495 (2003); 10.1063/1.1566511 An adaptable Boolean net trainable to control a computing robot AIP Conf. Proc. 465, 169 (1999); 10.1063/1.58263

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

CHAOS 25, 023104 (2015)

Controllability and observability of Boolean networks arising from biology Rui Li,1,2 Meng Yang,3 and Tianguang Chu2,a) 1

Key Laboratory of Systems and Control, Institute of Systems Science, Chinese Academy of Sciences, Beijing 100190, China 2 State Key Laboratory for Turbulence and Complex Systems, College of Engineering, Peking University, Beijing 100871, China 3 China Ship Development and Design Center, Wuhan 430064, China

(Received 23 November 2014; accepted 24 January 2015; published online 10 February 2015) Boolean networks are currently receiving considerable attention as a computational scheme for system level analysis and modeling of biological systems. Studying control-related problems in Boolean networks may reveal new insights into the intrinsic control in complex biological systems and enable us to develop strategies for manipulating biological systems using exogenous inputs. This paper considers controllability and observability of Boolean biological networks. We propose a new approach, which draws from the rich theory of symbolic computation, to solve the problems. Consequently, simple necessary and sufficient conditions for reachability, controllability, and observability are obtained, and algorithmic tests for controllability and observability which are based on the Gr€obner basis method are presented. As practical applications, we apply the proposed approach to several different biological systems, namely, the mammalian cell-cycle network, the T-cell activation network, the large granular lymphocyte survival signaling network, and the Drosophila segment polarity network, gaining novel insights into the control and/or monitoring of C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4907708] the specific biological systems. V

A salient problem in systems biology is to develop a control theory for complex and nonlinear biological systems. A number of mathematical models have been proposed to this end, among which perhaps the most attention has been given to the Boolean network (BN) model. Even though BNs seem to be oversimplified models of biological systems, they can provide good approximations to the nonlinear behavior appearing in many real biological networks. From a control theory viewpoint, a BN is controllable if it can be driven from any initial state to any final state within finite time by application of suitable binary inputs. This agrees with our intuitive notion of controlling biological systems, i.e., moving biological systems out of undesirable states (such as those associated with disease) and into desirable ones using appropriate intervention strategies. Controllability analysis is thereby a first and important step towards control synthesis for complex biological systems. On the other hand, an accurate description and control of a Boolean biological network is inherently limited by our ability to assess the network’s state. Although simultaneous measurement of all state variables furnishes a complete description of a biological system, in practice, the measurement is often limited to a subset or a function of state variables. Thus, one needs to infer the state of the whole system from the knowledge of inputs and available measurable outputs, motivating observability analysis of BNs. The present work is devoted to the study of controllability and observability of BNs. The discussion is based on the idea of representing a BN as a polynomial dynamical system. The results obtained provide algorithmic tests for a)

Author to whom correspondence should be addressed. Electronic mail: [email protected].

1054-1500/2015/25(2)/023104/15/$30.00

controllability and observability that can be easily implemented within the computer algebra environments, and are thus believed to scale well to some large biological networks. Using these results, one can readily discover the reachable states of a complex biological system and/ or identify necessary state variables whose measurements allow to reconstruct the system’s complete initial state. Given the fundamental roles, controllability and observability play in biological systems, this work offers an avenue for the development of a practical control theory for complex biological systems.

I. INTRODUCTION

A BN is a discrete-time dynamical system with binary state variables. The state variables interact with each other according to some logical rules, specified through a set of Boolean functions that determine the values of the state variables at the next time step, and thereby define the dynamics of the system. BNs are appealing to any situation, in which the states of system variables can be quantized to only two levels. For example, most biological phenomena are frequently described in a binary logical language such as “on and off,” “upregulated and downregulated,” and “responsive and nonresponsive,” even though they may manifest themselves in the continuous domain.1 Kauffman2 pioneered the use of BNs as simplified models for the complex interaction in the regulatory networks of living cells. Since Kauffman’s original work, BNs have found many applications in a wide spectrum of fields, including cell differentiation,3 immune response,4 evolution,5–7 neural networks,8,9 and gene regulation.10–14 These studies indicate that despite the fact of their

25, 023104-1

C 2015 AIP Publishing LLC V

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-2

Li, Yang, and Chu

conceptual simplicity, BNs retain in most cases meaningful information that can be used to make inferences regarding the systems they model and, hence, provide a nice framework that allows for qualitative and theoretical analysis. Many complex biological systems are known to possess exogenous inputs. One of the main objectives of BN modeling is to use the BN to develop control strategies for biological systems such as therapeutic intervention strategies that shift the state of a genetic regulatory network from an undesirable location to a desirable one. To be consistent with the binary nature of BNs, the control inputs will be allowed to take on only the binary values zero or one. From a biological perspective, a binary input may represent whether a certain treatment such as a medication intervention is applied or not at each time step. Owing to their great biological interest, over the years much effort has been devoted to the development of a control theory for BNs. In Ref. 15, a model checking based algorithm was presented for discovering control sequences that steer a BN to prescribed terminal states in a given time. A method to probabilistically model a BN was proposed in Ref. 16, and some verification and control problems were addressed using the probabilistic model. Cheng et al.17 developed a semi-tensor product approach to deal with BNs. The basic idea underlying this approach is to convert a BN into a discrete-time bilinear system whose state, input, and output variables are canonical vectors. Their work stimulated further research on the analysis and control of BNs; see, for instance, Refs. 18–28. Controllability and observability are among the key concepts of systems theory and are currently major focuses in studying BNs. Roughly speaking, controllability captures the ability to guide a system to a desired state through a suitable choice of input signals. Controllability analysis is thus a preliminary step towards intervention and regulation of biological systems modeled using BNs. On the other hand, observability describes the possibility of uniquely inferring the system’s initial state by measuring its output sequence. In many practical biological systems, experimental access is usually limited to a subset or a function of state variables. Hence, one needs to infer system states from experimentally accessible outputs. In Ref. 29, the problems of reachability, controllability, and observability in BNs were studied using the semi-tensor product tool. There, the observability analysis was addressed under the assumption that the BN is controllable. Using the input-state incidence matrix of BNs, Zhao et al.30 provided a necessary and sufficient condition for controllability and a sufficient condition for observability. In Ref. 31, necessary and sufficient conditions guaranteeing observability were obtained without any particular assumption concerning the controllability of the system. In Ref. 32, observability and reconstructibility properties for BNs were analyzed and two types of observers were proposed for reconstructible BNs. Laschov and Margaliot33 considered the problem of designing control sequences that steer a BN between two given states, while avoiding certain forbidden states, and used the result together with the Perron-Frobenius theory for nonnegative matrices to derive a controllability condition. Output-controllability for a class of

Chaos 25, 023104 (2015)

Boolean biological networks was discussed in Ref. 34. The controllability and observability analysis of BNs with timedelays was carried out in Refs. 35 and 36, respectively. The reachability, controllability, and observability problems for switched BNs were addressed in Refs. 37 and 38. More recently, the analysis was extended to asynchronous Boolean multiplex networks in Ref. 39, where controllability conditions were found. In this paper, we propose a new method to investigate the controllability and observability of BNs. The method is based on the well-known fact that every Boolean function can be represented by a polynomial with coefficients in {0, 1} (see, e.g., Ref. 40). Therefore, a BN with n state variables can be equivalently described as an n-dimensional polynomial system. This mathematical structure enables us to employ techniques and algorithms from symbolic computation. Symbolic computation approaches have proved quite useful for studying polynomial systems (see, e.g., Ref. 41). Recently, Veliz-Cuba42 applied symbolic computation theory to the reverse engineering of genetic regulatory networks and gave an algorithm for constructing BNs from experimental data (see also Ref. 40). We present here new necessary and sufficient conditions for the controllability and observability of a BN. The conditions naturally lead to algorithmic tests for controllability and observability, which exploit the Gr€obner basis algorithm43 and can be implemented in computer algebra systems such as Maple and Mathematica. Computational complexity is one of the current challenges that arise in the study of dynamics and control of BNs. In fact, even the determination of the existence of point attractors for BNs is, in general, strong NP-complete.44 The problems of determining whether a BN is controllable or observable have also been known to be NP-complete.45,46 High computational costs are thus unavoidable for controllability and observability tests unless P ¼ NP. Nevertheless, due to the fact that polynomial algebra computations do not depend so strongly on the size of logical models,47 it is reasonable to believe that the proposed method can scale well to large networks. As will be seen, we demonstrate by example that the method can indeed tackle some real-world biological systems that are almost impossible to deal with using existing methods (such as the semi-tensor product approach). The results thus definitely provide novel and valuable insights into control and monitoring of biological systems. Before concluding this introduction, we would like to mention very briefly some recent related advances in control and observation of continuous-valued networks. Controllability of continuous networks with linear dynamics was discussed in Ref. 48, where a methodology was developed to identify the minimal sets of driver nodes, whose control can guide the system’s entire dynamics. The role of each individual node in controlling the whole network was considered by the same authors in Refs. 49 and 50. An abrupt increase in the control success rate for linear time-invariant networks was identified in Ref. 51 when the number of driver nodes exceeds certain thresholds. Yang et al.52 characterized the emergence of macroscopic observable islands when the number of measurement nodes is increased. A graph-based

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-3

Li, Yang, and Chu

approach to predicting measurement nodes necessary for observability was proposed in Ref. 53. The control profile, which quantifies the different proportions of controlinducing structures present in a network, was studied in Ref. 54, where it was shown that existing synthetic network models fail to generate control profiles that are observed in real-world networks. In Ref. 55, an adaptation of the Barabasi-Albert model that permits tuning of the control profiles was presented. The control of continuous systems through compensatory perturbations was investigated in Refs. 56 and 57. It is worthwhile to point out that such control technique has been successfully translated to Boolean dynamical systems very recently; see Ref. 58. The remainder of this paper is organized as follows. Section II presents some results on polynomial rings and Boolean function rings that will be used later on, and reviews BN models. Sections III and IV analyze the controllability and observability of BNs, respectively. Algorithmic controllability and observability tests are obtained, and several different biological applications are discussed. Section V follows with concluding remarks. To improve readability of the paper, all the proofs are collected in the Appendix. II. PRELIMINARIES AND MODEL A. Mathematical preliminaries

Here, we introduce some notation and some basic algebraic facts needed later in the paper. We denote the ring of polynomials in x1,…, xn with coefficients in a field K by K[x1,…, xn]. Let I  K[x1,…, xn] be an ideal. The affine variety defined by I, denoted V(I), is the set n

fða1 ; …; an Þ 2 K : f ða1 ; …; an Þ ¼ 0 for all f 2 Ig: Given an affine variety V  Kn, we denote by I(V) the set ff 2 K½x1 ; …; xn  : f ða1 ; …; an Þ ¼ 0 for all ða1 ; …; an Þ 2 Vg: Note that I(V) is actually an ideal of K[x1,…,xn]. Given polynomials f1,…,fs 2 K[x1,…,xn], the ideal generated by f1,…,fs is denoted hf1 ; …; fs i. Fix a monomial order. A finite subset G ¼ {g1,…,g‘} of an ideal I is said to be a Gr€obner basis if hLTðg1 Þ; …; LTðg‘ Þi ¼ hLTðIÞi; where LT(gi) is the leading term of gi and hLTðIÞi is the ideal generated by the leading terms of all the elements in I. It is a known fact that43 every nonzero ideal I  K[x1,…, xn] has a Gr€ obner basis and that any Gr€obner basis for an ideal I is a basis of I, i.e., I ¼ hg1 ; …; g‘ i. A reduced Gr€obner basis for a polynomial ideal I is a Gr€obner basis G for I such that every polynomial g 2 G is monic and no monomial of g lies in hLTðG  fggÞi. Given a fixed monomial order on K[x1,…,xn], every ideal I  K[x1,…,xn] other than {0} possesses a unique reduced Gr€obner basis. This yields an ideal equality algorithm for deciding whether two sets of polynomials {f1,…, fs} and {g1,…,g‘} generate the same ideal: simply fix a monomial order and compute the reduced

Chaos 25, 023104 (2015)

Gr€obner basis for hf1 ; …; fs i and hg1 ; …; g‘ i. Then the ideals are equal if and only if the reduced Gr€obner bases are the same. Many computer algebra systems (including Maple and Mathematica) have a built-in command for computing a Gr€obner basis, whose elements are constant multiples of the elements in a reduced Gr€obner basis. Let F2 denote the field of integers modulo 2, and let Bn be the ring of Boolean functions u : Fn2 ! F2 . It is well known that Bn is isomorphic to F2 ½x1 ; …; xn =hx21 þ x1 ; …; x2n þxn i, the quotient ring of F2 ½x1 ; …; xn  by the ideal generated by the relations x21 þ x1 ; …; x2n þ xn (see, e.g., Ref. 59). Namely, every Boolean function u : Fn2 ! F2 can be represented by a polynomial f 2 F2 ½x1 ; …; xn  and the representation is unique provided that the degree of each variable in f is at most one. Thus, Bn is a finite ring. Given a polynomial f 2 F2 ½x1 ; …; xn , we let [f] denote the Boolean function in Bn represented by f. There is a close relation between ideals in Bn and ideals in F2 ½x1 ; …; xn . In fact, the Correspondence Theorem60 tells us that the ideals of F2 ½x1 ; …; xn  that contain hx21 þ x1 ; …; x2n þ xn i are in one-to-one correspondence with the ideals of Bn. If I is an ideal of F2 ½x1 ; …; xn  containing hx21 þ x1 ; …; x2n þ xn i, then this correspondence sends I to pðIÞ ¼ f½f  : f 2 Ig; if J is an ideal of Bn, then the ideal of F2 ½x1 ; …; xn  sent to J under this correspondence is p1 ðJÞ ¼ ff 2 F2 ½x1 ; …; xn  : ½f  2 Jg: Let I and J be ideals in Bn. The sum of I and J, denoted I þ J, is the set fu þ w : u 2 I; w 2 Jg. The product of I and J, denoted IJ, is the set of all finite sums fu1 w1 þ    þ us ws : s 2 N; ui 2 I; wi 2 Jg. Note that I þ J and IJ are ideals in Bn. If J1, J2,…, Js are ideals, we write J1 þ J 2 þ    þ Js ¼

s X

Ji ;

J1 J2    Js ¼

i¼1

s Y

Ji :

i¼1

Given Boolean functions u1 ; …; us 2 Bn , the ideal generated by u1 ; …; us in Bn is denoted hu1 ; …; us i. Since the ring Bn is finite, every ideal J  Bn is finitely generated. That is, there exist u1 ; …; us 2 Bn such that J ¼ hu1 ; …; us i. For any ideal J  Bn, we define VB ðJÞ ¼ fða1 ; …; an Þ 2 Fn2 : uða1 ; …; an Þ ¼ 0 for all u 2 Jg: For each subset W  Fn2 , we define IB ðWÞ ¼ fu 2 Bn : uða1 ; …; an Þ ¼ 0 for all ða1 ; …; an Þ 2 Wg: B. Boolean network model

By representing Boolean functions as polynomials over F2 , a BN with input u and output y can be described by the following equations:

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-4

Li, Yang, and Chu

Chaos 25, 023104 (2015)

0 2m Y c J1c ¼ Ivc ; Jsþ1 ¼@

x1 ðk þ 1Þ ¼ f1 ðx1 ðkÞ; …; xn ðkÞ; u1 ðkÞ; …; um ðkÞÞ; .. . xn ðk þ 1Þ ¼ fn ðx1 ðkÞ; …; xn ðkÞ; u1 ðkÞ; …; um ðkÞÞ; y1 ðkÞ ¼ h1 ðx1 ðkÞ; …; xn ðkÞÞ; .. . yp ðkÞ ¼ hp ðx1 ðkÞ; …; xn ðkÞÞ;

m X

(2)

(1)

Proposition 3.1: Let a; b 2 Fn2 . If ða; bÞ 2 VB ðJsc Þ for some s  1, then b is reachable from a. Proof: See the Appendix. 䊏 Consider the ideals Jsc for s  1. We observe that c for all s. Thus, we get the descending chain of ideals Jsc  Jsþ1 J1c  J2c  J3c     :

As u(k) ranges from (0, 0,…,0) to (1, 1,…,1), v(k) takes on all values from 1 to 2m. Given u ¼ ðc1 ; …; cm Þ 2 Fm 2 , let v be the decimal number corresponding to u, and let ðvÞ

fi ðx1 ; …; xn Þ ¼ fi ðx1 ; …; xn ; c1 ; …; cm Þ 2 F2 ½x1 ; …; xn  for i ¼ 1,…,n. Put ðvÞ

f ðvÞ ¼ ðf1 ; f2 ; …; fnðvÞ Þ: Then we can think of the BN (1) as a Boolean switched system switching between 2m subsystems corresponding to f (1), m f (2),…, f (2 ), with the value of the control u, or equivalently v, determining which subsystem is active at each time step. III. CONTROLLABILITY ANALYSIS OF BOOLEAN NETWORKS

Consider the BN (1). Let x(k) ¼ (x1(k),…, xn(k)) denote the n-dimensional state variable of the BN at time k, and let a; b 2 Fn2 . We say, that b is reachable from a if there exists a positive integer k and a control sequence {u(0), u(1),…,u(k  1)} that steers the BN from x(0) ¼ a to x(k) ¼ b. The BN is controllable if, for every choice of a; b 2 Fn2 , b is reachable from a.17 In other words, controllability means that it is possible to drive the network from an arbitrary initial state to an arbitrary finial state within finite time using external signals. We apply symbolic computation methods to test the controllability of BNs. We first require some notation. For v ¼ 1, 2,…,2m, let f (v) be as in Sec. II B. Given ðv1 Þ ðv2 Þ f ; f ; …; f ðvs Þ 2 ðF2 ½x1 ; …; xn Þn , we use the notation f ðv1 v2 …vs Þ to denote the composition f ðvs Þ 䊊    䊊 f ðv2 Þ 䊊 f ðv1 Þ , ðv v …v Þ and we write fi 1 2 s for the ith component of f ðv1 v2 …vs Þ , ðv v …v Þ ðv …v Þ ðv1 v2 …vs Þ ¼ ðf1 1 2 s ; …; fnðv1 v2 …vs Þ Þ. Consider f1 1 s i.e., f ðv1 …vs Þ þz1 ; …; fn þ zn as polynomials in F2 ½x1 ; …; xn ;z1 ;…;zn , ðv …v Þ and let ½fi 1 s þ zi  denote the Boolean function in B2n repðv …v Þ resented by fi 1 s þ zi . Put ðv …v Þ ðv …v Þ Ivc1 …vs ¼ h½f1 1 s þ z1 ;½f2 1 s þ z2 ;…;½fnðv1 …vs Þ þ zn i  B2n :

We define a sequence of ideals fJsc g1 s¼1 recursively as follows:

(3)

Since the ring B2n is finite, the descending chain (3) must stabilize. That is, there exists a positive integer N such that c c ¼ JNþ2 ¼ : JNc ¼ JNþ1

2mi ui ðkÞ:

i¼1

ðvÞ

Ivc1 …vsþ1 AJsc ; s ¼ 1; 2; 3; …:

v1 ;…;vsþ1 ¼1

v¼1

where xi ; ui ;yi 2 F2 ; fi 2 F2 ½x1 ;…;xn ;u1 ;…;um , and hi 2 F2 ½x1 ; …; xn . Note that at every time step k, the input u(k) ¼ (u1(k), u2(k),…, um(k)) can be equivalently represented using the decimal number vðkÞ ¼ 1 þ

1

2m Y

Write J c ¼ JNc . The following proposition asserts that the chain of ideals (3) is strictly decreasing until its stabilization, thus providing a direct method for determining the ideal Jc. c ¼ Jic for some i  1, then Proposition 3.2: If Jiþ1 c c J ¼ Ji . Proof: See the Appendix. 䊏 We can now derive a necessary and sufficient condition for controllability. Theorem 1. The BN (1) is controllable if and only if J c ¼ f½0g  B2n : Proof: See the Appendix. 䊏 Theorem 1 asserts that the controllability of a BN can be judged by calculating the ideals Jsc . Notice that directly using Eq. (2) to calculate Jsc involves computing the product of 2ms ideals Ivc1 …vs (1  v1,…,vs  2m). Such computational requirements may be formidable, in general. Our next result provides a computationally less expensive method for calculating Jsc , and hence a feasible controllability test can be obtained. Recall that for each integer s  1 the ideal Jsc is finitely generated. Proposition 3.3: Let s be a positive integer. Suppose that Jsc ¼ hu1 ; u2 ; …; u‘ i for some u1 ; u2 ; …; u‘ 2 B2n . For 2n v ¼ 1, 2,…, 2m, define av : F2n 2 ! F2 by av ða1 ; …; an ; b1 ; …; bn Þ ðvÞ

¼ ðf1 ða1 ; …; an Þ; …; fnðvÞ ða1 ; …; an Þ; b1 ; …; bn Þ; and let c I~v ¼ hu1 䊊 av ; u2 䊊 av ; …; u‘ 䊊 av i  B2n :

Then

c Jsþ1

¼

2m Y

! c I~v

J1c :

v¼1

Proof: See the Appendix. 䊏 We are now in a position to deduce the following controllability test algorithm.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-5

Li, Yang, and Chu

Theorem 2 (Controllability test algorithm). The controllability of the BN (1) can be tested in a finite number of steps by the following algorithm: (1) (a) Fix a monomial order. Let s ¼ 1 and compute J1c . (b) Compute the reduced Gr€ obner basis G1 ¼ fg11 ; …; g1‘1 g for p1 ðJ1c Þ  F2 ½x1 ; …; xn ; z1 ; …; zn . (2) (a) Let s ¼ s þ 1. (b) Compute ! 2m Y c Js ¼ h½gs1 1  䊊 av ;…;½gs1‘s1  䊊 av i h½g11 ; …; ½g1‘1 i:

Chaos 25, 023104 (2015)

0 h a ðJ c Þ ¼ @ sþ1 0 ¼@ 0

Gs ¼ fgs1 ; …; gs‘s g for p1 ðJsc Þ  F2 ½x1 ; …; xn ; z1 ; …; zn . (d) If Gs ¼ Gs1, define G ¼ Gs and go to step (3). If Gs 6¼ Gs1, go to step (2a). (3) (a) Find the reduced Gr€ obner basis G~ for hx21 þ x1 ; …; x2n 2 2 þxn ; z1 þ z1 ; …; zn þ zn i. ~ the BN is controllable. If G 6¼ G, ~ the BN is (b) If G ¼ G, not controllable. Proof: See the Appendix. 䊏 A state of a BN is said to be globally reachable if it is reachable from every initial state. Clearly, a BN is controllable if and only if all its states are globally reachable. In the remainder of this section, we consider the problem of determining whether a specific state of a BN is globally reachable. This problem is naturally closely related to the controllability problem and can also find applications in other control theoretic issues for BNs (see, e.g., Ref. 61). The symbolic computation approach is appropriate for studying this problem. Let a ¼ ða1 ; …; an Þ 2 Fn2 . The assignment f ðx1 ; …; xn ; z1 ; …; zn Þ 7! f ðx1 ; …; xn ; a1 ; …; an Þ defines an epimorphism of rings ha : F2 ½x1 ; …; xn ; z1 ; …; zn  ! F2 ½x1 ; …; xn : Since B2n ffi F2 ½x1 ; …; xn ; z1 ; …; zn =hx21 þ x1 ; …; x2n þ xn ; z21 þz1 ; …; z2n þ zn i and Bn ffi F2 ½x1 ; …; xn =hx21 þ x1 ; …; x2n þxn i, ha induces a ring epimorphism h a : B2n ! Bn given by uðx1 ; …; xn ; z1 ; …; zn Þ 7! uðx1 ; …; xn ; a1 ; …; an Þ: For any ideal J  B2n, we let h a ðJÞ denote the set f h a ðuÞ : u 2 B2n g. Then one sees easily that  h a ðJ1c Þ ¼

2m 2m Y Y ðvÞ  hh a ð½fi þ zi Þ : 1  i  ni h a ðIvc Þ ¼ v¼1

¼

2m Y v¼1

and

v¼1 ðvÞ

h½fi

þ ai  : 1  i  ni

h a ðIc A c v1 …vsþ1 Þ h a ðJs Þ

v1 ;…;vsþ1 ¼1 2m Y

1 ðv1 …vsþ1 Þ

hh a ð½fi

þzi Þ : 1  i  niAh a ðJsc Þ

v1 ;…;vsþ1 ¼1 2m Y

1

ðv …v Þ ¼@ h½fi 1 sþ1 þai  : 1  i  niAh a ðJsc Þ v1 ;…;vsþ1 ¼1

for s ¼ 1, 2, 3,…. Since Jsc ¼ J c for all s  N, it follows that the descending chain of ideals h a ðJ c Þ  h a ðJ c Þ  h a ðJ c Þ     1 2 3

v¼1

(c) Find the reduced Gr€ obner basis

1

2m Y

in Bn stabilizes with h a ðJ c Þ. Applying arguments similar to those used in the proofs of Theorems 1 and 2, we can prove the following results for global reachability. Proposition 3.4: Given a BN (1), a state a 2 Fn2 is globally reachable if and only if h a ðJ c Þ ¼ f½0g  Bn : Proposition 3.5: Given a BN (1) and a state a 2 Fn2 , the global reachability of a can be tested in a finite number of steps by the following algorithm: (1) (a) Fix a monomial order. Let s ¼ 1 and compute h a ðJ1c Þ. (b) Compute the reduced Gr€ obner basis G1 ¼ fg11 ; …; g1‘1 g for the ideal p1 ðh a ðJ1c ÞÞ  F2 ½x1 ; …; xn . (2) (a) Let s ¼ s þ 1. (b) Compute h a ðJ c Þ s

2m Y ¼ h½gs1 1 䊊 f ðvÞ ; …; ½gs1 ‘s1 䊊 f ðvÞ i

!

v¼1

 h½g11 ; …; ½g1‘1 i: (c) Find the reduced Gr€ obner basis Gs ¼ fgs1 ; …; gs‘s g for p1 ðh a ðJsc ÞÞ  F2 ½x1 ; …; xn . (d) If Gs ¼ Gs1, define G ¼ Gs and go to step (3). If Gs 6¼ Gs1, go to step (2a). (3) (a) Find the reduced Gr€ obner basis G~ for 2 2 hx1 þ x1 ; …; xn þ xn i. ~ a is globally reachable. If G 6¼ G, ~ a is not (b) If G ¼ G, globally reachable. Example 1 (Mammalian cell-cycle network). The cell cycle is a collection of highly ordered processes that result in the faithful duplication and transmission of genetic information from one cell generation to the next. In a normal dividing cell, the cell cycle is composed of four distinct phases: G1, S, G2, and M phases. During the S phase, cells replicate their DNA. This ultimately leads to the M phase, during which time cells are divided into two daughter cells. The G1 phase is the interval between the end of the previous M phase and

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-6

Li, Yang, and Chu

Chaos 25, 023104 (2015)

This network has 28 ¼ 256 states. Using the symbolic computation approach developed above, we can check the global reachability of these states. For this purpose we introduce x1 ¼ Rb;

x2 ¼ E2F;

x5 ¼ Cdc20;

x3 ¼ CycE;

x6 ¼ Cdh1;

x4 ¼ CycA;

x7 ¼ UbcH10;

x8 ¼ CycB;

u ¼ CycD; and we obtain the BN x1 ðk þ 1Þ ¼ 1 þ x3 ðkÞx4 ðkÞx8 ðkÞuðkÞ; x2 ðk þ 1Þ ¼ x1 ðkÞx4 ðkÞx8 ðkÞ; x3 ðk þ 1Þ ¼ 1 þ x1 ðkÞx2 ðkÞ; x4 ðk þ 1Þ ¼ 1 þ x1 ðkÞx5 ðkÞ þ x1 ðkÞx4 ðkÞx5 ðkÞ þ x1 ðkÞ

x2 ðkÞx4 ðkÞx5 ðkÞ þ x1 ðkÞx5 ðkÞx6 ðkÞx7 ðkÞ þ x1 ðkÞx4 ðkÞx5 ðkÞx6 ðkÞx7 ðkÞ FIG. 1. The mutated mammalian cell-cycle network. Pointed arrows correspond to positive regulatory interactions and blunt arrows to negative regulatory interactions. The control node is marked in red.

the beginning of DNA synthesis. The interval between the end of DNA synthesis and the beginning of mitosis is called the G2 phase. Cells in the G1 phase can exit the cell cycle to settle into a quiescent state, called the G0 phase, for a variable amount of time before reentering the cell cycle. The mammalian cell cycle is tightly regulated by a complex system of factors that either promote cell division cycle activity or arrest cells at G0 phase. Aberrant expression of cell cycle regulators has been implicated in the development of some cancers. Recently, Faryabi et al.62 have proposed a BN model that describes cell cycle progression in a cancerous scenario (see Fig. 1). The model postulates the mammalian cell cycle with a mutated phenotype, in which the cell might cycle in the absence of any extracellular growth factor. The model includes the following genes: CycD, Rb, E2F, CycE, CycA, Cdc20, Cdh1, UbcH10, and CycB. The expression of CycD reflects the state of the extracellular growth factor and thus represents the control input. The Boolean functions of the model are presented in Table I.

TABLE I. Boolean rules for the mutated mammalian cell-cycle network. Reproduced with permission from Faryabi et al., “Optimal constrained stationary intervention in gene regulatory networks,” EURASIP J. Bioinf. Syst. Biol. 2008, 620 767 (2008). Copyright 2008 SpringerOpen. Variable CycD Rb E2F CycE CycA Cdc20 Cdh1 UbcH10 CycB

Boolean rule Input CycD Ù CycE Ù CycA Ù CycB Rb Ù CycA Ù CycB E2F Ù Rb ðE2F Ù Rb Ù Cdc20 Ù ðCdh1 Ù UbcH10ÞÞ ÚðCycA Ù Rb Ù Cdc20 Ù ðCdh1 Ù UbcH10ÞÞ CycB ðCycA Ù CycBÞÚCdc20 ðCdh1 Ù UbcH10 Ù ðCdc20ÚCycAÚCycBÞÞÚCdh1 Cdc20 Ù Cdh1

þ x1 ðkÞx2 ðkÞx4 ðkÞx5 ðkÞx6 ðkÞx7 ðkÞ; x5 ðk þ 1Þ ¼ x8 ðkÞ; x6 ðk þ 1Þ ¼ 1 þ x5 ðkÞ þ x4 ðkÞx5 ðkÞx8 ðkÞ; x7 ðk þ 1Þ ¼ 1 þ x6 ðkÞ þ x6 ðkÞx7 ðkÞ þ x4 ðkÞx5 ðkÞx6 ðkÞx7 ðkÞx8 ðkÞ; x8 ðk þ 1Þ ¼ 1 þ x5 ðkÞ þ x5 ðkÞx6 ðkÞ: Now, we apply the methodology outlined in this section and we conclude that among all the 256 possible states exactly 19 of them are globally reachable. A list is presented in Table II for all the globally reachable states of this model. This means, in particular, that the system is not fully controllable if only the state of CycD is externally manipulated. TABLE II. The 19 globally reachable states of the mutated mammalian cellcycle network. The rows correspond to the states that can be reached from all the other states. Rb 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1

E2F

CycE

CycA

Cdc20

Cdh1

UbcH10

CycB

0 0 0 0 0 0 0 0 1 1 1 1 1 0 1 1 1 1 1

0 0 0 0 0 0 1 1 0 0 0 1 1 0 0 0 0 1 1

0 0 0 0 1 1 0 1 0 0 0 0 1 0 0 0 0 0 1

0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0

0 1 0 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1

0 0 1 1 1 1 1 0 0 1 1 0 0 0 0 1 1 0 0

0 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 1 0 0

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-7

Li, Yang, and Chu

Chaos 25, 023104 (2015)

IV. OBSERVABILITY ANALYSIS OF BOOLEAN NETWORKS

Consider the BN (1). For an initial condition xð0Þ 2 Fn2 , and a control sequence U ¼ {u(0), u(1),…}, with uðiÞ 2 Fm 2, let y(k; x(0), U) denote the output of the network at time k. We say, that two different states a; b 2 Fn2 are indistinguishable if y(k; a, U) ¼ y(k; b, U) for any positive integer k and any control sequence U ¼ {u(0), u(1),…, u(k  1)}. In this case, it is not possible to uniquely infer the network’s initial state by observing the output sequence. We say, that the BN is observable if every two distinct states are distinguishable.17 Note that, in contrast to Refs. 29, 31, and 32, our notion of observability is the more standard concept of multiple-experiment observability:63 every pair of initial states can be distinguished by some input sequence that may depend on the initial states. We will show below how it is possible to use the symbolic computation methods to test the observability of BNs. Again, we first need some notation. Let f ðv1 …vs Þ be as in Sec. III. For i ¼ 1, 2,…, p, we use ðv …v Þ hi 1 s to denote the composition hi 䊊 f ðv1 …vs Þ . Note that ðv …v Þ hi 1 s 2 F2 ½x1 ; …; xn  for each i. Put Ivo1 …vs

ðv …v Þ ðv …v Þ ¼h½h1 1 s ðx1 ; …; xn Þ þ h1 1 s ðz1 ; …; zn Þ; …; ½hpðv1 …vs Þ ðx1 ; …; xn Þ þ hpðv1 …vs Þ ðz1 ; …; zn Þi

 B2n :

We define a sequence of ideals in B2n recursively as follows:

observability directly from the ideal Jo, thereby saving the computation of the variety VB(Jo). Proposition 4.3: If J is an ideal in B2n such that VB ðJÞ ¼ fða; aÞ : a 2 Fn2 g, then J ¼ h½x1 þ z1 ; …; ½xn þ zn i. Proof: See the Appendix. 䊏 As an immediate consequence of Propositions 4.2 and 4.3, we obtain the following criterion for observability. Theorem 3. The BN (1) is observable if and only if J o ¼ h½x1 þ z1 ; …; ½xn þ zn i  B2n : The next result helps us to reduce the computational requirements when calculating the ideals Jso . Proposition 4.4: Let s  0 be an integer. Suppose that Jso ¼ hw1 ; w2 ; …; w‘ i for some w1, w2,…,w‘ 2 B2n. For v ¼ 1, 2n 2,…,2m, define bv : F2n 2 ! F2 by ðvÞ

bv ða1 ;…;an ;b1 ;…;bn Þ ¼ðf1 ða1 ;…;an Þ;…;fnðvÞ ða1 ;…;an Þ; ðvÞ

f1 ðb1 ;…;bn Þ;…;fnðvÞ ðb1 ;…;bn ÞÞ; and let o I~v ¼ hw1 bv ; w2 bv ; …; w‘ bv i  B2n :

Then o Jsþ1

¼

J0o

þ

2m X

o I~v :

v¼1

J0o ¼ h½h1 ðx1 ; …; xn Þ þ h1 ðz1 ; …; zn Þ; …;

(4)

Proof: See the Appendix. 䊏 Now, we can deduce the following method to check the observability of the BN (1). Theorem 4 (Observability test algorithm). The observability of the BN (1) can be tested in a finite number of steps by the following algorithm:

(5)

(1) (a) Fix a monomial order. Let s ¼ 0 and compute J0o . (b) Compute the reduced Gr€ obner basis

½hp ðx1 ; …; xn Þ þ hp ðz1 ; …; zn Þi; o Jsþ1

¼

Jso

þ

2m X

Ivo1 …vsþ1 ;

s ¼ 0; 1; 2; …:

v1 ;…;vsþ1 ¼1 o Observe that for all s  0 we have Jso  Jsþ1 , and so

J0o  J1o  J2o    

forms an ascending chain of ideals. Since B2n is finite, there exists an M  0 such that o o o ¼ JMþ1 ¼ JMþ2 ¼ : JM o . Write J o ¼ JM In fact, the chain (5) is strictly increasing until its stabilization, as the next proposition shows. o ¼ Jio for some i  0, then Proposition 4.1: If Jiþ1 o o J ¼ Ji . Proof: See the Appendix. 䊏 The following result tells us that the observability of a BN can be judged by computing the variety VB(Jo). Proposition 4.2: The BN (1) is observable if and only if VB ðJ o Þ ¼ fða; aÞ : a 2 Fn2 g. Proof: See the Appendix. 䊏 Observe that if J o ¼ h½x1 þ z1 ; …; ½xn þ zn i  B2n , then VB ðJ o Þ ¼ fða; aÞ : a 2 Fn2 g and, hence, the BN is observable by Proposition 4.2. We now show that J ¼ h½x1 þ z1 ; …; ½xn þ zn i is the unique ideal in B2n such that VB ðJÞ ¼ fða; aÞ : a 2 Fn2 g. This will allow us to check the

G0 ¼ fg01 ; …; g0‘0 g for p1 ðJ0o Þ  F2 ½x1 ; …; xn ; z1 ; …; zn . (2) (a) Let s ¼ s þ 1. (b) Compute Jso ¼ h½g01 ; …; ½g0‘0 i þ

2m X h½gs1 1  䊊 bv ; …; ½gs1 ‘s1  䊊 bv i: v¼1

(c) Find the reduced Gr€ obner basis Gs ¼ fgs1 ; …; gs‘s g for p1 ðJso Þ  F2 ½x1 ; …; xn ; z1 ; …; zn . (d) If Gs ¼ Gs1, define G ¼ Gs and go to step (3). If not, go to step (2a). (3) (a) Compute the reduced Gr€ obner basis G~ for 2 2 hx1 þ z1 ; …; xn þ zn ; z1 þ z1 ; …; zn þ zn i. ~ the BN is observable. If G 6¼ G, ~ the BN is (b) If G ¼ G, not observable. Proof: See the Appendix. 䊏

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-8

Li, Yang, and Chu

We now apply our test for observability to some Boolean biological networks. Example 2 (T-cell signaling network). T-cells are a type of white blood cells known as lymphocytes. These white blood cells play a central role in adaptive immunity and enable the immune system to mount specific immune responses. Once activated, T-lymphocytes may act as direct cytotoxic effectors or furnish help for other cells that are essential in the effector phase of the immune response: CD8þ T-cells destroy cells infected by viruses or tumor cells, and CD4þ T-cells assist the functions of other white blood cells in immunologic processes such as B-cells and monocytes. T-cells have the ability to recognize potentially dangerous agents and subsequently initiate an immune reaction against them. They do so by using T-cell receptors (TCRs) to detect foreign antigens bound to major histocompatibility complex (MHC) molecules, and then activate, through a signaling cascade, several transcription factors. These transcription factors, in turn, influence the cell’s fate such as proliferation. In the following, we discuss the application of our approach to two well-developed BN models concerning the signaling pathway machinery of T-cells—the T-cell activation network and the T-cell large granular lymphocyte (T-LGL) survival signaling network.

Chaos 25, 023104 (2015)

(1) T-cell activation network. A BN model for the activation of transcription factors involved in T-cell activation was identified in Ref. 64 (see Fig. 2). This BN includes 3 input variables and 37 state variables. The variables and Boolean functions of the model are given in Table III. The output of the model contains four transcription factors: AP1, CRE, NFAT, and NFjB. Here, to better illustrate the efficiency of the proposed algorithm, we address a more complicated observability problem. We assume that we can monitor a selected subset of state variables and we consider the problem of identifying a minimal subset of state variables such that the BN with these state variables as outputs is observable. In other words, observing the state variables in this subset allows to uniquely determine the network’s complete initial state. It is trivial to see that observability holds when taking all the 37 state variables as outputs. A next possible step is to take 36 out of the 37 state variables as outputs and check observability of the corresponding BN. Note that there are a total of 37 possible cases. Consider the case where the output includes all the state variables except for AP1. By using the procedure given in Theorem 4, we find that the BN is not observable, and so AP1 is needed to obtain observability. Proceeding in this way shows that in order to guarantee observability the output must contain the following set of 14 state variables:

FIG. 2. The T-cell activation network. Pointed arrows correspond to positive regulatory interactions and blunt arrows to negative regulatory interactions. The control nodes are marked in red. The minimal set of nodes, whose measurements allow us to reconstruct the initial states of all other variables, is shown in green.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-9

Li, Yang, and Chu

Chaos 25, 023104 (2015)

TABLE III. Boolean rules for the T-cell activation model. Adapted with permission from Klamt et al., “A methodology for the structural and functional analysis of signaling and regulatory networks,” BMC Bioinf. 7, 56 (2006). Copyright 2006 BioMed Central. Variable

Boolean rule

Variable

Boolean rule

Variable

Boolean rule

CD8 CD45 TCRlig AP1 Ca Calcin cCbl CRE CREB DAG ERK Fos Fyn

Input Input Input Fos Ù Jun IP3 Ca ZAP70 CREB Rsk PLCg(act) MEK ERK ðLck Ù CD45Þ ÚðTCRbind Ù CD45Þ

Gads Grb2Sos IKKbeta IP3 Itk IjB JNK Jun LAT Lck MEK NFAT NFjB PAGCsk

LAT LAT PKCth PLCg(act) SLP76 Ù ZAP70 IKKbeta SEK JNK ZAP70 PAGCsk Ù CD8 Ù CD45 Raf Calcin IjB FynÚTCRbind

PKCth PLCg(act)

DAG ðItk Ù PLCgðbindÞ Ù SLP76 Ù ZAP70Þ ÚðPLCgðbindÞ Ù Rlk Ù SLP76 Ù ZAP70Þ LAT Ras Grb2SosÚRasGRP1 DAG Ù PKCth Lck ERK PKCth Gads cCbl Ù TCRlig FynÚðLck Ù TCRbindÞ cCbl Ù Lck Ù TCRphos

S ¼ fAP1; CRE; Fos; Grb2Sos; Itk; Jun; NFAT; NFjB; PLCgðbindÞ; RasGRP1; Rlk; SLP76; TCRbind; TCRphosg: On the other hand, it is easy to check using Theorem 4 that observing the state variables in S suffices to guarantee observability. Thus, S is the unique minimal set with the desired property (Fig. 2). This implies, in particular, that we cannot have observability when the output includes only the transcription factors AP1, CRE, NFAT, and NFjB. (2) T-LGL survival signaling network. To further attest that the proposed algorithms scale well to large biological systems, we consider the 58-node BN introduced in Ref. 65. It is a Boolean model for the survival of cytotoxic Tlymphocytes in T-LGL leukemia. There are three deregulated nodes without upstream regulation: Stimuli, IL-15, and PDGF. As reported in Ref. 65, the model is capable of predicting all known signaling abnormalities in leukemic T-LGL when IL-15 and PDGF are constantly ON (i.e., taking the value of 1). Following this fact, we maintain the setting of IL-15 and PDGF as ON and regard Stimuli as the control input. By using our observability algorithm, we can efficiently identify the state variables that allow us to infer the state of the whole system. The result gives that aside from the trivial case of monitoring all state variables, observability holds only when the output includes all the state variables except for Apoptosis (a conceptual node serving as an indicator of cell fate, see Ref. 65). This finding suggests a system-theoretic perspective for understanding the intricacy of signaling activities in leukemic T-LGL: It is generally not possible to acquire the complete signaling information by measuring only a part of components in the signaling pathways. Example 3 (Drosophila segment polarity network). It is apparent that our observability results are easily adapted to the special case of a BN system without external inputs. We now consider an input-free Boolean system presented in Ref. 66 that attempts to model the segment polarity gene network in Drosophila (see Table IV). Like all other arthropods, the body of the fruit fly Drosophila melanogaster is

PLCg(bind) Raf Ras RasGRP1 Rlk Rsk SEK SLP76 TCRbind TCRphos ZAP70

composed of repeated units called segments. At the early stages, prior to morphological differentiation, these segments are marked out by a chemical blueprint in a process called determination. Segment determination is controlled by several related classes of genes organized in a hierarchical cascade, with one class being the segment polarity genes. The segment polarity genes represent the last step in the cascade and ultimately specify the pattern of each segment. While most of the preceding genes in the cascade act only transiently, the segment polarity genes are expressed throughout the life of the fly. The attractors of the segment polarity network capture the long-run behavior of the model and may be used to characterize the expression patterns of the segment polarity genes. From a biological point of view, one may be interested in observability at attractor states. The authors of

TABLE IV. Boolean rules for the non-spatial Drosophila model. Reproduced with permission from K. Willadsen and J. Wiles, “Robustness and state-space structure of Boolean gene regulatory models,” J. Theor. Biol. 249, 749–765 (2007). Copyright 2007 Elsevier. Variable

Boolean rule

SLP wg WG en EN hh HH ptc PTC PH SMO ci CI CIA CIR n WG n hh

SLP ðCIA Ù SLP Ù CIRÞÚððCIAÚSLPÞ Ù CIR Ù wgÞ wg n WG Ù SLP en EN Ù CIR hh CIA Ù EN Ù CIR ptcÚðPTC Ù n hhÞ PTC Ù n hh PTCÚn hh EN ci CI Ù ðSMOÚn hhÞ CI Ù ðSMOÚn hhÞ n WG n hh

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-10

Li, Yang, and Chu

Chaos 25, 023104 (2015)

TABLE V. Attractor states of the non-spatial Drosophila model. Reproduced with permission from K. Willadsen and J. Wiles, “Robustness and state-space structure of Boolean gene regulatory models,” J. Theor. Biol. 249, 749–765 (2007). Copyright 2007 Elsevier. Attractor 1 2 3 4 5 6 7 8 9 10

SLP

wg

WG

en

EN

hh

HH

ptc

PTC

PH

SMO

ci

CI

CIA

CIR

n WG

n hh

0 0 0 0 0 0 1 1 1 1

0 0 0 0 0 1 0 1 0 1

0 0 0 0 0 1 0 1 0 1

0 0 1 1 1 0 0 0 0 0

0 0 1 1 1 0 0 0 0 0

0 0 1 1 1 0 0 0 0 0

0 0 1 1 1 0 0 0 0 0

0 1 0 0 0 1 0 1 0 1

1 1 0 0 1 1 1 1 1 1

0 1 0 0 0 1 0 1 0 1

0 1 1 1 0 1 0 1 0 1

1 1 0 0 0 1 1 1 1 1

1 1 0 0 0 1 1 1 1 1

0 1 0 0 0 1 0 1 0 1

1 0 0 0 0 0 1 0 1 0

0 0 1 1 1 0 0 0 1 1

0 1 0 1 0 1 0 1 0 1

Ref. 66 enumerated the state space of the model and found that the system dynamics results in ten different attractors, all of which are point attractors (i.e., cycles of length one, see Table V). Using the methodology developed above, we have determined, for each attractor, the minimal subset of state variables such that the BN with these state variables as outputs is observable at the attractor. The results are summarized in Table VI. It can be seen that to achieve observability some state variables—WG, HH, and PH—always need to be measured, while some others—wg, en, EN, hh, PTC, ci, and CI—never need to be considered. The BN is observable for all the attractors when the output includes the following state variables: SLP, WG, HH, ptc, PH, SMO, CIA, CIR, n WG, and n hh (Fig. 3). Using the proposed observability test algorithm, it is easily shown that this output suffices to yield an observable BN. In other words, the output that guarantees observability for ten attractors achieves observability for all the 217 ¼ 131 072 states.

Typically, the polynomial algebra computations depend heavily on the complexity of the model’s polynomial form: the more terms the polynomial form has, the more time the Gr€obner basis calculation needs. The size of the model, in contrast, is not as much correlated to the calculation time.47 Although the complexity of computing a Gr€obner basis may be, in general, doubly exponential, it is fast, in practice, for polynomial ideals representing biological networks,67 thanks to the fact that most actual biological networks, especially, genetic regulatory networks are sparsely connected. The proposed method thus scales gracefully even to some largescale biological systems. Generally speaking, as is well known in control theory, controllability concerns the existence of a control to complete a desired task but does not need to tell about how to determine such a control. Consequently, a relevant topic

V. CONCLUSION

We have developed a new approach to the controllability and observability analysis of BNs. The approach is based on the idea of representing a BN in polynomial form. As a consequence, we have available to us the well-developed theory of commutative algebra and computational algebraic geometry, with a wide variety of implemented procedures. We have shown how it is possible to use these powerful tools to check controllability and observability of BNs. TABLE VI. Minimum state variables needed to be measured to achieve observability. Attractor 1 2 3 4 5 6 7 8 9 10

Variables SLP, WG, HH, ptc, PH, CIA, CIR WG, HH, PH, SMO WG, HH, PH, SMO, CIA, n hh WG, HH, PH, SMO, CIA, n hh WG, HH, ptc, PH, SMO, CIA SLP, WG, HH, PH, SMO SLP, WG, HH, ptc, PH, CIA, CIR, n WG SLP, WG, HH, PH, SMO, n WG WG, HH, ptc, PH, CIA, CIR, n WG WG, HH, PH, SMO, n WG

FIG. 3. The Drosophila segment polarity network model. Pointed arrows denote positive interactions and blunt arrows denote negative interactions. The minimum set of nodes, that guarantees observability for all the attractors, is shown in green.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-11

Li, Yang, and Chu

raised is that of control design. Analogously, observability guarantees that the outputs contain sufficient information to reconstruct the state of the whole system. In practice, to recover the system state from outputs, it often resorts to the design of state observers, an important and challenging subject in modern control theory. These topics, however, are beyond the scope of this paper, and will be pursued in the future. We finally remark that the dynamic updating scheme employed in this paper is the classical synchronous scheme: the values of all variables are updated simultaneously. This type of update assumes implicitly that the time scales of biological processes in the system are similar and state transitions of components are closely synchronized.68 On the other hand, in situations where time scales of biological events are different and vary over a wide range, asynchronous models offer a possible method of properly accounting for all such variations. Controllability and observability of asynchronous Boolean models (for example, the random order asynchronous models, the general asynchronous models, or the deterministic asynchronous models; see Ref. 69) may also be considered. Such extensions would be another promising direction for future research. ACKNOWLEDGMENTS

We would like to thank the anonymous reviewers for their many helpful comments and suggestions. This work was supported by the National Natural Science Foundation of China (Grant Nos. 61273111 and 61333001), by the National Basic Research Program of China (Grant No. 2012CB821200), and by the Guozhi Xu Posdoctoral Research Foundation.

Chaos 25, 023104 (2015)

Proof: If f 2 p1(J), then [f] ¼ [h1f1 þ    þ hsfs] for P some hi 2 F2 ½x1 ; …; xn , hence f  si¼1 hi fi 2 hx21 þ x1 ; …; x2n þ xn i. Thus, p1 ðJÞ  hf1 ; …; fs ; x21 þ x1 ; …; x2n þ xn i. The other inclusion is obvious. 䊏 Proof of Proposition 3.1: We consider the case s ¼ 1, from which the general case follows easily by induction on s. By Lemma 1(c), VB ðJ1c Þ ¼

(a) (b) (c) (d) (e)

I  J ) VB(J)  VB(I); J  IB(VB(J)); VB(IJ) ¼ VB(I) [ VB(J); VB(I þ J) ¼ VB(I) \ VB(J); if I ¼ hu1 ; …; us i and J ¼ hw1 ; …; wr i, then IJ ¼ hui wj : 1  i  s; 1  j  ri and I þ J ¼ hu1 ; …; us ; w1 ; …; wr i.

Proof: The results follow immediately from the definitions. 䊏 Lemma 2: If I ¼ hf1 ; …; fs i is an ideal in F2 ½x1 ; …; xn  containing hx21 þ x1 ; …; x2n þ xn i, then pðIÞ ¼ h½f1 ; …; ½fs i. Proof: One inclusion is obvious, and to prove the other inclusion pðIÞ  h½f1 ; …; ½fs i, note that for any f 2 I, we have f ¼ h1f1 þ    þ hsfs for some hi 2 F2 ½x1 ; …; xn . Hence, ½f  ¼ ½h1 ½f1  þ    þ ½hs ½fs  2 h½f1 ; …; ½fs i. 䊏 Lemma 3: If J ¼ hu1 ; …; us i  Bn , and if fi 2 F2 ½x1 ; …; xn  represents ui for i ¼ 1,…, s, then p1 ðJÞ ¼ hf1 ; …; fs ; x21 þ x1 ; …; x2n þ xn i:

VB ðIvc Þ:

v¼1

Hence, ða; bÞ 2 VB ðIvc Þ for some 1  v  2m. Write b ¼ (b1, b2,…,bn). Then we have ðvÞ

f1 ðaÞ þ b1 ¼ 0; ðvÞ

f2 ðaÞ þ b2 ¼ 0; .. . fnðvÞ ðaÞ þ bn ¼ 0; 䊏 so that b ¼ f (v)(a). Thus, b is reachable from a. Proof of Theorem 1: (Sufficiency). Since JNc ¼ J c ¼ f½0g, we see that VB ðJNc Þ 2n ¼ F2 . Hence, the BN is controllable by Proposition 3.1. (Necessity). It is clear that Jc  {[0]}. To establish the reverse inclusion, suppose that a; b 2 Fn2 , and write b ¼ (b1,…,bn). Since b is reachable from a, there exists k  1 and 1  v1,…,vk  2m such that f ðv1 …vk Þ ðaÞ ¼ b. Hence, ðv …v Þ fi 1 k ðaÞ þ bi ¼ 0 for i ¼ 1,…, n, so that ða; bÞ 2 VB ðIvc1 …vk Þ  VB ðJkc Þ;

APPENDIX: PROOFS OF MAIN RESULTS

We begin with some technical lemmas, which will be useful to prove the main results. Lemma 1: Let I, J be ideals in Bn. Then

2m [

c c by Lemma 1(c). Thus, VB ðJkc Þ ¼ F2n 2 . Since J  Jk , it fol2n c lows from Lemma 1(a) that VB ðJ Þ ¼ F2 . Hence, Lemma 1(b) implies

J c  IB ðVB ðJ c ÞÞ ¼ IB ðF2n 2 Þ ¼ f½0g: The theorem is proved. 䊏 Proof of Proposition 3.3: Repeated application of Eq. (2) gives Jsc ¼

Y

Ivc11

v11

 Y

  Y  Ivc21 v22    Ivcs1 …vss :

v21 ;v22

(A1)

vs1 ;…;vss

By Lemma 1(e), the right side of Eq. (A1) is equal to *

Y   Y  ðv Þ ðv …v Þ ½fiv 11 þ ziv11     ½fiv s1…vss ss þ zivs1 …vss  : v11

11

vs1 ;…;vss

s1

+

1  i1 ; …; i2m ; …; i1…1 ; …; i2m …2m  n : Then it follows easily that

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-12

c I~v ¼

*

Li, Yang, and Chu

Chaos 25, 023104 (2015)

Y   Y ! ðv Þ ðv …v Þ ½fiv 11 þziv11   ½fiv s1…vss ss þzivs1 …vss  䊊 11

v11

s1

vs1 ;…;vss

+ av :1i1 ;…;i2m ;…;i1…1 ;…;i2m …2m n * Y   Y  ðv Þ ðv …v Þ ½fiv 11 þziv11  䊊 av  ½fiv s1…vss ss þzivs1 …vss  䊊 av ¼ 11

v11

s1

vs1 ;…;vss

¼

c Þ must happen for some s. stabilizes, p1 ðJsc Þ ¼ p1 ðJs1 This implies that Gs ¼ Gs1, so that the algorithm must terminate after a finite number of steps. 䊏 Proof of Proposition 4.2: (Sufficiency). Suppose that there exist distinct states a; b 2 Fn2 which are indistinguishable. Then

+

hi ðaÞ ¼ hi ðbÞ;

:1i1 ;…;i2m ;…;i1…1 ;…;i2m …2m n *

J1c  J2c  J3c    

and for any s > 0

Y   Y  ðvv Þ ðvvs1 …vss Þ ½fiv 11 þziv11   ½fiv …v þz  : i v …v ss s1 ss 11

v11

vs1 ;…;vss

ðv …v Þ

+ (A2)

Hence, ða; bÞ 2 VB ðJ0o Þ;

Note that c Jsþ1 ¼

Ivc11

 Y

v11

¼J1c

ðv …v Þ

hi 1 s ðaÞ ¼ hi 1 s ðbÞ; v1 ; …; vs ¼ 1; 2; …; 2m ; i ¼ 1; 2; …; p:

s1

1i1 ;…;i2m ;…;i1…1 ;…;i2m …2m n :

Y



Ivc21 v22 

v21 ;v22

(Y Y v



Y

Ivcsþ11 …vsþ1sþ1

v22

!) c Ivv : sþ12 …vsþ1sþ1

Y

s ¼ 1; 2; 3; …;



vsþ11 ;…;vsþ1sþ1

  c  Ivv 22

i ¼ 1; 2; …; p;

vsþ12 ;…;vsþ1sþ1

ða; bÞ 2 VB ðIvo1 …vs Þ; v1 ; …; vs ¼ 1; 2; …; 2m :

Since 0 o VB ðJsþ1 Þ¼@

2m \

1 VB ðIvo1 …vsþ1 ÞA \ VB ðJso Þ

v1 ;…;vsþ1 ¼1

(A3) Using Lemma 1(e) and Eq. (A2), we see that for each 1  v  2m Y v22

  c  Ivv 22 *

¼ 

Y

c Ivv sþ1 2 …vsþ1 sþ1



vsþ1 2 ;…;vsþ1 sþ1

Y  ðvv Þ ½fiv 22 þ ziv22     v22

22

Y

vsþ1 2 ;…;vsþ1 sþ1

ðvvsþ1 2 …vsþ1 sþ1 Þ

½fiv

sþ1 2 …vsþ1 sþ1

þzivsþ1 2 …vsþ1 sþ1 



:

c 1  i1 ; …; i2m ; …; i1…1 ; …; i2m …2m  ni ¼ I~v :

(A4)

The result follows from Eqs. (A3) and (A4). 䊏 Proof of Proposition 3.2: It suffices to show that c c c ¼ Jiþ1 . Suppose that Jiþ1 ¼ Jic ¼ hu1 ; …; u‘ i for some Jiþ2 c u1 ; …; u‘ 2 B2n . Then Proposition 3.3 shows that Jiþ2 Q2m c c ¼ ð v¼1 hu1 䊊 av ; …; u‘ 䊊 av iÞJ1 ¼ Jiþ1 . 䊏 Proof of Theorem 2: To prove that the algorithm works, c ¼ h½gs1 1 ; …; first note that for every s  2, we have Js1 ½gs1 ‘s1 i by Lemma 2, and hence, we can compute Jsc as in step (2b). Next, note that if Gs ¼ Gs1 for some s  2, then c c p1 ðJsc Þ ¼ p1 ðJs1 Þ, hence Jsc ¼ Js1 , so that J c ¼ Jsc by Proposition 3.2. Thus, G is the reduced Gr€obner basis of p1(Jc). Since p1 ðf½0gÞ ¼ hx21 þ x1 ; …; x2n þ xn ; z21 þ z1 ; …; z2n þ zn i by Lemma 3, Theorem 1 implies that the BN is controllable ~ if and only if G ¼ G. It remains to show that the algorithm terminates. Since the descending chain of ideals

by Lemma 1(d), a simple induction argument shows that ða; bÞ 2 VB ðJso Þ for all s  0, so that ða; bÞ 2 VB ðJ o Þ ¼ fða; aÞ : a 2 Fn2 g. We have arrived at a contradiction. Therefore, every two distinct states are distinguishable, and hence, the BN is observable. (Necessity). To prove the inclusion VB ðJ o Þ  fða; aÞ : a 2 Fn2 g, suppose on the contrary that there are distinct a; b 2 Fn2 such that (a, b) 2 V B(Jo). Since the BN is observable, there exists a positive integer k and integers ðv …v Þ ðv …v Þ 1  v1,…,vk  2m such that hi 1 k ðaÞ 6¼ hi 1 k ðbÞ for some 1  i  p. Hence, ða; bÞ 62 VB ðIvo1 …vk Þ, so that ða; bÞ 62 VB ðJko Þ. However, since Jko  J o , it follows from Lemma 1(a) that ða; bÞ 2 VB ðJ o Þ  VB ðJko Þ, which is a contradiction. Therefore, VB ðJ o Þ  fða; aÞ : a 2 Fn2 g. The opposite inclusion is obvious, so that 䊏 VB ðJ o Þ ¼ fða; aÞ : a 2 Fn2 g. Proof of Proposition 4.3: We first show that IB ðVB ðh½x1 þ z1 ; …; ½xn þ zn iÞÞ ¼ h½x1 þ z1 ; …; ½xn þ zn i: (A5) The inclusion h½x1 þ z1 ; …; ½xn þ zn i  IB ðVB ðh½x1 þ z1 ; …; ½xn þ zn iÞÞ follows from Lemma 1(b). For the opposite inclusion, let u be an element of IB ðVB ðh½x1 þ z1 ; …; ½xn þ zn iÞÞ, represented by a polynomial g(x1,…,xn, z1,…,zn). Then g 2 Iðfða; aÞ : a 2 Fn2 gÞ: Write I ¼p1 ðh½x1 þz1 ;…;½xn þzn iÞF2 ½x1 ;…;xn ;z1 ;…;zn . By Lemma 3 we have

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.133.226.104 On: Mon, 22 Jun 2015 11:46:08

023104-13

Li, Yang, and Chu

Chaos 25, 023104 (2015)

I ¼ hx1 þ z1 ; …; xn þ zn ; x21 þ x1 ; …; x2n þ xn ; z21

þ

z1 ; …; z2n

þ zn i

¼ hx1 þ z1 ; …; xn þ zn ; z21 þ z1 ; …; z2n þ zn i: Using lexicographic order with

   zn , it is easy to verify that

is a Gr€ obner basis for I (see, e.g., Ref. 43). Then it follows that43 there is a unique r 2 F2 ½x1 ; …; xn ; z1 ; …; zn  with the following two properties:

(b)

No term of r is divisible by LTðx1 þ z1 Þ; …; LTðxn þ zn Þ; LTðz21 þ z1 Þ; …; LTðz2n þ zn Þ. There is q 2 I such that g ¼ q þ r.

Since LT(xi þ zi) ¼ xi and LTðz2i þ zi Þ ¼ z2i , property (a) shows that r is a polynomial only in z1,…,zn and can be written in the form r ¼ k0 þ

n X

q0 ¼ ðg1 þ 1Þ    ðg‘ þ 1Þ þ 1;

x 1    x n z1

G ¼ fx1 þ z1 ; …; xn þ zn ; z21 þ z1 ; …; z2n þ zn g

(a)

To establish the opposite inclusion, note that the Hilbert Basis Theorem implies that p1 ðJÞ ¼ hg1 ; …; g‘ i for some gj 2 F2 ½x1 ; …; xn ; z1 ; …; zn . Define

qi ¼ ðxi þ zi Þq0 ; Then q0 ða1 ; …; an ; b1 ; …; bn Þ ¼ 0

() gj ða1 ; …; an ; b1 ; …; bn Þ ¼ 0 for j ¼ 1; 2; …; ‘ () ða1 ; …; an ; b1 ; …; bn Þ 2 Vðp1 ðJÞÞ: It is easily verified that V (p1(J)) ¼ VB(J), so that q0 ða1 ; …; an ; b1 ; …; bn Þ ¼ 0 () ai ¼ bi for all 1  i  n: Fix 1  i  n. Since qi ða1 ; …; an ; b1 ; …; bn Þ ¼ ðai þ bi Þq0 ða1 ; …; an ; b1 ; …; bn Þ;

X

ki1 …is zi1    zis ; we see that

s¼1 1i1

Controllability and observability of Boolean networks arising from biology.

Boolean networks are currently receiving considerable attention as a computational scheme for system level analysis and modeling of biological systems...
1MB Sizes 0 Downloads 8 Views