Finite-state discrete-time Markov chain models of gene regulatory networks

In this study Markov chain models of gene regulatory networks (GRN) are developed. These models gives the ability to apply the well known theory and tools of Markov chains to GRN analysis. We introduce a new kind of the finite graph of the interactions called the combinatorial net that formally represent a GRN and the transition graphs constructed from interaction graphs. System dynamics are defined as a random walk on the transition graph that is some Markovian chain. A novel concurrent updating scheme (evolution rule) is developed to determine transitions in a transition graph. Our scheme is based on the firing of a random set of non-steady state vertices of a combinatorial net. We demonstrate that this novel scheme gives an advance in the modeling of the asynchronicity. Also we proof the theorem that the combinatorial nets with this updating scheme can asynchronously compute a maximal independent sets of graphs. As proof of concept, we present here a number of simple combinatorial models: a discrete model of auto-repression, a bi-stable switch, the Elowitz repressilator, a self-activation and show that this models exhibit well known properties.


Introduction
Efforts to study gene expression regulation networks has led to a detailed description many of them, and many more are to be identified in the near future. Therefore there is a need to develop methods of computational and theoretical analysis of GRNs. One of the most promising directions is to reduce the problem to the study of Markov chains, generated in some way from the GRN [1,2,3,4,5]. Usually Boolean networks [6] are used as a formal representation of GRN. Classification of process states, the study of long-term behavior [7], and development of optimal strategies for therapeutic intervention [8,9,10,7,11,12,13,14,15] provide good examples of this approach [16]. In contrast to the Boolean network, the Hopfield networks are defined using arithmetic operations [17]. It is a well-developed branch of science which deals with stochastic processes of asynchronous state switching as a result of interactions. As such they are similar to Boolean networks. A Hopfield like formalism also leads to the definition of the Markov chain. In the Hopfield network area essential results were obtained in the study of various update schemes [18], network oscillations [19], solving of combinatorial optimization problems [20,21,22,23,19,24,25,26], and estimating the rate of the convergence and many others. This makes it valuable to study the possibility of using Hopfield like networks for the construction of Markov chains from GRNs and other interactions graphs. We consider a GRN as a kind of interaction graph. Interaction (regulatory) graphs are emerged in various fields of the life science [27] . Nowadays, their transition graphs are often used to analyze properties of interactions (regulations). One promising way to understand the nature of the regulation or interactions represented by interaction graphs is to analyze the Markovian chain associated with their transition graphs.

Method
The proposed method may be viewed as a version of the Hopfield like network [17] where groups of randomly selected unstable units are updated in parallel [18].

The interactions graph and non steady state vertices
Let G = (V, E) be a directed graph, where V is set of vertices and E is set of edges. Let B = {0, 1} be a set of vertex states.We say that the vertex is active if the state of this vertex is equal to "1", otherwise we say that vertex is inactive. The map function M : V → B give us the state of each vertex. If some vertex v ∈ V , then M (v) is the state of the vertex v that correspond to map M . M (v) = 1 is equivalent to the vertex v is active under map M. M (v) = 0 is equivalent to the vertex v is inactive. The weight function W : E → R gives the value of each edge of the graph G,which represent the power of interactions. If e = (u, v) than we say that u is a direct ancestor of v and v is a direct descendant of u. The inf luence on v under the map M is defined as the sum of weights of edges from all direct active ancestors of vertex v. The influence on v under the map M is denoted by I(v, M ) (also called "the local field or the net input"). That is By definition, if forced state and current state of v are the same then state of the vertex v under the map M is steady. Next equation give us the formal definition:

The Random Set update rule
Now we consider a stochastic process{Y j , j = 1, 2, 3, . . . } that takes on the set of maps of some interaction graph G, where Y j denote the map of G at time period j. At each time period j for each non steady-state vertex v i under the map Y j we change the current vertex state to a forced state with the probability p i , and the current state stays unchanged with the probability 1 − p i . Let S = {v 1 , v 2 , . . . , v n } be a set of all non steady-state vertices at the time period j. Vertices chosen to change their state in one-step transition form a random set X ⊆ S. To make the new map that direct accessible from current map Y j all vertices from X simultaneously change their state whereas another vertices stay unchanged. Let P={p 1 , p 2 , . . . , p n } be some vector of numbers such that 0 ≤ p i ≤ 1, we refer to p i as the probability of the state changing (firing) of the non steady-state vertex v i . For any X ⊆ S let 1 X : S → B be indicator function Hence we assume that each v i acts independently to make random set X then the production of P X give us p X that is the firing probability of the random set X Evidently,

X⊆S v∈S
Now we can apply this definition of Random Set update rule and its probabilities to define the transitions graph of the combinatorial net models.

The transitions graph of the interactions graph
Let S be the set of non steady state vertices under map M = Y j . S represent a full set of vertices each of which can flip to a forced state at next j+1 time step. In the combinatorial model steady state vertices can not flip. To construct a transition graph we should define the full set of maps M 1 , M 2 , . . . that are direct reachable from the map M . Each pair (M, M i ) will correspond to one edge in transition graph. What is the set of maps M 1 , M 2 , . . ., which can be direct (by one-step transition) reached from map M ? To represent the independence and the asynchronicity we assume that any random set of non steady state vertices X ⊆ S can produce next map from current map. Let M X be map such that ∈ X, then we say that M X produced by X from M . That is,the random set X of non steady state vertices produce the map M X from a map M . The weight of edge from map M to a map M X are given by the probability defined by equation (4).

The random walk network dynamics
We suppose that whenever the process is in state M , there is a probability p X such that at next step it will be in state M X . This probability is defined for each random set X of non-steady state vertices of map M .

Random Set update rules are more general
It is well known that asynchronous and Random Set update rules are equivalent in the sense of global stable states [28]. But in the sense of the reachability of one state from another they are not equivalent. Figure 1 show the Mace combinatorial model that illustrate this fact .The vertex e provides constant level of repression for vertex c, that is equal -2. Let the vertex d of the Mace will be active at the start. Then it can activate both middle vertices a and b. Due to repression, the vertex c of the Mace can be activated only if both middle vertices a and b will be active simultaneously. Asynchronous (one at a time) updating exclude simultaneous activation of these vertices, but Random Set update rule do not. Synchronous update rule do not exclude simultaneous activation of a and b, but it make the system deterministic. The Random Set update rule is more general than both synchronous and asynchronous update rules, because it allow all possible paths of a system evolution. Therefore, transition graphs of both synchronous and asynchronous update rules are subgraphs of Random Set update graph.

Examples of combinatorial models
Next we use described above method to develop models of some important graphs of repressive interactions of self activating nodes and prove their main properties. Such models we call combinatorial models.

The combinatorial model of an Auto-repression
A negative auto-regulation or an auto-repression occurs when products of some gene represses its own gene. This form of a simple regulation serve as a basic building block of most important transcription networks [27,29]. An Auto-repression can produce oscillations. For example embryonic stem cells fluctuate between high and low Nanog expression and the Nanog activity is auto-repressive [30]. Our model of an auto-repression shown in Figure 3a and  [2], F (v 1 , M 0 ) = 0, then state of the vertex v 1 is non steady under the map M 0 . Let p 1 = 1/2 then the state of the vertex v 1 could be changed to active state with probability 1/2 and could stay unchanged with probability 1/2. The another state of the model are also non steady. Therefore, there is only non steady states in the model. Thus, it will oscillate infinitely between 0 and 1. Figure 3a show full transition graph of the auto-repression model.

The combinatorial model of a Bi-stable Switch
A Bi-stable Switch is bi-stable gene regulatory network that constructed from a two mutually repressive genes [31]. They are widely common in nature and most used in synthetic biology [32,33]. ODEs used to construct their mathematical models are convenient way for analyzing in detail some small circuits.
Here we develop techniques that can be used to construct models of large networks of bis-table switches and to prove some their important properties. Thus, we use probabilistic coarse-scale modeling approach [34] instead of fine-scale ODE modeling. Our model of Bis-table Switch shown in Fig 2b and Fig 3b  exhibit two steady states. Figure 2b present

