History-dependent dynamics in a generic model of ion channels – an analytic study

Recent experiments have demonstrated that the timescale of adaptation of single neurons and ion channel populations to stimuli slows down as the length of stimulation increases; in fact, no upper bound on temporal timescales seems to exist in such systems. Furthermore, patch clamp experiments on single ion channels have hinted at the existence of large, mostly unobservable, inactivation state spaces within a single ion channel. This raises the question of the relation between this multitude of inactivation states and the observed behavior. In this work we propose a minimal model for ion channel dynamics which does not assume any specific structure of the inactivation state space. The model is simple enough to render an analytical study possible. This leads to a clear and concise explanation of the experimentally observed exponential history-dependent relaxation in sodium channels in a voltage clamp setting, and shows that their recovery rate from slow inactivation must be voltage dependent. Furthermore, we predict that history-dependent relaxation cannot be created by overly sparse spiking activity. While the model was created with ion channel populations in mind, its simplicity and genericalness render it a good starting point for modeling similar effects in other systems, and for scaling up to higher levels such as single neurons which are also known to exhibit multiple time scales.

Introduction

Many recent experiments have demonstrated that the timescale of adaptation of a single neuron in response to periodic stimuli slows down as the period of stimulation increases (Fairhall et al., 2001; Lundstrom et al., 2008; Wark et al., 2009). At a sub-neuronal level, experiments on sodium (Toib et al., 1998; Melamed-Frank and Marom, 1999; Ellerkmann et al., 2001) and calcium (Uebachs et al., 2006) ion channel populations have shown that the timescale of the recovery from inactivation following a long duration of membrane depolarization increased with the length of the depolarization period. We refer to this type of behavior as history-dependence. Finally, patch clamp experiments on single ion channels have hinted at the existence of a large inactivation state space within a single ion channel (Liebovitch and Sullivan, 1987; Millhauser et al., 1988; Marom, 1998 and the references therein). These multi-level experimental findings raise several important questions. How are the behaviors observed at the different levels related (e.g., Lowen et al., 1999)? Specifically, is there a connection between the history-dependent timescale of adaptation in the neuron to the history-dependent behavior of ion channels? Does a multitude of inactivation states create the observed channel behavior? What is the functional significance of this history-dependent behavior (e.g., Wark et al., 2009)?

Although we do not address all these questions in this paper, we believe that in order to begin addressing them we first need to construct a simple working, and mathematically tractable, model of slow inactivation in ion channels. Such a model must reproduce the long term behavior in channel population experiments. Our main focus here is the experiment in Toib et al. (1998), which was performed on transfected oocytes and provides very clear empirical findings. In that experiment a membrane with a population of sodium channels of a single type was clamped at low voltage (−90 mV), then at high voltage (−10 mV), for varying lengths of time – from 10 ms to 5 min. During the high voltage stimulus, the sodium channels entered inactivation. Since the fraction of inactivated channels determines the membrane conductivity, by measuring the membrane current, it was possible to observe the dynamics of slow inactivation and recovery in the channel population. After stimulating the membrane with the high voltage clamp for a duration of t_stim seconds, the voltage was decreased and clamped back at the low value (−90 mV). At this low voltage level the channels recovered from inactivation. For short stimulations (t_stim < 1 s), the recovery was exponential with a single non-history-dependent timescale. After sufficiently long stimulations (t_stim > 1 s) the recovery was distinctly exponential and history dependent, the timescale of recovery monotonically increasing with the length of the inactivation period, as seen in Toib et al. (1998). Interestingly, early experiments on visual receptors already displayed similar behavior (Baylor et al. 1974; in particular see Figures 18A,B).

This history-dependence is generally thought to result from the large inactivation state space hinted at by the single channel patch clamp experiments, as suggested first by Toib et al. (1998). Previous modeling approaches, based on this idea, have already been suggested in the literature, but fall short in accurately reproducing this behavior. We present a comparative discussion in Section ‘Relation to Previous Work’. One difficulty in modeling channel behavior is that the nature of the protein conformation dynamics leading to the complex properties of ion channels at long timescales is currently ill-understood, precluding the construction of a full bottom-up biophysical model of ion channels. In fact, it is unclear whether such dauntingly complex low-level models would be useful in explaining phenomena at the level of current interest.

In this work we present a simple two-state generic mathematical model which requires very few assumptions on the nature of the inactivation state space, and which leads to concise explanations of observed experimental findings, and to concrete predictions for future experiments. We reproduce for the first time, to our knowledge, the main experimental finding from Toib et al. (1998), namely an exponential recovery process with a history-dependent timescale, as demonstrated in Figure 3. Using these results and other similar experiments (Ellerkmann et al., 2001; Hering et al., 2004; Uebachs et al., 2006), we narrow down the options for the model parameters at different voltages for several channel types, and explicitly address the issue of long-memory phenomena (Mercik and Weron, 2001). The model introduced here also provides many predictions. Qualitatively, we predict that temporally spaced spiking stimuli will have a significantly reduced effect on the timescale of channel recovery from inactivation, and that the rate of this recovery must be voltage dependent. Quantitatively, we derive a dynamic equation that fully defines an input–output relation between the membrane voltage and channel availability and solve it exactly in many important cases. Additionally, we develop expressions that describe all joint moments in the single channel and population.

We note that the potential contribution of this model goes beyond the specific system addressed in this work. As pointed out in Marom (2009, 2010), current models of channels and receptors (e.g., Faber et al., 2007) tend to suffer from an embarrassment of riches. In order to explain behaviors over an ever expanding range of timescales, these complex models often include multiple inactivation states. Since the number of states and their parameters are not directly observable, these models tend to be highly specific and are likely to suffer from over-fitting. Furthermore, such models always have an upper bound on their timescale. In this work, we introduce and thoroughly analyze, for the first time to our knowledge, a type of model that does not suffer from these limitations. Despite its simplicity, it provides a generalization of previous models, is based only on measurable quantities, does not possess an upper bound on its timescale and exhibits considerable analytical tractability. As such, it stands as an appealing alternative to previous approaches, and as a basic building block in the construction of higher level neuronal models. For example, the work of Lowen et al. (1999) clearly demonstrates the strong impact of a similar model at the channel level on the long term statistics of neuronal firing.

The outline of the article is as follows. We begin in Section ‘The Model’ by motivating the model and describing its structure, and then present several exact solutions for the channel dynamics for different types of voltage inputs in Section ‘Response of Channel to Different Types of Voltage Input’. In Section ‘Reproduction of Experimental Results’ we demonstrate how to reproduce experimental results from Toib et al. (1998), then we use these results to narrow down the possible parameter values, and make specific predictions. In Section ‘Temporal Correlations’ we discuss temporal correlations and long-memory, and in Section ‘Relation to Previous Work’ we explain the relation of this paper to previous theoretical work. Next, we summarize and discuss our results in Section ‘Discussion’, and present the mathematical details of the analytical methods used in Section ‘Methods’. In the Supplementary Material we address several technical issues and provide the simulation code.

Results

The Model

Background

The total conductivity, G(t), of a population of ion channels located on a piece of membrane is determined by A(t), the portion of channels which are available for conducting, through the equation G(t) = G_maxA(t), where G_max is the conductivity of the membrane when all channels are available. In the limit of many channels which are independent given the voltage, the law of large numbers implies that A(t) is approximately equal to p(t), the probability of a single channel to be available (see Supplementary Material). Thus, in order to understand the conductivity dynamics of a population of channels, we first calculate the probability of a single channel to be available, given the input voltage.

Ion channels are commonly modeled as a continuous-time Markov process with a finite, possibly large, number of states (Colquhoun and Hawkes, 1977). Each state is completely characterized by the probability density function of its residence time (the time required to leave that state) and the probabilities to shift to other states when this transition occurs. This division into ‘states’ implies that the dynamics of a single channel following a transition from one state to another is independent of the history prior to that transition. The term ‘Markov process’ implies that this ‘memoryless’ behavior is also maintained in each state, namely, the transition out of the state is independent of the time the channel already resided in that state. Slightly abusing notation, we refer to a state as ‘Markovian’ if the dynamics in that state are memoryless. For a Markovian state, the transition rates leading out of this state are constant, and the resulting residence time probability density function (RTPDF) is exponential.

In these common Markovian models, the states are divided into groups of ‘open’, ‘closed’ and ‘inactivated’ states. The channel may conduct ions only in an ‘open’ state, thereby actively participating in the generation of an action potential. As noted in Marom (1998) and Marom (2009), it is possible to lump together into a single Markovian ‘available’ state all the states from which the channel may change into ‘open’ quickly enough to participate in the creation of a single action potential. This is possible because transitions between the states that compose the available state are much faster than the transitions to and from non-available states (see Supplementary Material). In contrast, if we lump together all the remaining slow inactivation states into a single state, generally it is not a Markovian state. Therefore, we need to switch from the familiar concept of a Markov process to the more general ‘semi-Markov process’ in which some of the states may be non-Markovian (for exact definitions and properties, see Cinlar, 1975).

