#  The noise intensity of a Markov chain

Lukas Ramlow  Benjamin Lindner  Bernstein Center for Computational
Neuroscience Berlin, Philippstr. 13, Haus 2, 10115 Berlin, Germany  Physics
Department of Humboldt University Berlin, Newtonstr. 15, 12489 Berlin, Germany

###### Abstract

Stochastic transitions between discrete microscopic states play an important
role in many physical and biological systems. Often, these transitions lead to
fluctuations on a macroscopic scale. A classic example from neuroscience is
the stochastic opening and closing of ion channels and the resulting
fluctuations in membrane current. When the microscopic transitions are fast,
the macroscopic fluctuations are nearly uncorrelated and can be fully
characterized by their mean and noise intensity. We show how, for an arbitrary
Markov chain, the noise intensity can be determined from an algebraic
equation, based on the transition rate matrix. We demonstrate the validity of
the theory using an analytically tractable two-state Markovian dichotomous
noise, an eight-state model for a Calcium channel subunit (De Young-Keizer
model), and Markov models of the voltage-gated Sodium and Potassium channels
as they appear in a stochastic version of the Hodgkin-Huxley model.

##  I Introduction

For models of many fluctuation phenomena in physics, biology, chemistry, and
other fields, it is important to properly characterize the noise that drives a
given dynamical system. Examples include the firing of a neuron, driven by
channel noise WhiRub00 ; DorWhi05 ; SchFis10 ; FisSch12 ; MoeIan16 ; PuTho21
and by shot-noise-like input from other neurons HohBur01 ; RicGer06 ; WolLin10
; RicSwa10 ; BriDes15 ; DroLin17 , the fluctuations in the intensity of an
excitable laser HakSau67 ; DubKra99 ; LamGuz03 ; SchTia12 , or chemical
reactions in mesoscopically small volumes Gar85 ; Van92 ; Gil00 . In the past
decades, strongly simplified noise models, such as white Gaussian noise,
Poissonian shot noise, dichotomous noise, or an exponentially correlated
Ornstein-Uhlenbeck noise have often been used to describe the input noise in
these systems. As the models for the driving noise process become more
complex, one would like to use characteristics of the noise process that can
be used to fairly compare different noise models (and their effect on a
dynamical system). This fair comparison is already possible for simple noise
processes, such as different exponentially correlated processes like the
Gaussian Ornstein-Uhlenbeck process UhlOrn30 ; Ris84 ; HanJun95 , the
dichotomous telegraph noise HorLef83 ; Gar85 ; Ben06 ; DroLin14 , or the
exponentially distributed noise FarLin21  (which is, however, only
approximately exponentially correlated).

Simple characteristics of a stochastic (noise) process $x(t)$ are its
stationary mean and variance

$\displaystyle\mu=\left\langle
x(t)\right\rangle,\;\;\;\sigma^{2}=\left\langle\Delta
x^{2}\right\rangle=\left\langle x^{2}(t)\right\rangle-\left\langle
x(t)\right\rangle^{2},$ (1)

its correlation time

$\displaystyle\tau=\int_{0}^{\infty}dt^{\prime}\frac{\left\langle
x(t)x(t+t^{\prime})\right\rangle-\left\langle
x(t)\right\rangle^{2}}{\left\langle\Delta x^{2}\right\rangle},$ (2)

and its noise intensity

$\displaystyle D=\int_{0}^{\infty}dt^{\prime}\left\langle
x(t)x(t+t^{\prime})\right\rangle-\left\langle x(t)\right\rangle^{2}.$ (3)

The meaning of mean and variance are quite obvious. The correlation time (here
defined by an integral over the normalized autocorrelation function) provide
an order-of-magnitude estimate of the periods over which the process changes
significantly. Last but not least, the intensity (here defined by an integral
over the _un_ normalized autocorrelation function) captures how much of an
effect the process would have when driving a dynamical system. More
specifically, if $x(t)$ were the velocity of a Brownian particle, $D$ would
correspond to the diffusion coefficient, which is a reasonable measure of the
effect of the velocity noise on the position dynamics.

It is clear from the above definitions that correlation time, variance, and
intensity are connected by

$\displaystyle D=\sigma^{2}\tau,$ (4)

i.e. if we know two of the characteristics, we can easily compute the third
one. For processes described by a nonlinear Langevin equation, all four
characteristics (including also the mean value) can be expressed by
quadratures Ris89 ; see also refs. Lin07 ; Lin08 , which make the above-
mentioned connection between noise intensity and diffusion coefficient more
explicit.

For discrete-valued processes, $x(t)\in\\{x_{1},x_{2},\dots\\}$, governed by a
master equation, the mean and the variance can easily be calculated in all
cases where the stationary probability can be obtained. The calculation of the
noise intensity is more involved but has recently been worked out in our
previous paper and applied for a specific model RamFal23 ; RamFal23b . It
turned out that, in complete analogy to the procedure for the Langevin case
studied in Ris89 ; Lin07 , we can obtain closed-form expressions for these
characteristics. Multiple integrals in the continuous Langevin case correspond
to multiple sums in the Markov chain case. Here, we show that the theory is
not limited to a specific model, but can be applied to any random process
where the state transitions are governed by a master equation.

Our paper is structured as follows. We begin in sec. (II) with the general
framework for calculating the noise intensity. Then we discuss three examples.
In sec. (III) we illustrate the general result for the (simple and well-known)
case of Markovian dichotomous noise. In sec. (IV) we study the more
complicated case of an eight-state Markov model, as used in the De Young-
Keizer model to describe a subunit of calcium release channel. In sec. (V) we
study a stochastic version of the sodium and potassium currents as they appear
in the Hodgkin-Huxley model with channel noise. Finally, in sec. (VI), we
discuss further applications of our results, such as a general white-noise
approximation in cases where the microscopic transitions are much faster than
the dynamics of the driven system. However, we also point out some limitations
for systems where the noise is _not_ purely external and/or very fast but also
depends on the state of the driven system.