The combinatorial model of the Elowitz Repressilator
Elowitz repressilator consists of three genes [35]. Each of these genes are constitutively expressed. The first gene inhibits the transcription of the second gene, whose protein product in turn inhibits the expression of a third gene, and finally, third gene inhibits first gene expression, completing the cycle. Such a negative feedback loop lead to oscillations. The combinatorial model of Elowitz repressilator produce oscillations and consists of three vertices and three edge, which weights is equal to -1. Figure 2c present Table  1.

The combinatorial model of the self-activation
A constitutively expressed gene represent an example of a self-activation. Such gene do not require any interaction to be active. A combinatorial model of the self-activation consist of one vertex and no edges. In any case the influence on it equals 0 since there are no other vertices. Therefore a forced state of the vertex equals 1. Thus 1 is steady state and 0 is non steady state. Vertex starting in steady state will stay in it infinitely. Vertex starting in non steady state with probability p flip to steady state and with probability 1-p stay in non steady state. The amount of time periods T which the vertex spent in non steady state is the random variable. The distribution of this random variable is the shifted geometric distribution with parameter p. Figures 2d and 3d represents the combinatorial model of a self-activation.

The network of bi-stable switches
An independent set (IS) in a graph is a set of vertices no two of which are adjacent . An independent set is called maximal(MIS) if there is no independent set that it contains properly. Hopfield network whose stable states are exactly maximal independent sets was developed by Shrivastava [36]. An independent set in a graph is a clique in the complement graph and vice versa. Thus, cliques can be used to find or to enumerate MISs [20,21]. Finding independent sets (or cliques) has applications in various fields [37]. Combinatorial nets can be used to compute maximal independent sets of graphs in a distributed self-organization fashion. A stable states of a network of bi-stable switches derived from a graph are exactly maximal independent sets of its underlying graph. Now we consider simple graph H, that is, graph without directed, multiple, or weighted edges, and without self loops. Let C(H) denote the graph obtained from H by deleting an each undirected edge (u, v) of H and adding instead of this edge new two directed edges (u, v) and (v, u). The set of vertices of C(H) let be the same as the set of vertices of H. Let -1 be a weight of each edge of C(H). The Bi-stable Switch can be seen as combinatorial net derived from simple graph of order 2.
For example Figure 4b demonstrates the network derived from the graph shown in Figure 4a. The first switch is formed by the subgraph induced by the {1,2} set of vertices of the C(H) network . The second switch is formed by the {2,3} set of vertices. The vertex 2 is common one of these switches, therefore they interact by means of this vertex. Each edge of an underlying graph correspond to a switch in a derived network. If two incident edges share a common vertex then the corresponding switches interact because this vertex have the same state in both of them. We refer to C(H) as the derived network of bi-stable switches, and we refer to H as the underlying graph.  Proof. Assume by contradiction that u and v are adjacent active vertices, then I(u, M ) < 0 and the forced state of the u is inactive. If a state is steady and the forced state is inactive then the state of the vertex is inactive, whereas by conditions of the lemma u is active.
Lemma 2. If a map M of a combinatorial net that derived from a simple graph is steady and some vertex is inactive then there is at least one active vertex adjacent to it.
Proof. Assume by contradiction, that there are no active adjacent vertices of v, then influence I(v, M ) = 0 and then forced state of v is active. Hence the map M is steady state, forced state of vertex is equal to current state, then state of the vertex is active in contradiction to conditions of the lemma. Proof. Let us prove the independence first. Assume by contradiction that the set of active vertices of some steady state map C(H) is not an independent set of vertices of the graph H. Then there is a pair of adjacent vertices in H, which are simultaneously active under this steady state map M. But by lemma 1 there are no two adjacent active vertices under steady state map. This contradiction proves the independence of the set of active vertices under the steady state map. Now we consider maximality of the set of active vertices under a steady state map. Assume that there is a some inactive vertex which not adjacent to any active vertex. So that it may be added to this set to form the bigger independent set. But if such vertex exist then the map is not steady by lemma 2.

Conclusions
A similar approach to construct Markov chains for interaction graphs was developed in our earlier works for neural and gene regulatory networks [38,39,40,41,42]. Both approaches can be used to construct Markov chains of gene regulatory networks. Systems of mutually repressive elements are ubiquitous in nature. The network of bi-stable switches can be used to create models of their stable states and the self-evolution of such systems toward a stable states.