\[\def\bpi{\boldsymbol{\pi}} \def\bpit{\boldsymbol{\pi}^{\,T}} \def\bzero{\boldsymbol{0}} \def\be{\boldsymbol{e}}\]

Non-equilibrium steady states

The equilibrium formulation assumes detailed balance. However, when the state-transition diagram of a receptor model includes cycles, the steady-state probability distribution may not satisfy detailed balance.

Let us consider again the following state-transition diagram.

var('a12, a21, a13, a31, a23, a32, a34, a43')
d = {1: {2:a12, 3:a13}, 2: {1:a21, 3:a23}, 3: {2:a32, 1:a31, 4:a34}, 4: {3:a43}};
G = DiGraph(d,weighted=True)
vertex_positions = {1: (0, 0), 2: (1, 1.41), 3: (2, 0), 4: (4,0)}
G.plot(figsize=8,edge_labels=True,pos=vertex_positions,graph_border=True)

Generator matrix

The generator matrix \(Q\) for the Markov chain associated to \(G\) can be constructed from the weighted adjacency matrix \(A\), as follows

A = G.weighted_adjacency_matrix(sparse=False)
Q = A - diagonal_matrix(sum(A.T))
show(Q)
print('The rank of Q is',Q.rank())

\(\displaystyle \left(\begin{array}{rrrr} -a_{12} - a_{13} & a_{12} & a_{13} & 0 \\ a_{21} & -a_{21} - a_{23} & a_{23} & 0 \\ a_{31} & a_{32} & -a_{31} - a_{32} - a_{34} & a_{34} \\ 0 & 0 & a_{43} & -a_{43} \end{array}\right)\)

The rank of Q is 3

The generator of a Markov chain receptor model with \(n\) states has rank \(n-1\) due to conservation of probability. The four-state receptor model under consideration has a 4 states. The \(4 \times 4\) generator matrix \(Q\) has rank \(4-1=3\).

print('The rank of Q is',Q.rank())

The rank of Q is 3

The probability distribution solves a linear system of ordinary differential equations

The relationship between the adjacency matrix \(A\) and generator matrix \(Q\) can be written as

(6)\[\begin{equation} Q = A - \text{diag}(A\be) \end{equation}\]

where \(\be\) is a commensurate column vector of 1s. The following code confirms that \(Q \be = \bzero\), i.e., each row of \(Q\) sums to zero.

e = matrix([1,1,1,1]).T
show(Q*e)

\(\displaystyle \left(\begin{array}{r} 0 \\ 0 \\ 0 \\ 0 \end{array}\right)\)

The matrix \(Q\) is referred to as the generator matrix for the Markov chain because the probability distribution \(\bpit\) solves

(7)\[d\bpit/dt = \bpit Q \, . \]

Note

In the above expressions, the probability distribution \(\bpit\) is a row vector that multiplies \(Q\) on the left, while \(\be\) is a column vector that multiplies \(Q\) on the right. If one prefers to represent the steady-state probability distribution as a column vector, one can write \(\be^T Q^T = \bzero^T\) and \(d\bpi/dt = Q^T \bpi\).

Setting the left side of (7) to zero, it is evident that steady-state probability distribution \(\bpit\) solves \(\bpit Q = \bzero\) subject to \(\bpit \be = 1\). This expression is equivalent to \(\sum_i \pi_i = 1\) (conservation of probability).

Symbolic calcualtion of steady-state probability distribution

Using Sagemath we can symbolically solve \(\bpit Q = \bzero\) subject to \(\bpit \be = 1\).

To begin, we will unpack the four linear equations of \(\bpit Q = \bzero\), which is compact notation for four linear equations.

var('p1 p2 p3 p4')
p = vector([p1, p2, p3, p4])
pQ = p*Q
eq =[]
for lhs in pQ:
    show(lhs == 0)
    eq.append(lhs == 0)

\(\displaystyle -{\left(a_{12} + a_{13}\right)} p_{1} + a_{21} p_{2} + a_{31} p_{3} = 0\)

\(\displaystyle a_{12} p_{1} - {\left(a_{21} + a_{23}\right)} p_{2} + a_{32} p_{3} = 0\)

\(\displaystyle a_{13} p_{1} + a_{23} p_{2} - {\left(a_{31} + a_{32} + a_{34}\right)} p_{3} + a_{43} p_{4} = 0\)

\(\displaystyle a_{34} p_{3} - a_{43} p_{4} = 0\)

As discussed above, the generator matrix \(Q\) is rank 3. Thus, the fourth equation (a34*p3 - a43*p4 == 0) is superfluous. We replace this equation by the condition p1+p2+p3+p4 == 1 (the solution should be a normalized probability distribution).

eq[-1] = p1+p2+p3+p4 == 1
for q in eq:
    show(q)

\(\displaystyle -{\left(a_{12} + a_{13}\right)} p_{1} + a_{21} p_{2} + a_{31} p_{3} = 0\)

\(\displaystyle a_{12} p_{1} - {\left(a_{21} + a_{23}\right)} p_{2} + a_{32} p_{3} = 0\)

\(\displaystyle a_{13} p_{1} + a_{23} p_{2} - {\left(a_{31} + a_{32} + a_{34}\right)} p_{3} + a_{43} p_{4} = 0\)

\(\displaystyle p_{1} + p_{2} + p_{3} + p_{4} = 1\)

This system of four linear equations can now be solved to obtain the steady-state probability distribution.

z = solve(eq,list(p))
for i in range(4):
    f = z[0][i].rhs()
    print('p%s' % (i+1),'=',f.expand().factor())

p1 = (a21*a31 + a23*a31 + a21*a32)*a43/(a13*a21*a34 + a12*a23*a34 + a13*a23*a34 + a13*a21*a43 + a12*a23*a43 + a13*a23*a43 + a12*a31*a43 + a21*a31*a43 + a23*a31*a43 + a12*a32*a43 + a13*a32*a43 + a21*a32*a43)
p2 = (a12*a31 + a12*a32 + a13*a32)*a43/(a13*a21*a34 + a12*a23*a34 + a13*a23*a34 + a13*a21*a43 + a12*a23*a43 + a13*a23*a43 + a12*a31*a43 + a21*a31*a43 + a23*a31*a43 + a12*a32*a43 + a13*a32*a43 + a21*a32*a43)
p3 = (a13*a21 + a12*a23 + a13*a23)*a43/(a13*a21*a34 + a12*a23*a34 + a13*a23*a34 + a13*a21*a43 + a12*a23*a43 + a13*a23*a43 + a12*a31*a43 + a21*a31*a43 + a23*a31*a43 + a12*a32*a43 + a13*a32*a43 + a21*a32*a43)
p4 = (a13*a21 + a12*a23 + a13*a23)*a34/(a13*a21*a34 + a12*a23*a34 + a13*a23*a34 + a13*a21*a43 + a12*a23*a43 + a13*a23*a43 + a12*a31*a43 + a21*a31*a43 + a23*a31*a43 + a12*a32*a43 + a13*a32*a43 + a21*a32*a43)

This solution can be written more compactly as follows.

(8)\[\begin{split}p_1 & = \frac{z_1}{z_1+z_2+z_3+z_4}\\ p_2 & = \frac{z_2}{z_1+z_2+z_3+z_4}\\ p_3 & = \frac{z_3}{z_1+z_2+z_3+z_4}\\ p_4 & = \frac{z_4}{z_1+z_2+z_3+z_4}\end{split}\]

where

\[\begin{align*} z_1 & = (a_{21} a_{31} + a_{23} a_{31} + a_{21} a_{32} ) a_{43}\\ z_2 & = (a_{12} a_{31} + a_{12} a_{32} + a_{13} a_{32} ) a_{43}\\ z_3 & = (a_{13} a_{21} + a_{12} a_{23} + a_{13} a_{23} ) a_{43}\\ z_4 & = (a_{13} a_{21} + a_{12} a_{23} + a_{13} a_{23} ) a_{34} \, . \end{align*}\]

In general, this probability distribution is a non-equilibrium steady state. To see this, check to see if the distribution satisfies detailed balance, i.e., \(a_{ij} \pi_i = a_{ji} \pi_j\). Using the note above, we see that detailed balance implies \(a_{ij} z_i = a_{ji} z_j\), but \(a_{01} z_0 \neq a_{10} z_1\) in general.

Komolgorov’s criterion and equilibrium

If the product of rate constants around the cycle is the same in both directions, i.e., a12*a23*a31=a13*a32*a21 Komolgorov’s criterion is satisfied. In this case, the steady-state probability distribution is guaranteed to satisfy detailed balance. The code below usings this condition to repace a12 by a13*a32*a21/(a23*a31).

Show code cell content

Hide code cell content

def mysolve(p,Q):
    pQ = p*Q
    eq =[]
    for lhs in pQ:
       eq.append(lhs == 0)
    eq[-1] = p1+p2+p3+p4 == 1
    z = solve(eq,list(p))
    for i in range(4):
        f = z[0][i].rhs()
        print('p%s' % (i+1),'=',f.expand().factor())

The equilibrium steady-state probability distribution assuming the Komolgorov condition is

Q = Q.subs(a31=a13*a32*a21/(a12*a23))
mysolve(p,Q)

p1 = a21*a32*a43/(a12*a23*a34 + a12*a23*a43 + a12*a32*a43 + a21*a32*a43)
p2 = a12*a32*a43/(a12*a23*a34 + a12*a23*a43 + a12*a32*a43 + a21*a32*a43)
p3 = a12*a23*a43/(a12*a23*a34 + a12*a23*a43 + a12*a32*a43 + a21*a32*a43)
p4 = a12*a23*a34/(a12*a23*a34 + a12*a23*a43 + a12*a32*a43 + a21*a32*a43)

which can be written more compactly using () and

\[\begin{align*} z_1 & = a_{21} a_{32} a_{43} & z_2 & = a_{12} a_{32} a_{43} & z_3 & = a_{12} a_{23} a_{43} & z_4 & = a_{12} a_{23} a_{34} \, . \end{align*}\]

This distribution satisfies detailed balance, i.e., \(a_{ij} \pi_i = a_{ji} \pi_j\). As discussed previously, detailed balance implies that the steady-state probability distribution can be written in terms of equilibrium constants (as opposed to rate constants). Define \(\kappa_{j}=a_{ij}/a_{ji}\) for \(i < j\) whenever vertex \(i\) and \(j\) are adjacent. Dividing the numerator and denominator of each \(\pi_i\) by \(a_{12}a_{23}a_{34}\) yields the relative probabilities \(z_1 = 1\), \(z_2 = \kappa_{2}\), \(z_3 = \kappa_{2}\kappa_{3}\), and \(z_4 = \kappa_{2} \kappa_{3} \kappa_{4}\). Substituting these values into () gives the solution in terms of the equilibrium constants \(\kappa_{2}\), \(\kappa_{3}\), and \(\kappa_{4}\).

In the notation of the equilibrium formulation, this probability distribution,

(9)\[\begin{equation} [ \pi_1 : \pi_2 : \pi_3 : \pi_4] = [1 :\kappa_{2} : \kappa_{2} \kappa_{3} : \kappa_{2} \kappa_{3} \kappa_{4} ] \, , \end{equation}\]

corresponds to the following spanning tree rooted in state 1.

var('kappa_2, kappa_3, kappa_4')
d = {2: {1:kappa_2}, 3: {2:kappa_3}, 4: {3:kappa_4}};
G = DiGraph(d,weighted=True)
vertex_positions = {1: (0, 0), 2: (1, 1.41), 3: (2, 0), 4: (4,0)}
G.plot(figsize=8,edge_labels=True,pos=vertex_positions,graph_border=True)