Wavelet-based protocols for ion channel electrophysiology

Abstract

Background

Fluctuation-induced phenomena caused by both random and deterministic stimuli have been previously studied in a variety of contexts. They are based on the interplay between the spectro-temporal patterns of the signal and the kinetics of the system it is applied to. The aim of this study was to develop a method for designing fluctuating inputs into nonlinear system which would elicit the most desired system output and to implement the method to studies of ion channels.

Results

We describe an algorithm based on constructing the input as a superposition of wavelets and optimizing it according to a selected cost functional. The algorithm is applied to ion channel electrophysiology where the input is the fluctuating voltage delivered through a patch-clamp experimental apparatus and the output is the whole-cell ionic current. The algorithm is optimized to aid selection of Markov models of the gating kinetics of the voltage-gated Shaker K⁺ channel and tested by comparison of numerically obtained ionic currents predicted by different models with experimental data obtained from the Shaker K⁺ channels. Other applications and optimization criteria are also suggested.

Conclusion

The method described in this paper can be useful in development and testing of models of ion channel gating kinetics, developing voltage inputs that optimize certain nonequilibrium phenomena in ion channels, such as the kinetic focusing, and potentially has applications to other fields.

Background

Ion channels in cellular membranes are proteins that form gated pores to allow passive transport of ions down their electrochemical potential gradient [1]. They open and close in response to an appropriate gating stimulus, such as transmembrane voltage, ligand binding, or mechanical stress. In this paper we concentrate on voltage-gated ion channels. Their gating can be mathematically described as a Markov chain with voltage-dependent transition rates between a small set of discrete states. Such a model corresponds to a microscopic picture of the channel as a macromolecular assembly undergoing conformational changes. The gating is a generalized motion in a certain “energy landscape” and the discrete Markov states correspond to local energy minima. The transition rates reflect thermally-activated jumps over energy barriers separating the minima. In practice, both the topology (the number and connectivity of the discrete states) as well as the transition parameters of such models are developed to fit the experimental data that comes mostly in a form of electrophysiological recordings [2–9].

The technique of patch-clamping in ion channel electrophysiology allows measurements of currents flowing through the channels in cellular membranes [10]. When a glass micropipette containing a recording electrode touches the cell membrane it forms a tight seal (a gigaseal) around a patch of the membrane. Thus clamped patch can be ruptured, creating a direct electrical connection with the cytoplasm and allowing measurements of ionic currents through the entire cell surface (the whole-cell mode), or maintained. In the latter case the recording is of currents flowing through this small patch only, which contains very few, or even just a single channel (the single-channel mode). Typical recordings are performed under the voltage-clamp conditions, where the membrane potential is controlled by a patch-clamp amplifier and the resultant current measured. A specific form of the applied voltage depends on the details of the channel gating kinetics that the protocol is intended to probe. For instance, the activation protocols are based on a voltage step up from a hyperpolarized value, at which channels are closed, to various depolarized values. These protocols test the activation of channels at different voltages. In contrast, the tail protocols are designed to observe channel closing (deactivation) due to a re-polarizing voltage step. The channels are subjected to a depolarizing prepulse of a certain duration, the result of which is the opening of a number of channels. It is then followed by a re-polarizing voltage step and the channel deactivation (closing) can be observed.

These two protocols illustrate the current paradigm on which a vast majority of electrophysiological experiments are based: the voltage protocols are piecewise constant and consist of only few voltage steps at discrete times [2–9]. From a physical point of view this corresponds to observing an ensemble of ion channels under equilibrium or near-equilibrium conditions. For a constant voltage, the distribution of channels among discrete Markov states reaches an equilibrium form. Following a voltage step, this equilibrium is disturbed, but the ensemble simply relaxes to a new equilibrium distribution corresponding to the new voltage value. An entirely different approach, where the voltage fluctuates on a time scale comparable to the relaxation times of the channel kinetics has been discussed [10–17]. This approach forces the channels into nonequilibrium distributions which may lead to new phenomena not observable under equilibrium conditions. These types of fluctuation-induced effects in nonlinear systems are well known in various other areas of physics and in recent years there has been a growing interest in investigating such phenomena in biological systems, including ion channels. One of the most studied examples is the stochastic resonance [17–22], where a noise, either intrinsic or added to the system, improves the system response to weak time-dependent signals. Other examples include ratchets [23, 24], resonant activation [25–28], or the nonequilibrium kinetic focusing [29–31].

