Matrix Analysis of Synchronous Boolean Networks

The synchronous Boolean network (SBN) is a simple and powerful model for describing, analyzing, and simulating cellular biological networks. This paper seeks a complete understanding of the dynamics of such a model by employing a matrix method that relies on relating the network transition matrix to its function matrix via a self-inverse state matrix. A recursive ordering of the underlying basis vector leads to a simple recursive expression of this state matrix. Hence, the transition matrix is computed via multiplication of binary matrices over the simplest finite (Galois) field, namely the binary field GF(2), i.e., conventional matrix multiplication involving modulo-2 addition, or XOR addition. We demonstrate the conceptual simplicity and practical utility of our approach via an illustrative example, in which the transition matrix is readily obtained, and subsequently utilized (via its powers, characteristic equation, minimal equation, 1-eigenvectors, and 0-eigenvectors) to correctly predict both the transient behavior and the cyclic behavior of the network. Our matrix approach for computing the transition matrix is superior to the approach of scalar equations, which demands cumbersome manipulations and might fail to predict the exact network behavior. Our approach produces result that exactly replicate those obtained by methods employing the semi-tensor product (STP) of matrices, but achieves that without sophisticated ambiguity or unwarranted redundancy. Key wordsSynchronous Boolean networks, Transition matrix, Function matrix, Recursive ordering, Self-inverse state matrix, Galois field GF (2), Characteristic equation, Minimal equation, 1-eigenvectors, 0-eigenvectors.


