Potassium channels are integral proteins that enable rapid, selective transport of K⁺ ions across the cell membrane down their electrochemical gradient (Chung et al. ²⁰⁰⁷). In general, ion channels are not just simple pores with a constant permeability—they undergo conformational changes resulting in transitions between the conducting and non-conducting (open and closed) states. The probabilities of these states depend on channel-specific gating stimuli, of which the most common are membrane potential, ligand binding, and mechanical force (Hille ²⁰⁰¹).

Ion-channel proteins comprise several subunits with different functions. In particular, distinct structural domains are responsible for the stimulus sensing and for the channel activation. Thus, stimulus–response relationships require a cooperative interaction between the sensor domains and the gate. When a gating stimulus occurs at the sensor domain, the conformational change is conveyed to the activation gate and changes the transition probability between the open and closed states.

Voltage-dependent, calcium-activated potassium channels (BK) are characterized by a large single-channel conductance (~100–300 pS) (Latorre and Miller ¹⁹⁸³; Latorre et al. ¹⁹⁸⁹; Marty ¹⁹⁸³) compared with that of other potassium channels. These channels, which can be activated by membrane depolarization or by elevation of intracellular calcium concentration, can be found in many cells and tissues, for example neurons, chromaffin cells, the inner hair cells of cochlea, and muscles, and have an important function in several physiological processes (Cui et al. ²⁰⁰⁸). Irrespective of their unusually large conductance, BK channels remain highly selective for K⁺ ions over other cations (Cui et al. ²⁰⁰⁸).

In general, the putative structure of BK channels shares many similarities with that of voltage-dependent Kv channels (Ma et al. ²⁰⁰⁶; Liu et al. ^2008a, ^b) and ligand gated potassium channels, for example MthK (Latorre and Brauchi ²⁰⁰⁶; Jiang et al. ²⁰⁰¹, ²⁰⁰²; Shi et al. ²⁰⁰²; Xia et al. ²⁰⁰²; Tang et al. ²⁰⁰³; Yusifov et al. ²⁰⁰⁸; Kim et al. ²⁰⁰⁶, ²⁰⁰⁸). BK channels are tetramers of a channel protein encoded by the Slo1 gene, composed of seven transmembrane domains (S0–S6) and four cytosolic hydrophobic C-terminal domains (S7–S10) (Cox ²⁰⁰⁶; Cui et al. ²⁰⁰⁸; Latorre and Brauchi ²⁰⁰⁶). Functionally, segments S1–S4 form a voltage sensor domain (VSD), and segments S5–S6 form a pore-gate domain (PGD). The large cytoplasmic C-terminal domain serves as the primary ligand sensor. S0 is a segment that is absent from Kv and MthK channels and causes some problems in homology modelling, directing the N-terminus to the extracellular side of the membrane (Cui et al. ²⁰⁰⁸).

The simplest model of the channel’s gating kinetics can be described by a two-state system:

[Formula ID: Equ1]
1  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \begin{array}{*{20}c} C & {\begin{array}{*{20}c} {k_{\text{O}} } \\ \to \\ \leftarrow \\ {k_{\text{C}} } \\ \end{array} } & O \\ \end{array} $$\end{document}]

Invention of the patch-clamp technique (Sakmann and Neher ¹⁹⁹⁵) has enabled experimental observation of the ionic currents at a single-channel level, which gave the opportunity to obtain the probability densities of closed f_C(t) and open f_O(t) dwell time intervals. For many channels, the open dwell-time distribution can be reasonably approximated by a single exponential function, but the closed dwell-time distribution fits better to the sum of many exponentials of the form (Goychuk and Hänggi ²⁰⁰², ²⁰⁰³):

[Formula ID: Equ2]
2  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ f_{C} (t) = \sum\limits_{i = 1}^{N} {c_{i} \lambda_{i} \exp ( - \lambda_{i} t)} , $$\end{document}]

[Formula ID: Equ3]
3  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sum\limits_{i = 1}^{N} {c_{i} = 1} , $$\end{document}]

Quite often Eq. (2) can be replaced by a single stretched exponential (Millhauser et al. ¹⁹⁸⁸) or a power law function (Sansom et al. ¹⁹⁸⁹; Blatz and Magleby ¹⁹⁸⁶; Ring ¹⁹⁸⁶; Mercik and Weron ²⁰⁰¹; Läuger ¹⁹⁸⁸; Condat and Jäckle ¹⁹⁸⁹). When the dwell time distribution can be described by a sum of exponentials, Markov models with either a few or many open and closed states have been proposed. Such models assume the channel protein to be found in a number of discrete states, separated by relatively high potential barriers (Sansom et al. ¹⁹⁸⁹).

When multiple open and closed states are fitted to the model, a possibility arises for allosteric activation of the channel, as shown for BK channels as an example in Horrigan et al. (¹⁹⁹⁹), Horrigan and Aldrich (¹⁹⁹⁹), and Rothberg and Magleby (¹⁹⁹⁹). The HCA model proposed by Horrigan et al. (¹⁹⁹⁹) and Horrigan and Aldrich (¹⁹⁹⁹) is a classical, many-state Markov model of BK channel voltage-dependent gating. It reflects well the channel’s behaviour at different values of membrane potential in the absence of any calcium ions and agrees quantitatively with experimental results. A crucial feature of the HCA model is the allosteric mechanism of BK channel activation, which is different from the strict coupling proposed in the literature to describe gating of the Kv channel (an example of which is the Hodgkin–Huxley model) (Hille ²⁰⁰¹). Whereas in the strict coupling of Kv channel activation voltage sensor movement opens a channel, this is not so in allosteric activation. In the allosteric model (of BK channels) neither of the sensors directly opens the gate but all of them affect the gate, and each other, in an allosteric manner, and gate opening can even precede the sensor’s activation (Rothberg and Magleby ¹⁹⁹⁹).