Several of these effects can be very promising in investigating and controlling the kinetics of ion channels. For instance, one of the main goals of ion channel research is to develop a model of channel gating kinetics. Most commonly used are the discrete Markov chains, however even the basic features of these models are still disputed, e.g. cooperativity or Markovian character of gating. It has been suggested [10, 11, 16, 17] that this ambiguity results partially from a very incomplete set of experimental data and by expanding the data we can develop better and more unique models. In essence, if we have several models that adequately reproduce the existing set of experimental data, these models can be tested and disproved by comparing to new experimental data until one or more models fail to match them. One of the ways the available data sets can be expanded is by using fluctuating voltages in voltage-clamp electrophysiological experiments, contrary to a currently preferred method of applying piecewise constant voltages. The problem with this approach is in deciding what type of fluctuating voltage input would be most useful. In the previous studies the dichotomous noise has been used [10, 11] but the exact properties of the noise (amplitude, frequency, temporal asymmetry) have been chosen arbitrarily. What is needed is a systematic method of selecting a fluctuating voltage input that would generate maximally different responses from channel models, thus allowing us to efficiently select the model most compatible with all existing experimental data.

A related issue is that of controlling the gating of ion channels. It has an enormous practical importance in biology and medicine. So far the dominant approach is to use pharmacological agents (drugs and toxins) or global electric fields (defibrillator shocks) to modify channel behavior. A much more subtle approach, aimed at forcing ion channels into a specific conformational state has been also proposed [29]. Known as the nonequilibrium kinetic focusing, the method is based on applying fluctuating voltage to ion channels in order to enhance transitions into a selected state and suppress transitions out of that state. As a result, in the ensemble of ion channels, i.e. in a cell, most channels will occupy a selected state. This is also a nonequilibrium effect and is not achievable in standard electrophysiology protocols. The kinetic focusing has been studied numerically and analytically [29, 31] using the dichotomous noise stimulation but it is entirely unclear what type of stimulation would be optimal to achieve the maximum focusing in a desired state. These are two examples of problems in ion channel research where one needs a systematic method of selecting a fluctuating voltage for use in patch-clamping experiments which would optimize the experimental output, either by producing maximally divergent responses from kinetic models of gating, thus facilitating model selection, or by focusing ion channels into a specific conformational state.

In this paper we describe in detail such a method of designing an input to a system that maximizes a desired output. It is based on constructing a fluctuating signal as a superposition of wavelets in a dyadic wavelet basis and optimizing the wavelet coefficients for a specific system response. We developed the method for use with the experimental technique of patch-clamping but it could be applied in a variety of contexts where fluctuating system inputs are used. We concentrate in this paper on applications to ion channel electrophysiology, in particular since our method represents a radical departure from the current paradigm. We show how one can synthesize a signal with desired properties and implement it in a voltage-clamp experiment with ion channels.

We illustrate the method by applying to one of the two problems mentioned earlier in the introduction, i.e. by developing voltage protocols that maximize the difference between various Markov models for the same ion channels. The numerical results are compared to experimental data obtained from Shaker K⁺ channels. The idea of this method was mentioned in our previous work on the nonequilibrium response spectroscopy technique [17]. A version of this method was also used to analyze ionic currents in human heart sodium channels [16]. A related method, where a pulse composed of wavelets is designed to maximize a response of nonlinear systems has been also proposed recently in engineering [32].

Methods

Markov models of ion channels

Ion channel gating results from rearrangements of the tertiary structure of the channel proteins, i.e. transitions between certain meta-stable conformational states of the molecule, in response to changes in trans-membrane potential. These molecular states can be either conducting (open – O) or non-conducting (closed – C or inactivated – I). Examples of Markov models for the Shaker K⁺ channels are shown in Figure 1. Transitions between the states are thermally activated and, in case of voltage-gated channels, are voltage-dependent according to [33]:

(1)