##  II Noise intensity of a Markov chain

We consider a random process $x(t)$ with discrete states, where the
probability $p_{i}(t)$ of finding a state $i$ at time $t$ is determined by the
(homogeneous) master equation:

$\displaystyle\dot{\bm{p}}(t)=W\bm{p}(t)$ (5)

with the probability vector
$\bm{p}(t)=\begin{pmatrix}p_{1}(t)&p_{2}(t)&\dots\end{pmatrix}$ and the
transition rate matrix $W=(w_{ij})$. The entries $w_{ij}>0$ for $i\neq j$ are
the transition rates from a state $j$ to state a $i$ and $w_{jj}=-\sum_{i\neq
j}w_{ij}$. To fully characterize the process $x(t)$, each state $i$ is
assigned a specific value $x_{i}$, which are not necessarily different from
each other.

For such a random process, the calculation of the mean and the variance
follows standard procedures Gar85 ; Van92  and is based on the stationary
probability vector $\lim_{t\to\infty}\bm{p}(t)=\bm{p}$ (we indicate the
stationary state by omitting the time argument). This vector can be obtained
from the stationary master equation

$\displaystyle 0=W\bm{p}$ (6)

together with the normalization condition $\sum_{i}p_{i}=1$; the additional
condition is needed because of the rank deficiency of the matrix $W$.
Practically, the normalization can be incorporated by replacing an arbitrary
row of the matrix $W$ with ones and the corresponding entry in the zero vector
on the l.h.s. by a one. This leads to a linear system of equations solvable by
standard methods. Given the stationary probabilities $p_{i}$, the mean and the
variance can be calculated by

$\displaystyle\mu$ $\displaystyle=\sum x_{i}p_{i},$ $\displaystyle\sigma^{2}$
$\displaystyle=\sum_{i}(x_{i}-\mu)^{2}p_{i}.$ (7)

The calculation of the noise intensity is more advanced. We recall the
definition of the noise intensity by the integral over the autocorrelation
function eq. (3). In principle, the correlation function can be determined by
solving the time-dependent master eq. (5). However, it turns out that the
calculation of the time-dependent probability vector $\bm{p}(t)$ is not
necessary in order to calculate the noise intensity. Instead, taking advantage
of the fact that the integrated correlation function is of interest, an
algebraic equation can be found that determines the noise intensity and is not
much more complicated to solve than the equations that determine the mean or
variance.

To show this, we relate the noise intensity to the probabilities of the Markov
chain:

$\displaystyle D$ $\displaystyle=\int_{0}^{\infty}dt^{\prime}\,\langle
x(t+t^{\prime})x(t)\rangle-\langle x(t)\rangle^{2}$ (8)
$\displaystyle=\int_{0}^{\infty}dt^{\prime}\,\sum_{i,j}[x_{i}x_{j}p_{ij}(t^{\prime})p_{j}-x_{i}x_{j}p_{i}p_{j}]$
$\displaystyle=\sum_{i,j}x_{i}f_{ij}x_{j}p_{j}.$

where $p_{ij}(t^{\prime})=p(i,t+t^{\prime}|j,t)$ is the transition
probability, i.e. the probability of finding the state $i$ at $t+t^{\prime}$
given the state $j$ at time $t$. Since we are considering a homogeneous
process, this conditional probability does not depend on the absolute time
$t$, but only on the difference $t^{\prime}$ and is determined by the master
eq. (5) with the initial condition $p_{k}(t)=\delta_{kj}$. The auxiliary
function introduced in the last line of eq. (8),

$\displaystyle f_{ij}=\int_{0}^{\infty}dt^{\prime}\,p_{ij}(t^{\prime})-p_{i},$
(9)

is given by the integral over the difference between the transition and
stationary probabilities. Eq. (8) is of course just a reformulation of the
problem. However, it turns out that the auxiliary functions $f_{ij}$ for a
given $j$ can be calculated from a system of algebraic equations together with
an additional condition, a calculation that is very similar to that of the
stationary probabilities $p_{i}$. To see this, we formulate the master
equation where the state at some reference time $t$ has been specified
($p_{kj}(0)=\delta_{kj}$):

$\displaystyle\dot{p}_{kj}(t^{\prime})$
$\displaystyle=\sum_{i}w_{ki}p_{ij}(t^{\prime}),$ (10)
$\displaystyle\dot{p}_{kj}(t^{\prime})$
$\displaystyle=\sum_{i}w_{ki}[p_{ij}(t^{\prime})-p_{i}],$ $\displaystyle
p_{k}-\delta_{kj}$ $\displaystyle=\sum_{i}w_{ki}f_{ij}.$

To get from the first to the second line we have subtracted the stationary
master equation $0=\textstyle\sum_{i}w_{ki}p_{i}$. To get from the second to
the third line we have integrated over $t^{\prime}$, used eq. (10), and
exploited that
$\int_{0}^{\infty}dt^{\prime}\,\dot{p}_{kj}(t^{\prime})=p_{k}-\delta_{kj}$.
The last line in eq. (10) looks like an equation that uniquely determines
$f_{ij}$. However, because of the rank deficiency of $W$ we need an additional
condition that is obtained by observing that

$\displaystyle\sum_{i}f_{ij}=\int_{0}^{\infty}dt^{\prime}\,\sum_{i}p_{ij}(t^{\prime})-p_{i}=0.$
(11)

This condition is independent of $j$ and reflects that for any $t^{\prime}$
both the transition probability and the stationary probability are normalized
over the states $i$.

Finally, while eqs. (8) - (11) allow for the calculation of noise intensity,
they can be expressed more conveniently. For this purpose, we write eq. (8) in
matrix notation

$\displaystyle D=\bm{x}^{T}F\bm{y},$ (12)