The calcium activation of BK channels is characterized in the voltage-dependent Monod Wyman Changeux (VD-MWC) model (Cox et al. ¹⁹⁹⁷). Analogously to the previous case, because of incorporation of an allosteric mechanism, none of the Ca²⁺-activated sensor domains renders the channel conductive in a direct manner. Nevertheless, responding to the ligand binding on any of the multiple binding sites, the channel undergoes a conformational change that promotes the opening, and (extending the MWC formalism) indirectly alters the remaining binding sites by increasing the affinity for Ca²⁺. In this model, the open state has a higher Ca²⁺ affinity than a corresponding closed state with the same number of calcium ions bound. Thus, an open Ca²⁺-bound state is energetically preferred over a closed Ca²⁺-bound state. The voltage dependence of the gating is embedded in the equilibrium constants between the closed and open states.

Advanced approaches based on the HCA and MWC models have been proposed in the literature to account for both allosteric activation processes in a single model (Rothberg and Magleby ²⁰⁰⁰; Horrigan and Aldrich ²⁰⁰²). By extending the assumptions of its templates, these models reach a large number of states (e.g. 50-state two-tiered model); their purpose is to give an even more complete picture of the gating process over a wide range of voltages and calcium concentrations, assuming all the rate constants are correctly determined (which is not easy). What is important here it is that coupling between different functional parts of the channel may be strong, as for the sensor–gate interaction, or weak, as for the voltage sensor–calcium sensor interaction. An allosteric model allows the channel to open even in the absence of any voltage and calcium stimuli or to close when all the sensors are activated. In general, however, the greater the number of sensors activated the more energetically favoured an open conformation is. Both the voltage and the ligand sensors can switch from the resting to the activated state in either closed or open conformation of a channel.

The multi-exponential dwell-time distribution of the Markovian systems in some cases reduces to a single stretched exponential or a power law dependence. This observation leads to the introduction of the fractal (Liebovitch et al. ¹⁹⁸⁷) and diffusive models (Millhauser et al. ¹⁹⁸⁸; Läuger ¹⁹⁸⁸; Condat and Jäckle ¹⁹⁸⁹; Goychuk and Hänggi ²⁰⁰², ²⁰⁰³) to address the voltage-dependent channels (Kv). According to these models, the channel protein exists in a continuous conformational space with a rather smooth energy surface, in contrast with a few, well separated discrete conformations determined by the potential wells, as expected from the finite state Markovian models.

A classical model of channel gating by a Markovian diffusion process over a very (infinitely) large number of similar energy closed states was introduced by Millhauser et al. (¹⁹⁸⁸). This model produces the power law dependence in the dwell time distributions. It was generalised by Hänggi and Goychuk (²⁰⁰², ²⁰⁰³) who investigated the gating dynamics as a continuous diffusion process with part of the closed state potential being voltage dependent. The model explains how the closed state dwell time distribution of a voltage dependent channel can change its character from a power law to an exponential function, depending on the membrane polarization. Moreover, the model justifies the exponential dependence of the Hodgkin–Huxley model used to describe the experimental relationship between the voltage and the rate of opening of a single gate in the potassium channels.

Further generalisations of diffusion modelling lead to non-Markovian conformational subdiffusion. It has been shown (Goychuk and Hänggi ²⁰⁰⁸) that the subdiffusive generalization reproduces the main features of the closed time distribution and predicts the autocorrelation function of the conductance fluctuations. The non-Markovian nature of the channel currents suggested by Fuliński et al. (¹⁹⁹⁸) has also been addressed in different ways (Grzywna and Siwy ¹⁹⁹⁷; Grzywna et al. ¹⁹⁹⁹; Liebovitch and Sullivan ¹⁹⁸⁷; Liebovitch et al. ¹⁹⁸⁷). Interesting non-Markovian models were given by Grzywna and Siwy (¹⁹⁹⁷) and Grzywna et al. (¹⁹⁹⁹), who provide an alternative approach based on the assumption that the channel behaviour is governed by a nonlinear deterministic dynamics.

The models described capture many important features of the gating phenomena but not all of them. Among the missing features, a very interesting but not completely understood property of many channel species is the long-term memory, as seen by Hurst R/S analysis (Campos de Oliveira et al. ²⁰⁰⁶; Bandeira et al. ²⁰⁰⁸). The rescaled range analysis (R/S analysis, Hurst analysis; Hurst ¹⁹⁵¹; Hurst et al. ¹⁹⁶⁵; Feder ¹⁹⁸⁸; Borys ²⁰¹¹), when applied to the ionic current recordings, measures the correlation of adjacent dwell times of an ion channel (Varanda et al. ²⁰⁰⁰), commonly indicative of strong trend-reinforcing behaviour.