where α is a generic transition rate, V – is the membrane voltage, k – the Boltzmann constant, and T – the absolute temperature. q is the so-called gating charge characterizing the molecule’s sensitivity to external electric fields. Such a description is only a coarse approximation and more precisely the channel gating should be viewed as a motion of a “gating particle” in a certain energy landscape, subject to thermal fluctuations and governed by the Langevin or Fokker-Planck equation [29, 31, 34–36]. However, a picture of channel gating as a discrete Markov chain with the discrete states corresponding to the minima in this energy landscape has been very successful. It disregards the internal structure of each energy well, but there is very little evidence that this structure has any significant effect on macroscopically observable quantities in most applications. In practice, the main task in functional ion channel studies is to determine the topology of Markov models (i.e. the number, type and connectivity of discrete states) as well as the parameters for the state transition rates (the gating charges q and the rates at 0 mV membrane potential). The model topology and model parameters are determined by fitting model responses to various sets of experimental data. As we argued in the past [17], the data used for model development is usually incomplete, since it consists mostly of responses to piecewise-constant voltage stimulation. The result is that several distinct models can be developed for the same kinetic data.

Figure 1

Examples of Markov models for Shaker IR K ⁺ channels.

Wavelet decomposition of voltage inputs

The purpose of this paper is to describe and implement a method of constructing fluctuating voltage inputs for ion channel electrophysiology to obtain a certain type of response from ion channels. The method exploits the interplay between spectro-temporal patterns of the signal and the gating kinetics of the channel molecules. It requires simultaneous control over temporal and spectral properties of the voltage pulse. A tool that offers this kind of time-frequency localization is the wavelet analysis, developed in the last two decades [37, 38]. Wavelets are “localized waves”, i.e. oscillating functions ψ concentrated around a certain region in the argument space with zero mean:

(2)

and finite “energy”:

(3)

There are many different families of wavelets, some with compact support. In essence, the wavelet analysis is similar to Fourier analysis, however it has certain advantages. In Fourier transform the testing functions (trigonometric) have infinite support hence they are not very well suited to analysis of nonstationary signals. Fourier transform provides spectral information about the signal, however temporal information, although not lost, is hidden in the phases of the sines and cosines. By proper translation and dilation of the testing wavelet, the wavelet transform can provide both spectral and temporal information, i.e. it shows not only what frequencies are present, but also when they are present in a nonstationary signal.

Wavelet analysis begins with a mother wavelet, i.e. a function satisfying conditions (2)-(3). The simplest is the Haar wavelet:

(4)

which is also the first member of the Daubechies family of wavelets. By varying the dilation parameter a and the translation parameter b we obtain from the mother wavelet a set of normalized wavelets of the form:

(5)

The continuous wavelet transform CWT is the convolution of the signal x(t) with the members of this wavelet set:

(6)

The CWT coefficients T(a,b) can be plotted versus the dilation and translation parameters (a and b). More common are plots of the wavelet energy density E(a,b) = |T(a,b)|², known as the scalograms. The coefficients T(a,b) (or the energy E(a,b)) carry the information about how much the signal is compatible with the wavelet at scale a, at location b. One may also notice that since the wavelet is a function with oscillation within a certain, relatively narrow, frequency range, then moving to a different scale amounts to changing the dominant frequency of the wavelet. This can be expressed in terms of the mean wavelet frequency f ₀ [38]. In that sense a scalogram shows spectro-temporal patterns in the signal.

In principle, the CWT is an expansion of a signal x(t) into an uncountable set of wavelets {ψ _a,b }. Some wavelets, such as the Daubechies family, have been designed so that a countable subset of dilated and translated wavelets, of the form:

(7)

constitutes an orthonormal basis in L²(ℜ). In other words, every signal x(t) with finite “energy”, defined as

(8)

can be approximated arbitrarily well by a finite superposition of basis elements ψ _m,n . Such a basis is called the dyadic wavelet basis. Using wavelets (7) we can define the discrete wavelet transform (DWT) as:

(9)

Just like for the CWT, the DWT coefficients can be plotted versus the scale and translation indices (m and n), yielding the discrete transform plot. Since on the axes one has discrete indices, the plot has a block form. Knowing the DWT coefficients, a signal can be synthesized using:

$$x\left(t\right)={\displaystyle \sum _{m,n}{T}_{m,n}{\psi }_{m,n}}\left(t\right)$$

(10)

One of the applications is the so-called multiresolution analysis (MRA). It is an iterative procedure where a given signal is separated into a coarse approximation and the detail. The latter is expressed as a superposition of wavelets at a certain level m (with a given scale parameter m), while the former is expressed in terms of the companion function – the scaling function. Next the level-m approximation is separated into a coarser, level-m + 1 approximation and the detail. The latter is again a superposition of wavelets at scale m + 1. Ultimately, the signal can be written as a sum of the scaling function and a combination of wavelets at different scales. The orthogonality of wavelets assures that the signal details at different scales are uncorrelated. The DWT given by equation (9) is discrete in scale and position, but still is a continuous function of time. In practice signals are discretely sampled. In this case in each iteration of the MRA both the approximation and the details are subsampled by a factor of 2. The process can continue until the “coarsest” level (i.e. the mean value) is reached.

The technique we use to construct signals (voltage inputs) with desired properties is based on the inverse DWT (10) and is essentially the opposite of the MRA we just described. For a chosen wavelet type (we tested the Haar wavelets since they most closely resemble currently used voltage step protocols and the Daubechies 8 wavelets for their compact support and the degree of smoothness) we synthesized the signal using (10) and including a finite number of levels. This number of MRA levels was determined by physiological considerations and the bandwidth of the recording apparatus. The latter was typically of the order of 5–10 kHz, and from previous studies [10–17] we expected the effect of the fluctuating voltages to be concentrated in the 1–2 kHz region, hence we considered voltage waveforms composed of 8–10 MRA levels. Rather than imposing a finite energy requirement in the form (8) we put a constraint on the amplitude of the voltage fluctuations not to exceed a physiologically reasonable value. We chose 200 mV peak-to-peak since larger oscillations resulted in a quick loss of a gigaseal.

Channel model response and input optimization

In the previous paragraph we described how a signal with desired properties can be synthesized in a dyadic wavelet basis. We should also specify what these “desired” properties are. The voltage input is selected to elicit a certain response from ion channels. The method we are describing in this paper is flexible enough for a variety of purposes. It can be used to synthesize inputs that maximize differences in computed ionic currents between various models, hence aiding model development and testing, or inputs that maximize the parametric sensitivity of the model, as described in [13]. Another possible application to ion channels would be to test for the phenomenon called the nonequilibrium kinetic focusing, described in a numerical studies of ion channels [29, 31], and recently investigated experimentally [30]. The technique can be also used in applications to other fields [32].

In order to determine the model output, let us consider an n-state Markov model for an ion channel (see e.g. Figure 1). The current state of the channel can be described as an n—component normalized vector P (t) of probabilities {P _i (t)} of finding a channel in various Markov states. For a large ensemble of channels, e.g. in whole-cell mode patch-clamping experiments, these equal the state occupancies. The time evolution of the system is described by the master equation:

(11)

where the W [V(t)] is the transition matrix consisting of the transition rates (1) between various states. For a variable voltage input a formal solution to (11) can be written in a form:

(12)

where P (0) is the equilibrium probability distribution at the holding potential. We compute (12) numerically by iteration:

(13)

with a sufficiently small time step. P (t) represents the time evolution of the probability distribution. If the goal is to get a model output that is comparable with the experimentally measurable quantities, such as ionic currents, then the model current can be computed from the probability distribution using the Ohm’s law:

(14)

Here g ₀ g(V) O·P (t) is the total conductance, where O·P (t) is the projection of the probability vector onto the open state(s), g(V) - the nonlinear, voltage dependent part of the conductance, and g ₀ is the overall scaling factor dependent on the number of channels in the ensemble and the conductance of an open channel.

Genetic algorithm

The algorithm for voltage input synthesis is an optimization procedure. Depending on the goals we define the “cost functional” for the model output and then modify the inputs to optimize this functional. We use a version of a genetic algorithm applied in other contexts previously [10–17]. We begin with a first generation voltage pulse obtained from (10) with a set of randomly chosen wavelet coefficients T _m,n . From this set of coefficients we obtain a number, typically 10–20, of second generation sets by random perturbations according to:

(15)

where T is one of the wavelet coefficients, n numbers the generations of the algorithm, E ₀ and R are convergence parameters characterizing the initial range of the random perturbations (E ₀ ) and the rate at which the search narrows down with each generation (R), and X is a uniformly distributed random number from the interval (0,1). In practice X was generated using the pseudorandom number generator in Matlab, and for a run of 15,000 generations we used E ₀  = 0.5 and R = 0.0004. For each of the daughter sets we synthesize the input according to (10), compute the model output using (13), and evaluate the cost functional. The daughter set that optimizes the cost functional is chosen as the parent for the following generation. The procedure is iterated but as eq. (15) shows at each generation the range of random perturbations used to generate the daughter sets decreases. There are different possible variations of this scheme, where for instance the input is optimized level by level or all levels at once. Other methods are also possible for constructing the new daughter sets, for instance in [32] a differential evolution is described where the new generation inputs are obtained not by multiplicative random perturbation of the chosen “parent”, but according to a simple additive formula involving the “parent” and several randomly chosen intermediate sets. As mentioned in the previous section, there are practical constraints on this optimization procedure. For instance, the maximum amplitude of voltage fluctuations must be restricted to a certain range.

Channel electrophysiology

As an illustration of this numerical scheme we constructed voltage inputs designed to maximize the difference in model outputs between different Markov models for the Shaker K⁺ channels shown in Figure 1. The simulated data were compared to experimental currents recorded from the channels stably expressed in tSA201 cells. The details of experimental design were described previously [12–14]. Our electrophysiology apparatus consists of the Axopatch 200B amplifier with a CV 203 BU headstage (Axon Instruments Inc., Union City, CA), Nikon TS100 inverted microscope Nikon Instruments Inc., Lewisville, TX), Sutter MP-285 micromanipulator (Sutter Instrument Co., Novato, CA). The experiments were performed at 12°C, maintained by Physitemp temperature controller (Physitemp Instruments Inc., Clifton, NJ). Fluctuating voltage inputs were prepared as binary files with Matlab (Mathworks Inc., Natick, MA), read into the Pulse program (HEKA Electronik Gmbh, Lambrecht, Germany) and converted to analog form by ITC18 AD/DA converter (Instrutech Corp., Great Neck, NJ). The ionic currents were digitized at a rate between 50 and 200 kHz and stored on a hard disk. All analysis was performed with custom programs written in Matlab.

We studied a mutant Shaker channel (Shaker K Sk1), stably expressed in tsA201 cells (gift from D. Hanck). The cells were cultured in 35 mm Corning Petri dishes (Corning) in DMEM medium (ATCC, Manassas, VA) supplemented with 10% FBS (ATCC, Manassas, VA), 1% penicillin/streptomycin (Gibco BRL, Gaithersburg, MD) and 200 μg/ml Zeocin (Invitrogen, Carlsbad, CA) at 37°C in 5% CO₂ in a CO₂ incubator (Fisher Scientific, Pittsburgh, PA). Patch clamping pipettes were pulled from borosilicate capillary glass (Warner Instruments LLC, Hamden, CT) on a Sutter pipette puller (Sutter Instrument Co., Novato, CA).

Results

Models for Shaker K⁺channels

As an illustration of the proposed method we consider one of the main goals in ion channel studies – the development of kinetic models of channel gating. Such models are developed by fitting to various types of experimental data (mostly various ionic current recordings). If the set of experimental data is small, the number of models compatible with it is typically very large. Model refinement is achieved by introducing new types of data (e.g. gating currents) and by increasing the number of different voltage protocols. The paradigm of electrophysiology is to use voltage protocols consisting of static voltages changing only stepwise at few discrete times. Typical examples are the activation and tail protocols. Very few studies used randomly fluctuating voltages such as the dichotomous noise [10–12].

Several models have been developed for the Shaker K⁺ channels. We consider here the Bezanilla-Peroso-Stefani (BPS) [3], Schoppa-Sigworth (SS) [6], and the Zagotta-Hoshi-Aldrich (ZHA) [9] models (Figure 1). The BPS model is a simple linear chain of 8 states and was developed based on the activation and deactivation (tail) protocols, as well as “on” and “off” gating currents. The model contains 5 independent pairs of transition rates (forward α and backward β). According to (1) there are two model parameters (the rates at 0 mV and the gating charges) per transition. Zagotta et al. proposed several models in a series of 3 papers [7–9] describing a number of experiments based on stepped voltage protocols. The ZHA A and D models, which we analyze here, also take into account basic knowledge of the channel structure. Namely, it is believed that the Shaker channel consists of four symmetrical subunits, each of which undergoes an activation process. It is mostly postulated that the subunit activation is independent from each other and the channel becomes open when all four subunits are activated. The ZHA A model is the simplest possible such topology. The four subunits can exist in two states only (resting R, and activated A) and the activation of all of them is equivalent to the opening of the channel. The transition rates between the R and A states are the same for all subunits. This topology can also be represented in its expanded form where states C₁ to C₄ describe the channel conformations with different number of subunits in the A state (0 for C₁, 1 for C₂, 2 for C₃, 3 for C₄). The conformation with all four subunits activated is the open (O) state of the channel. The transition rates in the expanded topology follow directly from the rates in the abbreviated version hence the model involves only 1 pair of transition rates (4 model parameters). The ZHA D model introduces two complications. First, the subunit kinetics now involves two resting states, so it is itself a linear chain of three states. Second, there is another closed state (C_f) accessible only from the open state. Both corrections were introduced to improve the model fit to a variety of available experimental data. The expanded version of this topology is no longer linear and involves 16 distinct channel states, however the number of model parameters is the same as for the abbreviated version (three pairs of transition rates, i.e. 12 parameters). In that respect the BCS model is significantly more complex and computationally demanding, with 20 parameters. Finally, the SS model is the most complex [6]. Like the ZHA D model, it is rooted in the idea of 4 subunits undergoing a number (three in this case) independent transitions followed by two additional transitions that describe the concerted motion of the subunits. It also has three inactivated states reachable from the open state only. In the fully expanded form (not shown) this model has 38 states.

For all the models shown in Figure 1 we can write the transition matrices, according to (11). For instance the matrices for models ZHA A and BPS are as follows (for the ZHA D and SS models the matrices are 16 × 16 and 40 × 40, respectively, and are shown in the sample programs, see Additional files 1 and 2):

Model ZHA A:

(16)

Model BPS:

(17)

Voltage protocols optimized for maximal divergence of model responses

The primary criterion for a model selection is its ability to accurately reproduce various types of data. It is the choice of the type of data and the number of different data sets that determine how good and how unique the model is. The models described in the preceding sections were all developed using mostly different voltage-step protocols. As such they reproduce the experimental ionic current for any similar voltage-step protocol very well. Figure 2 shows model currents and the experimental current for the four models. Model parameters are listed in Table 1. Model ZHA A is a notable exception and it clearly fails even such simple test. Of the remaining three it is impossible to tell based on the typical activation and tail protocols which is the best approximation of the underlying molecular kinetics.

Figure 2

Comparison of model currents (triangles) and experimental ionic currents (solid line) for activation protocols for four Markov models of Shaker K⁺channels. The voltage protocol consisted of holding potential of −90 mV followed by a series of steps to potentials from −70 mV to +42 mV in 8 mV increments. Data is shown for models: A) ZHA A, B) ZHA D, C) BPS, D) SS.

Table 1 Model parameters used for current simulations using eq. ( 14 )

We used these models with the parameters optimized for the fit to standard activation and tail protocols. Our goal was to construct new fluctuating voltage waveforms that produce maximal differences in model responses (i.e. model currents) and thus could be used to select one of the models over the others. We used the numerical procedure described in Methods to construct the new voltage input as a superposition of wavelets. After computing the model response to this voltage input, i.e. the model ionic current, using eq. (14), we defined the difference in model response as the χ² error between the corresponding model currents (see Additional file 3). The voltage input was then iteratively optimized through the genetic algorithm to maximize model divergence. We have done this for different pairings of the models (ZHA A vs. SS, ZHA D vs. BPS, etc.) and for different number of wavelet decomposition levels. A sample MATLAB code for ZHA D vs. SS is shown in Additional file 4. We used the Daubechies 8 wavelets with 8, 9, and 10 levels of MRA, which corresponded to the frequency range commensurate with the recording bandwidth of the patch-clamping apparatus. Figure 3 shows three samples of such voltage protocols with different number of wavelet levels, i.e. different frequency content.

Figure 3

Sample voltage protocols optimized for model current difference. ZHA A vs. SS model with A) 8 wavelet levels, B) 9 wavelet levels, and C) 10 wavelet levels.

Comparison of model currents and experimental currents

Ionic currents were recorded from tSA201 cells stably transfected with the Shaker IR potassium channels, as described in Methods. We used the typical activation and tail protocols. The former were obtained for voltage steps from a holding potential of −90 mV to a series of values from −70 mV to 42 mV in 8 mV intervals. For the latter the protocol consisted of a 32 mV prepulse of 30 ms duration followed by a series of voltage steps from −120 mV to 48 mV in 12 mV intervals. The capacitative transients were removed from the experimental data using the standard P/4 method. The wavelet-based voltage pulses, described in the previous section, were implemented on our recording apparatus as described in the Methods. We monitored the bandwidth by obtaining the cell’s capacitance and the series resistance. Only cells with bandwidth exceeding 5 kHz were considered. For capacitative current transients we used the P/2 method, used also in our earlier paper [15], which has been effective in removing these transients even for continuously varying voltages.

Comparing model currents and ionic currents from electrophysiological experiments is the ultimate test for any model of channel gating. In most studies the fit is limited to the activation and tail ionic currents. Sometimes other voltage protocols are added but they typically consist of few voltage steps at discrete times. As noted earlier (see Figure 2) with the exception of ZHA A all models seem to adequately fit the standard sets of data for stepped-voltage protocols. We also computed responses of models ZHA A, ZHA D, BPS and SS to the wavelet-based voltage pulses using equation (14). In Figure 4, these model currents are plotted against the experimental ionic currents recorded for these same voltages waveforms. As expected the ZHA A model does a poor job approximating the experimental data. The new result is the BPS model which fit the activation and tail data very well but fails to fit the wavelet-based protocols. It suggests that the model has deficiencies and should be discarded in favor of models ZHA D or SS. As for the latter two this set of protocols still does not distinguish well between them, however, the algorithm can be used again to develop new set of pulses that would maximize differences in responses of these two remaining models. These protocols would then need to be implemented in another patch-clamping experiment and we one of them would prove superior to the other. Figure 5 shows the responses to a higher frequency voltage pulse constructed using 10 wavelet levels. While it still does not entirely distinguish between models ZHA D and SS, it shows that the former overestimates a little the rate of activation while the latter seems to underestimate steady-state currents.

Figure 4

Model currents for wavelet-based voltage pulses. Responses of the four models to voltage pulses from Figure 3 A) are shown. Model currents (solid line) are compared to the experimental ionic currents (dotted line). Data is shown for models: ZHA A, ZHA D, BPS, and SS.

Figure 5

Model currents for wavelet-based voltage pulses. Responses of the four models to voltage pulses from Figure 3C) are shown. Model currents (solid line) are compared to the experimental ionic currents (dotted line). Data is shown for models: ZHA A, ZHA D, BPS, and SS.

Discussion

The goal of this paper is to describe a new method for developing voltage protocols for ion channel patch-clamping experiments. It is expected that in order to obtain new information about channel gating kinetics or to explore new phenomena that ion channels may exhibit, one must go beyond the typical set of stepped voltage protocols. The method we develop allows tailoring the voltage inputs in a patch-clamping experiment to maximize a desired output. The design is based on the inverse discrete wavelet transform where the voltage pulse is constructed as a superposition of wavelets in a dyadic wavelet basis. Each such pulse is characterized by the corresponding set of wavelet coefficients and it can be synthesized using a variety of available numerical packages. The desired outcome of the patch-clamping experiment needs to be quantified through a suitably defined cost functional and then the design algorithm is essentially an optimization procedure where the set of wavelet coefficients is optimized to maximize (or minimize) the cost functional. The optimization procedure we used is a search of the parameter space through a random genetic algorithm. We used a method where the “daughter” parameter sets are generated by randomly perturbing the parameters of the “parent” set, with a range of random fluctuations decreasing through the generations. This ensures convergence but also avoids trapping in a local extremum of the cost functional. We used this method previously e.g. for fitting of kinetic models of channel gating. Other algorithms are possible, for instance [32] mentions generating daughter sets by an algebraic formula that combines randomly generated intermediate sets and the “parent” set. Also, in our application we used the Daubechies 8 wavelets, however other types, as well as wavelet packet could be used. An interesting case would be using the Haar wavelets as they would be a natural extension of the voltage-step protocols currently dominating ion channel experimental practice.

