Document Detail

From MEDLINE®/PubMed®, a database of the U.S. National Library of Medicine

 Full Text Journal Information Journal ID (nlm-ta): Eur Biophys J Journal ID (iso-abbrev): Eur. Biophys. J ISSN: 0175-7571 ISSN: 1432-1017 Publisher: Springer-Verlag, Berlin/Heidelberg Article Information Download PDF © The Author(s) 2012 Received Day: 6 Month: 9 Year: 2011 Revision Received Day: 6 Month: 3 Year: 2012 Accepted Day: 13 Month: 3 Year: 2012 Electronic publication date: Day: 7 Month: 4 Year: 2012 pmc-release publication date: Day: 7 Month: 4 Year: 2012 Print publication date: Month: 6 Year: 2012 Volume: 41 Issue: 6 First Page: 505 Last Page: 526 ID: 3359465 PubMed Id: 22484857 Publisher Id: 806 DOI: 10.1007/s00249-012-0806-8
 On the simple random-walk models of ion-channel gate dynamics reflecting long-term memory Agata Wawrzkiewicz1 Address: +48-32-2371067 +48-32-2371722 agata.wawrzkiewicz@gmail.com Krzysztof Pawelek1 Przemyslaw Borys1 Beata Dworakowska2 Zbigniew J. Grzywna1 1Department of Physical Chemistry and Technology of Polymers, Section of Physics and Applied Mathematics, Silesian University of Technology, Ks. M. Strzody 9, 44-100 Gliwice, Poland 2Division of Biophysics, Department of Physics, Warsaw University of Life Sciences—SGGW, Nowoursynowska 166, 02-787 Warsaw, Poland

Introduction
Theoretical background

Potassium channels are integral proteins that enable rapid, selective transport of K+ ions across the cell membrane down their electrochemical gradient (Chung et al. 2007). 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 2001).

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 1983; Latorre et al. 1989; Marty 1983) 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. 2008). Irrespective of their unusually large conductance, BK channels remain highly selective for K+ ions over other cations (Cui et al. 2008).

In general, the putative structure of BK channels shares many similarities with that of voltage-dependent Kv channels (Ma et al. 2006; Liu et al. 2008a, b) and ligand gated potassium channels, for example MthK (Latorre and Brauchi 2006; Jiang et al. 2001, 2002; Shi et al. 2002; Xia et al. 2002; Tang et al. 2003; Yusifov et al. 2008; Kim et al. 2006, 2008). 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 2006; Cui et al. 2008; Latorre and Brauchi 2006). 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. 2008).

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

[Formula ID: Equ1]
[\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}]
in which the states of the channel switch randomly between the open and closed conformations with gating-factor-dependent kinetic rate constants (kO, kC). This kind of description may belong to a Markov class of model, if the state transition probabilities depend only on the current state of the system (Fuliński et al. 1998).

Invention of the patch-clamp technique (Sakmann and Neher 1995) has enabled experimental observation of the ionic currents at a single-channel level, which gave the opportunity to obtain the probability densities of closed fC(t) and open fO(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 2002, 2003):

[Formula ID: Equ2]
[\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}]
with weight coefficients obeying:
[Formula ID: Equ3]
[\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. 1988) or a power law function (Sansom et al. 1989; Blatz and Magleby 1986; Ring 1986; Mercik and Weron 2001; Läuger 1988; Condat and Jäckle 1989). 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. 1989).

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. (1999), Horrigan and Aldrich (1999), and Rothberg and Magleby (1999). The HCA model proposed by Horrigan et al. (1999) and Horrigan and Aldrich (1999) 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 2001). 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 1999).

The calcium activation of BK channels is characterized in the voltage-dependent Monod Wyman Changeux (VD-MWC) model (Cox et al. 1997). Analogously to the previous case, because of incorporation of an allosteric mechanism, none of the Ca2+-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 Ca2+. In this model, the open state has a higher Ca2+ affinity than a corresponding closed state with the same number of calcium ions bound. Thus, an open Ca2+-bound state is energetically preferred over a closed Ca2+-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 2000; Horrigan and Aldrich 2002). 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. 1987) and diffusive models (Millhauser et al. 1988; Läuger 1988; Condat and Jäckle 1989; Goychuk and Hänggi 2002, 2003) 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. (1988). This model produces the power law dependence in the dwell time distributions. It was generalised by Hänggi and Goychuk (2002, 2003) 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 2008) 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. (1998) has also been addressed in different ways (Grzywna and Siwy 1997; Grzywna et al. 1999; Liebovitch and Sullivan 1987; Liebovitch et al. 1987). Interesting non-Markovian models were given by Grzywna and Siwy (1997) and Grzywna et al. (1999), 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. 2006; Bandeira et al. 2008). The rescaled range analysis (R/S analysis, Hurst analysis; Hurst 1951; Hurst et al. 1965; Feder 1988; Borys 2011), when applied to the ionic current recordings, measures the correlation of adjacent dwell times of an ion channel (Varanda et al. 2000), 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. 2000; Campos de Oliveira et al. 2006) 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. 2007; Varanda et al. 2000 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. 1988; Goychuk and Hänggi 2002, 2003). 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.

