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

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,

This influence is determined by the map function M and the weighting function W. Only if the weighting function W is assumed to be constant over time can it be said that the influence I(v, M) is the influence of map M on vertex v. If I(v, M) ≥ 0, it is said that the map M activates vertex v, otherwise that the map M represses (or inactivates) vertex v. Now the most important definition of a vertex steady state under map M can be provided. Let v be an active vertex under map M. If map M activates v, then the state of v is a steady state under map M, else the state of v is a non-steady state under map M. Now let v be an inactive vertex under map M. If map M inactivates v, then the state of v is a steady state under map M, else the state of v is a non-steady state under map M. If all vertices are in steady state under map M, then map M is called a steady-state map. The forced state of vertex v under map M is defined as follows:

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-steady-state 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 1 – p_i. Let S = {v₁, v₂, … , 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₁, p₂, … , 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 that 1_X(v_i) = 1 if v_i ∈ X, otherwise 1_X(v_i) = 0. Let P_X: S → [0.0, 1.0] denote the function such that P_X(v_i) = p_i if 1_X(v_i) = 1, otherwise P_X(v_i) = 1 − p_i. 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.

The transition graph of the interaction graph

Let S be the set of non-steady state vertices under map M = Y_j. S represents a full set of vertices, each of which can flip to a forced state at the next j + 1 time step. In the combinatorial model, steady-state vertices cannot flip. To construct a transition graph, the full set of maps M₁, M₂, … that are directly reachable from map M should be defined. Each pair (M, M_i) will correspond to one edge in the transition graph. What is the set of maps M₁, M₂, … that can be directly reached (by a one-step transition) from map M? To represent independence and asynchronicity, it is assumed that any random set of non-steady-state vertices X ⊆ S can produce the next map from the current map. Let M_X be a map such that M_X(v) = F(v, M) if v Ï X and M_X(v) = M(v) if v ∉ X; 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.

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³⁰. The model of auto-repression presented here and shown in Figure 3 and Figure 2a also exhibits stochastic oscillating behavior. Figure 2a presents an interaction graph G = (V, E) of the auto-repression model. The set V contains a single vertex v₁, and the set E consists of a single edge e1 = (v1, v1) from vertex v1 to itself. The weight of e1 is equal to -1. Let M₀ be a starting map, and let M₀(v₁) = 0, i.e., v₁ is inactive under M₀. Therefore, the influence of M₀ on vertex v₁ equals 0, I(v1, M₀) = 0. By Equation (2), F(v₁, M₀) = 1, and the state of vertex v₁ is non-steady under map M₀. Let p₁ = 1/2; then the state of vertex v₁ can change to an active state with probability 1/2 and can remain unchanged with probability 1/2. The other state of the model is also non-steady. 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.

Figure 2.

The interaction graph of: (a) The auto-repression model; (b) The bi-stable switch model; (c) The Elowitz repressilator model; (d) The self-activation model.

Figure 3. The transition graph of an auto-repression model.

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³¹. 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³⁴ 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₁, v₂} 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. Maps 00 and 11 are non-steady-state, whereas maps 01 and 10 are steady-state. To show that a map M is steady-state, it is necessary to show that each vertex v is steady-state under map M. To show that a vertex v is steady-state under map M, one must compute the value of the influence F(v, M) of map M on vertex v and compare it to the state of v under map M, i.e. M(v). Now it will be shown that map 01 is steady-state.

Figure 4. The transition graph of a bi-stable switch.

Combinatorial model of the Elowitz repressilator

The Elowitz repressilator consists of three genes³⁵, 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. Figure 2c present an interaction graph G = (V, E) of the Elowitz repressilator model, where the set V = {v₁, v₂, v₃} contains three vertices and the set E = {(v₁, v₂), (v₂, v₃), (v₃, v₁)} contains three edges. Oscillations are produced by circle {1,3,2,6,4,5} in Figure 5. Decimal, binary, and graphical representations of the state space of the Elowitz repressilator are shown in Table 1.

Figure 5. The transition graph of the Elowitz repressilator model.

Table 1. Decimal, binary, and graphical representations of the state space of the Elowitz repressilator.

Dec	Binary	Graphic	Non steady vertices	Direct reachable states
0	000			0, 1, 2, 3, 4, 5, 6, 7
1	001			1, 3
2	010			2, 6
3	011			3, 2
4	100			4, 5
5	101			5, 1
6	110			6, 4
7	111			0, 1, 2, 3, 4, 5, 6, 7

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 non-steady 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.

Figure 2d and Figure 6 represent the combinatorial model of self-activation.

Figure 6. The transition graph of a self-activation model.