Home | About | Journals | Submit | Contact Us | Français |

**|**Front Comput Neurosci**|**v.4; 2010**|**PMC2916672

Formats

Article sections

- Abstract
- Introduction
- Results
- Discussion
- Methods
- Conflict of Interest Statement
- Supplementary Material
- References

Authors

Related links

Front Comput Neurosci. 2010; 4: 3.

PMCID: PMC2916672

Edited by: David Hansel, University of Paris, France

Reviewed by: Nicolas Brunel, Centre National de la Recherche Scientifique, France; Carl van Vreeswijk, Centre National de la Recherche Scientifique, France

*Correspondence: Ron Meir, Department of Electrical Engineering, Technion, Haifa 32000, Israel.e-mail: li.ca.noinhcet.ee@riemr

Received 2009 December 14; Accepted 2010 March 2.

Copyright © 2010 Soudry and Meir.

This is an open-access article subject to an exclusive license agreement between the authors and the Frontiers Research Foundation, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are credited.

This article has been cited by other articles in PMC.

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.

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 (−90mV), then at high voltage (−10mV), for varying lengths of time – from 10ms to 5min. 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 (−90mV). At this low voltage level the channels recovered from inactivation. For short stimulations (*t*_{stim}<1s), the recovery was exponential with a single non-history-dependent timescale. After sufficiently long stimulations (*t*_{stim}>1s) 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 and Hodgkin 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 Figure3.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.

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*_{max}*A*(*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).

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 Figure11.

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

$${\psi}_{E}(t)=\gamma \text{exp}(-\gamma t),\text{\hspace{1em}}t\ge 0.$$

(1)

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

$${\psi}_{P}(t)=\frac{c/{t}_{0}}{{(1+t/{t}_{0})}^{c+1}},\text{\hspace{1em}}t\ge 0,$$

(2)