Description of the models and biological inspirations

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. 2007). Because of studies on KcsA and MthK channels (Jiang et al. 2002) 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. 2002; Long et al. 2005; Ding et al. 2005; Cui et al. 2008). 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 2004). 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 2006; Flynn and Zagotta 2001; Bruenning-Wright et al. 2007; Cui et al. 2008).

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. 2002). 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. 2005).

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 2002, 2003), 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. 1995). 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]
[\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 (DB is the boundary diffusion coefficient and DRC 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 1983):

[Formula ID: Equ5]
[\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]
[\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]
[\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, UB1, UTP, and UTP−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 2BMAX nodes, where BMAX 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 −BMAX/2 and BMAX/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 tRC was equal to 1 and the time step length to move the boundaries was set to tB = 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 (EH), the open state probability (Epop), and the open and closed state dwell time histograms fitting (Edw,open) and (Edw,closed) were less than 5 %, and their sum (E) was less than 10 %. The overall error was given by the formula:

[Formula ID: Equ8]
[\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]
[\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 Hexp is the Hurst exponent of the experimental series’ subsequent open and closed dwell times, Hsim is the Hurst exponent of a series of subsequent open and closed dwell times obtained from the simulated data, pop,exp is the open state probability calculated from the experimental data, pop,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 pi 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, UB1, UTP, and UTP−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 (pT), and to the boundaries (qT), are evaluated by use of Eqs. (5) and (6) (with pT = p and qT = q when RC < TP, and pT = q and qT = 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 UB1, 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.

Materials and methods
Experimental
Cell line and solutions

Human bronchial epithelial cells were cultured on Petri dishes in minimum essential Eagle’s medium (Sigma) with 10 % fetal calf serum, 100 units/ml penicillin, and 100 μg/ml (PAA) at 37 °C in 5 % CO2. For patch-clamp recordings, the culture medium was exchanged for an extracellular solution containing 4 mM CaCl2, 2 mM EGTA, 10 mM HEPES, and 135 mM potassium gluconate, pH 7.3.

Electrophysiology

Experimental results were recorded from outside-out patches of human bronchial epithelial cells. The measurements were performed at room-temperature (20–21 °C). In all experiments symmetrical solutions on either side of the cell membrane were used. Patch pipettes were pulled from borosilicate capillary tubes of 1.2 mm diameter (Clark Electromedical Instruments) on a Narishige micro puller. The tips of the pipettes were smoothed in a Narishige micro forge. The ion currents were recorded using an Axopatch 200B amplifier (Axon Instruments). The experimental data were low-pass filtered at 10 kHz and transferred to a computer at a sampling frequency of 20 kHz using Clampex 7 software (Axon Instruments). Single currents were initially analysed by use of pClamp 7 software (Axon Instruments). The single-channel recordings were analysed at a fixed pipette potential. The experiments were carried out at: −80, −60, −40, −20, 20, 40, 60, and 80 mV and five independent measurements were recorded for each voltage. Each experiment was performed on a different patch (so we obtained 8 × 5 = 40 different patch-clamp recordings in total). The channel current was measured at time intervals of Δt = 5 × 10−5 s. The ionic current measurement error was ΔI = 5 × 10−4 pA. Each of patch-clamp recordings lasted 300 s, thus the experimental time series comprised N = 6 × 106 current values at the applied time resolution of the measurement.

Data analysis
Event detection

Investigating experimental time series, we considered two modes of ion-channel conduction, here called the open (conducting) and closed (non-conducting) states. The threshold current value used to identify transitions between following states was evaluated by use of a procedure described elsewhere (Mercik et al. 1999). The ion current probability density function (PDF) approximated by the nonparametric kernel density estimate with the Epanechnikov kernel was plotted on a log–log scale. The resulting graph reflects the bimodal nature of the analysed data set. The estimate obtained may be regarded as a mixture of two unimodal densities that satisfy power laws. By means of linear regression we obtained the intervals in both component distributions where the power laws and the corresponding formulas expressing scaling relationships were valid. A point of intersection of the power law plots indicates a threshold current value, in relation to which we identify open and closed states.

It is well worth noticing that if one applies the double Gaussian method of threshold current determination, the difference between the values obtained and those estimated by the kernel density method is less than 5 % (for each patch-clamp time series obtained in our experiments). The nonparametric technique of the PDF approximation avoids the assumption that the analysed data belong to any particular distribution, so we prefer this approach.

The R/S Hurst analysis

The rescaled range analysis (R/S Hurst analysis) is an important statistical method used for testing whether a system under consideration has long-term memory (Hurst 1951; Varanda et al. 2000). The Hurst exponent (H) describes the time scaling of a range (R) in a random motion, normalized to zero mean increments, [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R \propto t^{H} .$$\end{document}] It can take values from 0 to 1, and provides a means of classification of temporal series in terms of predictability. A Hurst exponent of 0.5 is indicative of purely random behaviour of the system, and is an easily derived characteristic of the scaling of the standard deviation in Brownian motion (Risken 1996) (it is also characteristic of the scaling of R in the Brownian motion, as shown with somewhat larger effort by Borys 2011). In such a case there is no correlation among any element of the temporal series. If 0 < H < 0.5 the system is said to be anti-persistent, i.e. the range of corresponding anomalous Brownian motion grows slower than randomly. This means that a positive increment will tend to be followed by a negative one (and vice versa). The time series with a Hurst exponent value from 0.5 to 1 is persistent, or trending, i.e. the range of the corresponding Brownian motion grows faster than randomly. This means that a positive increment will tend to be followed by a positive one and a negative increment will tend to be followed by a negative one. The larger the H value, the stronger the trend, and the easier to predict the future behaviour of the system. In this paper we carry out a Hurst analysis for a sequence of adjacent open and closed times for both the BK channel recordings available, and the modelled data.

The procedure for evaluating the Hurst exponent can be found in many references (Hurst 1951; Hurst et al. 1965; Feder 1988; Varanda et al. 2000; Borys 2011). Details are also included in Appendix 1.

Results and discussion
Analysis of experimental data

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, po, for the single-channel recordings obtained at different voltages V at fixed calcium concentration [Ca2+] 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. 1965). The mean r2 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. (2000). Analogously, the calcium ion concentration does not affect values of the Hurst exponent significantly (Barbosa et al. 2007). This is consistent with our modelling approach in which we neglect the details of operation of the gating sensors.

Analysis of the generated time series

From both of our models we generated n = 5 time series with N = 6 × 106 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 (po) 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 (Hsh).

The mean r2 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 Hsh) and for those generated by Model 2—0.987 (0.995 for the corresponding Hsh).

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 Fd (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 po (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 (Hsh ≈ 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 2002, 2003). 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 (UTP − UTP−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.

Conclusion and outlook

In this work we have proposed two random-walk models of BK channel gate dynamics with long-term memory. Their construction was motivated by Hurst analysis applied to the experimental data measured in the voltage range from −80 to 80 mV and for a Ca2+ concentration of approximately 2 mM, supported by the results obtained by other authors under different experimental conditions (Bandeira et al. 2008; Barbosa et al. 2007; Varanda et al. 2000), which suggest that intrinsic long-term memory exists in the series of subsequent single channel state dwell-times.

According to the experiments, the mechanism responsible for the long-term memory seems not to be directly related to the voltage and calcium sensors’ activation (but may, still, be related to their activity). Taking into account the biological context of the BK channel system, it is reasonable to assume that thermal fluctuations of the activation gate machinery by itself could be responsible for the spontaneous transitions between the open and closed states of the channel under fixed conditions of the activating voltage and [Ca2+]. Nevertheless, the simple random fluctuations of an activation gate are not enough to produce the long-term memory, and a degree of synchronicity is required.

We have shown how to reproduce the long-term correlations in the gate dynamics by membrane thickness fluctuations, which affect the available size of conformational space and/or induce stress on the gating segments. The slow thermal fluctuation of the membrane density surrounding the activation gate was incorporated by introduction of the synchronous diffusion space boundary fluctuations in Model 1, and the drift force fluctuations in Model 2. The synchronicity is thought to originate from a channel reaction to squeezing or stretching. For example, if the cell membrane underlies squeezing, a part of the membrane (here: the channel’s surrounding) increases its density. This also affects the channel gate by narrowing the distribution of its conformational states.

The part of the channel protein responsible for channel gating has notably smaller mass than the whole channel protein immersed in the cell membrane. This difference manifests itself in the fluctuation time scales of the channel gate and its surrounding. To illustrate a connection between the physical and theoretical properties of the channel system consider the following reasoning:

Within a random walk framework, the diffusion coefficient D can be described by the formula:

[Formula ID: Equ17]
17  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D \propto v \cdot \delta = v^{2} \cdot \tau$$\end{document}]
where v is the velocity of a considered object, δ denotes the mean free path, and τ is the mean time between subsequent collisions with the environment.

The frequency of collisions is proportional to the object’s surface (S), so τ is inversely proportional to S. Taking this into account, and the proportionality between the mass and volume (m ~ V ~ S3/2), the following relationship is valid:

[Formula ID: Equ18]
18  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \propto m^{ - 2/3} .$$\end{document}]

According to the equipartition theorem:

[Formula ID: Equ19]
19  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$kT \propto \frac{1}{2}mv^{2} ,$$\end{document}]
where k is the Boltzmann constant and T is temperature.

By combining Eqs. (17)–(19), one obtains the relationship between the mass and diffusion coefficient:

[Formula ID: Equ20]
20  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \propto \left( \frac{2kT}{D} \right)^{3/5} .$$\end{document}]

Thus, the relationship between the diffusion coefficients (and time scales) used in our models corresponds to the ratio of the masses of the considered channel elements:

[Formula ID: Equ21]
21  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{m_{\text{S}} }}{{m_{\text{CG}} }} = \left( {\frac{{D_{\text{CG}} }}{{D_{\text{S}} }}} \right)^{3/5}$$\end{document}]