Whether the discrete Markovian models can produce a long-range memory was the subject of research (Varanda et al. ²⁰⁰⁰; Campos de Oliveira et al. ²⁰⁰⁶) on the 3, 4, and 11-state models. All the tests were negative. Whether Markovian models with more states can reproduce the long-range memory, as measured by the Hurst exponent, is, however, still an open question.

According to the experiments, the long-term memory in ion channel activity seems to be independent of external channel activating factors (Barbosa et al. ²⁰⁰⁷; Varanda et al. ²⁰⁰⁰ and the section “Analysis of experimental data” in this paper). It was observed that at each fixed voltage and calcium concentration (at a macro scale), and thus at each corresponding average sensor activation level (at a molecular scale), the long-term correlations in time series are of similar magnitude. It is, therefore, reasonable to link the memory with gate fluctuations under given conditions only. This point of view is also supported by the allosteric picture of BK channel gating dynamics, with gate fluctuations present even in the absence of gating sensor activation.

In this paper, we limit our considerations to gate fluctuations, ignoring the full HCA-MWC machinery; this has the advantage of showing the processes relevant to the long-term memory without obscuring them by the additional details of a gate–sensor coupling. For the gate dynamics we propose two simple random-walk models belonging to the diffusion class of models (Millhauser et al. ¹⁹⁸⁸; Goychuk and Hänggi ²⁰⁰², ²⁰⁰³). The models generate a dichotomous current time series as found in the experimental results (the patch-clamp recordings), and produce the long-term memory of the considered system. The open and closed-dwell time distributions evaluated for the simulated series are concurrent, at a qualitative level, with those observed experimentally. We also propose biophysical interpretations for both models.

It is stated in the “Theoretical background” that the long-term memory is independent of the sensor state. To check whether our models agree with such results, we need some simplified representation of biosystem activation. To propose a reasonable approach, we briefly review the ideas which are the basis of BK channel’s voltage and calcium sensor operation in relation to the gate.

The activation gate is a part of the channel protein responsible for blocking the ion flux through the pore (Chung et al. ²⁰⁰⁷). Because of studies on KcsA and MthK channels (Jiang et al. ²⁰⁰²) we have physical evidence of gate-open and gate-closed conformations of these potassium channels. A key difference between these is the bending of an inner helix segment at a glycine residue close to the selectivity filter, called the “Gly-hinge”.

In many potassium channels, the inner mouth of the pore (the channel’s entrance) is formed by the hydrophobic residues on the C-terminus of the inner helix (S6 in the canonical channel structure with six transmembrane segments). The residues are arranged in a bundle to block the potassium ion flux when the channel is closed. Because of the flexibility of the S6 domain at the Gly-hinge in the middle of S6 (sometimes assisted with the Pro-kinks several residues later), a gate-open conformation can be reached with inner helices well separated from each other (Jiang et al. ²⁰⁰²; Long et al. ²⁰⁰⁵; Ding et al. ²⁰⁰⁵; Cui et al. ²⁰⁰⁸). For the BK channels, it can be shown by mutations that the Gly-hinges are important but not mandatory structures to reach the open conformation, which contrasts with the Kv or MthK channels (Magidovich and Yifrach ²⁰⁰⁴). Experiments with quaternary ammonium ions show further that the BK channel’s selectivity filter entry is accessible in the gate-closed conformation, which suggests the existence of a secondary gate in the close vicinity of the selectivity filter, similarly to several other ligand-gated ion channels (Wilkens and Aldrich ²⁰⁰⁶; Flynn and Zagotta ²⁰⁰¹; Bruenning-Wright et al. ²⁰⁰⁷; Cui et al. ²⁰⁰⁸).

The voltage sensor of the BK channel also differs significantly from the Kv’s, i.e. the gating charge is not confined to S4, but rather spreads over the whole voltage-sensing domain (VSD), and is much smaller than in Kv channels. Furthermore, the ligand sensor is not exactly similar to MthK, because the gating ring of a BK channel has important differences from that of the MthK one (e.g. the gating ring of MthK is formed by eight identical RCK domains, whereas in the BK channel it is formed by four RCK1 and four RCK2 domains).

Ligand sensor operation relies on the interaction between a C-terminal intracellular RCK domain which forms a gating ring which expands on calcium binding, and the inner helix connected to the RCK via an S7-residue linker, whose mutations reveal a profound effect on the gating dynamics (Jiang et al. ²⁰⁰²). The voltage dependency of the BK channels is hypothesized on the basis of homology with the Kv channel gating, which depends on the motion of the charged S4 segment binding through the S4–S5 linker to the S6 segment at a Pro-X-Pro motif, triggering mechanical translations of S6 which result in gating, probably via Gly-hinge bending (Ding et al. ²⁰⁰⁵).

In our approach, the conformational dynamics of the activation gate is described by a discrete random walk process performed by the RC (reaction coordinate) over the conformational space. We could assign the RC a meaning of, e.g., the hinge angle but, because the actual gating blocks of a BK channel are not yet known, we treat RC as abstract coordinate in the conformation space. The effect of the sensors on the gate is allosteric in nature in the BK channels and is introduced by a drift force acting on the reaction coordinate. Such a drift force approach to sensor activation retains the allosteric picture of the system, in which sensor activation is not mandatory for gate opening. Furthermore, the drift force is a good candidate to mimic sensor action, because the sensor seems to operate by forcing the Gly-hinge. In general, this drift force may depend nonlinearly on the reaction coordinate, calcium concentration, and voltage, [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ F_{\text{D}} = F_{\text{D}} ([{\text{Ca}}^{2 + } ],V,{\text{RC}}). $$\end{document}] In our case however, the actual transmission of the gating stimuli is irrelevant (we have assumed that the long-term memory is independent of these stimuli), and we have only checked whether the activation bias modifies the Hurst effect. Having no hints on the RC dependency of the drift force, it was assumed to be constant, which corresponds to the linear ramp potential.

In the diffusive modelling we can assume that the open and closed states of the gate embrace two manifolds of energetically similar substates. The manifolds observed in a macroscopic state can be associated with a set of single protein structure vectors, nearly equal in energy, which retain the macroscopically observed property of the state. From such a viewpoint, the shifting between two different macro-states is realized by a diffusion process between adjacent conformations.

Because we approximate the conformational dynamics in terms of the one-dimensional reaction coordinate x(t) in a conformational potential U(x), as was shown by (Goychuk and Hänggi ²⁰⁰², ²⁰⁰³), the question which arises is, is this approximation valid for a large 3D system such as a protein? Research on protein folding suggests this may often not be true (Dill et al. ¹⁹⁹⁵). Nevertheless, if we consider the activation gate as a hinge, movement of which is effectively restricted to one degree of freedom (or the other possible conformations that can be projected on it), the idea of a one-dimensional reaction coordinate may be valid.

The activation gate’s conformational space is illustrated in Fig. 1. It is divided into two parts, which correspond to the open and closed states of the gate, separated by a threshold point TP. The open and closed states can have different energies, which reflect the tendency of a channel to maintain a conducting or non-conducting state under fixed external conditions. The reaction coordinate performs a random walk in the conformation space with positions before the threshold point (TP) indicating a closed state and those behind the threshold point indicating an open state, as shown in Fig. 1.

We propose that the long-term memory of the system may result from synchronized fluctuations of the activation gate conformational space boundaries (Model 1). The biophysics of such fluctuating boundaries may be explained by considering what happens after squeezing of the channel protein. This may originate from the local fluctuations in the membrane thickness near the channel location. Squeezing will increase the mass density of the channel protein (including the activation gate), and the atoms (being restricted by their Van der Waals’ volumes) have less “space” to move, reducing the number of energetically equivalent substates within a given macro-state. In other words, the boundaries of the conformational space move to reduce the maximum conformational difference between the substates in the closed and open conformations. The opposite may be observed in case of channel protein relaxation by the decrease in the local mass density. This behaviour is schematically depicted in Fig. 2.

Another way of introducing fluctuations of channel density to the gate dynamics is to add a synchronous drift force acting on the reaction coordinate (Model 2). In such picture, lateral compression or relaxation of the channel does not strictly limit the size of the conformational space, but rather produces a drift force acting on the reaction coordinate and favouring the occupancy of substates near the threshold (by relaxation) or near the boundaries (by compression). Such an approach reproduces the long-term memory of the channel system and enables the boundaries to be kept constant. It is interesting to notice that the effect of the compressive force is exactly opposite in Model 2 compared with Model 1. In Model 1, the compression favours rapid state switching whereas in Model 2 this regime is reached in a relaxed state. This may be important for experimental validation of these models. Model 2 drops the assumption of constant energy within the macro-state in the absence of gate–sensor coupling but energy variation among the substates is a smooth, slowly increasing or decreasing function of the RC (Fig. 5).

The conditions for Model 2 may occur if there is a protein segment, connected with gating, which has a large amount of space available for bending. In such case, reduction in the membrane thickness causes the segment to bend, i.e. a force is generated. If this bending can take place either toward or away from the open conformation, which is schematically shown in Fig. 3, then it can be described by the proposed model.

Model 1

The lattice variables of Model 1 are shown in Fig. 4. The conformational space of the gate was projected to one dimension and divided into two parts of equal length separated by a threshold point. The channel is closed if the reaction coordinate occupies the position below the threshold point, otherwise it is open. The transition between the open and the closed conformations requires overcoming of a transition potential barrier as shown on the energy landscape in Fig. 5a.

The range of possible reaction coordinate values in the conformational space is limited by the two moveable boundaries B1 and B2. The motion of B1 and B2 corresponds to thermal fluctuations in the membrane thickness and internal strains within the protein segments which restrict the conformational space of the closed and open states. The boundary fluctuations are synchronized in direction (e.g. when the conformational space of an open state shrinks, the same should happen to the conformational space of the closed state), and the probabilities of these increasing or reducing the accessible conformational space were set equal to 0.5.

Because the mass of the fluctuating membrane is much higher than the mass of the putative activation gate, the boundary fluctuations are realized on a larger time scale, i.e. if we assume that a change of the gate state is recognizable in one time step in simulation, then the change of boundaries may be performed after a given number of time steps according to the assumed relationship of the diffusion time scales:

[Formula ID: Equ4]
4  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ D_{\text{B}} = \frac{{D_{\text{RC}} }}{600}, $$\end{document}]

which results from fitting to the experimental data (D_B is the boundary diffusion coefficient and D_RC the reaction coordinate diffusion coefficient).

To test for the independence of the long-term correlations on sensor activation, a constant drift force from the channel sensors may be added to the system. The biassed random walk process can stand for the behaviour of a channel gate at a fixed value of the activating or inactivating stimulus. For example, in terms of a high value of the membrane potential (e.g. V = 80 mV) open substates are preferred, so the additional drift force will favour movement of the gate toward the open conformation. As a result, the energy structure of substates may change at different activation levels, as illustrated in Fig. 5b, c.

The probabilities of the value of the reaction coordinate (RC) decreasing (q) or increasing (p) were evaluated according to the formulas (Berg ¹⁹⁸³):

[Formula ID: Equ5]
5  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p = \frac{1}{2} - \frac{\Updelta U}{4kT}, $$\end{document}]

[Formula ID: Equ6]
6  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q = \frac{1}{2} + \frac{\Updelta U}{4kT}, $$\end{document}]

where k is the Boltzmann constant, T denotes absolute temperature and ΔU is a potential energy difference within a lattice step centred around the reaction coordinate. The potential U(x) is postulated to take the following form (depicted in Fig. 5):

[Formula ID: Equ7]
7  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {U(x) = (x - B1) \cdot A + U_{{B1}} ;\quad x \in \langle B1;{\text{TP}} - 1.5)} \\ {U(x) = (x - ({\text{TP}} - 1.5)) \cdot B + U_{{{\text{TP}} - 1.5}} ;\quad x \in \langle {\text{TP}} - 1.5;{\text{TP}})} \\ {\begin{array}{*{20}c} {U(x) = (x - {\text{TP}}) \cdot ( - B) + U_{\text{TP}} ;\quad x \in \langle {\text{TP}};{\text{TP}} + 1.5\rangle } \\ {U(x) = (x - ({\text{TP}} + 1.5)) \cdot A + U_{{{\text{TP}} - 1.5}} ;\quad x \in ({\text{TP}} + 1.5;B2\rangle } \\ {A = \frac{{U_{{{\text{TP}} - 1.5}} - U_{{B1}} }}{{\text {TP}} - 1.5 - B1}} \\ \end{array} } \\ \end{array} } \\ {B = \frac{{U_{\text{TP}} - U_{{{\text{TP}} - 1.5}} }}{{\text {TP}} - ({\text {TP}} - 1.5)}} \\ \end{array} } \right., $$\end{document}]

where B1 and B2 are the locations of the left-hand and right-hand boundaries, respectively, TP is the threshold point position, U_B1, U_TP, and U_TP−1.5 denote the potential values at given points, and A and B are the potential slopes between points B1 and TP − 1.5 and between points TP − 1.5 and TP, respectively. According to Eq. (7) the constant A represents the drift force toward the open or closed states of the system and the constant B describes the repulsive force originating from the potential barrier between the open and closed states.

The lattice size was set to 2B_MAX nodes, where B_MAX determines the maximum size of the open (and closed) state space. This space is further restricted by the reflecting boundary positions (B1, B2), which were initially set to −B_MAX/2 and B_MAX/2. The starting RC position was set to −1 (a closed substate nearest the threshold TR = 0). The time step length to move the reaction coordinate t_RC was equal to 1 and the time step length to move the boundaries was set to t_B = 600 (i.e. [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ D_{\text{B}} = \frac{{D_{\text{RC}} }}{600} $$\end{document}] ).

The model variables were estimated by the gradient optimization technique. The values were adjusted in the direction of decreasing error up to a point when the relative errors of the Hurst exponent (E_H), the open state probability (E_pop), and the open and closed state dwell time histograms fitting (E_dw,open) and (E_dw,closed) were less than 5 %, and their sum (E) was less than 10 %. The overall error was given by the formula:

[Formula ID: Equ8]
8  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ E = E_{\text{H}} + E_{\text{pop}} + E_{\text{dw,open}} + E_{\text{dw,closed}} $$\end{document}]

where:

[Formula ID: Equ9]
9  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ E_{\text{H}} = \frac{{\left| {H_{\exp } - H_{\text{sim}} } \right|}}{{H_{\exp } }}, $$\end{document}]

[Formula ID: Equ10]
10  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ E_{\text{pop}} = \frac{{\left| {p_{\text{op,exp}} - p_{\text{op,sim}} } \right|}}{{p_{\text{op,exp}} }}, $$\end{document}]

[Formula ID: Equ11]
11  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ E_{\text{dw,open}} = \frac{{\left| {\chi_{\text{exp,open}}^{2} - \chi_{\text{sim,open}}^{2} } \right|}}{{\chi_{{\exp ,{\text{open}}}}^{2} }}, $$\end{document}]

[Formula ID: Equ12]
12  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ E_{\text{dw,closed}} = \frac{{\left| {\chi_{{\exp ,{\text{closed}}}}^{2} - \chi_{\text{sim,closed}}^{2} } \right|}}{{\chi_{\text{exp,closed}}^{2} }}. $$\end{document}]

where H_exp is the Hurst exponent of the experimental series’ subsequent open and closed dwell times, H_sim is the Hurst exponent of a series of subsequent open and closed dwell times obtained from the simulated data, p_op,exp is the open state probability calculated from the experimental data, p_op,sim is the open state probability calculated from the simulated data, and χ_exp and χ_sim are obtained on the basis of an appropriate normalized dwell time histogram. χ_exp is calculated according to the formula:

[Formula ID: Equ13]
13  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \chi_{\exp }^{2} = \frac{1}{N}\sum\limits_{i = 1}^{N} {\delta_{i}^{2} } $$\end{document}]

where N is the number of the dwell time histogram bins and δ_i is the standard deviation of a given bin (the standard deviation corresponds to the five different patch-clamp traces obtained under the same conditions of voltage and calcium concentration).

χ_sim is described by the equation below:

[Formula ID: Equ14]
14  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \chi_{\text{sim}}^{2} = \frac{1}{N}\sum\limits_{i = 1}^{N} {(p_{i,\exp } - p_{{i,{\text{sim}}}} )^{2} } $$\end{document}]

where N is the number of the dwell time histogram bins and p_i is the probability corresponding to the ith histogram bin.

The sensitivity of the total error (E) to the terms of Model 1 is given in Appendix 3.

The optimum values of the terms in Model 1 are presented in Table 1.

The simulation details of Model 1 (and Model 2) are summarized in Appendix 2.

Model 2

Model 2 differs in the way the channel protein density fluctuations are introduced to the gate dynamics. The boundary positions (B1, B2) are kept constant, and additional, random drift force is added to the system which makes synchronized changes to the energy structure of the open and closed states, as presented in Fig. 6. This force causes the reaction coordinate to fluctuate near the threshold point (TP) or conversely, far away from it. The force fluctuations follow their own random walk, which is slower than the gating dynamics. The appropriate potential is depicted in Fig. 6 and is given by:

[Formula ID: Equ15]
15  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {U(x)= (x - B1) \cdot A + U_{B1} ;\quad x \in\langle B1;{\text{TP}} - 1.5)} \\ {U(x) = (x -({\text{TP}} - 1.5)) \cdot B + U_{{{\text{TP}} - 1.5}} ;\quad x \in \langle {\text{TP}} - 1.5;{\text{TP}})} \\{\begin{array}{*{20}c} {U(x) = (x - {\text{TP}}) \cdot ( - B)+ U_{\text{TP}} ;\quad x \in \langle {\text{TP}};{\text{TP}} +1.5\rangle } \\ {U(x) = (x - ({\text{TP}} + 1.5)) \cdot ( - A)+ U_{{{\text{TP}} - 1.5}} ;\quad x \in ({\text{TP}} +1.5;B2\rangle} \\ {A = \frac{{U_{{{\text{TP}} -1.5}} - U_{B1} }}{{{\rm TP} - 1.5 - B1}}} \\\end{array} } \\ \end{array} } \\ {B = \frac{{U_{\text{TP}} -U_{{{\text{TP}} - 1.5}} }}{{{\text{TP}} - ({\text{TP}} -1.5)}}} \\ \end{array} } \right., $$\end{document}]

where B1 and B2 are the locations of the left-hand and right-hand boundaries, respectively, TP is the threshold point position, U_B1, U_TP, and U_TP−1.5 denote the potential values in given points, and A and B are the potential slopes between points B1 and TP − 1.5, and between points TP − 1.5, and TP, respectively. The A constant represents the drift force which facilitates the RC locations near the threshold point or near the boundaries. The B constant describes the repulsive force originating from the potential barrier between open and closed states.

The effect of the sensors on gate operation is reflected by the threshold point location. The conformational space is divided into two parts not necessarily equal, and the TP position is chosen such that the probability of finding the gate in an open state resembles the channel open probability under the given experimental conditions (Fig. 7). Such a mechanism can be justified by considering the effect of the sensing domains (e.g. the S5–S6 linker) on the gate, where they move the part of protein responsible for gating to increase the number of accessible closed conformations and reduce the number of open ones, or the opposite. The allosteric picture of the channel’s activation is conserved within such approach.

The jump probabilities to the threshold point (p_T), and to the boundaries (q_T), are evaluated by use of Eqs. (5) and (6) (with p_T = p and q_T = q when RC < TP, and p_T = q and q_T = p when RC > TP), and the potential given by Eq. (15). The random walk performed by the RC is illustrated in Fig. 8.

The landscape potential slope (the drift force):

[Formula ID: Equ16]
16  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ F_{\text d} = \pm \frac{{U_{{{\text {TP}} - 1.5}} - U_{B1} }}{{\text {TP}} - 1.5 - {B1}} $$\end{document}]

is controlled by a separate unbiased random walk of U_B1, taking place in a time scale a factor of 1,200 slower than the time scale of the reaction coordinate. The potential slopes cannot exceed minimum and maximum values, as specified in Table 2, with the other simulation values.

The model terms are estimated, analogously as for Model 1, by the gradient optimization technique.

The samples of original single-channel patch-clamp recordings, and the corresponding current–voltage relationship, are shown in Fig. 9. As one can see in Fig. 9a, channel opening varies with the potential value in a way typical of voltage-dependent ion channels, i.e. membrane depolarisation leads to an increase in the frequency of occurrence of the open state. According to Fig 9b, the reversal potential was near 0 mV. Under the experimental conditions, the mean conductance was estimated to be 235.6 pS.

Table 3 shows the values of the Hurst exponent, H, and the channel opening probability, p_o, for the single-channel recordings obtained at different voltages V at fixed calcium concentration [Ca²⁺] equal to 2 mM. It reveals no significant effect of membrane depolarization on long-term correlations. There is just small variability around the mean of H, which is expected because H is a random variable with a standard deviation of 0.1 (Hurst et al. ¹⁹⁶⁵). The mean r² coefficient of linear regression (used by evaluating H) for all experimental dwell time series was equal to 0.992.

Concurrent results were reported by Varanda et al. (²⁰⁰⁰). Analogously, the calcium ion concentration does not affect values of the Hurst exponent significantly (Barbosa et al. ²⁰⁰⁷). This is consistent with our modelling approach in which we neglect the details of operation of the gating sensors.

From both of our models we generated n = 5 time series with N = 6 × 10⁶ samples for five different values of the variable responsible for sensor–gate coupling. For each of the time series, we evaluated the open state probability (p_o) and mean open and closed dwell times; we then performed Hurst analysis of the corresponding dwell-time series of subsequent openings and closings. To determine whether the long-range memory found for the simulated data is its intrinsic feature, the simulated data were shuffled, and the Hurst exponent was calculated for a randomized time series. Table 4 shows the values of the drift-regulating term, the mean values of open state probability, mean open and closed dwell times, the Hurst exponent corresponding to generated time series H, and the Hurst exponent obtained after shuffling (H_sh).

The mean r² coefficient of linear regression used to evaluate H for all dwell time series generated by Model 1 was equal to 0.994 (0.996 for the corresponding H_sh) and for those generated by Model 2—0.987 (0.995 for the corresponding H_sh).

The values obtained indicate that our models enable estimation of approximate ionic current characteristics in quite good agreement with those obtained experimentally under different external conditions. By tuning the regulatory variables F_d (Model 1) or TP (Model 2) one obtains different open probabilities in the generated time series, which mimic the behaviour of the experimental time series under different gate-activating conditions (e.g. membrane depolarization).

The mean open dwell times generated by Model 2 are in quantitative agreement with experimental values, and the opening probability grows with open interval length, as expected. This tendency is conserved at a qualitative level in the data generated by Model 1. The long-term memory in the simulated data is exhibited independently of p_o (as expected), and the values obtained for the Hurst exponent are close to those corresponding to the experimental open and closed dwell time series. After a shuffling procedure, the long-term memory has disappeared (H_sh ≈ 0.5), which means that it is an inherent property of the considered system. The open and closed state dwell-time distributions obtained for the time series generated by Models 1 and 2 were plotted on a log–log scale in Figs. 10 and 11, respectively, and compared with the experimental distribution.

As follows from Fig. 10 Model 1 generates a power law dependence in the closed dwell-time distribution. This result coincides qualitatively with the empirical data and the literature suggestions that this kind of dependence should be expected at least in some cases (Goychuk and Hänggi ²⁰⁰², ²⁰⁰³). In contrast, open dwell time distribution has an exponential tail. In both cases short dwell times are overestimated, and long dwell times are underestimated in comparison with experimental data.

Model 2 enables reconstruction of a correct open dwell time distribution, but the closed dwell time distribution is no longer power law, and is characterized by an exponential tail.

Values of the potential barrier energy separating open and closed states (U_TP − U_TP−1.5) equal to 1.0 kT (Model 1) and 0.2 kT (Model 2), and the diffusion coefficient of the density fluctuations equal to 600 (Model 1) and 1,200 (Model 2) times the diffusion coefficient of reaction coordinate were set to generate a time series which gave the current characteristics similar to the experimental characteristics and had long-term memory as observed in the measurements. The assumed height of the activation barrier is not very large, which concurs with the diffusive nature of the gate dynamics, in which a rather smooth transition from open to closed conformation is expected. Moreover, the energy structure of the conformational space expressed by a smooth function without significantly differing peaks looks quite reasonable from a biological viewpoint, taking into account the existence of a large number of similar protein allosteric states with similar energy.

This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Steps of the Hurst exponent evaluation procedure:

1. The given time series T₁, T₂,…,T_N of length N is divided into d subseries of n elements, such that [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ d \cdot n = N. $$\end{document}] The subseries obtained are denoted I_m (m = 1,…,d) and their elements as t_m,k (k = 1,…,n).
2. The mean value E_m and the standard deviation S_m for each subseries I_m are evaluated.
3. For each I_m an appropriate mean adjusted series X_m = x_m,1,…,x_m,d is found, in which the elements are defined as x_m,k = t_m,k − E_m.
4. The cumulative deviate series Y_m of elements [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ y_{m,k} = \sum\nolimits_{j = 1}^{k} {x_{m,k} } $$\end{document}] are calculated.
5. The ranges R_m described as: R_m = max(Y_1,m,…,Y_n,m) − min(Y_1,m,…,Y_n,m) are found.
6. The rescaled range [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \frac{{R_{m} }}{{S_{m} }} $$\end{document}] for each of the subseries is calculated. Then, the mean value [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left\langle {\frac{{R_{m} }}{{S_{m} }}} \right\rangle $$\end{document}] for the current division of the original time series is evaluated.

All aforementioned steps are then repeated for different lengths (n) of subseries. (What is often seen in the literature is that the lengths of subseries are chosen to be equal to the powers of 2 (n = 2^p, where p = 1, 2, 3,…)).

The resulting relationship between [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left\langle {\frac{{R_{m} }}{{S_{m} }}} \right\rangle $$\end{document}] and n is described by power law scaling of the form:

[Formula ID: Equ27]
27  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left\langle {\frac{{R_{m} }}{{S_{m} }}} \right\rangle_{n} = c \cdot n^{H} $$\end{document}]

Taking the logarithm from Eq. (27) one obtains:

[Formula ID: Equ28]
28  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \log \left( {\left\langle {\frac{{R_{m} }}{{S_{m} }}} \right\rangle_{n} } \right) = H \cdot \log (n) + \log (c) $$\end{document}]

The Hurst exponent may be evaluated by use of the least-squares linear regression method, and its value is equal to the slope of the plot of Eq. (28) as illustrated in Fig. 12.

When one time series is analysed the goodness of fit (R²) may be used as the error of the estimated Hurst exponent value. However, in the case when one has at least a few repetitions of the experiment, e.g. a series of independent measurements carried out under the same conditions (in our studies these are the patch-clamp time series for five patches on the same voltage values) it is better to use a standard deviation from the mean value obtained for the repetitions as a measure of the Hurst coefficient error. A fitting error for one series of finite length could not resemble the properties of a set of analogous series, because of the possibility of over or underestimation, owing to the poorly averaged terms in the fitting at large n values which are likely to perturb randomly the linearity of the data.

The Model 1 simulation procedure:

1. The lattice with the maximum number of nodes equal to 2B_MAX is set.
2. The initial positions of the reaction coordinate x and the fluctuation boundaries B1 and B2 are set. The threshold point (TP) is equal to 0. The B1 and B2 are symmetrically located around the TP.
3. The potential function (U(x)) described by Eq. (7) is associated with the lattice.
4. The position of the reaction coordinate is randomly changed for one step length in one time step, with the probability of movement to the right p and to the left q described by Eqs. (5) and (6), respectively.
5. The position of the reaction coordinate relative to the TP is checked. If the RC is at the right-hand side of the TP, the open state is recognized. Otherwise the closed state is stated.
6. If the RC reaches the B1 or B2 positions during its random walk, it is reflected to its previous position.
7. Steps 4, 5, 6 are repeated for a number of time steps, which is determined by Eq. (4).
8. The boundaries B1 and B2 are randomly and synchronously moved for one step length toward or away from the TP with equal probability. (If B1 reaches B_MIN or TP − 1 positions, it is reflected to its previous position. Analogously, if B2 reaches B_MAX or TP + 1 positions, it is reflected to its previous position.)
9. Steps 4–8 are repeated for a desired time series length.

The Model 2 simulation procedure:

1. The lattice with the maximum number of nodes equal to 2B2 and the position of the threshold point TP are set.
2. The initial position of the reaction coordinate x is set.
3. The potential function (U(x)) described by Eq. (15) is associated with the lattice. The initial “synchronous” drift is set to be equal to 0.
4. The position of the reaction coordinate is randomly changed for one step length in one time step with the probability of movement toward the TP (p_T) and away from it (q_T) described by Eqs. (5) and (6).
5. The position of the reaction coordinate relative to the TP is checked. If the RC is at the right-hand side of the TP, the open state is recognized. Otherwise the closed state is stated.
6. If the RC reaches the B1 or B2 positions during its random walk, it is reflected to its previous position.
7. Steps 4, 5, 6 are repeated for a number of time steps, which is determined by Eq. (4).
8. The potential slopes (and, as a consequence, the value of the synchronous drift force described by Eq. (16)) is randomly changed by a constant ±|Δk| with equal probability.
9. Steps 4–8 are repeated for a desired time series length.

The model variables are determined by use of the gradient optimization technique. Optimization was performed until the sum of relative errors is less than 10 %. The values chosen are provided in Tables 1 and 2.

To support the quality of the values found we provide the dependencies of the total error in the model variables when one of them is being changed and the rest are kept constant and equal to the optimum values (Figs. 13, 14, 15, 16, 17, 18, 19). The error plots, being realizations of a stochastic process (they are evaluated after a simulation run), reveal some variability between realizations of the simulation, which can be characterized by a standard deviation of ±2 %.

On each of the figures, the chosen value is denoted by a dot.

From Fig. 15, it could be roughly stated that the total errors in the range of potential barriers from 0.2 to 1.0 kT maintain quite similar values. In the simulations, we have chosen its value equal to 1.0, because for lower values the error component corresponding to the Hurst exponent (E_H) was relatively higher, as shown by the dashed line in Fig. 15. Minimizing the E_H is thought to be important, taking into account the main objectives of the model.

Considering the errors connected with |Δk| estimation, although the total error was lower for |Δk| = 0.010 than for |Δk| = 0.005 (Fig. 16), we have regarded the second value as optimum, because of the corresponding Hurst exponent errors, which were equal to 0.066 and 0.053 for |Δk| = 0.010 and 0.005, respectively.

Analysing the total error dependency on the time scales’ ratio (Fig. 17), one can state, that the ratios below 600 are characterized by relatively high total errors. In the range up to 1,800 one can observe errors of similar magnitude, which average to about 13 %. Then, for ratios larger than 1,800 the error grows again. Because the exact plot in range between 600 and 1,800 has variability depending on the simulation realization (it is a random variable), we have chosen the middle value of this interval, i.e. [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \frac{600 + 1,800}{2} = $$\end{document}] 1,200.

^arcu reaction coordinate unit

Errors are given as standard deviations

Hurst exponents: H for original, and H_sh for shuffled simulated time series. The aforesaid values were calculated as means from five time series (Model 1, Model 2) for each stimulus-regulatory value