History-Dependent Dynamics in a Generic Model of Ion Channels – An Analytic Study

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


INTRODUCTION
Many recent experiments have demonstrated that the timescale of adaptation of a single neuron in response to periodic stimuli slows down as the period of stimulation increases (Fairhall et al., 2001;Lundstrom et al., 2008;Wark et al., 2009).At a sub-neuronal level, experiments on sodium (Toib et al., 1998;Melamed-Frank and Marom, 1999;Ellerkmann et al., 2001) and calcium (Uebachs et al., 2006) ion channel populations have shown that the timescale of the recovery from inactivation following a long duration of membrane depolarization increased with the length of the depolarization period.We refer to this type of behavior as history-dependence.Finally, patch clamp experiments on single ion channels have hinted at the existence of a large inactivation state space within a single ion channel (Liebovitch and Sullivan, 1987;Millhauser et al., 1988;Marom, 1998 and the references therein).These multi-level experimental fi ndings raise several important questions.How are the behaviors observed at the different levels related (e.g., Lowen et al., 1999)?Specifi cally, 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 signifi cance 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 fi rst 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

Background
The total conductivity, G(t), of a population of ion channels located on a piece of membrane is determined by A(t), the portion of channels which are available for conducting, through the equation G(t) = G 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 fi rst 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 fi nite, 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 defi nitions and properties, see Cinlar, 1975).