Model description

In view of the above comments, motivated by Liebovitch (1989) and Marom (2009), we model the slow inactivation of a channel as a continuous-time semi-Markov process consisting only of two states, the Markovian ‘available’ state, and the non-Markovian ‘inactivated’ state, as shown in Figure 1.

FIGURE 1

Figure 1. Our two state model for slow inactivation and recovery of an ion channel is completely described by ψ(t), the residence time probability density function (RTPDF) of each state.

The RTPDF of the Markovian available state is exponential with parameter γ,

The inactivated state is non-Markovian, where we use the power-law RTPDF:

where γ and c are, in general, voltage-dependent, while t₀ is taken to be constant. The way the RTPDFs change with voltage is explained in Section ‘General Voltage Input’. From the normalization demand imposed on the RTPDF, γ, c and t₀ must all be strictly positive constants.

We comment that although the structure of the model is indeed very simple, it gives rise to an unexpected richness of behaviors as a function of the underlying parameters. This exact and detailed mathematical analysis, covering a broad range of parameter values, is presented here for the first time. Moreover, subtle mathematical issues arise in correctly characterizing the different regimes. None of these complexities appear when dealing with the more standard Markovian channel models.

Model motivation

Here we choose to model the, possibly complex, slow inactivation state space solely through a single non-Markovian state and its associated RTPDF. As stated above, the transition structure and rates by which the channel proteins change conformations between different slow inactivation states are still ill-understood. Moreover, the only quantity related to this state space that may be directly measured is the RTPDF. Previous single channel patch clamp experiments have already measured the RTPDFs of different channels. In many cases single channel patch clamp experiments have indeed fitted their measured RTPDF for the closed time with a similar power-law function as above. However, those experiments were performed on timescales of milliseconds-to-seconds. Thus, at this experimental stage, we can only speculate as to whether the RTPDF is still a power law on the timescale of seconds-to-minutes, which is the relevant timescale in Toib et al. (1998). Our reasons for using the particular two-parameter power-law RTPDF in Eq. 2 are the following.

Any smooth normalizable RTPDF ψ(t) must decay to zero as t → ∞. The asymptotic form of ψ(t) as t → ∞ is referred to as the ‘tail of the distribution’. In experiments done on channels at long timescales, such as Toib et al. (1998), this tail dominates the dynamics of the channel. We sought to investigate a simple RTPDF with a non-exponential tail behavior, allowing for a non-Markovian inactivated state, which may lead to a history-dependent behavior. The power-law RTPDF function in Eq. 2 is the simplest possible form. Such power-law RTPDFs frequently appear in models of disordered systems, e.g., spin-glasses (Bouchaud, 1992). Moreover, a similar power-law RTPDF was used in Lowen et al. (1999) to successfully simulate long term temporal behaviors at the level of a single neuron (in that model a specific voltage-dependent choice of c was made, and the available state had also a power-law RTPDF, rather than exponential).

The parametrization Eq. 2 includes two important elements – c as a parameter that determines the tail of the distribution, and t₀ as a lower boundary on the temporal resolution. For simplicity, we chose it to be constant. It is interesting to note that in limit c → ∞ ψ_P(t) approximates an exponential RTPDF corresponding to a Markovian state (see Supplementary Material).

Response of Channel to Different Types of Voltage Input

In this part we present several results on the behavior of the channel in response to different types of voltage input. All of the these results are derived and proved using analytical techniques (see Methods) and demonstrated numerically.

Relaxation to a steady state under a step voltage

First we investigate what happens to the channel when the voltage is maintained fixed, so that c and γ are constant. By ‘projecting’ the single inactivated state onto a continuum of Markovian states, we derive in Eq. 16 a dynamic equation for p(t), the probability of occupying the available state,

Here we used the initial condition p(0) = 1; an extension to general initial conditions is presented in Eq. 16. The intuition behind this equation is simple: the first term represents the loss of probability from the available state, while the second term represents the probability current that goes from availability to inactivation and back again, where we integrate over all the possible past inactivation times. Note that ψ_P(·) can generally be replaced in this equation by any RTPDF of the inactivated state, but we focus in our analysis on the case of the power-law RTPDF from Eq. 2.

We solve Eq. 3, and find that p(t) relaxes asymptotically in a power-law manner, to a steady state value,

where for 0 < c < 1, for c > 1, and p_∞ is the steady state value:

We can express p_∞ in the more familiar form of the steady state of a simple two-state Markov process, where is, as in the Markov scheme, the inverse of the mean inactivation time:

In Figure 2 we compare these asymptotic results with a numerical simulation. The agreement between the two results in this case indicates that the asymptotic results are accurate for a wide range of timescales.

FIGURE 2

Figure 2. Behavior of A(t) – comparison between numerical simulation and the analytical expression of the asymptotic behavior of A(t): (A) c = 0. 5, (B) c = 1.5. The remaining parameters are: N = 10⁶, γ = 1 Hz, t₀ = 1 s and a simulation step of Δt = 5 × 10⁻³ s. The agreement between the analytical and numerical results in this case demonstrates that the asymptotic analytical results are accurate for a wide range of timescales.

Notice that for c ≤ 1, for any choice of the other parameters, p(t) → 0. In other words, the channel eventually decays to complete inactivation, independently of the value of the other parameters (recall that γ > 0). Intuitively, in this case the residence time mean is infinite, and the rate of return to the available state is so slow, that ultimately it remains unoccupied. Assuming that real channels at rest voltage do not decay to complete inactivation, we conclude that c > 1 at rest voltage, namely c(V_rest) > 1.

History-dependent recovery timescale in response to voltage pulse

In Section ‘Relaxation to a steady state under a step voltage’ we studied the channel dynamics during a constant voltage step. The major feature of the model introduced here is its history-dependent dynamics, a notion which cannot be investigated for such a simple input. In order to quantify this notion, in the present section we consider a voltage pulse of constant amplitude and finite duration t, studying the recovery process immediately following the termination of the pulse (similarly to was done experimentally in Toib et al., 1998).

Consider a single channel at time t (the end of the voltage pulse). When the voltage is changed abruptly from one fixed value to another, c and γ also change. If during this voltage change the channel is in the available (Markovian) state, it does not ‘remember’ its history prior to the voltage change. If the channel is in the inactivated state during the voltage change, the subsequent dynamics depends on the prior history of the channel through the variable T – the duration of the channel’s sojourn in the inactivated state at time t. For the inactivated channel, the probability to recover at times between t and t + dt depends on T. This probability, divided by dt, is termed the time-dependent rate of recovery, and its inverse is the time-dependent timescale of recovery, τ. A simple derivation (Eq. 35) shows that τ depends linearly on T,

And so, immediately after a voltage change, the timescale of recovery is determined by Eq. 7. In this model the variable c changes instantly with voltage, while T is a continuous variable which retains the same value it had prior to the change. We note here that it is possible to model the effect of the voltage change on the time-dependent rate of recovery differently (see General Voltage Input).

For each channel, T is a random variable, and so in a population of channels, T has some distribution of values. We develop in Eq. 36 an exact asymptotic expression for this distribution, in the voltage pulse setting, where we set the initial condition so that all channels are initially available. We define 〈T〉, the mean duration of a channel in the inactivated state, and the coefficient of variation of the distribution of T, given by the ratio between the standard deviation and the mean of the distribution of T. The latter variable measures the dispersion of the distribution of T around its mean.

To approximate the timescale of recovery of the channel population after the voltage pulse, we can substitute T by 〈T〉 in Eq. 7. If CV_T ≪ 1 then this approximation is accurate, and the recovery after the pulse will follow an exponential form, to a good approximation. If CV_T ≫ 1 this approximation will not be accurate, since then the mean provides a poor representation of the distribution. Also, when CV_T ≫ 1 the recovery will be distinctly non-exponential.

In Section ‘Renewal Theory Approach’ we compute exactly the behavior of 〈T〉 and CV_T at the end of the voltage pulse, at time t. This behavior depends in an intricate way on the value of the parameter c during the pulse. And as we explained, the value of 〈T〉 and CV_T at the end of the pulse determines, through Eq. 7, the manner in which the channels recover immediately after the pulse:

• c < 1: 〈T〉 = (1 − c )t, The recovery timescale increase linearly with t, while maintaining a constant dispersion. Also, Since CV_T increases monotonically with c, the distribution of timescales becomes less dispersed for lower values of c, and hence, the recovery following the pulse becomes more exponential.

• 1 < c < 3: CV_T → ∞. Thus, the distribution of recovery timescales broadens constantly, which entails a non-exponential recovery (in this case, the mean, 〈T〉, provides a poor representation of the distribution).