where $f_{ij}$ is the entry in the $i$-th row and $j$-th column of the matrix
$F$ and the two vectors are given by
$\bm{x}=\begin{pmatrix}x_{1}&x_{2}&...\end{pmatrix}^{T}$ and
$\bm{y}=\begin{pmatrix}x_{1}p_{1}&x_{2}p_{2}&...\end{pmatrix}^{T}$. Similarly,
the set of linear eqs. (10) can be combined into a single matrix equation

$\displaystyle P-\mathbb{1}=WF$ (13)

where $P$ is a matrix in which each entry in the $i$-th row is given by the
stationary probability $p_{i}$ and $\mathbb{1}$ is the identity matrix. The
additional condition eq. (11) can be written as

$\displaystyle\begin{pmatrix}1&1&\dots&1\end{pmatrix}F=0,$ (14)

implying that each column of the matrix $F$ adds up to zero. Again,
practically, these conditions can be incorporated by replacing an arbitrary
row in $W$ by ones and the corresponding row of the matrix $P-\mathbb{1}$ by
zeroes.

In the following, we put the theory to the test for different models. The
first model is the well-known Markovian two-state model, the second is a more
involved eight-state model for a subunit of a Calcium channel (De Young-Keizer
model), the third example comprises the Sodium and Potassium channels in a
stochastic version of the Hodgkin-Huxley model of an excitable nerve membrane.

##  III A simple example:
Markovian dichotomous noise

As an introduction to the method, we consider a Markovian dichotomous noise
for which the noise intensity is known and can be calculated in several ways
HorLef83 . A dichotomous Markov noise $x(t)$ is a Markov process with two
levels $x_{1}$ and $x_{2}$, corresponding to two different states with
transition rates $\alpha$ and $\beta$ between them. A schematic representation
of the model is shown in Fig. 1.

Figure 1: State diagram of a dichotomous noise. The system consists of two
states $1$ and $2$ with corresponding levels $x_{1}$ and $x_{2}$. The
transitions between the two states occur at rates $\alpha$ and $\beta$.

To calculate the noise intensity according to eq. (12), we first determine the
stationary probabilities using the stationary master eq. (6)

$\displaystyle 0=\begin{pmatrix}-\alpha&\beta\\\
\alpha&-\beta\end{pmatrix}\begin{pmatrix}p_{1}\\\ p_{2}\end{pmatrix},$ (15)

together with the normalization condition $p_{1}+p_{2}=1$ and obtain
$p_{1}=\beta/(\alpha+\beta)$ and $p_{2}=\alpha/(\alpha+\beta)$. We can now
calculate $F$ using eq. (13)

$\displaystyle\begin{pmatrix}p_{1}&p_{1}\\\ p_{2}&p_{2}\\\ \end{pmatrix}$
$\displaystyle-\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}=$ (16)
$\displaystyle\begin{pmatrix}-\alpha&\beta\\\
\alpha&-\beta\end{pmatrix}\begin{pmatrix}f_{11}&f_{12}\\\ f_{21}&f_{22}\\\
\end{pmatrix}$

with the additional conditions that each column of $F$ sums to zero, i.e.
$f_{11}+f_{21}=0$ and $f_{12}+f_{22}=0$. This yields:

$\displaystyle F=\frac{1}{(\alpha+\beta)^{2}}\begin{pmatrix}\alpha&-\beta\\\
-\alpha&\beta\end{pmatrix}$ (17)

and allows to determine the noise intensity using eq. (12):

$\displaystyle D$
$\displaystyle=\frac{1}{(\alpha+\beta)^{3}}\begin{pmatrix}x_{1}&x_{2}\end{pmatrix}\begin{pmatrix}\alpha&-\beta\\\
-\alpha&\beta\end{pmatrix}\begin{pmatrix}x_{1}\beta\\\ x_{2}\alpha\\\
\end{pmatrix}$ (18)
$\displaystyle=\frac{\alpha\beta}{(\alpha+\beta)^{3}}(x_{1}-x_{2})^{2}$

an expression that is in agreement with the result presented in HorLef83 .

Figure 2: The noise intensity of a Markovian Dichotomous noise. Panels A and B
show the noise intensity as a function of the ratio between the transition
rates $\alpha/\beta$ or the correlation time $\tau=1/(\alpha+\beta)$. The
intensity has a maximum at $\alpha^{*}=\beta/2$ or $\tau^{*}=2/(3\beta)$.
Panels C and D show the noise intensity as a function of the ratio between the
values taken in the two states $x_{1}/x_{2}$ or the variance
$\sigma^{2}=\alpha\beta(x_{1}-x_{2})^{2}/(\alpha+\beta)^{2}$. The intensity
has a minimum at $x_{1}^{*}=x_{2}$ and scales linearly with the variance.
Parameters: $\beta=1$, $x_{2}=1$.

The noise intensity has a maximum as a function of the rate $\alpha$, keeping
the other rate $\beta$ fixed, at $\alpha^{*}=\beta/2$ (Fig. 2A). As a function
of the correlation time $\tau=1/(\alpha+\beta)$ has also a maximum at
$\tau^{*}=2/(3\beta)$ (Fig. 2B). If the intensity is plotted as a function of
the ratio $x_{1}/x_{2}$ (Fig. 2C), it has a minimum and vanishes for
$x_{1}^{*}=x_{2}$ because the variance vanishes in this case. When plotted as
a function of the variance $\sigma^{2}$ (Fig. 2D), the noise intensity
increases linearly according to eq. (4).

The calculation presented here serves only as a sanity check. The real
advantage of the method lies in the possibility of calculating the noise
intensity for more complicated transition rate matrices, as we will show in
the following.

##  IV A biophysical example:
Stochastic Ca2+ channel model

In this section, we consider a biophysical example and calculate the noise
intensity for a eight-state Markov model as illustrated in Fig. 3 and used in
the De Young-Keizer model to describe a single subunit of an inositol
trisphosphate (IP3) receptor YouKei92 . For such a model, no closed-form
expression for the noise intensity is known.