Model description
In view of the above comments, motivated by Liebovitch (1989) and Marom (2009), we model the slow inactivation of a channel as a continuous-time semi-Markov process consisting only of two states, the Markovian 'available' state, and the non-Markovian 'inactivated' state, as shown in Figure 1.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 fi ndings, and to concrete predictions for future experiments.We reproduce for the fi rst time, to our knowledge, the main experimental fi nding from Toib et al. (1998), namely an exponential recovery process with a history-dependent timescale, as demonstrated in Figure 3. Using these results and other similar experiments (Ellerkmann et al., 2001;Hering et al., 2004;Uebachs et al., 2006), we narrow down the options for the model parameters at different voltages for several channel types, and explicitly address the issue of long-memory phenomena (Mercik and Weron, 2001).The model introduced here also provides many predictions.Qualitatively, we predict that temporally spaced spiking stimuli will have a signifi cantly 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 defi nes 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 specifi c system addressed in this work.As pointed out in Marom (2009Marom ( , 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 specifi c and are likely to suffer from over-fi tting.Furthermore, such models always have an upper bound on their timescale.In this work, we introduce and thoroughly analyze, for the fi rst 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 fi ring.
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 specifi c 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, The RTPDF of the Markovian available state is exponential with parameter γ, The inactivated state is non-Markovian, where we use the powerlaw RTPDF: 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 fi rst time.Moreover, subtle mathematical issues arise in correctly characterizing the different regimes.None of these complexities appear when dealing with the more standard Markovian channel models.

Model motivation
Here we choose to model the, possibly complex, slow inactivation state space solely through a single non-Markovian state and its associated RTPDF.As stated above, the transition structure and rates by which the channel proteins change conformations between different slow inactivation states are still ill-understood.Moreover, the only quantity related to this state space that may be directly measured is the RTPDF.Previous single channel patch clamp experiments have already measured the RTPDFs of different channels.In many cases single channel patch clamp experiments have indeed fi tted 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 specifi c 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 (t) approximates an exponential RTPDF corresponding to a Markovian state (see Supplementary Material).

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

Relaxation to a steady state under a step voltage
First we investigate what happens to the channel when the voltage is maintained fi xed, so that c and γ are constant.By 'projecting' the single inactivated state onto a continuum of Markovian states, we derive in Eq. 16 a dynamic equation for p(t), the probability of occupying the available state, Here we used the initial condition p(0) = 1; an extension to general initial conditions is presented in Eq. 16.The intuition behind this equation is simple: the fi rst 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 fi nd that p(t) relaxes asymptotically in a power-law manner, to a steady state value, where q c t c = sin( )/( ) π πγ 0 for 0 < c < 1, q p t c = ∞ − 0 1 for c > 1, and p ∞ is the steady state value: We can express p ∞ in the more familiar form of the steady state of a simple two-state Markov process, p ∞ = + δ δ γ /( ) where δ is, as in the Markov scheme, the inverse of the mean inactivation time: In Figure 2 we compare these asymptotic results with a numerical simulation.The agreement between the two results in this case indicates that the asymptotic results are accurate for a wide range of timescales.
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 infi nite, and the rate of return to the available state is so slow, that ultimately it remains unoccupied.Assuming that real channels at rest voltage do not decay to complete inactivation, we conclude that c > 1 at rest voltage, namely c(V rest ) > 1.

History-dependent recovery timescale in response to voltage pulse
In Section 'Relaxation to a steady state under a step voltage' we studied the channel dynamics during a constant voltage step.The major feature of the model introduced here is its history-dependent dynamics, a notion which cannot be investigated for such a simple input.In order to quantify this notion, in the present section we consider a voltage pulse of constant amplitude and fi nite 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 fi xed value to another, c and γ also change.If during this voltage change the channel is in the available (Markovian) state, it does not 'remember' its history prior to the voltage change.If the channel is in the inactivated state during the voltage change, the subsequent dynamics depends on the prior history of the channel through the variable T -the duration of the channel's sojourn in the inactivated state at time t.For the inactivated channel, the probability to recover at times between t and t + dt depends on T. This probability, divided by dt, is termed the time-dependent rate of recovery, and its inverse is the time-dependent timescale of recovery, τ.A simple derivation (Eq. 35) shows that τ depends linearly on T, And so, immediately after a voltage change, the timescale of recovery is determined by Eq. 7. In this model the variable c changes instantly with voltage, while T is a continuous variable which retains the same value it had prior to the change.We note here that it is possible to model the effect of the voltage change on the time-dependent rate of recovery differently (see General Voltage Input).
For each channel, T is a random variable, and so in a population of channels, T has some distribution of values.We develop in Eq. 36 an exact asymptotic expression for this distribution, in the voltage pulse setting, where we set the initial condition so that all channels are initially available.We defi ne 〈T 〉, the mean duration of a channel in the inactivated state, and CV T T T σ / , the coeffi cient 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

Soudry and Meir
History-dependent dynamics the pulse will follow an exponential form, to a good approximation.
If CV T >> 1 this approximation will not be accurate, since then the mean provides a poor representation of the distribution.Also, when CV T >> 1 the recovery will be distinctly non-exponential.
In Section 'Renewal Theory Approach' we compute exactly the behavior of 〈T 〉 and CV T at the end of the voltage pulse, at time t.This behavior depends in an intricate way on the value of the parameter c during the pulse.And as we explained, the value of 〈T 〉 and CV T at the end of the pulse determines, through Eq. 7, the manner in which the channels recover immediately after the pulse: The recovery timescale increase linearly with t, while maintaining a constant dispersion.Also, Since CV T increases monotonically with c, the distribution of timescales becomes less dispersed for lower values of c, and hence, the recovery following the pulse becomes more exponential.
• 1 < c < 3: CV T → ∞.Thus, the distribution of recovery timescales broadens constantly, which entails a non-exponential recovery (in this case, the mean, 〈T 〉, provides a poor representation of the distribution).
3 Since both 〈T 〉 and CV T converge to a fi xed fi nite value, we get that the recovery following the pulse takes place at a single timescale.Higher values of c, lead to more exponential recovery.This implies that in this case the inactivated state is 'almost Markovian' .
The (technical) case of integer valued c is discussed in the Supplementary Material.Observe that the qualitative change of behavior noted above, results from the order of the fi rst infi nite moment of the RTPDF ψ P , which depends on the value of c.More specifi cally, the fi rst moment diverges for c ≤ 1, the second moment diverges for c ≤ 2, and so on.

Complex voltage input
So far, we have described the asymptotic channel dynamics when the voltage is constant, and also immediately after the voltage has jumped from one constant value to another.Here we briefl y discuss two other cases in which the channel dynamics may be solved analytically, and then discuss the general case of timevarying voltage.
The fi rst 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 c with the 'effective parameters' γ and ĉ given by the time-averaged values of γ and c, respectively, ˆ( ( )) ; ˆ( ( )) , where V t is the voltage at time t.Then we may again use the results derived in the case of constant voltage.Finally, when none of these approximations is valid, we show in Eq. 43 that it is possible to derive a closed form integro-differential equation generalizing Eq. 3 to the case of arbitrary input.

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

Exponential and history-dependent relaxation in response to long depolarizations
Using the model presented, we wish to reproduce the main experimental results in Toib et al. (1998).In this experiment, the membrane is clamped at a high voltage of −10 mV for varying length of time, t stim , and then clamped at low voltage of −90 mV -in which the recovery from inactivation occurred.This recovery was exponential and history dependent for every stimulus longer than 1 s.Analyzing these experimental observations using our model, we are able to restrict the possible values of the parameter c.We denote by c H , γ H and c L , γ L the model's parameter values during the high and low voltage phases, respectively.Based on the experimental results from Toib et al. (1998) we make the following claim.This claim is based on the characterization of the different regimes of c, provided in Section 'History-dependent recovery timescale in response to voltage pulse' .First, we infer that c L > 1.This follows from Eqs 4 and 5, since for c L < 1, the availability of the membrane declines to zero, a phenomenon which did not occur in the experiment.Second, we argue that c L is not in the range 1 < c L < 3, since if it were, the distribution of recovery timescales in the population would keep broadening with time, thus causing the recovery of A(t) to be non-exponential.Thus, we conclude that c L > 3.
Next, we argue that c H < 1.First, it is clear that c H < 3, since from Eq. 7 and the behavior of 〈T 〉 in this case, we know that during the inactivation period, 〈T 〉 is bounded by its steady state value 〈T 〉 = t 0 /(c H − 2).By substituting this value into Eq.7, and comparing with the result when 〈T 〉 = 0, we learn that in this case we should not expect to see large changes in the recovery timescales following different stimuli durations, contradicting the experimental results.Finally, c H cannot be in the range 1 < c L < 3 since, if it were, the distribution in the inactive state following the high voltage period would be very broad (a high value of CV T ), rendering an exponential relaxation impossible.
We thus conclude that in order to fi t the experimental results, we must set c H < 1 during the inactivation period and c L > 3 during the recovery phase.In this case, the timescale of recovery increases linearly with the stimulus duration, and the recovery is exponential as long as: 1. c H is small enough so that CV T << 1, since which is increasing in c H . 2. c L is large enough in the low voltage phase so that 〈T 〉 does not change considerably during the recovery.More formally, T >> τ = (T + t 0 )/c L , or more simply, c L is suffi ciently larger than 1. 3. γ L is small enough during the recovery period, so that during this period, the timescale of recovery does not change by 'freshly inactivated' channels; more formally, γ τ In this method of reproducing the experimental results the channel possesses 'infi nite memory' , since as long as the channel remains at the high voltage, 〈τ 〉 continues to increase linearly: We note here that it is possible to modify the model so that t 0 is the voltage-dependent parameter and c is fi xed at some value ( >1).In this case, the channel does not possess this 'infi nite memory' and the timescale of recovery is bounded.Therefore, in order to reproduce the experimental results in this modifi ed model, we have to assume that the increase in the timescale of recovery has some upper boundary, which was not yet reached in the experiment.Since it seems unnatural to assume the existence of such upper limit, this alternative model was not used.
In an additional experiment the recovery was examined under several voltages: −60, −90, −120 mV.The history-dependent behavior remained, as can be seen in Figure 6 in Toib et al. (1998).The recovery was similar for −90, −120 mV, indicating that c did not change much between these voltages.At −60 mV, the recovery was slower, and possibly less exponential, indicating perhaps that 1 < c L < 3, in that case.
A further experimental observation relates to the emergence of a history-dependent relaxation only for t stim ≥ 1 s.From this threshold point between the two modes of recovery and the results in Section 'History-dependent recovery timescale in response to voltage pulse' we get that in this experiment t 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.

Channel response to a spiking stimulus
Motivated by the results in Ellerkmann et al. (2001), Toib et al. (1998) and Uebachs et al. (2006), we study the model's dynamics when the voltage input is a periodic voltage spike train.As noted in Toib et al. (1998), such inputs are similar to fi ring patterns in neocortical neurons.This type of input may be important when considering the effects of a neuron's action potential on itself.We note that such an input does not present the entire picture, since synaptic inputs from other neurons are probably more realistically described as sums of continuous functions (Gerstner and Kistler, 2002).In any event, it is interesting to test the model's prediction in this setting, for which some experimental results are available (Toib et al., 1998;Ellerkmann et al., 2001;Uebachs et al., 2006).The voltage spikes are modeled here as a square wave -for T H seconds the voltage is set high, and for T L seconds FIGURE 3 | Simulation of experimental results.(A) Exponential and historydependent recovery in response to various stimulus lengths (R 2 > 0.99 for all exponential fi ts).Fraction of inactivation is defi ned as in Toib et al. (1998).Legend: different stimuli lengths, in seconds.(B) The increase in the timescale of recovery with stimulation length, in both the numerical simulation and the analytic expression (Eq.9) (with c H < 1).We see that the two results agree, especially for long stimulations.This demonstrates how, in this case, the asymptotic analytic result is a good approximation for the relevant timescales.Assuming that the input spike frequency is much higher than the transition timescales of the model, we use the effective parameter approximation.The behavior of the system may be inferred by replacing γ and c with the effective parameters, according to Eq. 8: It is important to notice in this case that if T H << T L (sparse spikes), then γ will be affected by γ H only if γ H >> γ L , and similarly for ˆ.
c Since we know from the experimental results that c H < c L , in this case we can approximate, ˆ.
c c L ≈ Using this approximation we address the dependence of γ on the voltage.From Toib et al. (1998) it is known that the sodium channel population goes into significant inactivation as a result of an action-potential-like stimulus.This is also the case in Ellerkmann et al. (2001).Since in this case ĉ c L ≈ (does not change much), then γ must be signifi cantly larger than γ L -so γ H must be signifi cantly higher than γ L (otherwise, the experimentally observed inactivation would not be reproduced).This means that γ(V) must increase with voltage, at least between the two voltage values used in this setup.
The distribution of recovery timescales in the channel population is affected mainly by the value of c (the dependence on γ is absent from the asymptotic form Eq. 36).If we wish to make the timescale of recovery change measurably in response to spike stimulation, ĉ must be signifi cantly different from c L .For example, if a neuron is affected solely by its own action potentials, then T H ≈ 1 ms, T H + T L > 10 ms, and therefore c c c L L > > ˆ. .0 9 So, in this case it would be hard to experimentally detect an increase of the recovery timescales due to the input (see Figure 5).This assertion will remain valid for any type of stimulus in which T H << T L -the channel may enter inactivation, but any history-dependent recovery will be mitigated.This complies with the results shown in Figure 5 of Ellerkmann et al. (2001), where we can see a relatively weak scaling of the recovery timescale with stimulus frequency -even though T H /(T H + T L ) is rather high in comparison with typical action potentials.
It is important to note that in the case ˆ, c > 3 the non-Markovian state can be approximated by a simple Markovian state, since in this case 〈T 〉 < t 0 , and thus τ = + ( ) drastically due to history-dependent effects.This may happen, for example, if c(V rest ) > 3, and the input to the neuron is suffi ciently 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 −55 mV depolarizations was observed, similar to what was measured in sodium channels (referred to as 'continuous' in Figure 8 in Uebachs et al., 2006).In contrast, after a spiking stimulus at a maximal voltage of 25 mV with T H /(T H + T L ) ≈ 0.056 ('mock APs' in Figure 8 in Uebachs et al., 2006), there was no observed increase in the timescale of recovery with the length of stimulation, in agreement with our model.Also, when the spiking stimulus was superimposed on the (continuous) step depolarization, the history-dependence of the timescale diminished ('continuous + mock APs' in Figure 8 in Uebachs et al., 2006).This may imply that c(V ) is not monotone in calcium channels -if c(V ) decreases between −95 and −55 mV and then increases between −55 and 25 mV, ĉ should change less during the stimulus, and therefore the recovery timescale should also change less.This result might explain why Hering et al. (2004) did not observe any change in the recovery timescale after different lengths of depolarizations, given at the higher voltage of −20 mV.

TEMPORAL CORRELATIONS
In Section 'The Joint Probability Distribution of the Availability' , we present a simple approach to calculating the joint moments and distributions to be in the available state at different times.Specifi cally, we derive Eq. 57, the asymptotic behavior of the stationary auto-covariance of a single channel, for c > 1: is not defi ned 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 . This implies that the process is characterized by a long-memory (see defi nition on page 42 of Beran, 1994).This result complies with the results of Mercik and Weron (2001) for potassium channels (on timescales <1 s), as also pointed out by Goychuk and Hanggi (2004).

RELATION TO PREVIOUS WORK
Several previous studies have attempted to construct a mathematical modeling framework which can account for long term temporal correlations, power-law behavior, and history-dependence of responses related to ion channels.The initial work along these lines was motivated by single channel experiments on time scales of milliseconds-to-seconds (see Liebovitch and Sullivan, 1987;Millhauser et al., 1988;Marom, 1998 and the references therein), and therefore does not directly imply the tail behavior of the RTPDF on the timescale of seconds-to-minutes, which was explored in Toib et al. (1998).In this brief comparative discussion we focus on models addressing channel dynamics.
As far as we are aware, the paper by Millhauser et al. (1988) was the fi rst 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 diffi culty.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 fi ndings, 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 fi nding from Toib et al. (1998), namely an exponential recovery process with a historydependent 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) (defi ned 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).Specifi cally, 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 fi nite 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 fi eld 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 powerlaw RTPDF with 1 < c < 2, which was used for the ion channels in that paper.

DISCUSSION
In this work we have thoroughly analyzed the dynamics of a generic two-state model of an ion channel consisting of a Markovian state and a non-Markovian state.We derived a general dynamic equation for the probability of the channel to be available, or equivalently, the fraction of available channels in a voltage clamped population.This dynamic equation, which is derived for both constant (Eq.3) and time-dependent (Eq.43) voltage, defi nes a direct input (voltage) to output (availability) relation.We derived explicit solutions to this equation in many important cases, and studied their properties.Specifi cally, 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 fl uctuations in the voltage (Eq.8).Also, we derived simple expressions for all joint moments of the availability and specifi cally, for the auto-covariance function of the channel (Eq.11).An interesting conclusion of this analysis is that the channel is characterized by four different 'modes' in which its behavior is rather different, depending on the value of c (see Table 1).
The main experimental fi nding 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 confi rmed 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 specifi c qualitative predictions are the following.The value of the parameter c at rest obeys c(V rest ) > 1 in all healthy channels (see Relaxation to a steady state under a step voltage).From Section 'History-dependent recovery timescale in response to voltage pulse' we conclude that the recovery from slow inactivation of NaII and NaIIA channels is voltage dependent, and for these channels, there exists some range of voltages for which 1 < c < 3, where the recovery becomes non-exponential.Finally, we predict that sparse spiking stimuli (e.g., the neuron's own action potentials) can induce slow inactivation but only minor changes in the timescale of recovery from inactivation (see Channel response to a spiking stimulus).
One of the more interesting questions pertains to the relation between the history-dependent channel dynamics and the generation of action potentials in the cell, and, more concretely, the possible functionality of such behavior (Lundstrom et al., 2008;Wark et al., 2009).In order address this issue clearly and fully, two further steps must be taken.
Experimentally, it is necessary to fi nd 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 longrange 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 fl exible, and may be used to describe other systems in which one can similarly defi ne separate available and unavailable states so that the transitions between

METHODS
In this section we establish several general analytical methods, which can be used in any two-state model in which one state is Markovian while the other state is non-Markovian, namely ψ E (t) is still exponential as in Eq. 1, but ψ P (t) may be replaced by some other general RTPDF, which we denote by ψ I (t).When establishing general results we use ψ I (t), and when specializing them to the power-law RTPDF we use ψ P (t).
We present two different formalisms -the Markovian-process based 'projection method' in Section ''Projection' of the Inactivated State' and the 'Renewal Theory' approach in Section 'Renewal Theory Approach' .This was done in order to explore different aspects of the channel dynamics in the case of constant voltage input.Moreover, we show how to map an existing Markov model to a two-state non-Markov model.In Section 'General Voltage Input' we use the presented methods to derive the dynamics in the case of general input.In Section 'Oscillating Voltage Input' we consider the case of rapidly oscillating input.Finally, in Section 'The Joint Probability Distribution of the Availability' we calculate the joint distributions of the channel.

'PROJECTION' OF THE INACTIVATED STATE
In this section we assume that the voltage input is constant, so that all the model parameters are also constant.We wish to make use of the well-established formalism of Markov processes to derive the dynamics of our non-Markovian model.To do so we needed to fi nd a way to replace the non-Markovian inactivated state with an equivalent Markov inactivation state space.For example, in the context of the specifi c channel model with the power-law RTPDF ψ P (t), there are many possible state spaces that may be used to produce it (Liebovitch, 1989).Some of them may be physiologically more accurate than others, but, mathematically, this is inconsequential to the channel dynamics.And so, we chose the simplest state space, in which all inactivation states are parallel -each state being connected only to the available state, as seen in Figure 6.Though the fi gure depicts a fi nite set of states, in the actual model we use a continuum of states, so that the equivalence between the models is exact.We use as the continuous 'index' of each activation state its recovery rate, δ.Since it is Markovian, the RTPDF function of each inactivation state δ is exponential, 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 specifi c inactivation state δ.This means that the inactivation rate from the available state into the inactivated states (δ,δ + dδ) is γf(δ)dδ, and so the total rate of inactivation is ∫ = ∞ 0 γ δ δ γ f d ( ) , as required.With this condition fulfi lled this multiple-state Markovian model is equivalent to our two state non-Markovian model if we ensure that the RTPDF of the aggregation of all the inactivated states is ψ I (t).To fulfi ll this condition, we use law of total probability and demand that:

Dynamic equations
In order to derive the dynamical equations of the non-Markovian model, we begin with the Markovian model described in Figure 6, present its underlying equations, and then take the continuum limit.
We denote π π π = ( , , , ), where p is probability of being in the available state, and π k is the probability of being in the inactivated state I k .The dynamic equation for the probability of this homogeneous Markov process is where , is the rate matrix.
Writing this explicitly, we get 0 yielding: Dividing the fi rst equation by dδ we obtain: with solution Inserting this into Eq.13 and switching the order of integration we get: ∞ ∫ Using Eq. 12, the expression can be re-written as: This result is quite intuitive: 1.
where the tilde denotes the Laplace transform.The fi nal term in the numerator results from the initial probability distribution in the inactivated states.If at t = 0 the channel is in the available state, namely p(0) = 1, then this term vanishes.

Projection of a power-law RTPDF
By using an inverse Laplace transform on the t variable in Eq. 12, we can easily show that in the case where ψ I (t) = ψ P (t), f(δ) is distributed according to the gamma distribution, where Γ(•) is the gamma function.Thus, the aggregation of all the inactivation states has the RTPDF: which equals ψ P (t) as required.

Asymptotic solution of the dynamic equations for power-law RTPDF
In this section we demonstrate how to derive an asymptotic solution of Eq. 17 in the case where ψ I (t) = ψ P (t).We assume throughout this section that c is non-integer; the limiting case where c is an integer is discussed in the Supplementary Material.In order to understand the behavior of p(t) for large values of t, we need to consider the small s limit of p s ( ).To do so, we fi rst calculate ψ P s ( ).In the Supplementary Material we prove that for non-integer values of c, Retaining the leading order in s, we have (see Supplementary Material): and substituting Eq. 20 into Eq.17, we fi nd: where we have introduced denoting the steady state value of p(t) -a result that can be confi rmed by using the fi nal value theorem on p s ( ).Also, Notice that in Eq. 22 we deliberately retained the s c−2 term, even though s c−2 = O(1) for c > 2. The reason for this will become clear soon.
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 defi ne m = ⎡c − 2⎤, the upper integer part of c − 2, and differentiate m times: Since that all the integer order terms with degree lower than m vanish after we differentiate for m time.Using the fact that: the linearity of the Laplace transform and the above Tauberian theorem, we get: where we used similar identities for Γ(•) as used previously, Eq. 23 and the fact that m is an integer.Notice that we can use this method to fi nd 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 L[f(t)] (assuming that L[f(t)] can be uniquely represented as a countable sum of integer and fractional powers of s).Because of this, any (minimal) linear time invariant system will behave asymptotically as a power law if it contains some component with a power-law kernel -which always introduces a fractional power of s into the system's transfer function in the Laplace domain.
In conclusion, we found the asymptotic behavior of p(t) as t → ∞, for all c, We derived this result for the case p(0) = 1, but it remains valid also when p( 0 )] δπ δ δ δ s d 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 fi nite supremum to the recovery timescales of all the inactivation states in the initial distribution.There is another way to see that the asymptotic solution in Eq. 24 is not too sensitive to the initial conditions.Suppose we set p(t i ) = 1 for some t i , instead of p( 0 and therefore Eq. 24 is changed by the addition of an asymptotically negligible term.This asymptotic insensitivity to initial conditions stands in contrast with the case where the inactivated state is Markovian, and the solution is exponential, In this case when we set the initial condition to t i , instead of 0, we get, Therefore, the solution is changed by a pre-factor, which is not asymptotically negligible as in the power-law case.

How to create a two-state non-Markov model from a general Markov model
So far, we used the Markov process formalism to fi nd 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?

Soudry and Meir
History-dependent dynamics 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 t dt t π π ( )/ ( ) , = Q where π( ) 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 π π π ( ) ( ( ), ( )).

t t t
A I = Similarly, we divide the rate matrix Q into block matrices corresponding to these sets of states: .
In order to calculate the RTPDFs of the inactivated states, we introduce a modifi ed Markov process with probability vector π π π ( ) ( ( ), ( )), t t t A I = and the rate matrix: The dynamics of the modifi ed 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 The solution of Eq. 25 can be written in the form, π π where the exponent of a matrix A is defi ned as exp( ) / !.
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 modifi ed process is in some available state at time t.This probability can be written as π A A T t u ( ) , where we defi ned u A T T = ( , , , ) , 1 1 1 … a column vector of ones with the same number of components as π A t ( ).Differentiating this expression, we get the RTPDF of the inactivated state, assuming that the channel has inactivated at t = 0: where we used Eqs 26 and 27, and π I ( ) 0 represents the probability distribution in the inactivated states immediately after the inactivation.Similarly, the RTPDF of the available state is where u I T is defi ned similarly to u A T .
Note that ψ I (t), ψ A (t) are always a sum of exponents, but after the transformation to the two state model is complete, we can approximate them by some simpler functional form.The form of these RTPDFs (in Eqs 28 and 29) pose a problem, through their dependence on π A ( ) 0 and π I ( ). 0 Assume that the channel inactivated at time t 1 , and recovered at time t 2 .Generally, π I t ( ) 1 may depend on the time that the channel was available before t 1 , and t 2 − t 1 may depend on π I t ( ).
1 Therefore, the time the channel spends in inactivation may depend on the previous duration of the channel in the available state, which may depend on the previous duration the channel in the inactivated state, and so on.This would prevent us from constructing a model with two independent non-Markov states.
There are some specifi c cases when this problem vanishes.One of them is when A can be approximated as single Markov state.This happens, as we mentioned in Section 'Background' , when the rates between the available states are much larger than the rates of inactivation.In this case, as we show in the Supplementary Material, where we defi ned 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.

RENEWAL THEORY APPROACH
Despite its merits, the use of the projection method may be limited in certain cases.The problem lies in the use of π δ (t), the probability distribution in the inactivation states.Knowledge of the probability to be in an 'inactivation state' implies that we have some additional information about the future recovery time of the channel.But, when a channel is inactivated, we know only for how long it has been inactivated, and not to which 'inactivation state' it has arrived.Since these inactivation states are not observable, this approach restricts the practical use of the initial conditions π δ (0) in Eq. 16.And so, in order to quantify the history-dependent behavior of the channel in terms of observable quantities, we turn to a different formalism that uses only the observable T, which is the time during which the channel has already been inactivated.In this section, as the in the last section, we again assume that the voltage input is constant, so all the model parameters are also constant.
For a channel that has been inactivated for time T, the probability of recovery in the next time interval Δt is where t R denotes the time it takes the channel to recover since its last inactivation.Therefore, the rate of recovery immediately after time T is and the time-dependent timescale of recovery is its inverse τ λ = Δ 1/ .
Next, we compute f t (T), the probability density function of T at time t, for a single channel, which well approximates the distribution of T throughout a suffi ciently large population.We did so using a 'Renewal Theory' formalism, which is described in Cox (1962).In renewal theory language T is called the 'backward-recurrence time' until the last inactivation event.
First, we denote by k n (t) the probability density function of the time of the n-th inactivation event, assuming that at t = 0 the channel was available (p(0) = 1).Clearly k n (t) is the probability density function of a random variable which is the sum of 2n In the Laplace domain, this translates to, Defi ning the rate of inactivation events: and in the Laplace domain: where we used the geometric series summation formula.Notice that h s 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 (T), conditioned on the fact that the channel is inactivated at time t.For a channel at time t to be inactivated for the last T seconds, there must be an inactivation event at time t − T and since that time the channel must not recover.Dividing this by (1 − p(t)), the probability that the channel is indeed inactivated, we get,

The case of power-law RTPDF
In the context of the specifi c model with power-law RTPDF the rate of recovery immediately after time T is and the time-dependent timescale of recovery is its inverse: The distribution of T in this case ψ Because of the multiplication in the last expression, it is generally diffi cult to write an expression for f t (T) in the Laplace domain.We therefore examine f t (T) in the limit t → ∞.We get, by using the same asymptotic derivation we used in deriving Eq. 24, that: Note that for 0 < c < 1, for any fi xed value of T, the function f t (T) vanishes for large values of t.This implies that the probability distribution shifts to infi nity as t increases; see Figure 7.
Calculating the fi rst two moments of T in each case, defi ning 〈T 〉 as the mean of T and σ T its standard deviation we obtain the following results by series expansion (see Supplementary Material).
We verifi ed numerically that all of the results for these moments and for f t (T) itself are valid for the relevant timescales, as is demonstrated in Figures S1 and S2 and in Supplementary Material.

GENERAL VOLTAGE INPUT
In the previous sections we used two different formalisms to analyze our model.One was the Projection Method, and in the other was the formalism of Renewal Theory.So far, these two approaches may seem unrelated, and dealt only with the case of constant parameters.In this section we shall use elements from both formalisms to derive a dynamic equation for p(t) for any time-varying parameters.For example, in the context of the spe-

Soudry and Meir
History-dependent dynamics cifi c channel model with power-law RTPDF, recall that the voltage input infl uences the channel through the variables c and γ, which are functions of the voltage.

Generalization of the time-dependent rate
One way to derive the general dynamic equation, is to generalize the defi nition of the time-dependent rate so that it also depends on the 'global time' t, where t is the current time, u is the time of the last transition between states, and t S is the time of the next switch in channel states.
The probability that the channel will switch back to the other state in the time interval (t, t′) can be calculated by dividing (t, t′) into n smaller time segments (t i , t i + Δt),i = 1,2,…,n.Breaking up the interval (t, t′) into these segments, we have, where we have defi ned The last equality can easily be established by taking the derivative of the fi nal term with respect to t' .
We have thus shown that: Next we introduce a 'time-dependent RTPDF' , can be thought of as an extension of the RTPDF of the type introduced in Section 'The Model' .This function is the probability density to switch back into the other state at time x: Now that we have defi ned the time-dependent RTPDFs in Eq. 42, we can generalize the dynamic equation developed in Eq. 16 (through the projection method).Setting p(0) = 1 for simplicity, we get, which allows us to calculate the behavior of p(t) for any timevarying input.Notice that the term ( ) ( ) t d 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 defi ned by π δ (0), which are unobservable.

Soudry and Meir
History-dependent dynamics

The case of power-law RTPDF
For a constant rate λ(t|u) = γ we get, as expected, the Markovian exponential RTPDF, while for a varying rate γ(t) we get, For the the non-Markovian state with a constant input λ(t|u) = c/(t − u + t 0 ), we get as expected, while for general input, if we assume that we obtain 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 , , and so we get the following dynamic equation, where we denoted And so, the initial conditions for the evolution of p(t) from time t 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.

Direct generalization of the RTPDF
Another way to generalize the channel model to include varying input is to directly generalize the defi nition of the RTPDF, instead of going through the time-dependent rate.This could be done if we assume that the two-state non-Markov model is based on some underlying many-states Markov model, as explained in Section 'How to create a two-state non-Markov model from a general Markov model' .In the case of time-varying rates, it is straightforward to generalize Eqs 30 and 31 and get where we defi ned γ π ( ) 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

Soudry and Meir
History-dependent dynamics 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 signifi cantly 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 fi nal choice between Eqs 45 and 52 or some other, more complicated expression, should be made through carefully designed experiments on channels.

OSCILLATING VOLTAGE INPUT
We examine the case where the input to the model oscillates with period T P .These oscillations may arise, for example, from periodic stimulation or some noise process (in which case T P is only approximate).For example, in the context of the specifi c channel model with power-law RTPDF the parameters γ, c may generally be voltage-dependent and vary with time, while t 0 is maintained constant.Intuitively, if T P is suffi ciently small in comparison with the transition timescales, then on suffi ciently large timescales we should not notice these fl uctuations in the parameters (even though they may be large in magnitude), and expect them to 'average-out' .
Here we prove this claim rigorously, and fi nd the exact conditions for this to happen.
For any time-dependent variable X S defi ne the mean over a period T P , Again, t denotes the current time, u is the time of the last transition between states, and t s is the time of the next switch in channel states.Using Eq. 41, we can write, for all θ > 0, where we denoted p(t i |t j ) as the probability of the channel to be in the available state at time t i assuming that at t j the channel was also in the available state.For arbitrary voltage input, each p(t i |t j ) term in this product is the solution of Eq. 43, with the initial condition p(t j ) = 1 (notice also that the lower limit of integration in Eq. 43 must be set to t j ).In the case that the voltage is constant, p(t i |t j ) = p(t i − t j ), where p(t) is the solution of Eq. 16 with the initial condition p(0) = 1.If we assume also, for simplicity, that p(0) = 1, we get, ) )

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