• c > 3: 〈T〉 = t₀/(c − 2), Since both 〈T〉 and CV_T converge to a fixed finite value, we get that the recovery following the pulse takes place at a single timescale. Higher values of c, lead to more exponential recovery. This implies that in this case the inactivated state is ‘almost Markovian’.

The (technical) case of integer valued c is discussed in the Supplementary Material. Observe that the qualitative change of behavior noted above, results from the order of the first infinite moment of the RTPDF ψ_P, which depends on the value of c. More specifically, the first moment diverges for c ≤ 1, the second moment diverges for c ≤ 2, and so on.

Complex voltage input

So far, we have described the asymptotic channel dynamics when the voltage is constant, and also immediately after the voltage has jumped from one constant value to another. Here we briefly discuss two other cases in which the channel dynamics may be solved analytically, and then discuss the general case of time-varying voltage.

The first case corresponds to the adiabatic limit, when the voltage changes on a timescale which is far slower than the timescales of the channel, namely the timescale of inactivation γ⁻¹, and the timescale of recovery τ. Notice that the timescales are themselves voltage dependent, and τ may even be history dependent, so it is not always trivial to determine whether we are indeed in the adiabatic limit. In this limit, we may assume that both p(t) and the distribution of T follow their steady states.

The second case occurs in the opposite limit, in which the voltage oscillates rapidly around some constant value with an effective period of T_P. By using the word ‘effective’ here we allow the oscillation to be stochastic (noise) with a period which is only approximate. In this case, we show in Eqs 55 and 56 that we can replace γ and c with the ‘effective parameters’ and given by the time-averaged values of γ and c, respectively,

where V_t is the voltage at time t. Then we may again use the results derived in the case of constant voltage.

Finally, when none of these approximations is valid, we show in Eq. 43 that it is possible to derive a closed form integro-differential equation generalizing Eq. 3 to the case of arbitrary input.

Reproduction of Experimental Results

As observed in Section ‘Response of Channel to Different Types of Voltage Input’, the channel dynamics depends intricately on the parameter c. In the present section we consider experimental results relating to channel dynamics, in order to test the predictions of our model, and, importantly, constrain the possible values of c.

Exponential and history-dependent relaxation in response to long depolarizations

Using the model presented, we wish to reproduce the main experimental results in Toib et al. (1998). In this experiment, the membrane is clamped at a high voltage of −10 mV for varying length of time, t_stim, and then clamped at low voltage of −90 mV – in which the recovery from inactivation occurred. This recovery was exponential and history dependent for every stimulus longer than 1 s. Analyzing these experimental observations using our model, we are able to restrict the possible values of the parameter c. We denote by c_H, γ_H and c_L, γ_L the model’s parameter values during the high and low voltage phases, respectively. Based on the experimental results from Toib et al. (1998) we make the following claim.

Claim: During inactivation c = c_H < 1, while during recovery c = c_L > 3.

This claim is based on the characterization of the different regimes of c, provided in Section ‘History-dependent recovery timescale in response to voltage pulse’. First, we infer that c_L > 1. This follows from Eqs 4 and 5, since for c_L < 1, the availability of the membrane declines to zero, a phenomenon which did not occur in the experiment. Second, we argue that c_L is not in the range 1 < c_L < 3, since if it were, the distribution of recovery timescales in the population would keep broadening with time, thus causing the recovery of A(t) to be non-exponential. Thus, we conclude that c_L > 3.

Next, we argue that c_H < 1. First, it is clear that c_H < 3, since from Eq. 7 and the behavior of 〈T〉 in this case, we know that during the inactivation period, 〈T〉 is bounded by its steady state value 〈T〉 = t₀/(c_H − 2). By substituting this value into Eq. 7, and comparing with the result when 〈T〉 = 0, we learn that in this case we should not expect to see large changes in the recovery timescales following different stimuli durations, contradicting the experimental results. Finally, c_H cannot be in the range 1 < c_L < 3 since, if it were, the distribution in the inactive state following the high voltage period would be very broad (a high value of CV_T), rendering an exponential relaxation impossible.

We thus conclude that in order to fit the experimental results, we must set c_H < 1 during the inactivation period and c_L > 3 during the recovery phase. In this case, the timescale of recovery increases linearly with the stimulus duration, and the recovery is exponential as long as:

1. c_H is small enough so that CV_T ≪ 1, since which is increasing in c_H.

2. c_L is large enough in the low voltage phase so that 〈T〉 does not change considerably during the recovery. More formally, T ≫ τ = (T + t₀)/c_L, or more simply, c_L is sufficiently larger than 1.

3. γ_L is small enough during the recovery period, so that during this period, the timescale of recovery does not change by ‘freshly inactivated’ channels; more formally,

In this method of reproducing the experimental results the channel possesses ‘infinite memory’, since as long as the channel remains at the high voltage, 〈τ〉 continues to increase linearly:

We note here that it is possible to modify the model so that t₀ is the voltage-dependent parameter and c is fixed at some value ( >1). In this case, the channel does not possess this ‘infinite memory’ and the timescale of recovery is bounded. Therefore, in order to reproduce the experimental results in this modified model, we have to assume that the increase in the timescale of recovery has some upper boundary, which was not yet reached in the experiment. Since it seems unnatural to assume the existence of such upper limit, this alternative model was not used.

In an additional experiment the recovery was examined under several voltages: −60, −90, −120 mV. The history-dependent behavior remained, as can be seen in Figure 6 in Toib et al. (1998). The recovery was similar for −90, −120 mV, indicating that c did not change much between these voltages. At −60 mV, the recovery was slower, and possibly less exponential, indicating perhaps that 1 < c_L < 3, in that case.

A further experimental observation relates to the emergence of a history-dependent relaxation only for t_stim ≥ 1 s. From this threshold point between the two modes of recovery and the results in Section ‘History-dependent recovery timescale in response to voltage pulse’ we get that in this experiment t₀ ∼ 1[s].

In any case, the resulting prediction of this model is that the rate of recovery from the inactivated state must decrease with voltage difference between the two values examined in the experiment. This is in accordance with the results of Figure 7 in Fleidervish et al. (1996) and Figure 4 in Ellerkmann et al. (2001), where it is observed that the timescale of recovery of slow inactivation in sodium channels increases monotonically with voltage. Moreover, if we assume, in accordance with these results, that c(V) is continuous for these sodium channels, a further prediction of our model is that there exists some range of voltages for which 1 < c_H < 3, where the recovery becomes non-exponential.

Channel response to a spiking stimulus

Motivated by the results in Ellerkmann et al. (2001), Toib et al. (1998) and Uebachs et al. (2006), we study the model’s dynamics when the voltage input is a periodic voltage spike train. As noted in Toib et al. (1998), such inputs are similar to firing patterns in neocortical neurons. This type of input may be important when considering the effects of a neuron’s action potential on itself. We note that such an input does not present the entire picture, since synaptic inputs from other neurons are probably more realistically described as sums of continuous functions (Gerstner and Kistler, 2002). In any event, it is interesting to test the model’s prediction in this setting, for which some experimental results are available (Toib et al., 1998; Ellerkmann et al., 2001; Uebachs et al., 2006). The voltage spikes are modeled here as a square wave – for T_H seconds the voltage is set high, and for T_L seconds the voltage is set low. As before, we denote by c_H, γ_H and c_L, γ_L the model parameters during high and low voltage, respectively (see Figure 4).

FIGURE 3

Figure 3. Simulation of experimental results. (A) Exponential and history-dependent recovery in response to various stimulus lengths (R² > 0.99 for all exponential fits). Fraction of inactivation is defined as in Toib et al. (1998). Legend: different stimuli lengths, in seconds. (B) The increase in the timescale of recovery with stimulation length, in both the numerical simulation and the analytic expression (Eq. 9) (with c_H < 1). We see that the two results agree, especially for long stimulations. This demonstrates how, in this case, the asymptotic analytic result is a good approximation for the relevant timescales. The linear relationship is expected, based on the analytic result, to persist for arbitrarily long stimuli. The parameters for the simulation were N = 10⁵, γ_L = 10⁻⁴ Hz, γ_H = 1 Hz, c_H = 0.2, c_L = 15, t₀ = 1 s and a simulation step of Δt = 10⁻³ s.

FIGURE 4

Figure 4. The voltage and model parameters, during a spiking stimulus. (A) Each stimulus period is divided into the time of the spike (T_H), and the inter-spike interval (T_L). (B) The parameter γ assumes the values γ_H during the high voltage, and γ_L during the low voltage. The value is the “effective parameter” – the average of the parameter – during a single period. As can be seen in the figure, in the case that T_H ≪ T_L, then can be significantly larger than γ_L if γ_H ≫ γ_L. This implies that the spikes may change the rate of inactivation. (C) Same as (B), except that the parameter c is displayed. As can be seen in the graph, in the case that T_H ≪ T_L, then . This implies that the spikes only marginally affect the rate of recovery.