Introduction
Many biological systems, such as cellular networks and genetic regulatory networks, are too complex to allow exact or nearly exact analysis. Therefore, these systems are usually studied in terms of synchronous Boolean networks (SBNs) as approximating models. In fact, an SBN is the simplest possible conceptual model that mimics or captures the essential features of the original systems, and that can be effectively used for describing, analyzing, and simulating these systems. An SBN is a set of nodes, each of which is either in state 1 (On) or state 0 (Off) at any given time t. Each node is updated at time (t+1) by inputs from any fixed subset of the set of nodes according to any desired logical rule. Since the total number of possible network states is finite ( 2 ) and the network changes states sequentially in discrete time steps, the network must necessarily return to a previous state in a finite time of at most 2 time points. All possible trajectories of the network are either cycles (loops or attractors) of any length from size 1 (a fixed point) to a maximum length of 2 , or transient states leading eventually to a cycle (attractor).
A total description of an SBN is typically achieved by solving a matrix equation, or by algorithms in which a matrix equation is implicit. A simpler description is possible when such a matrix equation is replaced by a scalar equation or a reduced scalar equation (Rushdi and Al-Otaibi, 2007;Rushdi, 2015). Unfortunately, the scalar-equation technique suffers from several shortcomings and limitations (Rushdi, 2015). Recently, the dynamics of SBNs have been explored by a novel matrix method utilizing a new matrix product, called the semi-tensor product (STP) of matrices (Cheng, et al., 2011;Li and Chu, 2012;Rushdi and Ghaleb, 2016;Cheng, 2019). The STP approach uses a certain matrix expression of logic for the numerical derivation of the transition matrix [ ] of the Boolean network. This matrix is then analyzed to deduce full information about the transient and cyclic behavior of the network. It is noted that the concepts of STP are quite involved, and though they lead to useful computational algorithms, they are not easy to grasp, and suffer from unwarranted duplications.
The purpose of this paper is to derive the transition matrix [ ], independently from any STP formulations or derivations, through the multiplication of binary matrices over GF(2), i.e., conventional multiplication associated with modulo-2 addition, or XOR addition. To attain this purpose, the paper utilizes well-known concepts relating the network transition matrix [ ] to its where [ ] is the state matrix of order n and dimensions 2 ×2 (Cull, 1971;Rushdi and Al-Otaibi, 2007 The organization of the remainder of this paper is as follows. Section 2 presents the methodology adopted herein. It starts by reviewing pertinent quantities and relations. Then, it presents our algorithm for computing the transition matrix [ ]. Finally, it explains how [ ] can be used to deduce the cyclic and transient behavior of the network without resorting to a construction of the network full state diagram. Section 3 presents our results in the form of an illustrative example of a synchronous Boolean network previously studied in the literature. The example not only demonstrates the binary matrix multiplications over GF(2) needed in constructing [ ] and its powers, but it also illustrates mathematical techniques employed in deriving necessary related entities. Section 4 discusses our findings via this example and asserts that our algorithm replicates in clear and succinct steps the results obtained by the STP methodology. Our transition-matrix predications of network transient and cyclic behavior are in full agreement with those obtained via a construction of the network state diagram. Section 5 concludes the paper.

Methodology
The methodology adopted herein is one of rigorous mathematical analysis and further demonstration by way of examples. In subsection 2.1, we present preliminary definitions and, in particular, identify four different forms of the basis vectors that are used in the representation of two-valued Boolean functions. These are the lexicographic basis (Cull, 1971), the recursive basis (Rushdi and Al-Otaibi, 2007), the conventional truth-table basis (Rushdi, 2018;Rushdi and Rushdi, 2018), and the STP basis (Cheng, 2019;Rushdi and Ghaleb, 2016). Subsection 2.2 outlines our algorithm for constructing the transition matrix [ ], while subsection 2.3 reviews how [ ] could be utilized in predicting network behavior via its powers, characteristic equation, minimal equation, 1-eigenvectors and 0-eigenvectors.

Preliminary Definitions
Functions used herein are defined over the simplest finite (Galois) field GF (2). The field has only two elements: the additive identity (0) and the multiplicative identity (1). The field addition (+) and multiplication (×) are defined as a modulo-2 addition and a 1-bit multiplication, which resemble the exclusive-OR (⊕) and the ANDing (∧) operations, respectively, in two-valued Boolean algebra. Any function of variables over GF(2) has a linear or Boolean-ring representation via a polynomial of 2 terms (called a Taylor, a Zhegalkin or a Reed-Müller polynomial), or via a vector of length 2 comprising binary coefficients of this polynomial. In the sequel, we interchangeably use two isomorphic representations: the one over GF(2) and that of the Boolean ring, and in particular, our (+) sign means modulo-2 addition or XOR operation, and our (×) sign (to be omitted and indicated simply by juxtaposition) depicts 1-bit multiplication or AND operation. A vector of 2 terms constitutes a basis for all functions of variables over GF(2) or equivalently all switching functions of variables, in the sense that any of these functions equals the (scalar) dot product of the vector of binary coefficients representing it with the afore-mentioned basis vector. Cull (1971) proposed the use of a lexicographic order for the basis vector: In fact, Cull replaced the leading 1 in (1) with a 0; an oversight since 0 is not a logical product at all, while 1 is the multiplicative identity (product of 0 variables). Rushdi and Al-Otaibi (2007) proposed a recursive order: which can be defined by the following recursive relation and boundary condition: An implication of the recursive relation in (3) is that none of the logical products (that are elements of ⃗⃗ ) is succeeded by products that are subsumed by it. We recall that a product subsumes another if the set of literals of the former is a superset of that of the latter (Rushdi and Rushdi, 2017). In the sequel, we will always implicitly use the recursive ordering in (3). We stress the existence of other ways of defining the basis vector that are in common use. A prominent one is the common truth-table or minterm order that can be given recursively as where ⨂ is the Kronecker product for two matrices. For = 3, this basis vector is: Clearly, if the elements of ⃗⃗⃗ are understood to depict binary numbers for 1 = 2 = 3 = 1, then the values of these numbers follow the natural order from 0 to (2 − 1). Another common basis used extensively in the literature is the one of expressing logical (switching) variables via semi-tensor products (STP) of matrices (Cheng, 2019;Rushdi and Ghaleb, 2016). Here, the basis is simply the STP of the pertinent variables: where the symbol ⋉ denotes the STP of matrices, and ⃗ = [ ̅ ] is a vector representing the variable . This corresponds to a recursive definition and an explicit one (for = 3) of the form: The values of the numbers corresponding to the elements of ⃗ ⃗⃗ follow a reverse order from (2 − 1) down to 0 . This reverse order results from defining ⃗ as [ ̅ ] and not as [ ̅ ] . Next, we review some pertinent terminology from Rushdi and Al-Otaibi (2007), and Rushdi (2015) concerning an SBN of nodes. Note that we use two distinct state vectors ⃗⃗ ( ) and ⃗ ( ). The STP community (see, e.g.,  typically uses only one state vector that is named ⃗ ( ) therein which corresponds to our ⃗⃗ ( ).

The Exact State Vector ⃗ ⃗⃗ ( ):
A binary column vector of dimension 2 , whose elements are all 0 except one element that is of value 1, which corresponds to the position where the product depicting the network state occurs in the basis vector ⃗⃗ given in (2) or (3).

The Cumulative State Vector ⃗⃗⃗ ( ):
A binary column vector of dimension 2 , which equals the basis vector ⃗⃗ evaluated at the network state at instant . Hence, ⃗ ( ) has 1 ′ in the positions representing any of the products of ⃗⃗ equated to 1, i.e., the products subsumed by the product depicting the network state. While the number of 1 ′ in the vector ⃗⃗ ( ) is 1, this number in the vector ⃗ ( ) is 2 , where = 0, 1, … , ( ⃗⃗ ( ) ≤ ⃗ ( )).

The Transition Matrix [ ]:
A 2 × 2 matrix with a single 1 and with 0's otherwise in each column (there may be more than one element of 1 in any row, and as many rows with all 0's).The transition matrix [ ] represents the state transition diagram of the network. The particular structure of [ ] asserts that the network is in only one state at a time instance and proceeds to a unique next state in the next instant. The transition matrix [ ] relates the exact next-state vector ⃗⃗ ( + 1) to the exact present-state one ⃗⃗ ( ) via: If the single 1 in column of [ ] is element , this means that state is the next state of state , abbreviated as: where ⃗ 2 is a column vector of length 2 whose elements are all 0 ′ with the exception of a single 1-element in position . The dynamics of an SBN is uniquely determined by the matrix equation (9), provided the initial state ⃗⃗ (0) is specified. The structure matrix [ ] in Cheng et al. (2011) is the same as [ ], but in disguise emanating from the differences in bases. The basis for [ ] in Cull (1971) is ⃗ given by (1) and the basis in Rushdi and AL-Otaibi (2007) is ⃗⃗ given by (2) or (3), while it is ⃗ ⃗⃗ given by (6-8) for [ ].

The Function Matrix [ ]:
A 2 × 2 matrix that has as its rows the vector representations of the 2 products of the scalar functions computed by the elements of the SBN (taken at a time (0 ≤ ≤ )). The function matrix relates two consecutive instances of the vector ⃗ as: where [0 −1 ] is the zero matrix of dimensions 2 −1 × 2 −1 . This recursive definition of [ ] is reminiscent of the transformation between the linear (Reed-Müller) representation and the minterm (truth-table) representation of a switching function (Rushdi and Ghaleb, 2016). The matrix [ ] is an involutory (self-inverse) matrix, i.e., where [ ] is the identity matrix of 2 × 2 elements. The state matrix [ ] is used to connect the transition matrix [ ] and the function matrix [ ] via the following similarity transformations (Cull, 1971;Rushdi, 2015): The similarity in (16) and (17)

Algorithm for Constructing the Transition Matrix [ ]
 Represent the current and next state values ( ) and ( + 1) of the ′ ℎ variable by and , respectively. Express next-state functions ( 1 , 2 , ⋯ , ), 1 ≤ ≤ , in linear form, i.e., in a Reed-Müller representation. Hence, construct their products, taken 0 at a time ( 0 = 1), taken one at a time ( , 1 ≤ ≤ ), taken two at a time ( , 1 ≤ < ≤ ), taken three at a time ( , 1 ≤ < < ≤ ), and so on till the overall product (∏ =1 ) is reached. Express each product in linear form.  Arrange the columns of the present states ( ) and the rows of the next-state functions ( ), both according to the recursive basis in (2) or (3), and construct the function matrix [ ] accordingly. Figure 1 illustrates such a construction for = 3.

Cyclic and Transient behavior
We review some of the mathematics utilized in the study of the cyclic and transient behaviors of an SBN. Most of the results cited herein were reported by Cull (1971), enhanced, elaborated, and proved by Rushdi and Al-Otaibi (2007), and utilized by . A state of the network ⃗⃗ ( ) is called cyclic if: and, otherwise, it is called transient. A set of cyclic states that map into one another are called a cycle (attractor). An example of a cycle of length is the set given by where is the total number of transient states and the subscripted ′ are the lengths of the various cycles. We note that the eigenvalues are the roots of (21), and hence the only possible eigenvalues in GF(2) are = 1 (corresponding to one-eigenvectors that represent cycles) and = 0 (corresponding to zero-eigenvectors that represent transient chains). The factorization of the binary characteristic equation is not always unique, so that one might not be able to specify all cycles exactly. Since matrices [ ] and [ ] are similar, they have the same characteristic equation: Here, 0 is the length of the longest transient chain and ℎ , ⋯ , are lengths of distinct cycles such that all other cycle lengths divide one of the ′ with no remainder. The cyclic behavior of an SBN can also be obtained from the ⃗ -type 1-eigenvectors of the function matrix [ ] or the ⃗⃗type 1-eigenvectors of the transition matrix [ ]. In fact, the cycles of the network are linearly independent 1-eigenvectors of [ ], and the 1-eigenvectors of [ ] are cycles or the unions of cycles (Cull, 1971;Rushdi and Al-Otaibi, 2007;. Likewise, the 0-eeigenvectors can be used in the study of the transient chains. There are as many independent 0-eigenvectors as there are chains. Each of these eigenvectors is of the form ⃗⃗ + ⃗⃗ , where ⃗⃗ ≠ ⃗⃗ . At least one of the states ⃗⃗ and ⃗⃗ is the last state of a transient chain, i.e., a state that is not on a cycle, but whose next state is on a cycle. Cheng et al. (2011) proved that the number 1 of fixed points for the network equals the trace (sum of diagonal elements) of [ ], i.e., 1 = ([ ]). They also asserted that the number of states on cycles of length which is a divisor of an integer is: and hence the number of cycles of length (1 ≤ ≤ 2 ) is given inductively by:
Note that relations (32) imply (31), but the converse is not true. Now, we note that the characteristic equation (29) To finalize this example, we indicate that all our findings can be easily verified. Figure 2 shows a Karnaugh-map representation for the three next-state functions 1 , 2 , 3 . For convenience, the current value of each map cell is shown in a small corner within the cell. From Figure 2, we can directly draw the complete state-diagram of the network in Figure 3. To check the number of is obviously 7 (the sum of all rows is 0, and there is no other linear dependency). This means that the number of cycles is (8 − 7) = 1.

Discussion
Starting from the scalar next-state functions of the network, we write its function matrix [ ] in terms of the various products of these functions. We then obtain the transition matrix [ ] utilizing its similarity with [ ], i.e., by pre-and post-multiplying [ ] with the state matrix [ ]. This step is self-checking, since it must produce a matrix [ ] whose column vectors are binary unit vectors. Since the implicit basis vector for the column vectors of [ ] represents the current states, while that of its row vectors represents the next state, the matrix [ ] can be used to automatically draw the full state diagram of the network. However, we have chosen in our example to utilize properties of [ ] to make direct deductions about network behavior without drawing the state diagram. We utilized the powers of [ ], or alternatively used its characteristic equation. The predictions made by the various methods are self-consistent and are in full agreement with results obtained earlier, notably via the STP machinery. These predictions also coincide exactly with those deduced from the full state diagram.

Conclusions
We utilized well-known results in the seminal work of Cull (1971), together with the improvement added by Rushdi and Al-Otaibi (2007), to present a method for constructing the transition matrix [ ] of an SBN. We then combined several approaches for inferring network behavior from mathematical properties of [ ] without drawing the full state diagram. Our work leads to a better understanding of the relation between the matrix and scalar approaches for studying SBNs. While our matrix method is algorithmic in nature and complete in scope, the scalar techniques use ad-hoc cumbersome heuristics, and might fail to predict the exact network behavior, especially its transient one. This paper is a continuation of earlier matrix methods, but it distinguishes itself by formulating a well-formed algorithm for constructing the matrix [ ], and by amalgamating various techniques for predicting network behavior from the mathematical properties of [ ]. We attempt to make the matrix analysis of SBNs independent of any particular implementation, as we deliberately avoid the use of semi-tensor products (STP) as our matrix approach. Many papers utilize STP as if it were the only tool for formulating matrix equations for Boolean networks. Admittedly, the STP is a powerful computational tool, but it is inefficient and not easy to follow. Further comparison between our matrix method and the STP one is warranted.