This method has many potential uses but we apply it to a specific example: testing of kinetic models of ion channel gating. It is known that the models are notoriously ambiguous. Different research groups propose very different models that match available data equally well, hence there is no way to decide which model (or maybe none) is the best approximation of reality. It is accepted that in order to improve the model design and testing process, we need to include more data and that new data should be substantially different from all existing data. Just adding more of the same will not improve the process. There are two approaches to this. One is to get new data that is of different physical nature. The best illustration would be to add gating currents to previously used ionic currents, or to use single channel data in addition to whole-cell data. This is not the focus of this manuscript. We consider the second approach where the new data is of the same physical nature as previously (i.e. the same physical quantity is measured) but the input to the system has been substantially altered. In recent years there has been interest in using fluctuating voltages as opposed to piecewise constant voltages [10–17, 29, 30], as well as using conductance hysteresis curves to differentiate between different models [39, 40]. The rationale is that this new type of stimulus causes a nonequilibrium response of the ion channel ensemble and some models that were built based on equilibrium (or near-equilibrium) experimental data will fail to reproduce the nonequilibrium data. The main issue with this approach is to decide what type of fluctuating voltage will work best. In the previous studies on this method (named the Nonequilibrium Response Spectroscopy) [10, 11, 17] the dichotomous noise was used but there was no systematic way to determine the properties of the noise that would allow the most efficient model evaluation. In fact, the choices were fairly arbitrary and amounted to scanning a range of noise amplitudes and frequencies and hoping to find the “right” range. The essential parameters are the range of fluctuations and the spectro-temporal properties, i.e. what frequencies are present in the signal and at what times (the signal need not be stationary). The method we described allows the selection of these properties and is essentially an optimization procedure. We have selected four models proposed by leading electrophysiology labs and defined the cost functional as the difference between the model currents generated by these models, measured as the χ² error. Then we constructed several pulses that maximize that cost functional, i.e. pulses which when input into these models generate maximally divergent outputs. These model currents were then compared to the experimental currents recorded in patch-clamping experiments, in response to the same voltage inputs. This facilitated the selection of the best models, and the process can be done iteratively. A preliminary version of this method was previously applied to models of a human heart sodium channel [16]. This is, however, only one possible application of this method. Another putative application would be to explore the possibility of kinetic focusing of ion channels. This term refers to a selective enhancement of a particular state of ion channels in response to a fluctuating stimulus in terms of the applied voltage. In the first paper describing the phenomenon [29] it was suggested that it could be generated by a dichotomous noise voltage and a preliminary experimental study [30] also used that same type of voltage stimulus. The method we describe here would make it possible to design rapidly fluctuating voltages that would maximize the occupancy of a particular conformational state of the ion channels. Further theoretical and experimental studies are in progress [31]. Finally, we mention that a related method has been proposed in engineering for optimizing outputs of nonlinear systems subject to energy constrains, and illustrated with examples of an aircraft model and of a bistable potential system [32].

Conclusions

In this article we describe a new method that represents a radical departure from the current paradigm in ion channel electrophysiology. We propose using fluctuating voltages through a patch-clamp apparatus. While the use of fluctuating inputs to nonlinear systems to elicit new types of system responses is well known and researched in different areas of physics, our method has the added advantage of allowing tailoring inputs to achieve a particular, most desired outcome. If we quantify the “desirability” of outcome through a suitable “cost functional”, through our method one can construct an input that optimizes this cost functional. The pulse design is based on the wavelet decomposition. In particular, we implement the method in ion channel electrophysiology where the input is the voltage applied in patch-clamp experiments, the output – the measured whole-cell ionic current and the “cost functional” – the fit of various Markov models to the experimental data. Our method can be used to aids selection and testing of models of channel gating kinetics. We also comment on other possible applications of this method to electrophysiology and other fields.