THE MARKOVIAN LIMIT OF THE NON-MARKOVIAN STATE
In Section 'Model Description' we defi ned the RTPDF of the non-Markovian state to be: If we naively take the limit that c → ∞, then ψ P (t) → δ (t), a Dirac delta function at zero, meaning that the channel will recover instantly.And so, in order for this limit to be meaningful, we need to 'zoom in' in the time axis to see what happens for very short times.To do this, we change the variables of the RTPDF -θ c t ⋅ , so that: where we used the limit lim x→∞ (1 + 1/x) x = e in the last line.In conclusion: in the limit of c → ∞ the non-Markovian state approximates a Markovian state, for very short times.

CALCULATION OF P s ( )
In this section we develop an expression for ψ P s ( ), the Laplace transform of ψ P (t).This will give Eq.19, which is the basis for the derivations in Section 'Asymptotic Solution of the Dynamic Equations for Power-Law .By the defi nition of the Laplace transform: The numbered steps are explained as follows: 1.This transform is properly defi ned for all s such as Re[s] ≥ 0.
2 2 2. We denoted z ≡ r(t + t 0 ). 3. We assumed here that c is not an integer, used Fubini's theorem (which we can use since <∞) to switch the order of summation and integration, and denoted F( ) ω to be some function of ω. 4. We used the fact that ψ( ) s must be piecewise analytic since for every n: which always exists ∀Re[s] > 0 -therefore F(ω) (t 0 r) c must be a function of s -and the only one possible is G (t 0 s) c , where G is some constant.
If we assume in step 3 that c is an integer, than we get instead: Also, since we can write: we can compare this derivation (which was done for complex s) with the one in Bender and Orszag (1978: 251-252) where it is shown that G = Γ (−c) < 0, in the case that c is not an integer, and if c is an integer then: where γ is Euler constant (do not confuse it with the inactivation rate).From the analyticity of ψ P s s ( ) Re[ ] ∀ >0 we come to the conclusion that G is identical also in our case.