(subscripts: S, “surrounding”; CG, channel gate).

The models predict a ratio between the gate’s diffusion coefficient and diffusion coefficient of the surroundings expressed by [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{D_{\text{S}} }}{{D_{\text{CG}} }} = \frac{{D_{\text{RC}} }}{{D_{\text{B}} }} = 600$$\end{document}] and [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{D_{\text{S}} }}{{D_{\text{CG}} }} = \frac{{D_{\text{RC}} }}{{D_{\text{DF}} }} =$$\end{document}] 1,200 for Model 1 and Model 2, respectively. The problem, however, is that such a gate diffusion coefficient may not necessarily be the “true” gate diffusion coefficient, and it may be an apparent one observed because of a low sampling rate in a self-similar gating process (Liebovitch and Toth 1990).

Considering the S6-helix of a BK channel of r = 2.5 Å radius, its diffusion coefficient in the cell membrane should take a value similar to D ≈ 4 × 10−8 cm2/s = 4 × 10−12 m2/s (Almeida and Vaz 1995). To estimate our apparent gating diffusion coefficient we need to approximate our reaction coordinate unit (rcu) by some physical quantity. It seems reasonable to assume that the S6-helix may function as the channel gate. The channel vestibule diameter is approximately equal to d ≈ 20 Å (Cui et al. 2008). From the symmetry of the BK channel system (which is a tetrameric protein) one S6-helix has a d/2 = 10 Å range of movements. In our simulation, this distance is projected on a lattice of length equal to 2BMAX length (Model 1). Thus:

[Formula ID: Equ22]
22  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1\,{\text{rcu}} \approx \frac{d}{{2 \times 2B_{{{\text{MAX}}}} }} = \delta = 0.25\,{\AA}$$\end{document}]
i.e. 1 rcu corresponds to 0.25 Å.

The simulation time step length is equal to τ = 5 × 10−5 s (because of the experimental sampling frequency of 20 kHz), which renders the gate diffusion coefficient DCG estimate:

[Formula ID: Equ23]
23  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{\text{CG}} = \frac{{\delta^{2} }}{6\tau } = 1.04 \times 10^{ - 17} \left[ {\frac{{m^{2} }}{s}} \right].$$\end{document}]

This value differs by six orders or magnitude from the expected value of DCG. This difference arises from the τ value restricted by the sampling frequency. Although more rapid channel gate fluctuations are not recorded, one cannot preclude their existence. According to Liebovitch and Toth (1990), the kinetic properties of patch-clamp recordings may have a self-similar nature. So, it is justified to consider the DCG as an apparent diffusion coefficient resulting from the observation time scale τ, which should be higher (lower τ, fixed δ in Eq. (23)) on the thermal fluctuation scale.

The diffusion coefficient of the channel’s surrounding DS was assumed in the simulation to satisfy Eq. (4) (Model 1), thus DS may be estimated as DCG/600 = 1.73 × 10−20[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left[ {\frac{{m^{2} }}{s}} \right].$$\end{document}] From Eq. (21), we can estimate the relationship between the actual mass of the fluctuating membrane (mm) and the channel gate (m), which is:

[Formula ID: Equ24]
24  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{m_{\text{m}} }}{m} = \left( {\frac{D}{{D_{\text{S}} }}} \right)^{3/5} = 1.04 \times 10^{5}$$\end{document}]

To check whether this large ratio is acceptable physically, one may estimate the volume of the membrane patch used in the patch-clamp experiment, which could act as a “fluctuating surrounding” of the assumed channel gate. The volume of an S6-helix V (channel gate) may be approximated as a cylindrical volume spanning the membrane, which is:

[Formula ID: Equ25]
25  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V = \pi \cdot r^{2} \cdot d_{\text{m}} \approx \pi \cdot (2.5)^{2} \cdot 70\,{\text{\AA}}^{ 3} = { 1,374}\,{\text{\AA}}^{ 3}$$\end{document}]
where dm = 70 Å, is the average membrane thickness.

The volume of the whole-cell membrane patch Vwm is described as:

[Formula ID: Equ26]
26  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{\text{wm}} = S_{\text{wm}} \cdot d_{\text{m}} = 10^{8} \cdot 70\, {\text{\AA}}^{ 3} = { 7} \times 10^{ 9} \,{\text{\AA}}^{ 3}$$\end{document}]
where Swm denotes the surface area of the whole cell membrane patch, Swm = 108 Å2 (Sakmann and Neher 1995).

The ratio of volumes yields Vwm/V = 5.09 × 106. By use of this result, and the ratio of masses mm/m = 1.04 × 105, one can infer that 2 % of the whole patch volume should be sufficient to act as the fluctuating “channel surrounding”. Such a result seems reasonable, and it may confirm, in some sense, the validity of our models.

The models can be further verified experimentally by measuring the membrane thickness fluctuations (possibly by changes in the membrane capacitance) simultaneously with the current recordings. In such an experiment, our model predicts a correlation between membrane thickness and local dwell time. A positive correlation between membrane thickness and dwell time would indicate Model 1, whereas a negative correlation would indicate Model 2.

Table 4 shows that the characteristics of the generated time series remain in a good agreement with those obtained experimentally (Table 3), in particular:

1. Our models enable generation of time series with different open state probabilities, and mean open dwell times which are in reasonable agreement with experimental values, recorded under different conditions of channel activating stimuli.
2. The Hurst analysis applied to the sequence of the opening and closing times in the simulated data furnishes mean values of 0.77 and 0.70 for Model 1, and Model 2, respectively. This result of the rescaled range analysis suggest that investigated data are long-term correlated, and the exponent remains in the allowed error range for the data (±0.1) (Hurst et al. 1965).

The Hurst exponents for the shuffled series show no long memory effects. In consequence, it is reasonable to assume that the trend-reinforcing behaviour is an inherent feature of the modelled systems.

Considering the experimental data obtained under the conditions described in the “Electrophysiology” section, it can be stated that the closed dwell-time distributions have power law dependence, and the open dwell-time distributions are characterized by an exponential tail. By use of Model 1 one can obtain the distributions which qualitatively reflect the empirical dependencies. By use of Model 2, qualitative and quantitative agreement in the open dwell time distributions is reached, but the experimental and simulated closed dwell time distribution types do not match exactly.

Details of the voltage and of the calcium sensors’ behaviour, and their effect on channel gate activity, were reduced in our models. The coupling between sensors and the gate was introduced in the simplest way possible: a constant drift acting on the reaction coordinate (Model 1) or a moveable threshold position (Model 2). As a consequence, our models provide only an approximate description of the dependence of gating on [Ca2+] and voltage, but the ideas presented can be introduced to more sophisticated models, for example MWC or HCA; this will, however, substantially increase their complexity.

We are grateful to Professor K. Dołowy from the Warsaw University of Life Sciences, Warsaw, Poland, for enabling us to carry out all patch-clamp experiments to obtain the original time series of ionic current recorded from BK potassium channels. The authors would like to thank The Ministry of Science and Higher Education for providing financial support under project N N508 409137.