The entire De Young-Keizer model describes the dynamics of the intracellular
calcium (Ca2+) concentration, which in many cells serves as a signaling
molecule to transmit information about extracellular stimuli (calcium
signaling) BerBoo98 ; Cla07 . The characteristic short periodic increases in
the intracellular Ca2+ concentration that carry the information can be caused
either by an influx of Ca2+ from the extracellular medium or by a release of
Ca2+ from an intracellular store, the endoplasmic reticulum (ER). In both
cases, stochastic transitions between discrete states of the ion channels give
rise to macroscopic fluctuations in the intracellular Ca2+ concentration. The
De Young-Keizer model covers the case where the Ca2+ signal is evoked by the
release of Ca2+ from the ER through the IP3 receptor channel. This receptor
channel in turn is assumed to consist of three independent and identical
subunits with three binding sites each: one for the second-messenger molecule
IP3, produced in the cell in response to an extracellular stimulus (IP3
pathway), and one activating Ca2+ binding site and one inhibiting Ca2+ binding
site. The Ca2+ current through a single IP3 receptor channel is given by

$\displaystyle
I_{\text{Ca}}=c_{1}x^{(1)}(t)x^{(2)}(t)x^{(3)}(t)([\text{Ca}\textsuperscript{2+}]_{\text{er}}-[\text{Ca}\textsuperscript{2+}]_{\text{i}}),$
(19)

where $[\text{Ca}\textsuperscript{2+}]_{\text{er}}$ and
$[\text{Ca}\textsuperscript{2+}]_{\text{i}}$ are the ER and intracellular
(cytosolic) Ca2+ concentrations, $c_{1}$ is the volume ratio between the ER
and the cytosol, and $x^{(n)}(t)$ are dichotomous (two-valued) stochastic
processes capturing the state of the three IP3 receptor subunits.

Figure 3: State diagram of a Ca2+ channel subunit YouKei92 . Panel A shows the eight-state model of a single IP3-receptor subunit. The states are denoted $ijk$, where each index represents one of the three binding sites for IP3 ($i$), activating Ca2+ ($j$), and inhibitory Ca2+ ($k$). An index is $1$ ($0$) if the binding site is occupied (unoccupied). The conducting state $110$ is highlighted in green. Panel B shows the transition rates on the front and back of the die. Panel C shows the transitions between the front and back faces. Binding rates are denoted $\alpha$ and depend linearly on the corresponding concentration ($\alpha_{i}=\hat{\alpha}_{i}[\text{IP}\textsubscript{3}]$ for $i=1,3$ and $\alpha_{i}=\hat{\alpha}_{i}[\text{Ca}\textsuperscript{2+}]_{\text{i}}$ for $i=2,4,5$), while unbinding rates are denoted $\beta$ and are constants. Parameters are according to Table 1. Table 1: Simulation parameters for a IP3-receptor subunit in the De Young-Keizer model YouKei92  Parameter | Value | Description
---|---|---
binding constants
$\hat{\alpha}_{1}$ / $\mu\text{M}^{-1}\text{s}^{-1}$ | 400 |  IP3
$\hat{\alpha}_{2}$ / $\mu\text{M}^{-1}\text{s}^{-1}$ | 0.2 |  Ca2+ inhibition
$\hat{\alpha}_{3}$ / $\mu\text{M}^{-1}\text{s}^{-1}$ | 400 |  IP3
$\hat{\alpha}_{4}$ / $\mu\text{M}^{-1}\text{s}^{-1}$ | 0.2 |  Ca2+ inhibition
$\hat{\alpha}_{5}$ / $\mu\text{M}^{-1}\text{s}^{-1}$ | 20 |  Ca2+ activation
dissociation constants $\gamma_{i}=\beta_{i}/\hat{\alpha}_{i}$
$\gamma_{1}$ / $\mu\text{M}$ | 0.13 |  IP3
$\gamma_{2}$ / $\mu\text{M}$ | 1.049 |  Ca2+ inhibition
$\gamma_{3}$ / nM | 943.4 |  IP3
$\gamma_{4}$ / nM | 144.5 |  Ca2+ inhibition
$\gamma_{5}$ / nM | 82.34 |  Ca2+ activation

The kinetics of a single subunit is described by the scheme shown in Fig. 3.
The eight possible states shown in Fig. 3A result from the fact that each of
the three binding sites can be in two possible states, occupied or unoccupied.
The subunit states are labeled $ijk$, where $i$, $j$, and $k$ indicate whether
the IP3, activating Ca2+, and inhibitory Ca2+ binding sites are occupied
($i,j,k=1$) or unoccupied ($i,j,k=0$). The entire channel is open when in all
subunits the IP3 and activating Ca2+ binding sites are occupied and the
inhibitory Ca2+ binding site is unoccupied, i.e. when all three subunits are
in the $110$ state, highlighted in green in Fig. 3A. Put differently, every
state is assigned a value according to
$x_{ijk}=\delta_{i1}\delta_{j1}\delta_{k0}$, i.e. the value of the conducting
state $110$ is $1$ and the value of every other state is $0$. The transition
rates between the states on the front and back faces of the cube are shown in
Fig. 3B, while the transition rates between the two faces are shown in Fig.
3C. Binding rates are denoted $\alpha$ and depend linearly on the IP3 or Ca2+
concentration according to the law of mass action
($\alpha_{i}=\hat{\alpha}_{i}[\text{IP}\textsubscript{3}]$ for $i=1,3$ and
$\alpha_{i}=\hat{\alpha}_{i}[\text{Ca}\textsuperscript{2+}]_{\text{i}}$ for
$i=2,4,5$), whereas unbinding rates are denoted $\beta$ and are constants (we
keep the original notation of De Young and Keizer in terms of $\hat{\alpha}$
and $\gamma=\beta/\hat{\alpha}$, see Table 1).

While De Young and Keizer considered the IP3 receptor and its subunits in the
thermodynamic limit, we calculate the variance $\sigma^{2}$, correlation time
$\tau$ and noise intensity $D$ for a single subunit $x^{(n)}(t)$. Although the
noise intensity for the Markov chain can be calculated analytically, the
expressions for the stationary probability vector $\bm{p}$ and the auxiliary
matrix $F$ are lengthy. Therefore, we determine these two statistics
numerically by inverting the transition rate matrix $W$. Since $W$ does not
have full rank, this requires some manipulation, that we mentioned already in
sec. (II) and are more detailed below.

Figure 4: Statistical measures of a stochastic Ca2+ channel subunit. Panels A,
B, and C show the variance $\sigma^{2}$, correlation time $\tau$, and noise
intensity $D$ as a function of the intracellular calcium concenntration
$[\text{Ca}\textsuperscript{2+}]_{\text{i}}$ for a stochastic IP3 receptor
subunit governed by the scheme illustrated in Fig. 3. Vertical lines indicate
the standard error calculated from ten simulations. The variance and noise
intensity are calculated according to eq. (7) and eq. (12), respectively. The
correlation time is determined as the ratio $\tau=D/\sigma^{2}$. Parameters
are according to Table 1.

To compute the stationary probability vector $\bm{p}$, we implement the
normalization condition by replacing all entries in an arbitrary row of $W$
with ones, and the corresponding entry in the zero vector on the l.h.s. of the
stationary master equation (eq. (6)) with a one. This removes a redundant row
in $W$, which can be obtained by linear combination of the other rows, and
replaces it with the normalization condition $\sum_{i}p_{i}=1$. The same trick
is used to compute the auxiliary matrix $F$, again replacing all entries in a
row of $W$ with ones and the corresponding row of the matrix
$P_{0}-\mathbb{1}$ on the l.h.s. of eq. (13) with zeros. This satisfies the
conditions $\sum_{i}f_{ij}=0$.

The results for three statistics are shown in Fig. 4. In all cases, the
numerically calculated values show excellent agreement with the theoretical
predictions, demonstrating that the method is applicable even when the
transition rate matrix is more complicated. Furthermore, the results show that
the variance alone is an insufficient measure to quantify the effect of a
random process on a driven variable. For example, while the variance is nearly
constant for low values of $[\text{Ca}\textsuperscript{2+}]_{\text{i}}$, the
noise intensity shows a pronounced maximum for an intermediate value.

##  V Voltage-gated channels:
Models of stochastic K+ and Na+ channels

As a third and final example, we calculate the noise intensity for two Markov
chains, as used in stochastic variants of the Hodgkin-Huxley model to describe
the gating of the potassium (K+) and sodium (Na+) channels FoxLu94 ; Fox97 ;
Koc99 . Similar to the example discussed in the previous section, microscopic
transitions between different discrete states of the ion channels lead to
stochastic ion currents and eventually to macroscopic fluctuations, here in
the voltage of an excitable membrane. We emphasize that we now consider the
characteristics of the current through an entire ion channel (K+ or Na+), in
contrast to the subunit activity addressed in the previous section.

The classical Hodgkin-Huxley model describes the dynamics of the membrane
potential $V$ and the generation of an action potential in a neuron by means
of a passive leak current, a voltage-dependent K+ current, and a voltage-
dependent Na+ current Izh07 ; GerKis14 . In a stochastic formulation the
latter two currents can be expressed by:

$\displaystyle I_{\text{K}}$
$\displaystyle=g_{\text{K}}n^{(1)}(t)n^{(2)}(t)n^{(3)}(t)n^{(4)}(t)(V-E_{\text{K}})$
(20) $\displaystyle I_{\text{Na}}$
$\displaystyle=g_{\text{Na}}m^{(1)}(t)m^{(2)}(t)m^{(3)}(t)h^{(1)}(t)(V-E_{\text{Na}})$

where $g_{\text{K}}$ and $g_{\text{Na}}$ are the maximal conductances and
$E_{\text{K}}$ and $E_{\text{Na}}$ are the reversal potentials. In our case,
the variables $n^{(i)}$, $m^{(j)}$, and $h^{(k)}$ are Markovian dichotomous
processes that capture the state of the subunits in a single K+ or Na+ channel
111The state is described as activated or deactivated for $n$ and $m$, and
inactivated or deinactivated for $h$ DayAbb01 . Accordingly to eq. (20), the
K+ channel consists of four activation gates of type $n$, whereas the Na+
channel consists of three activation gates of type $m$ and one inactivation
gate of type $h$. Only when all subunits are open, the channel is open.

In the original Hodgkin-Huxley model, the gating variables are deterministic
quantities bounded between zero and one and governed by the differential
equation:

$\displaystyle\tau_{x}(V)\dot{x}=\alpha_{x}(V)(1-x)-\beta_{x}(V)x,$ (21)

with $x=n,m,h$. In this case, $n$, $m$, and $h$ represent the fraction of open
subunits in a large ensemble. However, eq. (21) can also be interpreted as
$x(t)$ describing the probability of a single subunit taking a certain state
in a two-state system with voltage-dependent transition rates $\alpha_{x}(V)$
and $\beta_{x}(V)$ (similar to Fig. 1). This insight allows to formulate
stochastic variants of the Hodgkin-Huxley model consistent with the
deterministic model in the thermodynamic limit, where the state of each
subunit is represented by a Markovian dichotomous noise with a mean governed
by eq. (21) FoxLu94 ; Fox97 . In this formulation, the gating variables
correspond to the fluctuating fraction of open subunits in a finite ensemble.

Figure 5: State diagram of a K+ channel and a Na+ channel Koc99 . Panel A shows the five-state model of a stochastic K+ channel. The five states represent the number of activated $n$-type subunits (0 to 4) of the K+ channel. The transition rates are given above/below the arrows. All four subunits must be activated ($n_{4}$) for the K+ channel to open. Panel B shows the eight-state model of a stochastic Na+ channel. The eight states represent the number of activated $m$-type (0 to 3) and deinactivated $h$-type (0 to 1) subunits of the Na+ channel. All three subunits of type $m$ must be activated and the subunit of type $h$ must be deinactivated ($m_{3}h_{1}$) for the Na+ channel to open. Table 2: Simulation parameters for K+ and Na+ channels in the Hodgkin-Huxley model DayAbb01  Parameter | Value
---|---
$\alpha_{n}$ / $\text{s}^{-1}$ | $0.01(V+55)/[1-\exp(-0.1(V+55))]$
$\alpha_{m}$ / $\text{s}^{-1}$ | $0.1(V+40)/[1-\exp(-0.1(V+40))]$
$\alpha_{h}$ / $\text{s}^{-1}$ | $0.07\exp(-0.05(V+65))$
$\beta_{n}$ / $\text{s}^{-1}$ | $0.125\exp(-0.0125(V+65))$
$\beta_{m}$ / $\text{s}^{-1}$ | $4\exp(-0.0556(V+65))$
$\beta_{h}$ / $\text{s}^{-1}$ | $1/[1+\exp(-0.1(V+35))]$
$g_{\text{K}}$ / mS | 36
$g_{\text{Na}}$ / mS | 120
$E_{\text{K}}$ / mV | -77
$E_{\text{Na}}$ / mV | 50

In the following we calculate the noise intensity of the ion current through a
single K+ channel or a single Na+ channel. We have already emphasized that the
kinetics of a single subunit ($n^{(i)}(t)$, $m^{(j)}(t)$, and $h^{(k)}(t)$)
can be described by a Markovian dichotomous noise with a transition rate
matrix similar to the one used in sec. (III). One could be tempted to think
that the noise intensity of the product of a number of independent random
processes can be easily found from the intensities of the single processes.
However, we are not aware of a simple relation between the former and the
latter.

Figure 6: Statistical measures of a stochastic K+ and Na+ channel. Panels A
and B show the open probability, variance $\sigma^{2}$, correlation time
$\tau$, and noise intensity $D$ as a function of the membrane potential $V$
for the stochastic K+ and Na+ channels governed by the scheme illustrated in
Fig. 5 A and B, respectively. Vertical lines indicate the standard error
calculated from ten simulations. The variance and noise intensity are
calculated according to eq. (7) and eq. (12), respectively. The correlation
time is determined as the ratio $\tau=D/\sigma^{2}$. Parameters are according
to Table 2.

To calculate the noise intensity for the product, we need to formulate the
transition rate matrix for the random processes
$x(t)=n^{(1)}(t)n^{(2)}(t)n^{(3)}(t)n^{(4)}(t)$ for the K+ channel or
$x(t)=m^{(1)}(t)m^{(2)}(t)m^{(3)}(t)h^{(1)}(t)$ for the Na+ channel. Here, we
follow the formulation of Koc99  and use transition rate matrices
corresponding to the state diagrams in Fig. 5. The five possible states for
the K+ channel (Fig. 5A) result from the fact that the subunits are identical
and independent. Therefore, it is sufficient to describe the number of
subunits in the activated state. In this formulation, the transition rate from
the state $n_{3}$ (three activated $n$-type subunits) to $n_{4}$ is
$\alpha_{n}$, the rate at which the last deactivated gate is activated, and
the transition from $n_{3}$ to $n_{2}$ is $3\beta_{n}$, the rate at which one
out of three activated gates deactivates (see Koc99 ; DayAbb01 ). The entire
K+ channel is considered open when all gates are activated, i.e. the value of
the state $n_{4}$ is $1$ and the value of every other state is $0$. Similarly,
a reduced state diagram can be formulated for the Na+ channel (Fig. 5B). In
this case, the number of activated $m$-type subunits and number of inactivated
$h$-type subunit must be distinguished, resulting in eight different states.
The K+ channel is considered open, when all three $m$-type subunits are
activated and the $h$-type subunits deinactivated, i.e. the value of the state
$m_{3}h_{1}$ is $1$ and the value of every other state is $0$.

In Fig. 6 we compare simulation results and theoretical predictions of the
open probability, variance, correlation time, and noise intensity of the
stochastic K+ and Na+ currents according to eq. (20) where the gating
variables are described by the Markov schemes illustrated in Fig. 5A and 5B,
respectively. In both cases, the numerically calculations agree with the
theoretical predictions, demonstrating that the method is applicable to the
stochastic Hodgkin-Huxley model. However, we note that our method relies on
the assumption of a clamped voltage.

Regarding the interpretation of the obtained curves, we first note that the
two upper panels agree with the deterministic open probability of the
classical Hodgkin-Huxley model: A monotonically increasing function for the K+
channel (Fig. 6A1) and a non-monotonic function for the Na+ channel (Fig. 6B1)
due to the interplay between activation and inactivation. The latter maximum
implies maxima in the variance (Fig. 6B2) and the noise intensity (Fig. 6B4).
We note that there are also maxima in the characteristics of the K+ channel at
different voltage values (Fig. 6A2-A4). The maximum of the variance is
plausible because the open probability reaches zero and one in the limit of
extreme voltage values. This maximum then also entails a maximum of the noise
intensity.

##  VI Summary and discussion

In this paper, we have developed a general framework to characterize a noise
process that is described by a finite Markov chain, i.e. by a Master equation
with a finite number of states. More specifically, we demonstrated that the
calculation of the noise intensity and correlation time of the process is only
slightly more complicated than the computation of the steady state and its
mean and variance. We illustrated our general result by application to three
cases: (i) the dichotomous noise (for which all characteristics are, of
course, well known); (ii) a stochastic calcium channel subunit as it is used
in the De Young-Keizer model; (iii) sodium and potassium currents as used in
the Hodgkin-Huxley equation of action potential generation. In all these
cases, our comparison to stochastic simulations of the underlying discrete
dynamics agreed well with the analytical predictions of our formulas over a
wide range of tested parameters.