Assuming that the input spike frequency is much higher than the transition timescales of the model, we use the effective parameter approximation. The behavior of the system may be inferred by replacing γ and c with the effective parameters, according to Eq. 8:

It is important to notice in this case that if T_H ≪ T_L (sparse spikes), then will be affected by γ_H only if γ_H ≫ γ_L, and similarly for Since we know from the experimental results that c_H < c_L, in this case we can approximate, Using this approximation we address the dependence of γ on the voltage. From Toib et al. (1998) it is known that the sodium channel population goes into significant inactivation as a result of an action-potential-like stimulus. This is also the case in Ellerkmann et al. (2001). Since in this case (does not change much), then must be significantly larger than γ_L – so γ_H must be significantly higher than γ_L (otherwise, the experimentally observed inactivation would not be reproduced). This means that γ(V) must increase with voltage, at least between the two voltage values used in this setup.

The distribution of recovery timescales in the channel population is affected mainly by the value of c (the dependence on γ is absent from the asymptotic form Eq. 36). If we wish to make the timescale of recovery change measurably in response to spike stimulation, must be significantly different from c_L. For example, if a neuron is affected solely by its own action potentials, then T_H ≈ 1 ms, T_H + T_L > 10 ms, and therefore So, in this case it would be hard to experimentally detect an increase of the recovery timescales due to the input (see Figure 5). This assertion will remain valid for any type of stimulus in which T_H ≪ T_L – the channel may enter inactivation, but any history-dependent recovery will be mitigated. This complies with the results shown in Figure 5 of Ellerkmann et al. (2001), where we can see a relatively weak scaling of the recovery timescale with stimulus frequency – even though T_H/(T_H + T_L) is rather high in comparison with typical action potentials.

FIGURE 5

Figure 5. (A) Entry into inactivation for different stimulation frequencies. (B) A closer look at the first 0.5 s of the entry into inactivation, on a logarithmic y-scale: each jump is caused by a single spike, while the overall trend displays an exponential relaxation, with a timescale linear in the frequency of stimulation (through which is determined by Eq. 10) (C) Recovery timescales at different lengths of spiking stimulation and for different frequencies. (D) Recovery timescales at different frequencies for different lengths of spiking stimulations. Parameters: T_H = 2 ms, N = 10⁵, γ _H = 2 Hz, γ_L = 10⁻⁴ Hz, c_H = 0.2, c_L = 5, t₀ = 3 s and a time step of Δt = 10⁻³ s. Notice that for f = 500 Hz we have a continuous stimulation, and the timescale of recovery increases linearly with stimulation length, as seen in Figure 3. For spiking stimuli we observe a weak scaling with frequency (as in Ellerkmann et al., 2001), but observe no long term scaling with stimulation time for spiking stimuli. In graphs (C) and (D) we omitted some of the continuous stimulation data points (where the increase in timescale is very large) in order to facilitate a comparison between the other data points.

It is important to note that in the case the non-Markovian state can be approximated by a simple Markovian state, since in this case 〈T〉 < t₀, and thus does not change drastically due to history-dependent effects. This may happen, for example, if c(V_rest) > 3, and the input to the neuron is sufficiently sparse (so that it does not cause continuous depolarizations for long times).

In Uebachs et al. (2006), a related stimulation experiment was performed on three types of calcium channels, each of which possesses different inactivation properties (which may be described in our model using different values of γ, c and t₀). A slow-down of the timescale of recovery after a long −55 mV depolarizations was observed, similar to what was measured in sodium channels (referred to as ‘continuous’ in Figure 8 in Uebachs et al., 2006). In contrast, after a spiking stimulus at a maximal voltage of 25 mV with T_H/(T_H + T_L) ≈ 0.056 (‘mock APs’ in Figure 8 in Uebachs et al., 2006), there was no observed increase in the timescale of recovery with the length of stimulation, in agreement with our model. Also, when the spiking stimulus was superimposed on the (continuous) step depolarization, the history-dependence of the timescale diminished (‘continuous + mock APs’ in Figure 8 in Uebachs et al., 2006). This may imply that c(V) is not monotone in calcium channels – if c(V) decreases between −95 and −55 mV and then increases between −55 and 25 mV, should change less during the stimulus, and therefore the recovery timescale should also change less. This result might explain why Hering et al. (2004) did not observe any change in the recovery timescale after different lengths of depolarizations, given at the higher voltage of −20 mV.

Temporal Correlations

In Section ‘The Joint Probability Distribution of the Availability’, we present a simple approach to calculating the joint moments and distributions to be in the available state at different times. Specifically, we derive Eq. 57, the asymptotic behavior of the stationary auto-covariance of a single channel, for c > 1:

For 0 < c < 1, k(t) is not defined since in that case the channel is not stationary and decays to complete inactivation. Also, in the Supplementary Material we show that the auto-covariance function of the population availability under voltage clamp is simply k(t)/N, where N is the size of the channel population.

Notice that for 1 < c < 2, we get that This implies that the process is characterized by a long-memory (see definition on page 42 of Beran, 1994). This result complies with the results of Mercik and Weron (2001) for potassium channels (on timescales <1 s), as also pointed out by Goychuk and Hanggi (2004).

Relation to Previous Work

Several previous studies have attempted to construct a mathematical modeling framework which can account for long term temporal correlations, power-law behavior, and history-dependence of responses related to ion channels. The initial work along these lines was motivated by single channel experiments on time scales of milliseconds-to-seconds (see Liebovitch and Sullivan, 1987; Millhauser et al., 1988; Marom, 1998 and the references therein), and therefore does not directly imply the tail behavior of the RTPDF on the timescale of seconds-to-minutes, which was explored in Toib et al. (1998). In this brief comparative discussion we focus on models addressing channel dynamics.

As far as we are aware, the paper by Millhauser et al. (1988) was the first to propose a microscopic framework for generating RTPDFs with power-law temporal dependence, based on a large set of inactivation Markovian states. This work considered RTPDFs of the form t^−α for 1/2 ≤ α ≤ 3/2 which correspond to −1/2 ≤ c ≤ 1/2 in our model. As was mentioned in that work, α ≤ 1 (or c ≤ 0) contradicts the normalization of the RTPDF, and therefore cannot be correct for long timescales. However, for intermediate time scales, similar to those observed in single channel experiments, such behavior may be possible. A main focus of that work was proving that when the transitions between these inactivation states resemble a diffusion process, it leads to a RTPDF of the form t^−3/2, which corresponds to c = 1/2. As stated in Section ‘Relaxation to a steady state under a step voltage’, such a value of c always leads to complete inactivation of the channel in the long run, and therefore cannot be correct for longer timescales. All later models based on this diffusion model suffer from the same type of difficulty. A related line of work was pursued around the same time in Liebovitch and Sullivan (1987). This approach, as well as that of Lowen and Teich (1995), proposed several types of RTPDF, which reduce to a power law in special cases. However, we note that the main concern of Liebovitch and Sullivan (1987), Lowen and Teich (1995) and Millhauser et al. (1988) was the establishment of RTPDFs consistent with single channel experimental findings, rather than on providing an analysis of the long term channel dynamics that is the main focus of this paper.

A more recent line of work (Goychuk and Hanggi, 2004), leaning on earlier work in the statistical physics community, formulated channel dynamics in the form of the so-called generalized Master equation. This work allows for general asymptotic power-law dependence of the RTPDF, and was shown to lead to experimentally measured (Mercik and Weron, 2001) power-law decay of single channel correlations, as we obtain in Eq. 11. However, this work does not directly address the main experimental finding from Toib et al. (1998), namely an exponential recovery process with a history-dependent time scale. Conceptually, two main features distinguish this work from ours. First, input (voltage) dependent parameters, an essential feature of our work, cannot be dealt with through an approach based on the kernels used in Goychuk and Hanggi (2004) (defined through their Laplace transforms). Second, our methods additionally enable the explicit calculation of the channel probability to be available, and of its history-dependence and joint moments (to any order), in many important cases.

The model of Millhauser et al. (1988) was recently used in Gilboa et al. (2005) in order to directly explain the experimental results in Toib et al. (1998). They were able to show, through numerical simulation and approximate analytic solutions of an equivalent diffusion model, that multiple time scales history-dependent behavior (of the type described in Results) is indeed reproduced within the diffusion model. However, the exponential nature of the recovery was not addressed. More recently Friedlander and Brenner (2009) analyzed history-dependent phenomena in the case of RTPDFs with power-law behavior corresponding to 0 < c < 1 in our model, using the framework developed in Goychuk and Hanggi (2004). Specifically, they focused on the mapping between the input (voltage) and the output (availability) for step inputs, demonstrating power-law recovery. Our work corroborates these results, and provides exact analytic expressions for the recovery rates. Moreover, the input–output view is extended to include arbitrary time-dependent inputs and input dependent recovery.

A different line of work was recently suggested in Marom (2009), whereby the complex history-dependent channel dynamics is reduced to a single local in time logistic like equation, where the recovery rate depends on the level of activation in a power-law fashion. This approach can be viewed as a zero order approximation of the full dynamics given in Eq. 43, whereby the entire history is replaced by a single reporter, which is a function of the current availability. More generally, one can envisage replacing the complete dynamics by a finite set of differential equations containing a truncated list of moments of the distribution of recovery time scales. Such an approach, while introducing a further approximation step, falls within the widely studied field of Markovian population dynamics, and offers, due to its mathematical simplicity, the potential of being smoothly incorporated into higher level models of single neurons, viewed as a population of channels.

Along similar lines, the work of Lowen et al. (1999) uses a population of semi-Markovian with power-law RTPDF, coupled to an equation describing the voltage dynamics, to numerically reproduce many of the long term behaviors of single neurons. It is quite reasonable to assume that the long-range correlations reported in Lowen et al. (1999) are the direct result of the power-law RTPDF with 1 < c < 2, which was used for the ion channels in that paper.

Discussion

In this work we have thoroughly analyzed the dynamics of a generic two-state model of an ion channel consisting of a Markovian state and a non-Markovian state. We derived a general dynamic equation for the probability of the channel to be available, or equivalently, the fraction of available channels in a voltage clamped population. This dynamic equation, which is derived for both constant (Eq. 3) and time-dependent (Eq. 43) voltage, defines a direct input (voltage) to output (availability) relation. We derived explicit solutions to this equation in many important cases, and studied their properties. Specifically, we considered the asymptotic power-law approach of the channels to steady state under constant voltage (Eq. 4), the distribution of recovery timescales in the population following an abrupt change in voltage (Eq. 36), and the time-averaging of rapid fluctuations in the voltage (Eq. 8). Also, we derived simple expressions for all joint moments of the availability and specifically, for the auto-covariance function of the channel (Eq. 11). An interesting conclusion of this analysis is that the channel is characterized by four different ‘modes’ in which its behavior is rather different, depending on the value of c (see Table 1).

TABLE 1

Table 1. Modes of behavior for different ranges of c.

The main experimental finding from Toib et al. (1998), an exponential recovery process with a history-dependent timescale, was fully reproduced. We constrained, for several channel types and different voltages, the model parameter space (voltage-dependent c and γ and constant t₀) by using Toib et al. (1998) as well as other channel population experiments (Ellerkmann et al., 2001; Hering et al., 2004; Uebachs et al., 2006). Also, we addressed the issue of long-memory phenomena (Mercik and Weron, 2001).

Our model presents many quantitative predictions on the dynamics of ion channel populations, for different kinds of voltage inputs, especially for long times. These predictions can be confirmed by experiments similar to Toib et al. (1998), but in which voltage stimuli are more variable – in both voltage values and temporal shape. Some of the specific qualitative predictions are the following. The value of the parameter c at rest obeys c(V_rest) > 1 in all healthy channels (see Relaxation to a steady state under a step voltage). From Section ‘History-dependent recovery timescale in response to voltage pulse’ we conclude that the recovery from slow inactivation of NaII and NaIIA channels is voltage dependent, and for these channels, there exists some range of voltages for which 1 < c < 3, where the recovery becomes non-exponential. Finally, we predict that sparse spiking stimuli (e.g., the neuron’s own action potentials) can induce slow inactivation but only minor changes in the timescale of recovery from inactivation (see Channel response to a spiking stimulus).

One of the more interesting questions pertains to the relation between the history-dependent channel dynamics and the generation of action potentials in the cell, and, more concretely, the possible functionality of such behavior (Lundstrom et al., 2008; Wark et al., 2009). In order address this issue clearly and fully, two further steps must be taken.

Experimentally, it is necessary to find the correct model parameters (voltage-dependent c and γ and constant t₀) for different voltages and different types of channels. In particular, the value of c at V_rest and its average value during physiologically realistic cellular stimuli are critical for determining whether history-dependent relaxation and long term temporal correlations actually occur during normal cell activity. If, for example, c > 3 effectively for all channels in the cell, then we should not expect that these phenomena affect action potential generation – since in this case the inactivated state is approximately Markovian.

Theoretically, it is necessary to construct a model of a neuron that incorporates the type of channel studied here, as well as other types of channels which occur in the cell membrane. Using such a model we could determine the nature of the feedback interaction between channel activity and the membrane voltage, and the impact of this interaction on the action potential generation probability (e.g., Lowen et al., 1999; Gilboa et al., 2005). For example, it is quite reasonable to assume that the history-dependent relaxation exhibited at the single neuron level (Fairhall et al., 2001; Lundstrom et al., 2008; Wark et al., 2009) may be caused by the channel mechanism discussed here – especially since slow inactivation in sodium channels is known to have a strong effect on neuronal adaptation at these timescales (Fleidervish et al., 1996; Powers et al., 1999). If this is indeed the same mechanism described by our model, then replacing the continuous stimulation by spike stimulation should greatly reduce such behavior. The persistence of history-dependent relaxation in this setting, would imply that other processes (e.g., multiple interacting channel types) are in place. Another issue that can be explored using this method, is whether possible long term temporal correlations in the channel are the cause of similar long-memory phenomena at the cellular level (Soen and Braun, 2000). In fact, the work of Lowen et al. (1999) argues persuasively along these lines, as many of the long-range temporal behaviors of neurons are well replicated within a simple model for the membrane potential incorporating non-Markovian channel dynamics. Indeed, extending the mathematical tractability of the present approach to the higher level of a neuron is a major theoretical challenge.

Finally, it is important to observe that most of the mathematical results presented in this work are general, and can be extended to other models in which the RTPDF of the unavailable state is not a power law. The mathematical framework that was developed here to model ion channels is quite flexible, and may be used to describe other systems in which one can similarly define separate available and unavailable states so that the transitions between these states are random and depend only on the residence time of the state and the external input (e.g., voltage or ligand concentration). For example, inactivation of cellular receptors could be modeled similarly, as pointed out by Friedlander and Brenner (2009). The method for mapping a general Markov model to our two-state non-Markov model is described in detail in Section ‘How to create a two-state non-Markov model from a general Markov model’.

Methods

In this section we establish several general analytical methods, which can be used in any two-state model in which one state is Markovian while the other state is non-Markovian, namely ψ_E(t) is still exponential as in Eq. 1, but ψ_P(t) may be replaced by some other general RTPDF, which we denote by ψ_I(t). When establishing general results we use ψ_I(t), and when specializing them to the power-law RTPDF we use ψ_P(t).

We present two different formalisms – the Markovian-process based ‘projection method’ in Section ‘‘Projection’ of the Inactivated State’ and the ‘Renewal Theory’ approach in Section ‘Renewal Theory Approach’. This was done in order to explore different aspects of the channel dynamics in the case of constant voltage input. Moreover, we show how to map an existing Markov model to a two-state non-Markov model. In Section ‘General Voltage Input’ we use the presented methods to derive the dynamics in the case of general input. In Section ‘Oscillating Voltage Input’ we consider the case of rapidly oscillating input. Finally, in Section ‘The Joint Probability Distribution of the Availability’ we calculate the joint distributions of the channel.

‘Projection’ of the Inactivated State

In this section we assume that the voltage input is constant, so that all the model parameters are also constant. We wish to make use of the well-established formalism of Markov processes to derive the dynamics of our non-Markovian model. To do so we needed to find a way to replace the non-Markovian inactivated state with an equivalent Markov inactivation state space. For example, in the context of the specific channel model with the power-law RTPDF ψ_P(t), there are many possible state spaces that may be used to produce it (Liebovitch, 1989). Some of them may be physiologically more accurate than others, but, mathematically, this is inconsequential to the channel dynamics. And so, we chose the simplest state space, in which all inactivation states are parallel – each state being connected only to the available state, as seen in Figure 6. Though the figure depicts a finite set of states, in the actual model we use a continuum of states, so that the equivalence between the models is exact. We use as the continuous ‘index’ of each activation state its recovery rate, δ. Since it is Markovian, the RTPDF function of each inactivation state δ is exponential,

FIGURE 6

Figure 6. Approximation of a semi-Markovian channel by a number of Markovian inactivation states.

Each time the channel is inactivated, it goes into one of the inactivation states. We denote by f(δ) the probability density function to go from the available state into a specific inactivation state δ. This means that the inactivation rate from the available state into the inactivated states (δ,δ + dδ) is γf(δ)dδ, and so the total rate of inactivation is as required. With this condition fulfilled this multiple-state Markovian model is equivalent to our two state non-Markovian model if we ensure that the RTPDF of the aggregation of all the inactivated states is ψ_I(t).To fulfill this condition, we use law of total probability and demand that:

Dynamic equations

In order to derive the dynamical equations of the non-Markovian model, we begin with the Markovian model described in Figure 6, present its underlying equations, and then take the continuum limit.

We denote where p is probability of being in the available state, and π_k is the probability of being in the inactivated state I_k. The dynamic equation for the probability of this homogeneous Markov process is

where

is the rate matrix.

Writing this explicitly, we get

In the exact continuum limit we substitute yielding:

Dividing the first equation by dδ we obtain:

with solution

Inserting this into Eq. 13 and switching the order of integration we get:

Using Eq. 12, the expression can be re-written as:

This result is quite intuitive:

1. The term −γp(t) represents the probability current out of the available state.

2. The term represents the probability current that goes from availability to inactivation and back again – where we sum over all the possibilities to inactivate at time t′ and then recover at time t.

3. The term represents the probability current from the initial inactivation states – the probability to be at each inactivation state at time t = 0, times the probability density to recover exactly at time t.

Such equations are commonly dealt with in the Laplace domain. Using a Laplace transform on Eq. 16, we get:

where the tilde denotes the Laplace transform. The final term in the numerator results from the initial probability distribution in the inactivated states. If at t = 0 the channel is in the available state, namely p(0) = 1, then this term vanishes.

Projection of a power-law RTPDF

By using an inverse Laplace transform on the t variable in Eq. 12, we can easily show that in the case where ψ_I(t) = ψ_P(t), f(δ) is distributed according to the gamma distribution,

where Γ(·) is the gamma function. Thus, the aggregation of all the inactivation states has the RTPDF:

which equals ψ_P(t) as required.

Asymptotic solution of the dynamic equations for power-law RTPDF

In this section we demonstrate how to derive an asymptotic solution of Eq. 17 in the case where ψ_I(t) = ψ_P(t). We assume throughout this section that c is non-integer; the limiting case where c is an integer is discussed in the Supplementary Material. In order to understand the behavior of p(t) for large values of t, we need to consider the small s limit of To do so, we first calculate In the Supplementary Material we prove that for non-integer values of c,

Retaining the leading order in s, we have (see Supplementary Material):

and substituting Eq. 20 into Eq. 17, we find:

If c < 1,

while if c > 1,

where we have introduced

denoting the steady state value of p(t) – a result that can be confirmed by using the final value theorem on Also, Notice that in Eq. 22 we deliberately retained the s^c−2 term, even though s^c−2 = O(1) for c > 2. The reason for this will become clear soon.

We wish to find the asymptotic behavior of as t → ∞. To do this we use a Tauberian theorem (Theorem 4 on page 446 in Feller, 1971), stating that f(t) ∼ [K/Γ(α)]t^α−1 as t → ∞ if, and only if, ℒ[f(t)(t) ∼ Ks^−∞ as s → 0, where K is some constant and α > 0.

For 0 < c < 1, a direct application of the theorem to Eq. 21 gives:

where we used the result Γ(c)Γ(−c) = −π/c sin(πc) (page 3 in Erdeyi, 1953, vol. I).

For 1 < c < 2, we first subtract from Eq. 22 the Laplace transform of the steady state value, p_∞s⁻¹, and then use the Tauberian theorem. From the linearity of the Laplace transform we get:

where we used the result Γ(−c)/Γ(2 − c) = −1/(c(1 − c)) (page 3 in Erdeyi, 1953, vol. I), and Eq. 23.

For c > 2, we can again subtract the steady state term p_∞s⁻¹ from Eq. 22, but now we are left with only positive powers of s, and therefore we cannot directly apply the Tauberian theorem. To circumvent this problem, we define m = ⌈c − 2⌉, the upper integer part of c − 2, and differentiate m times:

Since that all the integer order terms with degree lower than m vanish after we differentiate for m time. Using the fact that:

the linearity of the Laplace transform and the above Tauberian theorem, we get:

where we used similar identities for Γ(·) as used previously, Eq. 23 and the fact that m is an integer. Notice that we can use this method to find higher orders of the asymptotic behavior of p(t) from the higher fractional order terms in s.

An interesting implication of the above derivation is that the asymptotic behavior of f(t) is completely determined by the minimal fractional power of s in ℒ[f(t)] (assuming that ℒ[f(t)] can be uniquely represented as a countable sum of integer and fractional powers of s). Because of this, any (minimal) linear time invariant system will behave asymptotically as a power law if it contains some component with a power-law kernel – which always introduces a fractional power of s into the system’s transfer function in the Laplace domain.

In conclusion, we found the asymptotic behavior of p(t) as t → ∞, for all c,

We derived this result for the case p(0) = 1, but it remains valid also when p(0) ≠ 1 if the expression is analytic at the origin. For example, this always is true if for all δ < δ₀, π_δ(0) = 0 for some δ₀ > 0 – meaning that there exists a finite supremum to the recovery timescales of all the inactivation states in the initial distribution. There is another way to see that the asymptotic solution in Eq. 24 is not too sensitive to the initial conditions. Suppose we set p(t_i) = 1 for some t_i, instead of p(0) = 1. The asymptotic solution is p(t − t_i). Since

and therefore Eq. 24 is changed by the addition of an asymptotically negligible term. This asymptotic insensitivity to initial conditions stands in contrast with the case where the inactivated state is Markovian, and the solution is exponential,

In this case when we set the initial condition to t_i, instead of 0, we get,

Therefore, the solution is changed by a pre-factor, which is not asymptotically negligible as in the power-law case.

How to create a two-state non-Markov model from a general Markov model

So far, we used the Markov process formalism to find a dynamic equation for a non-Markov two state model. Now consider the converse problem: suppose we have an existing Markov model for an ion channel, how can we transform it to a two-state non-Markov model?

To answer this question we use a method based on the ideas presented in Colquhoun and Hawkes (1977). Recall that the dynamics of a homogeneous Markov process are governed by the equation where is the probability vector to be at each state, and Q is the rate Matrix. We divide the Markov state space into two disjoint sets of states, according to their conductivity – A, the available states and I, the inactivated states. These two sets will the two states in the non-Markov model. We divide the probability vector into two vectors according to these sets Similarly, we divide the rate matrix Q into block matrices corresponding to these sets of states:

In order to calculate the RTPDFs of the inactivated states, we introduce a modified Markov process with probability vector and the rate matrix:

The dynamics of the modified process are similar to the original process, except that now the available states have become absorbing – each time the process reaches these states, it remains there forever. From we get,

The solution of Eq. 25 can be written in the form,

where the exponent of a matrix A is defined as

Now assume that meaning that both processes at t = 0 are distributed identically, over the inactivated states. The probability that the original process has jumped to some available state before time t is equal to probability that the modified process is in some available state at time t. This probability can be written as where we defined a column vector of ones with the same number of components as Differentiating this expression, we get the RTPDF of the inactivated state, assuming that the channel has inactivated at t = 0:

where we used Eqs 26 and 27, and represents the probability distribution in the inactivated states immediately after the inactivation. Similarly, the RTPDF of the available state is

where is defined similarly to

Note that ψ_I(t), ψ_A(t) are always a sum of exponents, but after the transformation to the two state model is complete, we can approximate them by some simpler functional form. The form of these RTPDFs (in Eqs 28 and 29) pose a problem, through their dependence on and Assume that the channel inactivated at time t₁, and recovered at time t₂. Generally, may depend on the time that the channel was available before t₁, and t₂ − t₁ may depend on Therefore, the time the channel spends in inactivation may depend on the previous duration of the channel in the available state, which may depend on the previous duration the channel in the inactivated state, and so on? This would prevent us from constructing a model with two independent non-Markov states.

There are some specific cases when this problem vanishes. One of them is when A can be approximated as single Markov state. This happens, as we mentioned in Section ‘Background’, when the rates between the available states are much larger than the rates of inactivation. In this case, as we show in the Supplementary Material,

where we defined and is the stationary distribution in the available set of states, when all the transitions to inactivation are forbidden. This is the two-state model which we aimed for: a Markovian available state, and a non-Markovian inactivated state.

Renewal Theory Approach

Despite its merits, the use of the projection method may be limited in certain cases. The problem lies in the use of π_δ(t), the probability distribution in the inactivation states. Knowledge of the probability to be in an ‘inactivation state’ implies that we have some additional information about the future recovery time of the channel. But, when a channel is inactivated, we know only for how long it has been inactivated, and not to which ‘inactivation state’ it has arrived. Since these inactivation states are not observable, this approach restricts the practical use of the initial conditions π_δ(0) in Eq. 16. And so, in order to quantify the history-dependent behavior of the channel in terms of observable quantities, we turn to a different formalism that uses only the observable T, which is the time during which the channel has already been inactivated. In this section, as the in the last section, we again assume that the voltage input is constant, so all the model parameters are also constant.

For a channel that has been inactivated for time T, the probability of recovery in the next time interval Δt is

where t_R denotes the time it takes the channel to recover since its last inactivation.

Therefore, the rate of recovery immediately after time T is

and the time-dependent timescale of recovery is its inverse

Next, we compute f_t(T), the probability density function of T at time t, for a single channel, which well approximates the distribution of T throughout a sufficiently large population. We did so using a ‘Renewal Theory’ formalism, which is described in Cox (1962). In renewal theory language T is called the ‘backward-recurrence time’ until the last inactivation event.

First, we denote by k_n(t) the probability density function of the time of the n-th inactivation event, assuming that at t = 0 the channel was available (p(0) = 1). Clearly k_n(t) is the probability density function of a random variable which is the sum of 2n − 1 independent random variables – corresponding to all the residence times in the available state, and in the inactivated state, that occurred until the n-th inactivation event. Since the RTPDF of X_m, Y_m are, respectively ψ_E(t) and ψ_I(t), and the probability density function of a sum of independent random variables is a convolution of the probability density functions of all the random variables, we get,

In the Laplace domain, this translates to,

Defining the rate of inactivation events:

we have and in the Laplace domain:

where we used the geometric series summation formula. Notice that and so h(t) = γp(t), which is not surprising, since the rate of inactivation events is expected to be the probability to be in the available state, times the rate of inactivation in that state. We can now write an expression for f_t(T), conditioned on the fact that the channel is inactivated at time t. For a channel at time t to be inactivated for the last T seconds, there must be an inactivation event at time t − T and since that time the channel must not recover. Dividing this by (1 − p(t)), the probability that the channel is indeed inactivated, we get,

The case of power-law RTPDF

In the context of the specific model with power-law RTPDF the rate of recovery immediately after time T is

and the time-dependent timescale of recovery is its inverse:

The distribution of T in this case ψ_I(t) = ψ_P(t) is f_t(T) = (1 – 9(t))^–1γp(t – T)(1 + T/t₀)^–c), (0 ≤ T ≤ t).

Because of the multiplication in the last expression, it is generally difficult to write an expression for f_t(T) in the Laplace domain. We therefore examine f_t(T) in the limit t → ∞. We get, by using the same asymptotic derivation we used in deriving Eq. 24, that:

Note that for 0 < c < 1, for any fixed value of T, the function f_t(T) vanishes for large values of t. This implies that the probability distribution shifts to infinity as t increases; see Figure 7.

FIGURE 7

Figure 7. A sample of f_t(T) at different values of t and c: (A) c = 0. 5, (B) c = 1.5. The circles represent results of the numerical simulation, while the lines are the analytical asymptotic expressions (Eq. 36). The numerical results are based on the same simulation as in Figure 2. Since T is proportional to the recovery timescale, through Eq. 7, f_t(T) is equivalent to the distribution of recovery timescales. Note that for c = 0.5 the distribution shifts rightward as a function of time, as explained following Eq. 36.

Calculating the first two moments of T in each case, defining 〈T〉 as the mean of T and σ_T its standard deviation we obtain the following results by series expansion (see Supplementary Material).

We verified numerically that all of the results for these moments and for f_t(T) itself are valid for the relevant timescales, as is demonstrated in Figures S1 and S2 and in Supplementary Material.

General Voltage Input

In the previous sections we used two different formalisms to analyze our model. One was the Projection Method, and in the other was the formalism of Renewal Theory. So far, these two approaches may seem unrelated, and dealt only with the case of constant parameters. In this section we shall use elements from both formalisms to derive a dynamic equation for p(t) for any time-varying parameters. For example, in the context of the specific channel model with power-law RTPDF, recall that the voltage input influences the channel through the variables c and γ, which are functions of the voltage.

Generalization of the time-dependent rate

One way to derive the general dynamic equation, is to generalize the definition of the time-dependent rate so that it also depends on the ‘global time’ t,

where t is the current time, u is the time of the last transition between states, and t_S is the time of the next switch in channel states. The probability that the channel will switch back to the other state in the time interval (t, t′) can be calculated by dividing (t, t′) into n smaller time segments (t_i, t_i + Δt),i = 1,2,…,n. Breaking up the interval (t, t′) into these segments, we have,

where we have defined The last equality can easily be established by taking the derivative of the final term with respect to t’.

We have thus shown that:

Next we introduce a ‘time-dependent RTPDF’, which can be thought of as an extension of the RTPDF of the type introduced in Section ‘The Model’. This function is the probability density to switch back into the other state at time x:

Now that we have defined the time-dependent RTPDFs in Eq. 42, we can generalize the dynamic equation developed in Eq. 16 (through the projection method). Setting p(0) = 1 for simplicity, we get,

which allows us to calculate the behavior of p(t) for any time-varying input. Notice that the term that appeared in Eq. 16, has been removed because we chose p(0) = 1. Since the input voltage is now general, choosing p(0) = 1 is less of a constraint than before, as we can simply choose t = 0 to be some time at which the channel was available. By doing so we got rid of the initial conditions defined by π_δ(0), which are unobservable.

The case of power-law RTPDF

For a constant rate λ(t|u) = γ we get, as expected, the Markovian exponential RTPDF,

while for a varying rate γ(t) we get,

For the the non-Markovian state with a constant input λ(t|u) = c/(t − u + t₀), we get as expected,

while for general input, if we assume that

we obtain

For example, if c(t) = c₁ for t < t₁ and c(t) = c₂ for t ≥ t_1, then we can write, for all t ≥ t₁,

and so we get the following dynamic equation,

where we denoted

And so, the initial conditions for the evolution of p(t) from time t₁ onward, are completely determined by the last term in Eq. 47. Also, it is important to note that the expression on Eq. 47 is easier to calculate numerically than Eq. 43, since all the integrals can be written in a convolution form.

Direct generalization of the RTPDF

Another way to generalize the channel model to include varying input is to directly generalize the definition of the RTPDF, instead of going through the time-dependent rate. This could be done if we assume that the two-state non-Markov model is based on some underlying many-states Markov model, as explained in Section ‘How to create a two-state non-Markov model from a general Markov model’. In the case of time-varying rates, it is straightforward to generalize Eqs 30 and 31 and get

where we defined and is the stationary distribution in the available set of states at time t : when all the transitions to inactivation are forbidden, and we assume that the input changes slower than time it takes to reach the stationary distribution over the available states. We can substitute Eq. 49 into Eq. 43 to get the dynamic equation for p(t), and derive an expression for the generalized time-dependent rate (similarly to the derivation in Eqs 32 and 33), yielding

where the rightmost expression allows to calculate λ_I(t|u) in a causal way, without knowing future inputs. Notice that here λ_I(t|u) may generally depend on the entire history of the input since time u.

For example, suppose we use this method on the same state space described in the case of projection method (see ‘Projection’ of the Inactivated State) and assume that only the rates of inactivation change with voltage. Taking the continuum limit of Eq. 49, it is straightforward to get the following generalization of Eq. 12, in the power-law case:

where is obvious generalization of Eq. 18. Using Eq. 50 we get,

Notice that Eq. 52 differs from Eq. 45, and also that Eq. 51 differs from Eq. 46. Recall that Eq. 45. was an assumption, which we used to derive the RTPDF in 46. In contrast, Eq. 52 was derived from the RTPDF at Eq. 51, which was itself derived from the assumed Markovian model of the projection method. The different assumptions led to different results. Specifically, this shows that when we take into account time-varying inputs, the Markov model used in the projection method is no longer equivalent to all other models. In this work we chose to model the channel behavior through Eq. 45, which is easier to implement numerically. Using Eq. 52 instead should not significantly change the results for the various inputs examined in this work (constant input, step input, oscillating input and slowly changing input) though it should have an effect for other types of inputs. The final choice between Eqs 45 and 52 or some other, more complicated expression, should be made through carefully designed experiments on channels.

Oscillating Voltage Input

We examine the case where the input to the model oscillates with period T_P. These oscillations may arise, for example, from periodic stimulation or some noise process (in which case T_P is only approximate). For example, in the context of the specific channel model with power-law RTPDF the parameters γ, c may generally be voltage-dependent and vary with time, while t₀ is maintained constant. Intuitively, if T_P is sufficiently small in comparison with the transition timescales, then on sufficiently large timescales we should not notice these fluctuations in the parameters (even though they may be large in magnitude), and expect them to ‘average-out’. Here we prove this claim rigorously, and find the exact conditions for this to happen.

For any time-dependent variable X_S define the mean over a period T_P,

Again, t denotes the current time, u is the time of the last transition between states, and t_s is the time of the next switch in channel states. Using Eq. 41, we can write, for all θ > 0,

where, in order to obtain the final line, we replaced λ(z|u) by in the penultimate line. From the above derivation, if we can choose θ so that and we can safely approximate,

Then, for every timescale above θ we can replace λ(t|u) by

The case of power-law RTPDF

In the context of our model, this means that if, for all t, T_P ≪ γ⁻¹(t), we can approximate,

where t_i is the time of inactivation, and the index t was suppressed in since it is assumed to be, at least approximately, constant.

Also, if T_P ≪ min_s[t₀/c(s)] and T_P ≪ t₀ then we can approximate, using Eq. 45,

where again, t_R is the time of recovery, and the t index was suppressed in

In conclusion, for rapidly oscillating input, we can replace γ, c by their time-averaged expressions .

The Joint Probability Distribution of the Availability

So far we have dealt only with the marginal probability to be in the available state – p(t). To complete the description of the process, we discuss here joint probability distributions.

The available state is Markovian, which makes it is easy to derive p(t₁,t₂,…,t_k), the joint probability distribution to be in the available state at some arbitrary times t₁,…,t_k,

where we denoted p(t_i|t_j) as the probability of the channel to be in the available state at time t_i assuming that at t_j the channel was also in the available state. For arbitrary voltage input, each p(t_i|t_j) term in this product is the solution of Eq. 43, with the initial condition p(t_j) = 1 (notice also that the lower limit of integration in Eq. 43 must be set to t_j). In the case that the voltage is constant, p(t_i|t_j) = p(t_i − t_j), where p(t) is the solution of Eq. 16 with the initial condition p(0) = 1. If we assume also, for simplicity, that p(0) = 1, we get,

where p(t) is the solution of Eq. 16.

Next, we compute the joint moments and the auto-covariance function. Define S(t) to be the channel state at time t and introduce the indicator function I(t) to equal 1 if, and only if, S(t) = A, and zero otherwise, then,

where the angular brackets indicate an average. Thus the joint distributions are equal to the joint moments of the available state.

Using the above results, we can easily calculate k(t₁, t₂), the availability auto-covariance function, where we assume that t₂ > t₁,

The case of power-law RTPDF

For example, in the context of our channel model, and for constant voltage input, we use Eq. 4 and asymptotically obtain, for c > 1,

where k(t₂ − t₁) = k(t₁, t₂) is the stationary auto-covariance.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors are grateful to Erez Braun, Naama Brenner, Asaf Gal and Avner Wallach for many insightful discussions. Special thanks to Yuval Elhanati, Yariv Kafri and Shimon Marom for the considerable effort they invested in improving the quality and focus of this manuscript. The work of R. Meir is partially supported by grant No. 665/08 from the Israel Science Foundation.

Supplementary Material

The Supplementary Material for this article can be found online at http://www.frontiersin.org/computational neuroscience/paper/10.3389/fncom.2010.00003

References

Baylor, D. A., Hodgkin, A. L., and Lamb, T. D. (1974). The electrical response of turtle cones to flashes and steps of light. J. Physiol. 242, 685–727.

Beran, J. (1994). Statistics for Long-Memory Processes. Chapman and Hall/CRC: New York.

Bouchaud, J. P. (1992). Weak ergodicity breaking and aging in disordered systems. J. Phys. I France 2, 1705–1713.

Cinlar, E. (1975). Introduction to Stochastic Processes. Prentice Hall: New Jersey.

Colquhoun, D., and Hawkes, A. G. (1977). Relaxation and fluctuations of membrane currents that flow through drug-operated channels. Proc. R. Soc. Lond., B, Biol. Sci. 199, 231–262.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cox, D. R. (1962). Renewal Theory. London, Methuen.

Ellerkmann, R. K., Riazanski, V., Elger, C. E., Urban, B. W., and Beck, H. (2001). Slow recovery from inactivation regulates the availability of voltage-dependent Na+ channels in hippocampal granule cells, hilar neurons and basket cells. J. Physiol. (Lond.) 532, 385–397.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Erdeyi, A. ed. (1953). Higher Transcendental Functions. New York: McGraw-Hill.

Faber, G. M., Silva, J., Livshitz, L., and Rudy, Y. (2007). Kinetic properties of the cardiac L-type Ca2+ channel and its role in myocyte electrophysiology: a theoretical investigation. Biophys. J. 92, 1522–1543.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fairhall, A. L., Lewen, G. D., Bialek, W., and de Ruyter van Steveninck, R. R. (2001). Efficiency and ambiguity in an adaptive neural code. Nature 412, 787–792.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Feller, W. (1971). An Introduction to Probability Theory and Its Applications, 2nd Edn, Vol. 2. Wiley: New York.

Fleidervish, I. A., Friedman, A., and Gutnick, M. J. (1996). Slow inactivation of Na current and slow cumulative spike adaptation in mouse and guinea-pig neocortical neurones in slices. J. Physiol. (Lond.) 493, 83–97.

Pubmed Abstract | Pubmed Full Text

Friedlander, T., and Brenner, N. (2009). Adaptive response by state-dependent inactivation. Proc. Natl. Acad Sci. U.S.A. 106, 22558.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gerstner, W., and Kistler, W. M. (2002). Spiking Neuron Models. Cambridge University Press: Cambridge.

Gilboa, G., Chen, R., and Brenner, N. (2005). History-dependent multiple-time-scale dynamics in a single-neuron model. J. Neurosci. 25, 6479–6489.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Goychuk, I., and Hanggi, P. (2004). Fractional diffusion modeling of ion channel gating. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 70, 051915.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hering, J., Feltz, A., and Lambert, R. C. (2004). Slow inactivation of the CaV3. 1 isotype of T-type calcium channels. J. Physiol. 555, 331–344.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liebovitch, L. S. (1989). Analysis of fractal ion channel gating kinetics: kinetic rates, energy levels, and activation energies. Math. Biosci. 93, 97–115.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liebovitch, L. S., and Sullivan, J. M. (1987). Fractal analysis of a voltage-dependent potassium channel from cultured mouse hippocampal neurons. Biophys. J. 52, 979–988.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lowen, S. B., Liebovitch, L. S., and White, J. A. (1999). Fractal ion-channel behavior generates fractal firing patterns in neuronal models. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 59, 5970–5980.

CrossRef Full Text

Lowen, S. B., and Teich, M. C. (1995). Estimation and simulation of fractal stochastic point processes. Fractals 3, 183–210.

CrossRef Full Text

Lundstrom, B. N., Higgs, M. H., Spain, W. J., and Fairhall, A. L. (2008). Fractional differentiation by neocortical pyramidal neurons. Nat. Neurosci. 11, 1335–1342.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Marom, S. (1998). Slow changes in the availability of voltage-gated ion channels: effects on the dynamics of excitable membranes. J. Membr. Biol. 161, 105–113.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Marom, S. (2010). Neural timescales or lack thereof. Prog. Neurobiol. 90, 16–28.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Marom, S. (2009). Adaptive transition rates in excitable membranes. Front. Comput. Neurosci. 3:2. doi: 10.3389/neuro.10.002.2009.

CrossRef Full Text

Melamed-Frank, M., and Marom, S. (1999). A global defect in scaling relationship between electrical activity and availability of muscle sodium channels in hyperkalemic periodic paralysis. Pflugers Arch. 438, 213–217.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mercik, S., and Weron, K. (2001). Stochastic origins of the long-range correlations of ionic current fluctuations in membrane channels. Phy. Rev. E Stat. Nonlin. Soft Matter Phys. 63, 51910.

Millhauser, G. L., Salpeter, E. E., and Oswald, R. E. (1988). Diffusion models of ion-channel gating and the origin of power-law distributions from single-channel recording. Proc. Natl. Acad. Sci. U.S.A. 85, 1503–1507.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Powers, R. K., Sawczuk, A., Musick, J. R., and Binder, M. D. (1999). Multiple mechanisms of spike-frequency adaptation in motoneurones. J. Physiol. (Paris) 93, 101–114.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Soen, Y., and Braun, E. (2000). Scale-invariant fluctuations at different levels of organization in developing heart cell networks. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 61, 2216–2219.

CrossRef Full Text

Toib, A., Lyakhov, V., and Marom, S. (1998). Interaction between duration of activity and time course of recovery from slow inactivation in mammalian brain Na+ channels. J. Neurosci. 18, 1893–1903.

Pubmed Abstract | Pubmed Full Text

Uebachs, M., Schaub, C., Perez-Reyes, E., and Beck, H. (2006). T-type Ca2+ channels encode prior neuronal activity as modulated recovery rates. J. Physiol. (Lond.) 571, 519–536.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wark, B., Fairhall, A., and Rieke, F. (2009). Timescales of inference in visual adaptation. Neuron 61, 750–761.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text