Open Access

References
 Almeida PFF,Vaz WLC. Lipowsky R,Sackmann ELateral diffusion in membranesStructure and dynamics of membranes from cells to vesiclesYear: 1995AmsterdamElsevier305359 Bandeira HT,Barbosa CTF,Campos De Oliveira RA,Aguiar JF,Nogueira RA. Chaotic model and memory in single calcium-activated potassium channel kineticsChaosYear: 20081803313619045474 Barbosa CTF,Campos de Oliveira RA,Nogueira RA. Calcium does not change memory in single calcium-activated potassium channel kineticsJ Intell Fuzzy SystYear: 20075477484 Berg HC. Random walks in biologyYear: 1983PrincetonPrinceton University Press Blatz AL,Magleby KL. Quantitative description of three modes of activity of fast chloride channels from rat skeletal muscleJ PhysiolYear: 19863781411742432249 Borys P. On the card tricks, Nile floods and the Hurst exponentFotonYear: 2011113422 Bruenning-Wright A,Lee WS,Adelman JP,Maylie J. Evidence for a deep pore activation gate in small conductance Ca2+-activated K+ channelsJ Gen PhysiolYear: 200713060161017998394 Campos de Oliveira RA,Barbosa CTF,Consoni LHA,Rodrigues ARA,Varanda WA,Nogueira RA. Long-term correlation in single calcium-activated potassium channel kineticsPhys AYear: 20063641322 Chung S-H,Andersen OS,Krishnamurthy V. Biological membrane ion channels. Dynamics, structure and applicationsYear: 2007New YorkSpringer Condat CA,Jäckle J. Closed-time distribution of ionic channels. Analytical solution to a one-dimensional defect-diffusion modelBiophys JYear: 1989559159252470430 Cox DH. Chung S-H,Andersen OS,Krishnamorthy VBKCa-channel structure and functionBiological membrane ion channels. Dynamics, structure and applicationYear: 2006New YorkSpringer171210 Cox DH,Cui J,Aldrich RW. Allosteric gating of a large conductance Ca-activated K+ channelsJ Gen PhysiolYear: 19971102572819276753 Cui J,Yang H,Lee US. Molecular mechanism of BK channel activationCell Mol Life SciYear: 20086685287519099186 Dill KA,Bromberg S,Yue K,Fiebig KM,Yee DP,Thomas PD,Chan HS. Principles of protein folding—a perspective from exact modelsProtein SciYear: 199545616027613459 Ding S,Ingleby L,Ahern CA,Horn R. Investigating the putative glycine hinge in Shaker potassium channelJ Gen PhysiolYear: 200512621322616103276 Feder J. FractalsYear: 1988New YorkPlenum Press Flynn GE,Zagotta WN. Conformational changes in S6 coupled to the opening of cyclic nucleotide-gated channelsNeuronYear: 20013068969811430803 Fuliński A,Grzywna ZJ,Mellor I,Siwy Z,Usherwood PNR. Non-Markovian character of ionic current fluctuations in membrane channelsPhys Rev EYear: 199858919924 Goychuk I,Hänggi P. Ion channel gating: a first-passage time analysis of the Kramers typePNASYear: 20029963552355611891285 Goychuk I,Hänggi P. The role of conformational diffusion in ion-channel gatingPhys AYear: 2003325918 Goychuk I,Hänggi P. Fractional diffusion modeling of ion-channel gatingPhys Rev EYear: 200870051915 Grzywna ZJ,Siwy Z. Chaos in ionic transport through membranesInt J Bifurcation ChaosYear: 1997511151123 Grzywna ZJ,Siwy Z,Fuliński A,Mellor I,Usherwood PNR. Chaos in the potassium current through channels of locust muscle membraneCell Mol Biol LettYear: 199943754 Hille B. Ion Channels of Excitable membranesYear: 2001Sunderland, MASinauer Horrigan FT,Aldrich RW. Allosteric voltage gating of potassium channels II. Mslo channel gating charge movement in the absence of Ca2+J Gen PhysiolYear: 199911430533610436004 Horrigan FT,Aldrich RW. Coupling between voltage sensor activation, Ca2+ binding and channel opening in large conductance (BK) potassium channelsJ Gen PhysiolYear: 200212026730512198087 Horrigan FT,Cui J,Aldrich RW. Allosteric voltage gating of potassium channels I. Mslo ionic currents in the absence of Ca2+J Gen PhysiolYear: 199911427730410436003 Hurst HE. Long-term storage capacity of reservoirsTrans Am Soc Civ EngYear: 1951116770808 Hurst HE,Black RP,Simaika YM. Long-Term StorageYear: 1965Constable, LondonAn Experimental Study Jiang Y,Pico A,Cadene M,Chait BT,MacKinnon R. Structure if RCK domain from the E. coli K+ channel and demonstration of its presence in the human BK channelNeuronYear: 20012959360111301020 Jiang Y,Lee A,Chen J,Cadene M,Chait BT,MacKinnon R. Crystal structure and mechanism of a calcium-gated potassium channelNatureYear: 200241751552212037559 Kim HJ,Lim HH,Rho SH,Eom SH,Park CS. Hydrophobic interface between two regulators of K+ conductance domains critical for calcium-dependent activation of large-conductance Ca2+-activated K+ channelsJ Biol ChemYear: 20062813857317040919 Kim HJ, Lim HH, Rho SH, Bao L, Lee JH, Cox DH, Kim do H, Park CS (2008) Modulation of the conductance-voltage relationship of the BKCa channel by mutations at the putative flexible interface between two RCK domains. Biophys J 94:446-456 Latorre R,Brauchi S. Large conductance Ca2+-activated K+ (BK) channel: Activation by Ca2+ and voltageBiol ResYear: 20063938540117106573 Latorre R,Miller C. Conduction and selectivity in potassium channelsJ Membr BiolYear: 19837111306300405 Latorre R,Oberhauser A,Labarca P,Alvarez O. Varieties of calcium-activated potassium channelsAnnu Rev PhysiolYear: 1989513853992653189 Läuger P. Internal motions in proteins and gating kinetics of ionic channelsBiophys JYear: 1988538778842456104 Liebovitch LS,Sullivan JM. Fractal analysis of a voltage-dependent potassium channel from cultured mouse hippocampal neuronsBiophys JYear: 1987529799882447974 Liebovitch LS,Toth TI. Using Fractals to Understand the Opening and Closing of Ion ChannelsAnn Biomed EngYear: 1990181771941693478 Liebovitch LS,Fishbarg J,Koniarek JP. Ion channel kinetics: a model based on fractal scaling rather than multistate Markov processesMath BiosciYear: 1987843768 Liu G,Zakharov SI,Yang L,Deng SX,Landry DW,Karlin A,Marx SO. Position and role of the BK channel alpha subunit S0 helix inferred from disulfide cross-linkingJ Gen PhysiolYear: 200813153754818474637 Liu G,Zakharov SI,Yang L,Wu RS,Deng SX,Landry DW,Karlin A,Marx SO. Locations of the beta1 transmembrane helices in the BK potassium channelProc Natl Acad Sci USAYear: 2008105107271073218669652 Long SB,Campbell EB,MacKinnon R. Crystal structure of a mammalian voltage-dependent Shaker family K+ channelScienceYear: 200530989790316002581 Ma Z,Lou XJ,Horrigan FT. Role of charged residues in the S1–S4 voltage sensor of BK channelsJ Gen PhysiolYear: 200612730932816505150 Magidovich E,Yifrach O. Conserved gating hinge in ligand-and voltage-dependent K+ channelsBiochemistryYear: 200443132421324715491131 Marty A. Ca2+-dependent K+ channels with large unitary conductanceTrends NeuroscYear: 19836262265 Mercik S,Weron K. Stochastic origins of the long-range correlations of ionic current fluctuations in membrane channelsPhys Rev EYear: 200163051910 Mercik S,Weron K,Siwy Z. Statistical analysis of ionic current fluctuations in membrane channelsPhys Rev EYear: 19996073437348 Millhauser GL,Salpeter EE,Oswald RE. Diffusion models of ion-channel gating and the origin of power-law distributions from single-channel recordingProc Natl Acad Sci USAYear: 198885150315072449693 Ring A. Brief closures of gramicidin A channels in lipid bilayer membranesBiochim Biophys ActaYear: 19868566466532421773 Risken H. The Fokker-Planck equation: methods of solutions and applicationsYear: 19962BerlinSpringer Rothberg BS,Magleby KL. Gating kinetics of single large-conductance Ca2+-activated K+ channels in high Ca2+ suggest a two-tiered allosteric gating mechanismJ Gen PhysiolYear: 19991149312410398695 Rothberg BS,Magleby KL. Voltage and Ca2+ activation of single large-conductance Ca2+-activated K+ channels described by a two-tiered allosteric gating mechanismJ Gen PhysiolYear: 2000116759910871641 Sakmann B,Neher N. Single-Channel RecordingYear: 19952New YorkEds Plenum Sansom MSP,Ball FG,Kerry CJ,McGee R,Ramsey RL,Usherwood PNR. Markov, fractal, diffusion, and related models of ion-channel gating. A comparison with experimental data from two ion channelsBiophys J Year: 198956122912432482085 Shi J,Krishnamoorthy G,Yang Y,Hu L,Chaturvedi N,Harilal D,Qin J,Cui J. Mechanism of magnesium activation of calcium-activated potassium channelsNatureYear: 200241887688012192410 Tang XD,Xu R,Reynolds MF,Garcia ML,Heinemann SH,Hoshi T. Haem can bind to and inhibit mammalian calcium-dependent Slo1 BK channelsNatureYear: 200342553153514523450 Varanda WA,Liebovitch LS,Figueiroa JN,Nogueira RA. Hurst analysis applied to the study of single calcium-activated potassium channel kineticsJ Theor BiolYear: 200020634335310988020 Wilkens CM,Aldrich RW. State-independent block of BK channels by an intracellular quaternary ammoniumJ Gen PhysiolYear: 200612834736416940557 Xia XM,Zeng X,Lingle CJ. Multiple regulatory sides in large-conductance calcium-activated potassium channelsNatureYear: 200241888088412192411 Yusifov T,Savalli N,Gandhi CS,Ottolia M,Olcese R. The RCK2 domain of the human BKCa channel is a calcium sensorProc Natl Acad Sci USAYear: 200810537638118162557
Appendix 1

Steps of the Hurst exponent evaluation procedure:

1. The given time series T1, T2,…,TN 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 Im (m = 1,…,d) and their elements as tm,k (k = 1,…,n).
2. The mean value Em and the standard deviation Sm for each subseries Im are evaluated.
3. For each Im an appropriate mean adjusted series Xm = xm,1,…,xm,d is found, in which the elements are defined as xm,k = tm,k − Em.
4. The cumulative deviate series Ym 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 Rm described as: Rm = max(Y1,m,…,Yn,m min(Y1,m,…,Yn,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 = 2p, 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}]
where H denotes the Hurst exponent and c is a constant.

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

Appendix 2

The Model 1 simulation procedure:

1. The lattice with the maximum number of nodes equal to 2BMAX 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 BMIN or TP − 1 positions, it is reflected to its previous position. Analogously, if B2 reaches BMAX 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 (pT) and away from it (qT) 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.

Appendix 3

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 (EH) was relatively higher, as shown by the dashed line in Fig. 15. Minimizing the EH 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.

Figures
 [Figure ID: Fig1] Fig. 1  Schematic illustration of a channel gate’s possible state transitions. One can observe openings (C → O), closings (O → C), and changes among open (O → O) or closed (C → C) states. TP threshold point separating open and closed states, x reaction coordinate [Figure ID: Fig2] Fig. 2  Schematic illustration of the changes of local membrane mass density and its effect on channel protein. The range of channel gate’s accessible conformations can be reduced or increased as a result of membrane squeezing and its relaxation, respectively, as denoted by the frame on the lower “diamond” in the figure [Figure ID: Fig3] Fig. 3  Squeezing of the membrane will generate a force acting on the channel segment, which is related to gating, and causing its bending. Bending motion may lead either to the gate’s closed or open conformations (a). Stretching of the membrane results in occurrence of a force which enables relaxation of the bent segment. The force strength depends on current thickness of membrane (and consequently, degree of bending) (b) [Figure ID: Fig4] Fig. 4  Schematic representation of Model 1. a The one-dimensional lattice has 20 nodes, among which first 10 represent closed states of the channel whereas the other nodes correspond to open states. Boundaries (B1, B) can move, simultaneously reducing (b) or increasing (c) the space accessible to the reaction coordinate (RC). Without any external force, no direction of motion is preferred (d), otherwise a drift in the direction down the potential energy gradient occurs (e), which changes the probability of finding the system in a given state [Figure ID: Fig5] Fig. 5  Schematic representation of the potential function associated with Model 1. a No force arising from sensor activity is acting, in particular, on the direction of the gate. The corresponding conformational potential energy U(x) of the reaction coordinate is “flat” and symmetric around the threshold point separating open and closed substates. If a nonzero drift force is present, it is represented by the conformational potential energy U(x) sketched on b or c. In b external potential field is expressed by a decreasing ramp function. As a result, open states will be preferred. Conversely, in a linearly increasing potential field closed states will be preferred (c) [Figure ID: Fig6] Fig. 6  Schematic representation of the potential function associated with Model 2. a The random drift force facilitates movement toward boundaries. As a result, the conformational potential energy U(x) of the reaction coordinate increases from the boundary B1 to the threshold point TP and decreases from TP to the second boundary B2. In the opposite case, because of the effect of random force, positions around the threshold are preferred. The appropriate potential energy U(x) is sketched in b [Figure ID: Fig7] Fig. 7  Schematic representation of the diffusive space used in Model 2. Boundaries (B1, B2) and the threshold position (TP) are fixed during the simulation. a If no gating stimulus is present in the system, the manifolds of closed and open states are of the same length. b The different numbers of possible open and closed substates account for preference of the channel for a particular macroscopic state. Among 20 nodes of the lattice used as the diffusive space first six indicate the closed substates of the gate, the other 14 nodes indicate open substates [Figure ID: Fig8] Fig. 8  Schematic representation of Model 2. The reaction coordinate (RC) shifts either toward the threshold position (TP) with probability pT represented as an arrow pointing toward TP (a, b, pT > qT) or away from it (arrow pointing toward B1 or B2) (c, d, qT > pT). The potential slope and, consequently, the values of pT and qT evolve randomly during the simulation (at a slower rate than evolution of the reaction coordinate) [Figure ID: Fig9] Fig. 9  a Samples of the original ion current signal recorded from a single BK channel over the range of membrane potentials in fixed calcium concentration ([Ca2+] = 2 mM). Dashed line indicates closed state of the channel. b Current–voltage relationship for the BK channel recordings presented in a [Figure ID: Fig10] Fig. 10  Closed (a) and open (b) residence time distribution for simulated time by use of Model 1, and recorded experimentally at [Ca2+] = 2 mM and V = 80 mV. Simulation data are given in Table 1, and the drift coefficient was chosen appropriately to obtain the open state probability possibly close to experimental data. The time unit for modelled data was rescaled to the experimental value. The fitting procedure was performed for experimental data, first, in linear coordinates; the results were then transformed to the logarithmic scale [Figure ID: Fig11] Fig. 11  Closed (a) and open (b) residence time distribution for simulated time by use of Model 2, and recorded experimentally at [Ca2+] = 2 mM and V = 80 mV. Simulation data are given in Table 2, and the threshold point was chosen appropriately to obtain the open state probability possibly close to experimental data. The time unit for modelled data was rescaled to the experimental value. The fitting procedure was performed for experimental data, first, in linear coordinates; the results were then transformed to the logarithmic scale [Figure ID: Fig12] Fig. 12  The Hurst exponent value (H) may be obtained by plotting Eq. (28) and evaluating its slope. In the figure, H is equal to 0.77 [Figure ID: Fig13] Fig. 13  The total error as a function of BMAX (Model 1). The values of all other variables are given in Table 1 [Figure ID: Fig14] Fig. 14  The total error as a function of the boundaries and reaction coordinate time scales ratio (Model 1). The values of all other variables are given in Table 1 [Figure ID: Fig15] Fig. 15  The total error (E) (values on the leftY axis) as a function the height of the potential barrier (Model 1). The appropriate error components corresponding to the Hurst exponent (EH) are represented by the dashed line (values on the rightY axis). The values of all other variables are given in Table 1 [Figure ID: Fig16] Fig. 16  The total error as a function of the drift force’s increment |Δk| (Model 2). The values of all other variables are given in Table 2 [Figure ID: Fig17] Fig. 17  The total error as a function of the boundaries and reaction coordinate time scales’ ratio (Model 2). The values of all other variables are given in Table 2 [Figure ID: Fig18] Fig. 18  The total error as a function the height of the potential barrier (Model 2). The values of all other variables are given in Table 2 [Figure ID: Fig19] Fig. 19  The total error as a function the boundary B1 location (Model 2). The values of all other variables are given in Table 2

Tables
[TableWrap ID: Tab1] Table 1

Values of the terms used in the Model 1 simulation

Term Value
BMAX 14 rcu
BMIN −14 rcu
TP 0 rcu
Reaction coordinate diffusion coefficient DRC 104 rcu2/s
Barriers diffusion coefficient DB DRC/600
UTP − UTP−1.5 1.0 kT
Initial position of RC −1 rcu
Initial position of B1, B2 Floor (BMIN/2), Floor (BMAX/2) rcua

arcu reaction coordinate unit

[TableWrap ID: Tab2] Table 2

Values of the terms used in simulation of Model 2

Term Value
B1 −18 rcu
B2 18 rcu
Reaction coordinate diffusion coefficient DRC 104 rcu2/s
Drift force diffusion coefficient DDF DRC/1,200
UTP − UTP−1.5 0.2 kT
Initial position of RC −1 rcu
Initial drift force 0.00 kT/rcu
Drift force’s increment |Δk| 0.005 kT/rcu
Max value of Fd 0.20 kT/rcu
Min value of Fd −0.20 kT/rcu

[TableWrap ID: Tab3] Table 3

Mean experimental values of the Hurst exponent (H), the open state probability (po), and mean open and closed dwell time (Τopen and Τclosed) obtained at [Ca2+] = 2 mM under different voltage conditions