The computation of noise intensity and correlation time has particular
importance in the context of the so-called diffusion approximation. In several
situations of interest the discrete fluctuations described by the Master
equation can be well approximated by a white Gaussian noise. This stochastic
process is fully characterized by its mean value and its noise intensity, for
which we derived a simple expression above. Once this approximation has been
made, the apparatus of nonlinear diffusion processes, in particular, the
Fokker-Planck equation for the evolution of the probability density, can be
used. In order to learn whether this approximation is really justified for a
specific system, it is crucial to know the correlation time of the noise and
to test whether it is much shorter than all other time scales in the system -
only if this is the case, we are permitted to neglect the temporal
correlations of noise entirely. This may also apply in the more complicated
situation in which both the mean and intensity depend on the dynamical
variable(s) of the driven system, i.e. when there is a feedback between the
dynamical variable(s) and the noise statistics.

Let us revisit the dynamics of the calcium subunit, for which the correlation
time was shown in Fig. 4B as a function of the (clamped) calcium
concentration. The maximum correlation time is below 3 seconds in this case.
If we now take into account that calcium is _not_ clamped but in fact obeys a
dynamics on the time scale of tens of seconds to multiple minutes BouMar08 ;
ThuSmi11 , we may justify to approximate the stochastic activity of the single
subunit by a Gaussian white noise with a calcium-dependent mean value and a
calcium-dependent noise intensity. This is true when the Ca2+ concentration is
below some spiking threshold, and it does not include the dynamics that is
responsible for the spike shape. For another calcium channel (cluster) model,
this has been carried out in a integrate-and-fire type model of intracellular
calcium excitability and thoroughly tested by us and a collaborator RamFal23 ;
RamFal23b .

For the gating variables of the Hodgkin-Huxley model a similar argument may be
possible. A naive version of the approximation is not justified to describe
the generation of the action potential. It is exactly the nonlinear interplay
between the voltage and gating dynamics that gives the action potential its
characteristic shape, i.e. the upstroke of the spike caused by the positive
feedback of sodium-channel opening upon an initial (and sufficiently strong)
depolarization and the downstroke due to the slower inactivation of sodium
channels and the opening of potassium channels. The membrane time constant in
the original Hodgkin-Huxley model (roughly, the time scale of the voltage
dynamics) is of the order of 3ms HodHux52  and thus comparable to the
correlation time of the potassium channel fluctuations (according to Fig. 6b
around 2.5ms for voltage values around the resting potential). Hence, in this
case it is recommendable to abstain from a white-noise approximation. Indeed,
different approximation schemes that are based on having a large number of
channels, have been divised, see e.g. the classical studies by Fox et al.
FoxLu94 ; Fox97 , who approximated the gating dynamics by chemical Langevin
equations, and more recent contributions which use stochastic shielding to
obtain numerically efficient descriptions of the inherent stochasticity
PuTho21 . We note that our results are still useful, because in experiment the
voltage _can and is routinely_ clamped to a prescribed value and currents
through specific channels can be isolated (methods for this are for instance
discussed in the textbook by Izhikevich Izh07 ), and in this situation our
formulas give exact results for the characteristics of the respective current
fluctuations. The same method can be applied to more complicated kinetic
schemes of channel states, see e.g. the review on the many models of sodium
channels Pat91 .

Of course, the two above cases are more involved in the sense that not always
the dynamics of the Markov chain itself is affected by the variable it drives.
When the output of the Markov chain acts as an external noise on a system, no
white-noise approximation has to be made and our results then simply provide
the most important noise characteristics of this stochastic process, making it
comparable to simpler noise models (white or low-pass filtered Gaussian noise,
white Poissonian noise, or colored dichotomous noise).

## References

  * (1) J. A. White, J. T. Rubinstein, and A. R. Kay.  Channel noise in neurons.  Trends Neurosci., page 131, 2000.
  * (2) A. D. Dorval and J. A. White.  Channel noise is essential for perithreshold oscillations in entorhinal stellate neurons.  J. Neurosci., 25:10025, 2005.
  * (3) T. Schwalger, K. Fisch, J. Benda, and B. Lindner.  How noisy adaptation of neurons shapes interspike interval histograms and correlations.  PLoS Comp. Biol., 6:e1001026, 2010.
  * (4) K. Fisch, T. Schwalger, B. Lindner, A. Herz, and J. Benda.  Channel noise from both slow adaptation currents and fast currents is required to explain spike-response variability in a sensory neuron.  J. Neurosci., 32:17332, 2012.
  * (5) B. Moezzi, N. Iannella, and M. D. McDonnell.  Ion channel noise can explain firing correlation in auditory nerves.  J. Comput. Neurosci., 41:193, 2016.
  * (6) S. Pu and P. J. Thomas.  Resolving molecular contributions of ion channel noise to interspike interval variability through stochastic shielding.  Biol. Cybern., 115(3):267–302, 2021.
  * (7) N. Hohn and A. N. Burkitt.  Shot noise in the leaky integrate-and-fire neuron.  Phys. Rev. E., 63:031902, 2001.
  * (8) M. J. E. Richardson and W. Gerstner.  Statistics of subthreshold neuronal voltage fluctuations due to conductance-based synaptic shot noise.  Chaos, 16:026106, 2006.
  * (9) L. Wolff and B. Lindner.  Mean, variance, and autocorrelation of subthreshold potential fluctuations driven by filtered conductance shot-noise.  Neural Comput., 22:94, 2010.
  * (10) M. J. E. Richardson and R. Swarbrick.  Firing-rate response of a neuron receiving excitatory and inhibitory synaptic shot noise.  Phys. Rev. Lett., 105:178102, 2010.
  * (11) M. Brigham and A. Destexhe.  Nonstationary filtered shot-noise processes and applications to neuronal membranes.  Phys. Rev. E, 91:062102, 2015.
  * (12) F. Droste and B. Lindner.  Exact analytical results for integrate-and-fire neurons driven by excitatory shot noise.  J. Comp. Neurosci., 43:81, 2017.
  * (13) H. Haken, H. Sauermann, C. Schmid, and H. D. Vollmer.  Theory of laser noise in the phase locking region.  Z. Physik, 206:369, 1967.
  * (14) J. L. A. Dubbeldam, B. Krauskopf, and D. Lenstra.  Excitability and coherence resonance in lasers with saturable absorber.  Phys. Rev. E., 60:6580, 1999.
  * (15) W. S. Lam, P. N. Guzdar, and R. Roy.  Effect of spontaneous emission noise and modulation on semiconductor lasers near threshold with optical feedback.  Int J Mod Phys B, 17:4123, 2003.
  * (16) T. Schwalger, J. Tiana-Alsina, M. C. Torrent, J. Garcia-Ojalvo, and B. Lindner.  Interspike-interval correlations induced by two-state switching in an excitable system.  Epl-Europhys. Lett., 99:10004, 2012.
  * (17) C. W. Gardiner.  Handbook of Stochastic Methods.  Springer-Verlag, Berlin, 1985.
  * (18) N. G. van Kampen.  Stochastic Processes in Physics and Chemistry.  North-Holland, Amsterdam, 1992.
  * (19) D. T. Gillespie.  The chemical langevin equation.  J. Chem. Phys., 113:297, 2000.
  * (20) G. E. Uhlenbeck and L. S. Ornstein.  On the theory of the Brownian motion.  Phys. Rev., 36:823, 1930.
  * (21) H. Risken.  The Fokker-Planck Equation.  Springer, Berlin, 1984.
  * (22) P. Hänggi and P. Jung.  Colored noise in dynamical-systems.  Adv. Chem. Phys., 89:239, 1995.
  * (23) W. Horsthemke and R. Lefever.  Noise-Induced Transitions: Theory and Applications in Physics, Chemistry, and Biology.  Springer, Berlin, 1983.
  * (24) I. Bena.  Dichotomous markov noise: Exact results for out-of-equilibrium systems.  Int. J. Mod. Phys. B., 20:2825, 2006.
  * (25) F. Droste and B. Lindner.  Integrate-and-fire neurons driven by asymmetric dichotomous noise.  Biol. Cybern., 108:825, 2014.
  * (26) G. N. Farah and B. Lindner.  Exponentially distributed noise - its correlation function and its effect on nonlinear dynamics.  J. Phys. A: Math. Theor., 54:035003, 2021.
  * (27) H. Risken.  The Fokker-Planck Equation (second edition).  Springer, Berlin, 1989.
  * (28) B. Lindner.  The diffusion coefficient of nonlinear Brownian motion.  New J. Phys., 9:136, 2007.
  * (29) B. Lindner.  Diffusion coefficient of a Brownian particle with a friction function given by a power law.  J. Stat. Phys., 130:523, 2008.
  * (30) L. Ramlow, M. Falcke, and B. Lindner.  An integrate-and-fire approach to Ca2+ signaling. Part I: Renewal model.  Biophys. J., 122:713, 2023.
  * (31) L. Ramlow, M. Falcke, and B. Lindner.  An integrate-and-fire approach to Ca2+ signaling. Part II: Cumulative refractoriness.  Biophys. J., 122:4710, 2023.
  * (32) G. W. De Young and J. Keizer.  A single-pool inositol 1,4,5-trisphosphate-receptor-based model for agonist-stimulated oscillations in Ca2+ concentration.  PNAS, 89(20):9895–9899, 1992.
  * (33) M. J. Berridge, M. D. Bootman, and P. Lipp.  Calcium - a life and death signal.  Nature, 395:645, 1998.
  * (34) D. E. Clapham.  Calcium signaling.  Cell, 131(6):1047–1058, 2007.
  * (35) R. F. Fox and Y. N. LU.  Emergent collective behavior in large numbers of globally coupled independently stochastic ion channels.  Phys. Rev. E., 49:3421, 1994.
  * (36) R. F. Fox.  Stochastic versions of the Hodgkin-Huxley equations.  Biophys. J., 72:2068, 1997.
  * (37) C. Koch.  Biophysics of Computation - Information Processing in Single Neurons.  Oxford University Press, New York, Oxford, 1999.
  * (38) E. M. Izhikevich.  Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting.  The MIT Press, Cambridge, London, 2007.
  * (39) W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski.  Neuronal Dynamics From single neurons to networks and models of cognition.  Cambridge University Press, Cambridge, 2014.
  * (40) The state is described as activated or deactivated for $n$ and $m$, and inactivated or deinactivated for $h$ DayAbb01 .
  * (41) P. Dayan and L. F. Abbott.  Theoretical Neuroscience.  MIT Press, Cambridge MA, 2001.
  * (42) M.J. Boulware and J.S. Marchant.  Timing in cellular Ca2+ signaling.  Curr. Biol., 18(17):R769–R776, 2008.
  * (43) K. Thurley, I. F. Smith, S. C. Tovey, C. W. Taylor, I. Parker, and M. Falcke.  Timescales of ip3-evoked Ca2+ spikes emerge from Ca2+ puffs only at the cellular level.  Biophys. J., 101(11):2638–2644, 2011.
  * (44) A. L. Hodgkin and A. F. Huxley.  A quantitative description of membrane current and its application to conduction and excitation in nerve.  J. Physiol. (Lond.), 10:500, 1952.
  * (45) J. Patlak.  Molecular Kinetics of Voltage-Dependent Na+ Channels.  Physiol. Rev., 71:1047, 1991.