ASYMPTOTIC BEHAVIOR OF P s ( )
In this section we examine the expressions for ψ P s ( ) which were derived in Eqs Sl and S2 for various values of c and retain the lowest orders in s.For non-integer values of c this will give Eq.20, which appears in Section 'Asymptotic Solution of the Dynamic Equations for Power-Law RTPDF' .
When c is not an integer: Soudry and Meir History-dependent dynamics For c = 2,

ASYMPTOTIC BEHAVIOR OF p s ( ) FOR INTEGER VALUES OF c
In Section 'Asymptotic Solution of the Dynamic Equations for Power-Law RTPDF' we derive an asymptotic solution for p s ( ) in the case that c is not an integer.In this section we show a similar derivation for the integer case.
Assuming p(0) = 1, we use the Eq.17: For c = 2 we get: is, as in the non-integer case, the steady state of p(t) obtained from the from the fi nal value theorem.Unfortunately, since the Tauberian theorem does not apply to logarithmic expressions we cannot apply it to infer the asymptotic behavior of p(t).From numerical simulations we see that the present case (integer valued c) is not pathological, and we can simply take it as the limit of the non-integer case.

THE MOMENTS OF f t (T)
In this section we derive Eqs 37-40.
To do this we need calculate the moments of Eq. 36: First we examine the case 1 < c The mean: The behavior of the last expression in the limit t → ∞ is determined by the value of c.
The second moment: ) ) The behavior of the fi nal expression in the limit t → ∞ is again determined by the value of c.
Using the above expression we get: and, For the case of integer valued c, again we observed by numerical simulations that it can be well approximated by the limiting case of the non-integer case.
In Figures S1 and S2 we show a comparison between numerical simulations and the asymptotic exact results we derived in this section, for different values of c: c = 0.5,1.5,2.5,3.5 (from top to bottom).The close agreement again demonstrates the validity of the asymptotic results in Eqs 37-40, for a wide range of timescales.

WHY IS THE AVAILABLE CLUSTER OF STATE MARKOVIAN?
In Section 'Background' we mentioned that transitions between the states that compose the available state are much faster than the transitions to and from the non-available states, and therefore the available cluster of states is, to a good approximation, Markovian.
In this section we use the methods presented in Section 'How to Create a Two-State Non-Markov Model from a General Markov Model' to show why this reasoning is correct.Recall Eqs 28 and 29: for the RTPDFs of the 'lumped' states.When the above conditions hold, we will show that ψ A (t) = γ exp (−γt) for some non-negative γ, and fi nd an explicit form for ψ I (t).
We denote by M the number of available states, and by L the number of inactivated states.Recall that the diagonal element of any Markovian rate matrix equals the negative sum of all non-diagonal elements.Since the transitions within the available states are much faster than transitions outside of this state, we can write: where Q AA is the rate matrix for the transitions inside the available states, with all the transitions to inactivation removed from the diagonal, and εD is a diagonal matrix that contains all the diagonal components of inactivation removed from Q AA .Formally, we can write ( ) ( ) .εD Since the rates of transitions to inactivation are much smaller than transitions inside the available state ε l − a small parameter.Using perturbation theory (Landau and Lifshitz, 1981, chapter VI) ε The last equation holds for all ε, therefore we can equate every order in ε separately.For the zeroth order in ε, Note that this pertubative expansion holds only in the case where the eigenvalue E m is non-degenerate.This applies to E 0 0 = , if we assume that Q AA is indecomposable, so that it has only a single stationary distribution π A s : this stationary distribution is φ 0 , the eigenvector of Q AA that corresponds to the eigenvalue 0. Since φ 0 is the only eigenvector with eigenvalue 0, and since Q AA A T u = 0 (because Q AA is a rate matrix), we get, from Eq. S6:  This the fi rst order correction to the 0 eigenvalue of π A s or, more simply put, the approximate eigenvalue of π A s as an approximate eigenvector of Q aa .The corrections to the other eigenvalues of Q AA are negligible compared to their original, non-zero, values.Additionally, we note that we can write: Thus, we obtained an exponential distribution as in Eq. 30, and the cluster of available states is Markovian.
Next, we calculate the RTPDF of the inactivated state in this case.Assume that the channel was in the available state from time 0 to t and then inactivated.From the law of total probability:

THE RENEWAL EQUATION
In this section we show how can we re-write Eq. 43 so it resembles the form of a Renewal Equation (see Cinlar, 1975).This expression helps gain a better understanding of the dynamic behavior.This time dependent generalization of the Renewal Equation, has a clear intuitive meaning.This fi rst term represents the probability that the channel has not inactivated since time t = 0.The second is a sum over all the probabilities for the channel to inactivate at time u, recover at time v, and remain at the available state until time t.
Note that in the constant input case the system becomes time invariant, and the integral with non-stationary kernels ψ Ψ

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

FIGURE 2 |
FIGURE 2 | Behavior of A(t) -comparison between numerical simulation and the analytical expression of the asymptotic behavior of A(t):(A) c = 0.5, (B) c = 1.5.The remaining parameters are: N = 10 6 , γ = 1 Hz, t 0 = 1 s and a simulation step of Δt = 5 × 10 −3 s.The agreement between the analytical and numerical results in this case demonstrates that the asymptotic analytical results are accurate for a wide range of timescales.
Claim: During inactivation c = c H < 1, while during recovery c = c L > 3.
The linear relationship is expected, based on the analytic result, to persist for arbitrarily long stimuli.The parameters for the simulation were N = 10 5 , γ L = 10 −4 Hz, γ H = 1 Hz, c H = 0.2, c L = 15, t 0 = 1 s and a simulation step of Δt = 10 −3 s.April 2010 | Volume 4 | Article 3 | 7 Soudry and Meir History-dependent dynamics the voltage is set low.As before, we denote by c H , γ H and c L , γ L the model parameters during high and low voltage, respectively (see Figure 4).

FIGURE 4 |
FIGURE 4 | The voltage and model parameters, during a spiking stimulus.(A) Each stimulus period is divided into the time of the spike (T H ), and the interspike interval (T L ). (B) The parameter γ assumes the values γ H during the high voltage, and γ L during the low voltage.The value γ is the "effective parameter"the average of the parameter -during a single period.As can be seen in the

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

FIGURE 6 |Frontiers
FIGURE 6 | Approximation of a semi-Markovian channel by a number of Markovian inactivation states.

FIGURE 7 |
FIGURE 7 | A sample of f t (T) at different values of t and c: (A) c = 0.5, (B) c = 1.5.The circles represent results of the numerical simulation, while the lines are the analytical asymptotic expressions (Eq.36).The numerical results are based on the same simulation as in Figure 2. Since T is and eigenvectors of Q AA and by { } eigenvectors of Q AA , we get:

FIGURE S1 |
FIGURE S1 | Comparison between numerical simulations and asymptotic analytical expressions for (T ), the mean of current duration of inactivated channels in the inactivated state, for different values of c.The parameters used are: a channel population of N = 10 6 , γ = l Hz, t 0 = 1 s and a simulation step of dt = 5 × 10 −3 s.

FIGURE S2 |
FIGURE S2 | Comparison between and numerical simulations and asymptotic analytical expressions for σ T , the standard deviation of current duration of inactivated channels in the inactivated state, for different values of c.The parameters used are: a channel population of N = 10 6 , γ = l Hz, t 0 = 1 s and a simulation step of dt = 5 × 10 −3 s.
are constants and we recall that for every m, E m < 0. Since for any nonzero value of m, probability distribution (and this only occurs if a 0 = 1).Therefore for t ≥ 0, defi ne the survival probability function in the Markovian state: , 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.Now we can re-write the last equation and get:

Table 1 | Modes of behavior for different ranges of c.
Friedlander and Brenner (2009)nt dynamics 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 byFriedlander 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'.
(2, 3)The channel has a stationary steady state in which it is partially available and recovery is non-exponential.(3,∞)Thechannel has a stationary steady state in which it is partially available, the recovery timescale distribution has a fi nite mean and dispersion and is near-exponential.The inactivated state is 'almost Markovian'.April 2010 | Volume 4 | Article 3 | 11 Such equations are commonly dealt with in the Laplace domain.Using a Laplace transform on Eq. 16, we get: ) = 1.The asymptotic solution is p(t − t i ).Since the n-th inactivation event.Since the RTPDF of X m , Y m are, respectively ψ E (t) and ψ I (t), and the probability density function of a sum of independent random variables is a convolution of the probability density functions of all the random variables, we get,

The case of power-law RTPDF In
the context of our model, this means that if, for all t, T P << γ −1 (t), we can approximate, is the time of inactivation, and the index t was suppressed in ˆ, γ t since it is assumed to be, at least approximately, constant.Also, if T P << min s [t 0 /c(s)] and T P << t 0 then we can approximate, using Eq.45,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 t 1 ,…,t k , If we multiply the general dynamic equation for the probability: Note that the Laplace domain representation of the equation is limited, since we cannot generalize it to the nonstationary, varying input case, as we did in the time domain.