where γ and *c* are, in general, voltage-dependent, while *t*_{0} 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*_{0} 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.

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*_{0} 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}*(

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.

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,

$$\frac{d}{dt}p(t)=-\gamma p(t)+{\displaystyle \underset{0}{\overset{t}{\int}}\gamma p\left({t}^{\prime}\right){\psi}_{P}\left(t-{t}^{\prime}\right)d{t}^{\prime}}.$$

(3)

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,

$$p(t)={p}_{\infty}+\left(1-{p}_{\infty}\right)q{t}^{-|1-c|},$$

(4)

where $q=\text{sin}(\pi c)/(\pi \gamma {t}_{0}^{c})$ for 0<*c*<1, $q={p}_{\infty}{t}_{0}^{c-1}$ for *c*>1, and *p*_{∞} is the steady state value:

$${p}_{\infty}=\{\begin{array}{ll}{\scriptscriptstyle \frac{c-1}{\gamma {t}_{0}+c-1}}\hfill & \text{if\hspace{1em}}c\ge 1,\hfill \\ 0\hfill & \text{if\hspace{1em}}0<c<1.\hfill \end{array}$$

(5)

We can express *p*_{∞} in the more familiar form of the steady state of a simple two-state Markov process, ${p}_{\infty}=\tilde{\delta}/(\tilde{\delta}+\gamma )$ where $\tilde{\delta}$ is, as in the Markov scheme, the inverse of the mean inactivation time:

$$\tilde{\delta}={\left({\displaystyle \underset{0}{\overset{\infty}{\int}}t{\psi}_{P}(t)\text{\hspace{0.17em}}dt}\right)}^{-1}=\{\begin{array}{ll}0\hfill & \text{if\hspace{1em}}c\le 1,\hfill \\ {\scriptscriptstyle \frac{c-1}{{t}_{0}}}\hfill & \text{if\hspace{1em}}c>1.\hfill \end{array}$$

(6)

In Figure Figure22 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.

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.

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*,

$$\tau =\frac{T+{t}_{0}}{c}.$$

(7)

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 $C{V}_{T}\triangleq {\sigma}_{T}/\langle T\rangle ,$ 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

In Section ‘Renewal Theory Approach’ we compute exactly the behavior of *T* and *CV _{T}* at the end of the voltage pulse, at time

*c*<1:*T*=(1−*c*)*t*, $C{V}_{T}=\sqrt{c/(2(1-c)).}$ The recovery timescale increase linearly with*t*, while maintaining a constant dispersion. Also, Since*CV*increases monotonically with_{T}*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*→∞. Thus, the distribution of recovery timescales broadens constantly, which entails a non-exponential recovery (in this case, the mean,_{T}*T*, provides a poor representation of the distribution). *c*>3:*T*=*t*_{0}/(*c*−2), $C{V}_{T}=\sqrt{(c-1)/(c-3)}.$ Since both*T*and*CV*converge to a fixed finite value, we get that the recovery following the pulse takes place at a single timescale. Higher values of_{T}*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

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 γ^{−1}, 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

$$\widehat{\gamma}=\frac{1}{{T}_{P}}{\displaystyle \underset{0}{\overset{{T}_{P}}{\int}}\gamma (V(t))\text{\hspace{0.17em}}dt}\text{\hspace{0.17em}};\text{\hspace{0.17em}\hspace{1em}}\widehat{c}=\frac{1}{{T}_{P}}{\displaystyle \underset{0}{\overset{{T}_{P}}{\int}}c(V(t))\text{\hspace{0.17em}}dt},$$

(8)

where *V _{t}* is the voltage at time

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.

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*.

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 −10mV for varying length of time, *t*_{stim}, and then clamped at low voltage of −90mV – in which the recovery from inactivation occurred. This recovery was exponential and history dependent for every stimulus longer than 1s. Analyzing these experimental observations using our model, we are able to restrict the possible values of the parameter *c*. We denote by *c _{H}*, γ

Claim: During inactivation *c*=*c _{H}*<1, while during recovery

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

Next, we argue that *c _{H}*<1. First, it is clear that

We thus conclude that in order to fit the experimental results, we must set *c _{H}*<1 during the inactivation period and

*c*is small enough so that_{H}*CV*1, since_{T}*CV*$C{V}_{T}=\sqrt{{c}_{H}/(2(1-{c}_{H}))},$ which is increasing in_{T}*c*._{H}*c*is large enough in the low voltage phase so that_{L}*T*does not change considerably during the recovery. More formally,*T*τ=(*T*+*t*_{0})/*c*, or more simply,_{L}*c*is sufficiently larger than 1._{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, ${\gamma}_{L}^{-1}\gg \tau =(T+{t}_{0})/{c}_{L}.$_{L}

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:

$$\langle \tau \rangle =\frac{\left(1-{c}_{H}\right)\cdot {t}_{\text{stim}}+{t}_{0}}{{c}_{L}}$$

(9)

We note here that it is possible to modify the model so that *t*_{0} 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, −120mV. The history-dependent behavior remained, as can be seen in Figure 6 in Toib et al. (1998). The recovery was similar for −90, −120mV, indicating that *c* did not change much between these voltages. At −60mV, 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}≥1s. 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*_{0}~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.

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

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:

$$\widehat{\gamma}=\frac{{T}_{H}{\gamma}_{H}+{T}_{L}{\gamma}_{L}}{{T}_{H}+{T}_{L}};\text{\hspace{1em}}\widehat{c}=\frac{{T}_{H}{c}_{H}+{T}_{L}{c}_{L}}{{T}_{H}+{T}_{L}}.$$

(10)

It is important to notice in this case that if *T _{H}*

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, $\widehat{c}$ must be significantly different from *c _{L}*. For example, if a neuron is affected solely by its own action potentials, then

It is important to note that in the case $\widehat{c}>3,$ the non-Markovian state can be approximated by a simple Markovian state, since in this case *T*<*t*_{0}, and thus $\tau =\left(T+{t}_{0}\right)/\widehat{c}$ 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*_{0}). A slow-down of the timescale of recovery after a long −55mV depolarizations was observed, similar to what was measured in sodium channels (referred to as ‘continuous’ in Figure8 in Uebachs et al., 2006). In contrast, after a spiking stimulus at a maximal voltage of 25mV with *T _{H}*/(

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:

$$k(t)=\left[{p}_{\infty}^{2}\left(1-{p}_{\infty}\right){t}_{0}^{c-1}\right]{t}^{-c+1}.$$

(11)

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 ${\int}_{0}^{\infty}k(t)\text{\hspace{0.17em}}dt=\infty .$ 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 <1s), as also pointed out by Goychuk and Hanggi (2004).

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.

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 Table11).

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*_{0}) 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 Na*II* and Na*IIA* 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*_{0}) 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’.

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}*(

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.

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}*(

$${\psi}_{\delta}(t)=\delta \mathrm{exp}(-\delta t),\text{\hspace{1em}}t\ge 0\text{\hspace{0.17em}\hspace{0.17em}}(0\le \delta <\infty ).$$

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 ${\int}_{0}^{\infty}\gamma f(\delta )\text{\hspace{0.17em}}d\delta =\gamma ,$ 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}*(

$${\psi}_{I}(t)={\displaystyle \underset{0}{\overset{\infty}{\int}}f(\delta ){\psi}_{\delta}(t)\text{\hspace{0.17em}}d\delta}.$$

(12)

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

We denote $\overline{\pi}=(p,{\pi}_{1},\dots ,{\pi}_{N}),$ where *p* is probability of being in the available state, and π* _{k}* is the probability of being in the inactivated state

$$\frac{d}{dt}\overline{\pi}(t)=\overline{\pi}(t)\hspace{0.17em}\mathbf{\text{Q}},$$

where

$$\mathbf{\text{Q}}=\left(\begin{array}{ccccc}-{\displaystyle {\sum}_{k=1}^{N}{\gamma}_{k}}& {\gamma}_{1}& {\gamma}_{2}& \cdots & {\gamma}_{N}\\ {\delta}_{1}& -{\delta}_{1}& 0& 0& 0\\ {\delta}_{2}& 0& -{\delta}_{2}& 0& 0\\ \vdots & 0& 0& \ddots & 0\\ {\delta}_{N}& 0& 0& 0& -{\delta}_{N}\end{array}\right)\text{\hspace{0.17em}},$$

is the rate matrix.

Writing this explicitly, we get

$$\frac{d}{dt}{\pi}_{k}\left(t\right)={\gamma}_{k}p\left(t\right)-{\delta}_{k}{\pi}_{k}\left(t\right),k=1,2,\dots ,N$$

$$\frac{d}{dt}p(t)=-{\displaystyle \sum _{k=1}^{N}{\gamma}_{k}p(t)}+{\displaystyle \sum _{k=1}^{N}{\pi}_{k}(t){\delta}_{k}}$$

In the exact continuum limit we substitute ${\gamma}_{k}\to \gamma f(\delta )d\delta ,\text{\hspace{0.17em}\hspace{0.17em}}{\pi}_{k}(t)\to {\pi}_{\delta}(t)d\delta ,\text{\hspace{0.17em}\hspace{0.17em}}{\sum}_{k=1}^{N}\to {\int}_{0}^{\infty},$ yielding:

$$\frac{d}{dt}{\pi}_{\delta}(t)\text{\hspace{0.17em}}d\delta =p(t)\text{\hspace{0.17em}}\gamma f(\delta )\text{\hspace{0.17em}}d\delta -\delta {\pi}_{\delta}(t)\text{\hspace{0.17em}}d\delta ,$$

$$\frac{d}{dt}p(t)=-\gamma p(t)+{\displaystyle \underset{0}{\overset{\infty}{\int}}{\pi}_{\delta}(t)\text{\hspace{0.17em}}\delta \text{\hspace{0.17em}}d\delta}.$$

(13)

Dividing the first equation by *d*δ we obtain:

$$\frac{d}{dt}{\pi}_{\delta}(t)=\gamma f(\delta )p(t)-\delta {\pi}_{\delta}(t)$$

(14)

with solution

$${\pi}_{\delta}(t)={\pi}_{\delta}(0){e}^{-\delta t}+\gamma f(\delta ){\displaystyle \underset{0}{\overset{t}{\int}}p\left({t}^{\prime}\right){e}^{-\delta \left(t-{t}^{\prime}\right)}d{t}^{\prime}}$$

(15)

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

$$\frac{d}{dt}p(t)=-\gamma p(t)+{\displaystyle \underset{0}{\overset{\infty}{\int}}{\pi}_{\delta}(0){\text{e}}^{-\delta t}d\delta}+\gamma {\displaystyle \underset{0}{\overset{t}{\int}}d{t}^{\prime}p\left({t}^{\prime}\right)}{\displaystyle \underset{0}{\overset{\infty}{\int}}f(\delta ){\text{e}}^{-\delta \left(t-{t}^{\prime}\right)}\delta \text{\hspace{0.17em}}d\delta}$$

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

$$\frac{d}{dt}p(t)=-\gamma p(t)+\gamma {\displaystyle \underset{0}{\overset{t}{\int}}p\left({t}^{\prime}\right){\psi}_{I}\left(t-{t}^{\prime}\right)d{t}^{\prime}}+{\displaystyle \underset{0}{\overset{\infty}{\int}}{\pi}_{\delta}(0){\psi}_{\delta}(t)\text{\hspace{0.17em}}d\delta}$$

(16)

This result is quite intuitive:

- The term −γ
*p*(*t*) represents the probability current out of the available state. - The term $\gamma {\int}_{0}^{t}p({t}^{\prime}){\psi}_{I}(t-{t}^{\prime})\text{\hspace{0.17em}}d{t}^{\prime}$ 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*. - The term ${\int}_{0}^{\infty}{\pi}_{\delta}(0){\psi}_{\delta}(t)\text{\hspace{0.17em}}d\delta $ 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:

$$s\tilde{p}(s)-p(0)=-\gamma \tilde{p}(s)+\gamma \tilde{p}(s){\tilde{\psi}}_{I}(s)+{\displaystyle \underset{0}{\overset{\infty}{\int}}\frac{\delta {\pi}_{\delta}(0)}{s+\delta}d\delta ,}$$

$$\tilde{p}(s)=\frac{p(0)+{\displaystyle {\int}_{0}^{\infty}{\scriptscriptstyle \frac{\delta {\pi}_{\delta}\left(0\right)}{s+\delta}}}\text{\hspace{0.17em}}d\delta}{s+\gamma \left(1-{\tilde{\psi}}_{I}(s)\right)},$$

(17)

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.

By using an inverse Laplace transform on the *t* variable in Eq. 12, we can easily show that in the case where ψ* _{I}*(

$$f(\delta )=\frac{{t}_{0}^{c}}{\Gamma (c)}{\delta}^{c-1}\mathrm{exp}\left(-\delta {t}_{0}\right),\text{\hspace{1em}}\delta >0,$$

(18)

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

$$\underset{0}{\overset{\infty}{\int}}f(\delta ){\psi}_{\delta}(t)\text{\hspace{0.17em}}d\delta}=\frac{{t}_{0}^{c}}{\Gamma (c)}{\displaystyle \underset{0}{\overset{\infty}{\int}}{\delta}^{c}\mathrm{exp}\left(-\delta \left(t+{t}_{0}\right)\right)d\delta$$

which equals ψ* _{P}*(

In this section we demonstrate how to derive an asymptotic solution of Eq. 17 in the case where ψ* _{I}*(

$${\tilde{\psi}}_{P}(s)=c{e}^{s{t}_{0}}\left[\Gamma (-c)\cdot {\left({t}_{0}s\right)}^{c}-{\displaystyle \sum _{k=0}^{\infty}\frac{{\left(-s{t}_{0}\right)}^{k}}{k!\left(k-c\right)}}\right]$$

(19)

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

$${\tilde{\psi}}_{P}(s)=1+\left(\frac{{t}_{0}}{1-c}\right)s+\frac{{t}_{0}^{2}}{(c-2)(c-1)}{s}^{2}+c\Gamma (-c){t}_{0}^{c}\cdot {s}^{c}+O\left({s}^{\mathrm{min}\left(3,c+1\right)}\right)$$

(20)

$$\tilde{p}(s)=\frac{1}{s-{\scriptscriptstyle \frac{\gamma {t}_{0}}{1-c}}s-c\gamma \Gamma (-c){t}_{0}^{c}{s}^{c}+O\left({s}^{\mathrm{min}\left(2,c+1\right)}\right)}.$$

If *c*<1,

$$\tilde{p}(s)=\frac{-1}{c\gamma \Gamma (-c){t}_{0}^{c}}{s}^{-c}+O\left({s}^{1-2c}\right).$$

(21)

while if *c*>1,

$$\tilde{p}(s)={p}_{\infty}\cdot {s}^{-1}+c\gamma \Gamma (-c){p}_{\infty}^{2}{t}_{0}^{c}{s}^{c-2}+O(1)$$

(22)

where we have introduced

$${p}_{\infty}=\{\begin{array}{ll}{\scriptscriptstyle \frac{c-1}{\gamma {t}_{0}+c-1}}\hfill & \text{if\hspace{1em}}c\ge 1,\hfill \\ 0\hfill & \text{if\hspace{1em}}0<c<1,\hfill \end{array}$$

(23)

denoting the steady state value of *p*(*t*) – a result that can be confirmed by using the final value theorem on $\tilde{p}(s).$ Also, Notice that in Eq. 22 we deliberately retained the *s ^{c−2}* term, even though

We wish to find the asymptotic behavior of $p(t)=\mathcal{L}{\text{\hspace{0.05em}}}^{-1}[\tilde{p}(s)]$ as $t\to \infty $. To do this we use a Tauberian theorem (Theorem 4 on page 446 in Feller, 1971), stating that *f*(*t*)~[*K*/Γ(α)]*t ^{α−1}* as if, and only if, [

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

$$p(t)=\frac{-1}{c\gamma \Gamma (-c)\Gamma (c){t}_{0}^{c}}{t}^{c-1}=\frac{\mathrm{sin}(\pi c)}{\gamma \pi {t}_{0}^{c}}{t}^{c-1},$$

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*^{−1}, and then use the Tauberian theorem. From the linearity of the Laplace transform we get:

$$\begin{array}{c}p(t)-{p}^{\ast}=c\gamma \frac{\Gamma (-c)}{\Gamma (2-c)}{p}_{\infty}^{2}{t}_{0}^{c}{t}^{1-c},\\ =\left(1-{p}_{\infty}\right){p}_{\infty}{t}_{0}^{c-1}{t}^{1-c}\end{array}$$

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*^{−1} 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:

$$\frac{{d}^{m}}{d{s}^{m}}\left(\tilde{p}(s)-{p}_{\infty}{s}^{-1}\right)=c\gamma \Gamma (-c){p}_{\infty}^{2}{t}_{0}^{c}\frac{\Gamma (c-1)}{\Gamma (c-1-m)}{s}^{c-2-m}+O(1)$$

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

$$L\left[{t}^{m}f(t)\right]={(-1)}^{m}\frac{{d}^{m}}{d{s}^{m}}L[f(t)],$$

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

$$\begin{array}{c}{\left(-1\right)}^{m}{t}^{m}\cdot \left(p(t)-{p}_{\infty}\right)=c\gamma {p}_{\infty}^{2}\frac{\Gamma (-c)\Gamma (c-1)}{\Gamma (c-1-m)\text{\hspace{0.17em}}\Gamma (m+2-c)}{t}_{0}^{c}{t}^{1-c+m},\\ p(t)-{p}_{\infty}=\left(1-{p}_{\infty}\right){p}_{\infty}\frac{{(-1)}^{m}\mathrm{sin}\left((c-2-m)\pi \right)}{\mathrm{sin}(c\pi )}{t}_{0}^{c-1}{t}^{1-c}\\ =\left(1-{p}_{\infty}\right){p}_{\infty}{t}_{0}^{c-1}{t}^{1-c}\end{array}$$

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*,

$$p(t)=\{\begin{array}{ll}{p}_{\infty}+\left(1-{p}_{\infty}\right){p}_{\infty}{t}_{0}^{c-1}{t}^{1-c}\hfill & \text{if\hspace{1em}}1<c\hfill \\ {\scriptscriptstyle \frac{\mathrm{sin}(\pi c)}{\pi \gamma {t}_{0}^{c}}}{t}^{c-1}\hfill & \text{if\hspace{1em}}0<c<1\hfill \end{array}$$

(24)

We derived this result for the case *p*(0)=1, but it remains valid also when *p*(0)≠1 if the expression ${\int}_{0}^{\infty}[\delta {\pi}_{\delta}(0)/(s+\delta )]\text{\hspace{0.17em}}d\delta $ is analytic at the origin. For example, this always is true if for all δ<δ_{0}, π_{δ}(0)=0 for some δ_{0}>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

$$\begin{array}{c}{\left(t-{t}_{i}\right)}^{a}={t}^{a}{\left(1-\frac{{t}_{i}}{t}\right)}^{a}={t}^{a}\left(1-\frac{a{t}_{i}}{t}+O\left({t}^{-2}\right)\right)\\ ={t}^{a}+O\left({t}^{a-1}\right),\end{array}$$

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,

$$p(t)={p}_{\infty}+\left(p(0)-{p}_{\infty}\right){\text{e}}^{-at}.$$

In this case when we set the initial condition to *t _{i}*, instead of 0, we get,

$$\begin{array}{c}p(t)={p}_{\infty}+\left(p\left({t}_{i}\right)-{p}_{\infty}\right){e}^{-a\left(t-{t}_{i}\right)}\\ ={p}_{\infty}+{\text{e}}^{a{t}_{i}}\left(p\left(0\right)-{p}_{\infty}\right){\text{e}}^{-at}.\end{array}$$

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

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 $d\overline{\pi}(t)/dt=\overline{\pi}(t)\hspace{0.17em}\mathbf{\text{Q}},$ where $\overline{\pi}(t)$ 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 $\overline{\pi}(t)=({\overline{\pi}}_{A}(t),\text{\hspace{0.17em}}{\overline{\pi}}_{I}(t)).$ Similarly, we divide the rate matrix **Q** into block matrices corresponding to these sets of states:

$$\mathbf{\text{Q\hspace{0.17em}}}=\left(\begin{array}{cc}{\mathbf{\text{Q}}}_{AA}& {\mathbf{\text{Q}}}_{AI}\\ {\mathbf{\text{Q}}}_{IA}& {\mathbf{\text{Q}}}_{II}\end{array}\right).$$

In order to calculate the RTPDFs of the inactivated states, we introduce a modified Markov process with probability vector $\tilde{\pi}(t)=({\tilde{\pi}}_{A}(t),\text{\hspace{0.17em}}{\tilde{\pi}}_{I}(t)),$ and the rate matrix:

$$\tilde{\mathbf{\text{Q}}}\text{\hspace{0.17em}}=\left(\begin{array}{cc}0& 0\\ {\mathbf{\text{Q}}}_{IA}& {\mathbf{\text{Q}}}_{II}\end{array}\right).$$

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 $d\tilde{\pi}(t)/dt=\tilde{\pi}(t)\tilde{\mathbf{\text{Q}}}$ we get,

$$\frac{d}{dt}{\tilde{\pi}}_{I}(t)={\tilde{\pi}}_{I}(t){\mathbf{\text{Q}}}_{II}$$

(25)

$$\frac{d}{dt}{\tilde{\pi}}_{A}(t)={\tilde{\pi}}_{I}(t){\mathbf{\text{Q}}}_{IA}$$

(26)

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

$${\tilde{\pi}}_{I}(t)={\tilde{\pi}}_{I}(0)\mathrm{exp}\left({\mathbf{\text{Q}}}_{II}t\right),\text{\hspace{0.17em}}t\ge 0,$$

(27)

where the exponent of a matrix **A** is defined as $\text{exp}(\mathbf{\text{A}})={\sum}_{k=0}^{\infty}{\mathbf{\text{A}}}^{k}/k!.$

Now assume that $\overline{\pi}(0)=\tilde{\pi}(0)=(0,\text{\hspace{0.17em}}{\overline{\pi}}_{I}(0)),$ 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 ${\tilde{\pi}}_{A}(t){\overline{u}}_{A}^{T},$ where we defined ${\overline{u}}_{A}^{T}={(1,1,\dots ,1)}^{T},$ a column vector of ones with the same number of components as ${\overline{\pi}}_{A}(t).$ Differentiating this expression, we get the RTPDF of the inactivated state, assuming that the channel has inactivated at *t*=0:

$${\psi}_{I}(t)=\frac{d}{dt}\left({\tilde{\pi}}_{A}(t){\overline{u}}_{A}^{T}\right)={\overline{\pi}}_{I}(0)\mathrm{exp}({\mathbf{\text{Q}}}_{II}t){\mathbf{\text{Q}}}_{IA}{\overline{u}}_{A}^{T},\text{\hspace{1em}}t\ge 0,$$

(28)

where we used Eqs 26 and 27, and ${\overline{\pi}}_{I}(0)$ represents the probability distribution in the inactivated states immediately after the inactivation. Similarly, the RTPDF of the available state is

$${\psi}_{A}(t)={\overline{\pi}}_{A}(0)\mathrm{exp}\left({\mathbf{\text{Q}}}_{AA}t\right){\mathbf{\text{Q}}}_{AI}{\overline{u}}_{I}^{T},\text{\hspace{1em}}t\ge 0.$$

(29)

where ${\overline{u}}_{I}^{T}$ is defined similarly to ${\overline{u}}_{A}^{T}.$

Note that ψ* _{I}*(

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,

$${\psi}_{A}(t)={\psi}_{E}(t)=\gamma \mathrm{exp}(-\gamma t),\text{\hspace{1em}}t\ge 0,$$

(30)

$${\psi}_{I}(t)={\overline{\pi}}_{A}^{s}{G}_{AI}\mathrm{exp}\left({\mathbf{\text{Q}}}_{II}t\right){\mathbf{\text{Q}}}_{IA}{\overline{u}}_{A}^{T},\text{\hspace{1em}}t\ge 0,$$

(31)

where we defined $\gamma ={\overline{\pi}}_{A}^{s}{\mathbf{\text{Q}}}_{AI}{\overline{u}}_{I}^{T},$ ${({G}_{AI})}_{mn}={({\mathbf{\text{Q}}}_{AI})}_{mn}/[{\sum}_{l=1}^{L}{({\mathbf{\text{Q}}}_{AI})}_{ml}]$ and ${\overline{\pi}}_{A}^{s}$ 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.

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

$$\begin{array}{c}P\left\{{t}_{R}<T+\Delta t|{t}_{R}>T\right\}=\frac{P\left\{T<{t}_{R}<T+\Delta t\right\}}{P\left({t}_{R}>T\right)},\\ =\frac{{\psi}_{I}(T)\Delta t+O\left(\Delta {t}^{2}\right)}{{\displaystyle {\int}_{T}^{\infty}{\psi}_{I}(z)\text{\hspace{0.17em}}dz}},\end{array}$$

(32)

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

$$\begin{array}{c}\lambda \text{\hspace{0.17em}}\stackrel{\Delta}{=}\text{\hspace{0.17em}}\underset{\Delta t\to 0}{\mathrm{lim}}\frac{1}{\Delta t}P\left({t}_{R}<T+\Delta t|{t}_{R}>T\right)\\ =-\frac{d}{dT}\mathrm{log}\left({\displaystyle \underset{T}{\overset{\infty}{\int}}{\psi}_{I}(z)\text{\hspace{0.17em}}dz}\right),\text{\hspace{1em}}(T\ge 0)\end{array}$$

(33)

and the time-dependent timescale of recovery is its inverse $\tau \text{\hspace{0.17em}}\stackrel{\Delta}{=}\text{\hspace{0.17em}}1/\lambda .$

Next, we compute *f _{t}*(

First, we denote by *k _{n}*(

$${k}_{n}(t)={\psi}_{E}(t)*\underset{\text{n}-1\text{\hspace{0.17em}times}}{\underbrace{\left({\psi}_{I}(t)*{\psi}_{E}(t)\right)*\left({\psi}_{I}(t)*{\psi}_{E}(t)\right)*\dots}}$$

In the Laplace domain, this translates to,

$${\tilde{k}}_{n}(s)\triangleq L\left[{k}_{n}(t)\right]=\left(\frac{\gamma}{\gamma +s}\right){\left({\tilde{\psi}}_{I}(s)\cdot \frac{\gamma}{\gamma +s}\right)}^{n-1}.$$

(34)

Defining the rate of inactivation events:

$$h(t)\triangleq \underset{\Delta t\to 0}{\mathrm{lim}}\frac{1}{\Delta t}P\text{\hspace{1em}}(\text{some\hspace{0.17em}inactivation\hspace{0.17em}event\hspace{0.17em}occured\hspace{0.17em}in\hspace{0.17em}}(t,t+\Delta t)),$$

we have $h(t)={\sum}_{n=1}^{\infty}{k}_{n}(t)$ and in the Laplace domain:

$$\tilde{h}(s)={\displaystyle \underset{n=1}{\overset{\infty}{\int}}{\tilde{k}}_{n}(s)}=\frac{\gamma}{s+\gamma \left(1-{\tilde{\psi}}_{I}(s)\right)}$$

where we used the geometric series summation formula. Notice that $\tilde{h}(s)=\gamma \tilde{p}(s)$ 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}*(

$$\begin{array}{l}{f}_{t}\left(T\right)={\left(1-p\left(t\right)\right)}^{-1}h\left(t-T\right){\displaystyle \underset{T}{\overset{\infty}{\int}}{\text{\psi}}_{I}\left(u\right)du,\text{\hspace{0.17em}\hspace{0.17em}}\left(0\le T\le t\right)}\\ ={\left(1-p\left(t\right)\right)}^{-1}\text{\gamma}\left(t-T\right){\displaystyle \underset{T}{\overset{\infty}{\int}}{\text{\psi}}_{I}\left(u\right)du,\text{\hspace{0.17em}\hspace{0.17em}}\left(0\le T\le t\right)}\end{array}$$

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

$$\lambda =-\frac{d}{dT}\mathrm{log}\left({\displaystyle \underset{T}{\overset{\infty}{\int}}{\psi}_{P}(u)\text{\hspace{0.17em}}du}\right)=\frac{c}{T+{t}_{0}},\text{\hspace{1em}}(T\ge 0)$$

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

$$\tau \stackrel{\Delta}{=}\frac{1}{\lambda}=\frac{T+{t}_{0}}{c},\text{\hspace{1em}}(T\ge 0).$$

(35)

The distribution of *T* in this case ψ* _{I}*(

Because of the multiplication in the last expression, it is generally difficult to write an expression for *f _{t}*(

$${f}_{t}(T)=\{\begin{array}{ll}{\scriptscriptstyle \frac{c-1}{{t}_{0}}}{\left(1+{\scriptscriptstyle \frac{T}{{t}_{0}}}\right)}^{-c}\hfill & \text{if\hspace{1em}}1<c\hfill \\ {\scriptscriptstyle \frac{\mathrm{sin}(\pi c)}{\pi {t}_{0}^{c}}}{(t-T)}^{c-1}{\left(1+{\scriptscriptstyle \frac{T}{{t}_{0}}}\right)}^{-c}\hfill & \text{if\hspace{1em}}0<c<1\hfill \end{array},\text{\hspace{1em}}(0\le T\le t)$$

(36)

Note that for 0<*c*<1, for any fixed value of *T*, the function *f _{t}*(

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

$$0<c<1:\{\begin{array}{l}\langle T\rangle =(1-c)t+O(1)\hfill \\ {\sigma}_{T}=\sqrt{{\scriptscriptstyle \frac{c(1-c)}{2}}}\text{\hspace{0.17em}}t+O(1)\hfill \end{array}$$

(37)

$$1<c<2:\{\begin{array}{l}\langle T\rangle ={\scriptscriptstyle \frac{c-1}{2-c}}{t}_{0}^{c-1}{t}^{2-c}+O(1)\hfill \\ {\sigma}_{T}=\sqrt{{\scriptscriptstyle \frac{c-1}{3-c}}}\text{\hspace{0.17em}}{t}_{0}^{\left(c-1\right)/2}{t}^{\left(3-c\right)/2}+O\left({t}^{2-c}\right)\hfill \end{array}$$

(38)

$$2<c<3:\{\begin{array}{l}\langle T\rangle ={\scriptscriptstyle \frac{{t}_{0}}{c-2}}+O\left({t}^{2-c}\right)\hfill \\ {\sigma}_{T}=\sqrt{{\scriptscriptstyle \frac{1-c}{c-3}}}\text{\hspace{0.17em}}{t}_{0}^{\left(c-1\right)/2}\cdot {t}^{\left(3-c\right)/2}+O(1)\hfill \end{array}$$

(39)

$$3<c:\{\begin{array}{l}\langle T\rangle ={\scriptscriptstyle \frac{{t}_{0}}{c-2}}+O\left({t}^{2-c}\right)\hfill \\ {\sigma}_{T}={\scriptscriptstyle \frac{{t}_{0}}{c-2}}\sqrt{{\scriptscriptstyle \frac{c-1}{c-3}}}+O\left({t}^{\left(3-c\right)/2}\right)\hfill \end{array}$$

(40)

We verified numerically that all of the results for these moments and for *f _{t}*(

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.

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*,

$$\lambda \left(t|u\right)\triangleq \underset{\Delta t\to 0}{\mathrm{lim}}\frac{1}{\Delta t}P\left({t}_{S}\le t+\Delta t|{t}_{S}\ge t\ge u\right),$$

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 (

$$\begin{array}{l}P\left({t}_{S}\le {t}^{\prime}|{t}_{S}\ge t\ge u\right)\\ =\underset{\Delta t\to 0}{\mathrm{lim}}\text{\hspace{0.17em}}{\displaystyle \sum _{i=1}^{n}\left({\displaystyle \prod _{k=1}^{i-1}\left(1-\lambda \left({t}_{k}|u\right)\Delta t+O\left({\left(\lambda \left({t}_{k}|u\right)\Delta t\right)}^{2}\right)\right)}\right)}\text{\hspace{0.17em}}\lambda \left({t}_{i}|u\right)\Delta t\\ =\underset{\Delta t\to 0}{\mathrm{lim}}{\displaystyle \sum _{i=1}^{n}\mathrm{exp}\left({\displaystyle \sum _{k=1}^{i-1}\mathrm{ln}\left(1-\lambda \left({t}_{k}|u\right)\Delta t+O\left({\left(\lambda \left({t}_{k}|u\right)\Delta t\right)}^{2}\right)\right)}\right)\lambda \left({t}_{i}|u\right)\Delta t}\\ =\underset{\Delta t\to 0}{\mathrm{lim}}{\displaystyle \sum _{i=1}^{n}\mathrm{exp}\left(-{\displaystyle \sum _{k=1}^{i-1}\left(\lambda \left({t}_{k}|u\right)\Delta t+O\left({\left(\lambda \left({t}_{k}|u\right)\Delta t\right)}^{2}\right)\right)}\right)\lambda \left({t}_{i}|u\right)\Delta t}\\ ={\displaystyle \underset{t}{\overset{{t}^{\prime}}{\int}}\mathrm{exp}\left(-{\displaystyle \underset{t}{\overset{v}{\int}}\lambda \left(z|u\right)dz}\right)\lambda \left(v|u\right)dv}\\ =1-\mathrm{exp}\left(-{\displaystyle \underset{t}{\overset{{t}^{\prime}}{\int}}\lambda \left(z|u\right)dz}\right),\end{array}$$

where we have defined ${\Pi}_{k=1}^{0}=1,$ ${\Sigma}_{k=1}^{0}=0.$ The last equality can easily be established by taking the derivative of the final term with respect to *t*’.

We have thus shown that:

$$P\left({t}_{S}\le {t}^{\prime}|{t}_{S}\ge t\ge u\right)=1-\mathrm{exp}\left(-{\displaystyle \underset{t}{\overset{{t}^{\prime}}{\int}}\lambda \left(z|u\right)dz}\right)$$

(41)

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*:

$$\begin{array}{c}\psi \left(x|y\right)\triangleq \frac{d}{dx}\left[P\left({t}_{S}\ge x|{t}_{S}\ge t=u=y\right)\right]\\ =\mathrm{exp}\left(-{\displaystyle \underset{y}{\overset{x}{\int}}\lambda \left(z|y\right)dz}\right)\lambda \left(x|y\right).\end{array}$$

(42)

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,

$$\frac{d}{dt}p(t)=-\gamma (t)p(t)+{\displaystyle \underset{0}{\overset{t}{\int}}\gamma (u)p(u){\psi}_{I}\left(t|u\right)du},$$

(43)

which allows us to calculate the behavior of *p*(*t*) for any time-varying input. Notice that the term ${\int}_{0}^{\infty}{\pi}_{\delta}(0){\psi}_{\delta}(t)\text{\hspace{0.17em}}d\delta $ 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.

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

$${\text{\psi}}_{E}\left(t|u\right)=\text{\gamma}\mathrm{exp}\left(-\text{\gamma}\left(t-u\right)\right)\left(t-u\right)$$

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

$${\psi}_{E}\left(t|u\right)=\gamma (t)\mathrm{exp}\left(-{\displaystyle \underset{u}{\overset{t}{\int}}\gamma (z)\text{\hspace{0.17em}}dz}\right)\text{\hspace{0.17em}\hspace{0.17em}}(t>u).$$

(44)

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

$$\begin{array}{c}{\psi}_{P}\left(t|u\right)=\mathrm{exp}\left(-{\displaystyle \underset{u}{\overset{t}{\int}}\frac{c}{z-u+{t}_{0}}\text{\hspace{0.17em}}dz}\right)\frac{c}{t-u+{t}_{0}}\\ =\frac{c/{t}_{0}}{{(1+(t-u)/{t}_{0})}^{c+1}}\text{\hspace{1em}}(t>u),\end{array}$$

while for general input, if we assume that

$$\lambda \left(t|u\right)=\frac{c(t)}{t-u+{t}_{0}},$$

(45)

we obtain

$${\psi}_{P}\left(t|u\right)=\mathrm{exp}\left(-{\displaystyle \underset{u}{\overset{t}{\int}}\frac{c(z)}{z-u+{t}_{0}}dz}\right)\frac{c(t)}{t-u+{t}_{0}}\text{\hspace{0.17em}}(t>u),$$

(46)

For example, if *c*(*t*)=*c*_{1} for *t*<*t*_{1} and *c*(*t*)=*c*_{2} for *t*≥*t*_{1,} then we can write, for all *t*≥*t*_{1},

$$\begin{array}{c}{\psi}_{P}\left(t|u\right)=\mathrm{exp}\left(-{\displaystyle \underset{u}{\overset{t}{\int}}\frac{c(z)}{z-u+{t}_{0}}dz}\right)\frac{c(t)}{t-u+{t}_{0}}\\ =\{\begin{array}{cc}\frac{{c}_{\text{2}}}{{t}_{\text{0}}}\frac{{\left(\text{1+(}t\text{}-\text{}u\text{)/}{t}_{0}\right)}^{{c}_{2}-{c}_{1}}}{{\left(\text{1+(}t\text{}-\text{}u\text{)/}{t}_{0}\right)}^{{c}_{2}+1}},& u{t}_{1}\\ \frac{{c}_{\text{2}}}{{t}_{\text{0}}}\frac{1}{{\left(\text{1+(}t\text{}-\text{}u\text{)/}{t}_{0}\right)}^{{c}_{2}+1}},& u\ge {t}_{1}\end{array},\end{array}$$

and so we get the following dynamic equation,

$$\begin{array}{c}\frac{d}{dt}p(t)=-\gamma (t)p(t)+{\displaystyle \underset{{t}_{1}}{\overset{t}{\int}}\gamma (u)p(u){\psi}_{P}(t-u)\text{\hspace{0.17em}}du}\\ \text{\hspace{0.17em}\hspace{0.17em}}+{\displaystyle \underset{0}{\overset{{t}_{1}}{\int}}\gamma (u)p(u){\left(1+\left({t}_{1}-u\right)/{t}_{0}\right)}^{{c}_{2}-{c}_{1}}{\psi}_{P}(t-u)\text{\hspace{0.17em}}du,}\end{array}$$

(47)

where we denoted ${\psi}_{P}(t)=({c}_{2}/{t}_{0}){(1+t/{t}_{0})}^{-{c}_{2}-1}.$

And so, the initial conditions for the evolution of *p*(*t*) from time *t*_{1} 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.

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

$${\psi}_{A}\left(t|u\right)=\gamma (t)\mathrm{exp}\left(-{\displaystyle \underset{u}{\overset{t}{\int}}\gamma (z)\text{\hspace{0.17em}}dz}\right),\text{\hspace{1em}}t\ge 0,$$

(48)

$${\psi}_{I}\left(t|u\right)={\overline{\pi}}_{A}^{s}(u){G}_{AI}(u)\mathrm{exp}\left({\displaystyle \underset{u}{\overset{t}{\int}}{\mathbf{\text{Q}}}_{II}(z)\text{\hspace{0.17em}}dz}\right){\mathbf{\text{Q}}}_{IA}(t){\overline{u}}_{A}^{T},\text{\hspace{1em}}t\ge 0,$$

(49)

where we defined $\gamma (t)={\overline{\pi}}_{A}^{s}(t){\mathbf{\text{Q}}}_{AI}(t){\overline{u}}_{I}^{T},$ ${({G}_{AI}(t))}_{mn}={({\mathbf{\text{Q}}}_{AI}(t))}_{mn}/({\sum}_{l=1}^{L}{({\mathbf{\text{Q}}}_{AI}(t))}_{ml})$ and ${\overline{\pi}}_{A}^{s}(t)$ 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

$${\lambda}_{I}\left(t|u\right)=\frac{{\psi}_{I}\left(t|u\right)}{{\displaystyle {\int}_{t}^{\infty}{\psi}_{I}\left(z|u\right)\text{\hspace{0.17em}}dz}}=\frac{{\psi}_{I}\left(t|u\right)}{1-{\displaystyle {\int}_{u}^{t}{\psi}_{I}\left(z|u\right)dz}}$$

(50)

where the rightmost expression allows to calculate λ* _{I}*(

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:

$$\begin{array}{c}{\lambda}_{I}\left(t|u\right)=\frac{{\psi}_{I}\left(t|u\right)}{{\displaystyle {\int}_{t}^{\infty}{\psi}_{I}\left(z|u\right)\text{\hspace{0.17em}}dz}}=\frac{{\psi}_{I}\left(t|u\right)}{1-{\displaystyle {\int}_{u}^{t}{\psi}_{I}\left(z|u\right)dz}}\\ =\frac{c(u)/{t}_{0}}{{(\text{1+(}t-u\text{)}/{t}_{0})}^{c(u)+1}}\end{array}$$

(51)

where ${f}_{u}(\delta )=({t}_{0}^{c(u)}/\Gamma (c(u))){\delta}^{c(u)-1}\mathrm{exp}(-\delta t),\text{\hspace{0.17em}}\delta >0,$ is obvious generalization of Eq. 18. Using Eq. 50 we get,

$${\lambda}_{P}\left(t|u\right)=\frac{c(u)}{t-u+{t}_{0}}$$

(52)

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.

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

For any time-dependent variable X* _{S}* define the mean over a period

$${\widehat{X}}_{t}\triangleq \frac{1}{{T}_{P}}{\displaystyle \underset{t}{\overset{t+{T}_{P}}{\int}}{X}_{z}\text{\hspace{0.17em}}dz}.$$

(53)

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,

$$\begin{array}{l}P\left({t}_{S}\le t+\theta |{t}_{S}\ge t\ge u\right)\\ =1-\mathrm{exp}\left(-{\displaystyle \underset{t}{\overset{t+\theta}{\int}}\lambda \left(z|u\right)dz}\right)\\ ={\displaystyle \underset{t}{\overset{t+\theta}{\int}}\lambda (z|u)\text{\hspace{0.17em}}dz}+O\left({\left({\displaystyle \underset{t}{\overset{t+\theta}{\int}}\lambda \left(z|u\right)dz}\right)}^{2}\right)\\ =\widehat{\lambda}(t|u)\theta +O\left({\displaystyle \underset{t}{\overset{t+\theta}{\int}}\left(\lambda \left(z|u\right)-\widehat{\lambda}(t|u)\right)dz}\right)+O\left({\left(\widehat{\lambda}(t|u)\theta \right)}^{2}\right),\end{array}$$

where, in order to obtain the final line, we replaced λ(*z*|*u*) by $\widehat{\lambda}(t|u)+(\lambda (z|u)-\widehat{\lambda}(t|u))$ in the penultimate line. From the above derivation, if we can choose θ so that $\widehat{\lambda}(t|u)\text{\hspace{0.17em}}\theta \ll 1,$ and ${\int}_{t}^{t+\theta}|\lambda (z|u)-\widehat{\lambda}(t|u)|du\ll \widehat{\lambda}(t|u)\text{\hspace{0.17em}}\theta ,$ we can safely approximate,

$$P\left({t}_{S}\le t+\theta |{t}_{S}\ge t\ge u\right)\approx \widehat{\lambda}(t|u)\text{\hspace{0.17em}}\theta .$$

(54)

Then, for every timescale above θ we can replace λ(*t*|*u*) by $\widehat{\lambda}(t|u).$

In the context of our model, this means that if, for all *t*, *T _{P}*γ

$$P\left({t}_{I}\le t+\theta |{t}_{I}\ge t\ge u\right)\approx \widehat{\gamma}\text{\hspace{0.17em}}\theta ,$$

(55)

where *t _{I}* is the time of inactivation, and the index

Also, if *T _{P}*min

$$P\left({t}_{R}\le t+\theta |{t}_{R}\ge t\ge u\right)\approx \frac{\widehat{c}}{t+\theta -u+{t}_{0}}\theta ,$$

(56)

where again, *t _{R}* is the time of recovery, and the

In conclusion, for rapidly oscillating input, we can replace γ, *c* by their time-averaged expressions $\widehat{\gamma},\text{\hspace{0.17em}}\widehat{c}$.

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*_{1},*t*_{2},…,*t _{k}*), the joint probability distribution to be in the available state at some arbitrary times

$$p\left({t}_{1},{t}_{2},\dots ,{t}_{k}\right)=p\left({t}_{1}\right)p\left({t}_{2}|{t}_{1}\right)\cdots p\left({t}_{k}|{t}_{k-1}\right),$$

where we denoted *p*(*t _{i}*|

$$p\left({t}_{1},{t}_{2},\dots ,{t}_{k}\right)=p\left({t}_{1}\right)p\left({t}_{2}-{t}_{1}\right)\cdots p\left({t}_{k}-{t}_{k-1}\right)$$

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 (*t*) to equal 1 if, and only if, *S*(*t*)=*A*, and zero otherwise, then,

$$\langle \mathcal{I}\left({t}_{1}\right)\mathcal{I}\left({t}_{2}\right)\cdots \mathcal{I}\left({t}_{k}\right)\rangle =p\left({t}_{1},{t}_{2},\dots ,{t}_{k}\right),$$

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*_{1},*t*_{2}), the availability auto-covariance function, where we assume that *t*_{2}>*t*_{1},

$$\begin{array}{c}k\left({t}_{1},{t}_{2}\right)\triangleq \langle \mathcal{I}\left({t}_{1}\right)\rangle \langle \mathcal{I}\left({t}_{2}\right)\rangle -\langle \mathcal{I}\left({t}_{1}\right)\rangle \langle \mathcal{I}\left({t}_{2}\right)\rangle ,\\ =p\left({t}_{1}\right)\left(p\left({t}_{2}|{t}_{1}\right)-p\left({t}_{2}\right)\right).\end{array}$$

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

$$k(t)={p}_{\infty}^{2}\left(1-{p}_{\infty}\right){t}_{0}^{c-1}\cdot {t}^{-|1-c|}$$

(57)

where *k*(*t*_{2}−*t*_{1})=*k*(*t*_{1}, *t*_{2}) is the stationary auto-covariance.

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.

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

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 ScienceFoundation.

- Baylor D. A., Hodgkin A. L., Lamb T. D. (1974). The electrical response of turtle cones to flashes and steps of light. J. Physiol. 242, 685–727 [PubMed]
- 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., 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–26210.1098/rspb.1977.0137 [PubMed] [Cross Ref]
- Cox D. R. (1962). Renewal Theory. London, Methuen
- Ellerkmann R. K., Riazanski V., Elger C. E., Urban B. W., 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–39710.1111/j.1469-7793.2001.0385f.x [PubMed] [Cross Ref]
- Erdeyi A. ed. (1953). Higher Transcendental Functions. New York: McGraw-Hill
- Faber G. M., Silva J., Livshitz L., 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–154310.1529/biophysj.106.088807 [PMC free article] [PubMed] [Cross Ref]
- Fairhall A. L., Lewen G. D., Bialek W., de Ruyter van Steveninck R. R. (2001). Efficiency and ambiguity in an adaptive neural code. Nature 412, 787–79210.1038/35090500 [PubMed] [Cross Ref]
- Feller W. (1971). An Introduction to Probability Theory and Its Applications, 2nd Edn, Vol. 2 Wiley: New York
- Fleidervish I. A., Friedman A., 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]
- Friedlander T., Brenner N. (2009). Adaptive response by state-dependent inactivation. Proc. Natl. Acad Sci. U.S.A. 106, 22558.10.1073/pnas.0902146106 [PubMed] [Cross Ref]
- Gerstner W., Kistler W. M. (2002). Spiking Neuron Models. Cambridge University Press: Cambridge
- Gilboa G., Chen R., Brenner N. (2005). History-dependent multiple-time-scale dynamics in a single-neuron model. J. Neurosci. 25, 6479–648910.1523/JNEUROSCI.0763-05.2005 [PubMed] [Cross Ref]
- Goychuk I., Hanggi P. (2004). Fractional diffusion modeling of ion channel gating. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 70, 051915.10.1103/PhysRevE.70.051915 [PubMed] [Cross Ref]
- Hering J., Feltz A., Lambert R. C. (2004). Slow inactivation of the CaV3. 1 isotype of T-type calcium channels. J. Physiol. 555, 331–34410.1113/jphysiol.2003.054361 [PubMed] [Cross Ref]
- Liebovitch L. S. (1989). Analysis of fractal ion channel gating kinetics: kinetic rates, energy levels, and activation energies. Math. Biosci. 93, 97–11510.1016/0025-5564(89)90015-1 [PubMed] [Cross Ref]
- Liebovitch L. S., Sullivan J. M. (1987). Fractal analysis of a voltage-dependent potassium channel from cultured mouse hippocampal neurons. Biophys. J. 52, 979–98810.1016/S0006-3495(87)83290-3 [PubMed] [Cross Ref]
- Lowen S. B., Liebovitch L. S., 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–598010.1103/PhysRevE.59.5970 [PubMed] [Cross Ref]
- Lowen S. B., Teich M. C. (1995). Estimation and simulation of fractal stochastic point processes. Fractals 3, 183–21010.1142/S0218348X95000151 [Cross Ref]
- Lundstrom B. N., Higgs M. H., Spain W. J., Fairhall A. L. (2008). Fractional differentiation by neocortical pyramidal neurons. Nat. Neurosci. 11, 1335–134210.1038/nn.2212 [PMC free article] [PubMed] [Cross Ref]
- 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–11310.1007/s002329900318 [PubMed] [Cross Ref]
- Marom S. (2010). Neural timescales or lack thereof. Prog. Neurobiol. 90, 16–2810.1016/j.pneurobio.2009.10.003 [PubMed] [Cross Ref]
- Marom S. (2009). Adaptive transition rates in excitable membranes. Front. Comput. Neurosci. 3:2.10.3389/neuro.10.002.2009 [PMC free article] [PubMed] [Cross Ref]
- Melamed-Frank M., 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–21710.1007/s004240050900 [PubMed] [Cross Ref]
- Mercik S., 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 [PubMed]
- Millhauser G. L., Salpeter E. E., 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–150710.1073/pnas.85.5.1503 [PubMed] [Cross Ref]
- Powers R. K., Sawczuk A., Musick J. R., Binder M. D. (1999). Multiple mechanisms of spike-frequency adaptation in motoneurones. J. Physiol. (Paris) 93, 101–11410.1016/S0928-4257(99)80141-7 [PubMed] [Cross Ref]
- Soen Y., 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–221910.1103/PhysRevE.61.R2216 [PubMed] [Cross Ref]
- Toib A., Lyakhov V., 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]
- Uebachs M., Schaub C., Perez-Reyes E., Beck H. (2006). T-type Ca2+ channels encode prior neuronal activity as modulated recovery rates. J. Physiol. (Lond.) 571, 519–53610.1113/jphysiol.2005.103614 [PubMed] [Cross Ref]
- Wark B., Fairhall A., Rieke F. (2009). Timescales of inference in visual adaptation. Neuron 61, 750–76110.1016/j.neuron.2009.01.019 [PMC free article] [PubMed] [Cross Ref]

Articles from Frontiers in Computational Neuroscience are provided here courtesy of **Frontiers Media SA**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |