Non-equilibrium steady states

\[\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)

A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.2.6 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.

If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.

Traceback (most recent call last):  File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/usr/lib/python3/dist-packages/sage/repl/ipython_kernel/__main__.py", line 3, in <module>
    IPKernelApp.launch_instance(kernel_class=SageKernel)
  File "/usr/lib/python3/dist-packages/traitlets/config/application.py", line 846, in launch_instance
    app.start()
  File "/usr/lib/python3/dist-packages/ipykernel/kernelapp.py", line 677, in start
    self.io_loop.start()
  File "/usr/lib/python3/dist-packages/tornado/platform/asyncio.py", line 199, in start
    self.asyncio_loop.run_forever()
  File "/usr/lib/python3.10/asyncio/base_events.py", line 603, in run_forever
    self._run_once()
  File "/usr/lib/python3.10/asyncio/base_events.py", line 1909, in _run_once
    handle._run()
  File "/usr/lib/python3.10/asyncio/events.py", line 80, in _run
    self._context.run(self._callback, *self._args)
  File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 461, in dispatch_queue
    await self.process_one()
  File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 450, in process_one
    await dispatch(*args)
  File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 357, in dispatch_shell
    await result
  File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 652, in execute_request
    reply_content = await reply_content
  File "/usr/lib/python3/dist-packages/ipykernel/ipkernel.py", line 353, in do_execute
    res = shell.run_cell(code, store_history=store_history, silent=silent)
  File "/usr/lib/python3/dist-packages/ipykernel/zmqshell.py", line 532, in run_cell
    return super().run_cell(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2914, in run_cell
    result = self._run_cell(
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2960, in _run_cell
    return runner(coro)
  File "/usr/lib/python3/dist-packages/IPython/core/async_helpers.py", line 78, in _pseudo_sync_runner
    coro.send(None)
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3185, in run_cell_async
    has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3377, in run_ast_nodes
    if (await self.run_code(code, result,  async_=asy)):
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3457, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "/tmp/ipykernel_12038/3336743817.py", line 5, in <module>
    G.plot(figsize=Integer(8),edge_labels=True,pos=vertex_positions,graph_border=True)
  File "/usr/lib/python3/dist-packages/IPython/core/displayhook.py", line 262, in __call__
    format_dict, md_dict = self.compute_format_data(result)
  File "/usr/lib/python3/dist-packages/IPython/core/displayhook.py", line 151, in compute_format_data
    return self.shell.display_formatter.format(result)
  File "/usr/lib/python3/dist-packages/sage/repl/display/formatter.py", line 181, in format
    sage_format, sage_metadata = self.dm.displayhook(obj)
  File "/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py", line 825, in displayhook
    plain_text, rich_output = self._rich_output_formatter(obj, dict())
  File "/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py", line 643, in _rich_output_formatter
    rich_output = self._call_rich_repr(obj, rich_repr_kwds)
  File "/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py", line 603, in _call_rich_repr
    return obj._rich_repr_(self)
  File "/usr/lib/python3/dist-packages/sage/plot/graphics.py", line 1000, in _rich_repr_
    return display_manager.graphics_from_save(
  File "/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py", line 731, in graphics_from_save
    save_function(filename, **kwds)
  File "/usr/lib/python3/dist-packages/sage/misc/decorators.py", line 410, in wrapper
    return func(*args, **kwds)
  File "/usr/lib/python3/dist-packages/sage/plot/graphics.py", line 3296, in save
    from matplotlib import rcParams
  File "/usr/lib/python3/dist-packages/matplotlib/__init__.py", line 109, in <module>
    from . import _api, _version, cbook, docstring, rcsetup
  File "/usr/lib/python3/dist-packages/matplotlib/rcsetup.py", line 27, in <module>
    from matplotlib.colors import Colormap, is_color_like
  File "/usr/lib/python3/dist-packages/matplotlib/colors.py", line 56, in <module>
    from matplotlib import _api, cbook, scale
  File "/usr/lib/python3/dist-packages/matplotlib/scale.py", line 23, in <module>
    from matplotlib.ticker import (
  File "/usr/lib/python3/dist-packages/matplotlib/ticker.py", line 136, in <module>
    from matplotlib import transforms as mtransforms
  File "/usr/lib/python3/dist-packages/matplotlib/transforms.py", line 46, in <module>
    from matplotlib._path import (

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
AttributeError: _ARRAY_API not found

/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py:608: RichReprWarning: Exception in _rich_repr_ while displaying object: numpy.core.multiarray failed to import
  warnings.warn(

A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.2.6 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.

If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.

Traceback (most recent call last):  File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/usr/lib/python3/dist-packages/sage/repl/ipython_kernel/__main__.py", line 3, in <module>
    IPKernelApp.launch_instance(kernel_class=SageKernel)
  File "/usr/lib/python3/dist-packages/traitlets/config/application.py", line 846, in launch_instance
    app.start()
  File "/usr/lib/python3/dist-packages/ipykernel/kernelapp.py", line 677, in start
    self.io_loop.start()
  File "/usr/lib/python3/dist-packages/tornado/platform/asyncio.py", line 199, in start
    self.asyncio_loop.run_forever()
  File "/usr/lib/python3.10/asyncio/base_events.py", line 603, in run_forever
    self._run_once()
  File "/usr/lib/python3.10/asyncio/base_events.py", line 1909, in _run_once
    handle._run()
  File "/usr/lib/python3.10/asyncio/events.py", line 80, in _run
    self._context.run(self._callback, *self._args)
  File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 461, in dispatch_queue
    await self.process_one()
  File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 450, in process_one
    await dispatch(*args)
  File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 357, in dispatch_shell
    await result
  File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 652, in execute_request
    reply_content = await reply_content
  File "/usr/lib/python3/dist-packages/ipykernel/ipkernel.py", line 353, in do_execute
    res = shell.run_cell(code, store_history=store_history, silent=silent)
  File "/usr/lib/python3/dist-packages/ipykernel/zmqshell.py", line 532, in run_cell
    return super().run_cell(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2914, in run_cell
    result = self._run_cell(
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2960, in _run_cell
    return runner(coro)
  File "/usr/lib/python3/dist-packages/IPython/core/async_helpers.py", line 78, in _pseudo_sync_runner
    coro.send(None)
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3185, in run_cell_async
    has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3377, in run_ast_nodes
    if (await self.run_code(code, result,  async_=asy)):
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3457, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "/tmp/ipykernel_12038/3336743817.py", line 5, in <module>
    G.plot(figsize=Integer(8),edge_labels=True,pos=vertex_positions,graph_border=True)
  File "/usr/lib/python3/dist-packages/IPython/core/displayhook.py", line 262, in __call__
    format_dict, md_dict = self.compute_format_data(result)
  File "/usr/lib/python3/dist-packages/IPython/core/displayhook.py", line 151, in compute_format_data
    return self.shell.display_formatter.format(result)
  File "/usr/lib/python3/dist-packages/sage/repl/display/formatter.py", line 186, in format
    if (not isinstance(obj, (IPYTHON_NATIVE_TYPES, Figure)) and
  File "/usr/lib/python3/dist-packages/matplotlib/__init__.py", line 109, in <module>
    from . import _api, _version, cbook, docstring, rcsetup
  File "/usr/lib/python3/dist-packages/matplotlib/rcsetup.py", line 27, in <module>
    from matplotlib.colors import Colormap, is_color_like
  File "/usr/lib/python3/dist-packages/matplotlib/colors.py", line 56, in <module>
    from matplotlib import _api, cbook, scale
  File "/usr/lib/python3/dist-packages/matplotlib/scale.py", line 23, in <module>
    from matplotlib.ticker import (
  File "/usr/lib/python3/dist-packages/matplotlib/ticker.py", line 136, in <module>
    from matplotlib import transforms as mtransforms
  File "/usr/lib/python3/dist-packages/matplotlib/transforms.py", line 46, in <module>
    from matplotlib._path import (

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
AttributeError: _ARRAY_API not found

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
/tmp/ipykernel_12038/3336743817.py in <module>
      3 G = DiGraph(d,weighted=True)
      4 vertex_positions = {Integer(1): (Integer(0), Integer(0)), Integer(2): (Integer(1), RealNumber('1.41')), Integer(3): (Integer(2), Integer(0)), Integer(4): (Integer(4),Integer(0))}
----> 5 G.plot(figsize=Integer(8),edge_labels=True,pos=vertex_positions,graph_border=True)

/usr/lib/python3/dist-packages/IPython/core/displayhook.py in __call__(self, result)
    260             self.start_displayhook()
    261             self.write_output_prompt()
--> 262             format_dict, md_dict = self.compute_format_data(result)
    263             self.update_user_ns(result)
    264             self.fill_exec_result(result)

/usr/lib/python3/dist-packages/IPython/core/displayhook.py in compute_format_data(self, result)
    149 
    150         """
--> 151         return self.shell.display_formatter.format(result)
    152 
    153     # This can be set to True by the write_output_prompt method in a subclass

/usr/lib/python3/dist-packages/sage/repl/display/formatter.py in format(self, obj, include, exclude)
    184         # use Sage rich output for any except those native to IPython, but only
    185         # if it is not plain and dull
--> 186         if (not isinstance(obj, (IPYTHON_NATIVE_TYPES, Figure)) and
    187             not set(sage_format.keys()).issubset([PLAIN_TEXT])):
    188             return sage_format, sage_metadata

/usr/lib/python3/dist-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport.__instancecheck__ (build/cythonized/sage/misc/lazy_import.c:7695)()
    914             True
    915         """
--> 916         return isinstance(x, self.get_object())
    917 
    918     def __subclasscheck__(self, x):

/usr/lib/python3/dist-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport.get_object (build/cythonized/sage/misc/lazy_import.c:2612)()
    215         if likely(self._object is not None):
    216             return self._object
--> 217         return self._get_object()
    218 
    219     cpdef _get_object(self):

/usr/lib/python3/dist-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport._get_object (build/cythonized/sage/misc/lazy_import.c:3073)()
    255             if self._feature:
    256                 raise FeatureNotPresentError(self._feature, reason=f'Importing {self._name} failed: {e}')
--> 257             raise
    258 
    259         name = self._as_name

/usr/lib/python3/dist-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport._get_object (build/cythonized/sage/misc/lazy_import.c:2935)()
    251 
    252         try:
--> 253             self._object = getattr(__import__(self._module, {}, {}, [self._name]), self._name)
    254         except ImportError as e:
    255             if self._feature:

/usr/lib/python3/dist-packages/matplotlib/__init__.py in <module>
    107 # cbook must import matplotlib only within function
    108 # definitions, so it is safe to import from it here.
--> 109 from . import _api, _version, cbook, docstring, rcsetup
    110 from matplotlib.cbook import MatplotlibDeprecationWarning, sanitize_sequence
    111 from matplotlib.cbook import mplDeprecation  # deprecated

/usr/lib/python3/dist-packages/matplotlib/rcsetup.py in <module>
     25 from matplotlib import _api, cbook
     26 from matplotlib.cbook import ls_mapper
---> 27 from matplotlib.colors import Colormap, is_color_like
     28 from matplotlib.fontconfig_pattern import parse_fontconfig_pattern
     29 from matplotlib._enums import JoinStyle, CapStyle

/usr/lib/python3/dist-packages/matplotlib/colors.py in <module>
     54 import matplotlib as mpl
     55 import numpy as np
---> 56 from matplotlib import _api, cbook, scale
     57 from ._color_data import BASE_COLORS, TABLEAU_COLORS, CSS4_COLORS, XKCD_COLORS
     58 

/usr/lib/python3/dist-packages/matplotlib/scale.py in <module>
     21 import matplotlib as mpl
     22 from matplotlib import _api, docstring
---> 23 from matplotlib.ticker import (
     24     NullFormatter, ScalarFormatter, LogFormatterSciNotation, LogitFormatter,
     25     NullLocator, LogLocator, AutoLocator, AutoMinorLocator,

/usr/lib/python3/dist-packages/matplotlib/ticker.py in <module>
    134 import matplotlib as mpl
    135 from matplotlib import _api, cbook
--> 136 from matplotlib import transforms as mtransforms
    137 
    138 _log = logging.getLogger(__name__)

/usr/lib/python3/dist-packages/matplotlib/transforms.py in <module>
     44 
     45 from matplotlib import _api
---> 46 from matplotlib._path import (
     47     affine_transform, count_bboxes_overlapping_bbox, update_path_extents)
     48 from .path import Path

ImportError: numpy.core.multiarray failed to import

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())

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 4 states. The \(4 \times 4\) generator matrix \(Q\) has rank \(4-1=3\).

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

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)

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 calculation 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)

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)

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())

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_{12} z_1 \neq a_{21} z_2\) in general.

Kolmogorov’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 Kolmogorov’s criterion is satisfied. In this case, the steady-state probability distribution is guaranteed to satisfy detailed balance. The code below uses this condition to replace a12 by a13*a32*a21/(a23*a31).

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 Kolmogorov condition is

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

which can be written more compactly using (8) 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 (8) 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)