V (mV) H ± ΔH po ± Δpo Τopen ± ΔΤopen (ms) Τclosed ± ΔΤclosed (ms)
−80 0.70 ± 0.05 0.17 ± 0.02 0.70 ± 0.10 0.46 ± 0.10
−60 0.76 ± 0.02 0.33 ± 0.07 1.00 ± 0.30 0.60 ± 0.20
−40 0.63 ± 0.03 0.63 ± 0.03 1.30 ± 0.15 0.80 ± 0.10
−20 0.73 ± 0.04 0.78 ± 0.07 1.20 ± 0.20 0.25 ± 0.10
20 0.82 ± 0.05 0.86 ± 0.07 1.90 ± 0.35 0.25 ± 0.10
40 0.63 ± 0.04 0.92 ± 0.01 2.80 ± 0.15 0.20 ± 0.05
60 0.75 ± 0.08 0.94 ± 0.04 2.85 ± 0.85 0.65 ± 0.35
80 0.77 ± 0.04 0.95 ± 0.06 3.10 ± 0.40 0.25 ± 0.10

Errors are given as standard deviations

[TableWrap ID: Tab4] Table 4

The values of characteristics arising from the gating stimulus, and the mean values and standard deviations (SD) of open state probability (po)

Fd (kT/rcu) po ± SD H ± SD Τopen ± ΔΤopen (ms) Τclosed ± ΔΤclosed (ms) Hsh ± SD
Model 1
0.40 0.15 ± 0.01 0.74 ± 0.02 0.32 ± 0.01 1.75 ± 0.15 0.51 ± 0.01
0.20 0.25 ± 0.01 0.79 ± 0.01 0.46 ± 0.02 1.38 ± 0.11 0.52 ± 0.01
0.00 0.50 ± 0.01 0.82 ± 0.01 0.74 ± 0.07 0.74 ± 0.07 0.52 ± 0.01
−0.20 0.74 ± 0.01 0.79 ± 0.01 1.25 ± 0.09 0.43 ± 0.20 0.51 ± 0.01
−0.40 0.85 ± 0.01 0.73 ± 0.01 1.62 ± 0.04 0.31 ± 0.01 0.53 ± 0.01
TP po ± SD H ± SD Τopen ± ΔΤopen (ms) Τclosed ± ΔΤclosed (ms) Hsh ± SD
Model 2
14 0.16 ± 0.02 0.69 ± 0.01 0.63 ± 0.01 3.25 ± 0.44 0.51 ± 0.01
7 0.32 ± 0.02 0.71 ± 0.01 1.35 ± 0.08 2.93 ± 0.46 0.53 ± 0.01
0 0.50 ± 0.01 0.72 ± 0.02 2.16 ± 0.28 2.13 ± 0.25 0.52 ± 0.01
−7 0.68 ± 0.02 0.71 ± 0.01 2.87 ± 0.45 1.32 ± 0.08 0.52 ± 0.01
−14 0.85 ± 0.02 0.68 ± 0.01 3.79 ± 0.66 0.64 ± 0.02 0.51 ± 0.01

Hurst exponents: H for original, and Hsh 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

 Article Categories:Original Paper Keywords: Keywords Random walk process, Conformational diffusion, Hurst analysis, BK channels, Activation gate.

Previous Document:  Ionizable side chains at catalytic active sites of enzymes.
Next Document:  Importance of observation interval in two-dimensional video analysis of individual diatom cells.