Finite-state discrete-time Markov chain models of gene regulatory networks [version 1; peer review: peer review discontinued]

In this study, Markov chain models of gene regulatory networks (GRN) are developed. These models make it possible to apply the well-known theory and tools of Markov chains to GRN analysis. A new kind of finite interaction graph called a combinatorial net is introduced to represent formally a GRN and its transition graphs constructed from interaction graphs. The system dynamics are defined as a random walk on the transition graph, which is a Markov chain. A novel concurrent updating scheme (evolution rule) is developed to determine transitions in a transition graph. The proposed scheme is based on the firing of a random set of non-steady-state vertices in a combinatorial net. It is demonstrated that this novel scheme represents an advance in asynchronicity modeling. The theorem that combinatorial nets with this updating scheme can asynchronously compute a maximal independent set of graphs is also proved. As proof of concept, a number of simple combinatorial models are presented here: a discrete auto-regression model, a bistable switch, an Elowitz repressilator, and a self-activation model, and it is shown that these models exhibit well-known properties.


Introduction
Efforts to study gene-expression regulation networks has led to detailed descriptions of many such networks, and many more can be expected to be identified in the near future. Therefore, there is a need to develop methods of computational and theoretical analysis of gene regulatory networks (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 a GRN. Classification of process states, studies of long-term behavior 7 , and development of optimal strategies for therapeutic intervention [7][8][9][10][11][12][13][14][15] provide good examples of this approach 16 . In contrast to Boolean networks, Hopfield networks are defined using arithmetic operations 17 . They are a well-developed branch of science dealing 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 a Markov chain. In the Hopfield network field, essential results have been obtained in the study of various update schemes 18 , network oscillations 19 , solutions of combinatorial optimization problems [19][20][21][22][23][24][25][26] , estimation of convergence rates, and many other problems. This makes it valuable to study the possibility of using Hopfield like-networks to construct Markov chains from GRNs and other interaction graphs. Here, a GRN is considered to be a kind of interaction graph. Interaction (regulatory) graphs have emerged in various fields of the life sciences 27 . In recent years, their transition graphs have often been used to analyze the properties of interactions (regulations). One promising way to understand the nature of the regulations or interactions represented by interaction graphs is to analyze the Markov chain associated with their transition graphs.

Method
The proposed method may be viewed as a version of the Hopfieldlike network 17 where groups of randomly selected unstable units are updated in parallel 18 .
Interaction graphs and non-steady-state vertices Let G = (V, E) be a directed graph, where V is a set of vertices and E is a set of edges. Let B = {0, 1} be a set of vertex states. It is said that the vertex is active if the state of this vertex is equal to "1", otherwise it is said that the vertex is inactive. The mapping function M: V → B gives the state of each vertex. For a given vertex v ∈ V, M(v) is the state of vertex v that corresponds to map M. M(v) = 1 is equivalent to saying that vertex v is active under map M. M(v) = 0 is equivalent to saying that vertex v is inactive. The weighting function W: E → R gives value of each edge of graph G, which represents the power of interactions. If e = (u, v), e ∈ E, then u is said to be a direct ancestor of v and v to be a direct descendant of u. The influence 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 the forced state and the current state of v are the same, then the current state of vertex v under map M is steady. The following equation provides a formal definition: The random set update rule Now 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 denotes the map of G at time period j. At each time period j for each non-steadystate vertex v i under map Y j , the current vertex state is changed to a forced state with probability p i , and the current state remains unchanged with probability 1p i . Let S = {v 1 , v 2 , … , v n } be a set of all non-steady-state vertices at time period j. The vertices chosen to change their state in a one-step transition constitute a random set X ⊆ S. To create a new map that is directly accessible from the current map Y j , all vertices in X simultaneously change their state, whereas other vertices remain unchanged. Let P={p 1 , p 2 , … , p n } be some vector of numbers such that 0 ≤ p i ≤ 1, where p i is referred to as the probability of a state change (firing) of the non-steady-state vertex v i . For any X ⊆ S, let 1 X : S → B be an indicator function such Hence, it can be assumed that each v i acts independently to create a random set X. Then the production of P X gives p X , which is the firing probability of the random set X: Evidently, Now this definition of the random set update rule and its probabilities can be used to define the transition graph of the combinatorial net model.

Examples of combinatorial models
The method described above will now be used to develop models of some important graphs of repressive interactions of self-activating nodes and to prove their main properties. Such models are called combinatorial models. Each combinatorial model consists of the interaction graph (combinatorial net) and the corresponding transition graph.
The combinatorial model of an auto-repression A negative auto-regulation or an auto-repression occurs when the products of a certain gene represse their own gene. This form of simple regulation serves as a basic building block for most important transcription networks 27,29 . Auto-repression can produce oscillations. For example, embryonic stem cells fluctuate between high and low Nanog expression, and Nanog activity is auto-repressive 30 .
The model of auto-repression presented here and shown in Figure 3 and Figure  This means that there are only non-steady states in the model. Therefore, it will oscillate infinitely between 0 and 1. Figure 3 shows the full transition graph of the auto-repression model.
then M X can be said to be produced by X from M. In other words, the random set X of non-steady-state vertices produces the directly reachable map M X from a map M. The weights of the edges from map M to map M X are given by the probability defined by Equation (4).

Random-walk network dynamics
Assume that whenever the process is in state M, there is a probability p X that at the 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.

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

Combinatorial model of a bi-stable switch
A bi-stable switch is a bi-stable gene regulatory network that is constructed from two mutually repressive genes 31 . These are very common in nature and extensively used in synthetic biology 32,33 . The ordinary differential equations (ODEs) used to construct their mathematical models are a convenient way to analyze some small circuits in detail. In this research, techniques have been developed that can be used to construct models of large networks of bi-stable switches and to prove some of their important properties. For this purpose, a probabilistic coarse-scale modeling approach 34 has been used here instead of fine-scale ODE modeling. The proposed model of a bi-stable switch illustrated in Figure 2b and Figure 4 exhibits two steady maps. Figure 2b presents an interaction graph G = (V, E) of the bi-stable switch model. The set V = {v 1 , v 2 } contains two vertices, and the set E contains two edges with weights of -1. Let the probability of firing a non-steady-state vertex be 1/2. Figure 4 presents the transition graph of the model.

Combinatorial model of the Elowitz repressilator
The Elowitz repressilator consists of three genes 35 , each of which is 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, the third gene inhibits the first gene's expression, completing the cycle. Such a negative feedback loop leads to oscillations. The combinatorial model of an Elowitz repressilator produces oscillations and consists of three vertices and three edges with weights equal to -1.  Table 1.
Combinatorial model of self-activation A constitutively expressed gene is an example of self-activation. Such genes do not require any interaction to be active. A combinatorial model of self-activation consists of one vertex and no edges. In any case, the influence on it equals 0 because there are no other vertices. Therefore, a forced state of the vertex equals 1. Therefore, 1 is a steady state and 0 is a non-steady state. A vertex starting in steady state will stay in it infinitely. A vertex starting in a nonsteady state will flip to steady state with probability p and stay in non-steady state with probability 1-p. The amount of time T which the vertex spends in non-steady state is the random variable. The distribution of this random variable is a shifted geometric distribution with parameter p.

A 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. A 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. Therefore, cliques can be used to find or to enumerate MISs 20,21 . Finding independent sets (or cliques) has applications in various fields 37 Figure 7b illustrates the combinatorial network derived from the graph shown in Figure 7a. 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. Vertex 2 is a common member of these switches, and therefore they can interact by means of this vertex. Each edge of an underlying graph corresponds to a switch in a derived network. If two incident edges share a common vertex, then the corresponding switches interact because this vertex has the same state in both switches. C(H) is referred to as the derived network of a bi-stable switch, and H is referred to as the underlying graph.

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