article: The voltage trace of neuronal activities can follow multiple timescale dynamics that arise from correlated membrane conductances. Such processes can result in power-law behavior in which the membrane voltage cannot be characterized with a single time constant. The emergent effect of these membrane correlations is a non-Markovian process that can be modeled with a fractional derivative. A fractional derivative is a non-local process in which the value of the variable is determined by integrating a temporal weighted voltage trace, also called the memory trace. Here we developed and analyzed a fractional leaky integrate-and-fire model in which the exponent of the fractional derivative can vary from 0 to 1, with 1 representing the normal derivative. As the exponent of the fractional derivative decreases, the weights of the voltage trace increase. Thus, the value of the voltage is increasingly correlated with the trajectory of the voltage in the past. By varying only the fractional exponent, our model can reproduce upward and downward spike adaptations found experimentally in neocortical pyramidal cells and tectal neurons in vitro. The model also produces spikes with longer first-spike latency and high inter-spike variability with power-law distribution. We further analyze spike adaptation and the responses to noisy and oscillatory input. The fractional model generates reliable spike patterns in response to noisy input. Overall, the spiking activity of the fractional leaky integrate-and-fire model deviates from the spiking activity of the Markovian model and reflects the temporal accumulated intrinsic membrane dynamics that affect the response of the neuron to external stimulation.
  The leaky integrator properties of a neuron are determined by the membrane resistance and capacitance which define a single time constant for the membrane voltage dynamics [1]–[3]. However, the voltage trace of real neurons can follow multiple timescale dynamics [4]–[6] that arise from the interaction of multiple active membrane conductances [7]–[11]. Such processes can result in power-law behavior in which the membrane voltage cannot be characterized with a single time constant [5], [12]–[15]. Since power-law dynamics can span all the scales of interest of neuronal behavior [16]–[20], it is necessary to develop a framework to study such processes and their effect on the electrical and computational capacities of neurons. In the classical leaky integrate-and-fire model the temporal evolution of the voltage is local [21], [22]. The value of the voltage at a given time depends only on the value of the voltage in the immediate previous time step. Such a process is called Markovian. However, coupling of active conductances does not allow the value of the voltage to be memoryless [11], [17], [23]–[26]. Instead, long time correlations affect the membrane voltage for hundreds of milliseconds. The emergent effect of these membrane correlations is a non-Markovian process that can be modeled with a fractional derivative [27]–[31]. A fractional derivative represents a non-local process [32]–[34] in which the value of the variable is determined by integrating a temporal weighted voltage trace, also called the memory trace. Although fractional derivatives and integrals are almost as old as traditional calculus [32], [35], [36], they have not been widely used due to limited computer power. In the fractional integrate-and-fire model the exponent of the fractional derivative goes from 0 to 1, with 1 representing the normal derivative. As the exponent of the fractional derivative decreases, the weights of the voltage trace increase. Thus, the value of the voltage is increasingly correlated with the trajectory of the voltage in the past. We developed and analyzed a fractional leaky integrate-and-fire model. The only parameters of the model are the conductance, capacitance, and the fractional exponent. By varying the fractional exponent our model reproduces the upward and downward spike-frequency adaptations found experimentally in pyramidal neurons [37], [38], tectal neurons [39] and fast-spiking cells of layer IV in the rat barrel cortex [40]. Furthermore, the model replicates not only the adapting firing rate but also the long first-spike latency seen in pyramidal neurons in layer V [38]. The model also produces spikes with longer first-spike latency and high inter-spike variability with power-law distribution, which cannot be reproduced by the classical integrate-and-fire model. We further analyze spike adaptation and the responses to noisy and oscillatory input. Overall, the spiking activity of the fractional integrate-and-fire model deviates from the spiking activity of the Markovian model and reflects the temporal accumulated intrinsic membrane dynamics that affect the response of the neuron to external stimulation.
  We define the fractional leaky integrate-and-fire model as (1) along with the fire-and-reset condition (2) where is the membrane potential, and is the order (exponent) of the fractional derivative, with. In the case of, the fractional model is the same as the classical leaky integrate-and-fire model. When the membrane potential reaches a threshold (), a spike is generated and is reset to for a refractory period. The passive membrane time constant is. Parameter values are given in Table 1 (see Methods). For, the fractional derivative of the voltage (in Eq. 1) can be defined with the Caputo [45] fractional derivative (3) By numerically integrating the above fractional derivative (Eq. 3) using the L1 scheme [46], we approximate the fractional derivative of order, where, (4) where, and is the value of time such that. For all simulations, we use the time step ms. By combining the right sides of Eqs. 1 and 4, and solving for at time () that depends on all past values of (from to), we obtain (5) where we define the Markov term weighted by the gamma function as (6) and the voltage-memory trace as (7) The voltage-memory trace (Eq. 7) can be further divided into the differentiation of past voltage () (Eq. 8) weighted by a function that depends on (Eq. 9): (8) (9) where is the time counter for past events and the positive integer corresponds to the value of time (=) at which the voltage is integrated (Eq. 5). The voltage-memory trace contains information of all the previous voltage activity of the neuron. Clearly, this is a computationally intensive problem due to the expanding matrix over time. We integrate this equation using our recently developed Fractional Integration Toolbox [47]. Since the fractional integration in Eq. 5 needs at least two inputs, is first integrated using the classical leaky integrate-and-fire model. In this section we show the spiking properties of the fractional leaky integrate-and-fire model and compare our results to experimental data, mainly from cortical pyramidal neurons. The voltage memory provides non-local dynamics that affects the spiking activity of the cell. The voltage-memory trace decays over time and is dependent on the value of. For the process is identical to the classical leaky integrator, while for values of the past trajectory activity increasingly contributes to the present value of the voltage. Since the weight is always positive, the sign of the voltage-memory trace depends on. With positive applied current the fractional model generates action potentials (Fig. 7A), with positive until the cell fires a spike and is reset (Fig. 7B). However, when there is a spike and voltage is reset, becomes negative. After the voltage escapes from the refractory period, the voltage-memory trace is positive until the next spike (Fig. 7C). As opposed to the classical leaky integrate-and-fire model, the memory trace accumulates over multiple spiking events, changes its dynamics, and thus it affects the ISI (Fig. 7D). The weight of the voltage-memory trace is determined by the fractional order. The weight is 0 for and it is always positive for. Fig. 7E shows the results of a simulation with standard parameters for 100 time steps (). The x-axis corresponds to the 1–100 temporal weights at t = N. The y-axis corresponds to the value of each weight. The increase in weight can be interpreted as the influence of the past state on the future state of the voltage. In a classical leaky integrator any past value is forgotten as a function of the time constant. In the fractional leaky integrator all the past values could continue to influence the future behavior of the system, particularly, for low values of. Fig. 7F illustrates this point by showing the dynamics of the weight of the initial condition () as a function of. The other important term that comes from the fractional derivative is the fractional coefficient which is a weight factor for the Markov process (Eq. 6). When becomes smaller, this function grows rapidly and it makes the effect of the input current on the voltage dynamics stronger (Fig. 7G). It is the combination of the weighted Markov process and the opposite effect of the memory trace that contribute to the long term spiking dynamics of the fractional model. Most of our results suggest that the value of has to be low in order to reproduce the spike timing adaptation observed experimentally. Although the fractional model has a continuous dependence on the power-law dynamics cause the effects to be nonlinear. For close to 1 the effects of the Markov term weighted by the gamma function dominate the dynamics. It is only when decreases that the voltage-memory trace can slow down the time evolution of the voltage. This is illustrated by plotting the value of the Markov term versus the memory trace for simulations in which we apply a step current (Fig. 8). When the memory trace is zero and the voltage only moves along the Markov term axis. As the value of decreases the voltage trajectory is deflected, taking longer to depolarize. The power-law dependency results that when the value of in that the memory trace dominates in the initial moments of the depolarization and slows down the dynamics (Fig. 8). When the input current is constant and, the sub-threshold voltage dynamics (Eq. 1) have an analytic solution which is the same as the solution of the classical integrate-and-fire model [73]. Similarly, for constant input current Langlands et al. [74] derived the analytic solution of the sub-threshold voltage for from the fractional cable equation. In that model the integration of the memory trace is restarted after every spike, thus wiping out the memory trace (see Methods). Such a system is capable of reproducing the delay to first-spike to constant input but then produces regular spiking (Fig. 9A). Our fractional model replicates the same result when we reset the memory trace after every spike (Fig. 9B). However, our model greatly differs from the analytical solution when taking into account the cumulative effect of the memory trace across multiple spiking cycles (Fig. 9C–D). Although both the analytic and full model with memory reset can capture short term memory, they do not produce spike adaptation. Hence, the full fractional model without any memory reset captures the multi-scale processes that spans the spiking activity of neurons.
  The fractional model can produce different degrees of adapting electrical activities by modifying the fractional exponent. A fractional order derivative captures the long-range correlations of the a system models that results in non-local dynamics. Fractional differential equations have been used in biological systems to capture the long-term memory effects of the dynamics [33], [74], [76]–[81]. For example, fractional order derivatives have been observed in the vestibular-ocular system [27], [82] and in the gating dynamics of ion channels [44], [83]. Power law statistical distributions of a neuronal response also exhibit fractional order dynamics [37]. The voltage dynamics in the fractional order model depend on both the Markov term (immediate past) and the voltage-memory trace that integrates all past voltage values. The behavior of the voltage-memory trace is similar to the behavior of the adaptive filter in the work of Pozzorini et al. [59], although in their work this filter is described as the sum of the spike-triggered current and a moving threshold. The voltage-memory trace also corresponds to the adaptation integral used in other works [18]. In this context our fractional order leaky integrate-and-fire model is a unified mathematical and computational framework that can be used to describe power-law dynamics and long-range correlations in neuronal activity. The fractional derivative can capture relationships between the distribution of conductances that can be complicated to model using explicit techniques. In biophysical systems, the voltage-memory trace might represent spike-triggered mechanisms that cause adaptation. For example, the voltage-memory trace might represent the slowly inactivating potassium-like current that induces the upward spike adaptation shown in Layer V pyramidal neurons of primary motor cortex [38], suggesting that these neurons are fractional differentiators. Alternatively, the voltage-memory trace can correspond to other adaptation currents such as calcium-activated after-hyperpolarization currents [5], slow sodium-channel inactivation currents [84], [85], or a combination of several adaptation currents. Many studies have used models with slow adaptation currents and exponential functions to analyze spike time adaptation [5], [54], [84]. Some of these models can produce similar properties found in fractional dynamics. For example, the Hodgkin-Huxley model with slow after hyper-polarization currents can produce multiple timescale adaptation processes [37]. The Generalized leaky integrate-and-fire model with an adaptive filter (GLIF-) produces spike adaptation with power-law dynamics [59]. Both the power-law dynamics and history-dependent properties of the GLIF- model correspond to that of the fractional leaky integrate-and-fire model. However, the fractional model provides a general way by simplifying the complicated details shown in other models. The fractional model exhibits spike adaptation with power-law dynamics by integrating all the past voltage values without any additional adaptation currents. Power law functions generalize the mechanism underlying exponential processes and are better alternatives to describe scale invariant spike adaptation [17], [18], [42], [59], [86], [87]. The fractional model shows new directions for studying spike adaptation using fractional derivatives and power-law dynamics instead of classical derivatives and exponential functions. The key parameter in our fractional model is. Experimental results have suggested that can be as small as 0. 15 for neocortical pyramidal neurons [37]. Our fitting of these data also resulted in a value of. Fits to the response of Layer 5 pyramidal motor cortex neurons resulted in a value of [38]. Thus, the experimental results and our modeling analysis suggest that the biophysically important values of are when. However, it is clear that if is much closer to 0, the system takes a longer time to generate spikes or never fires spikes at all, depending on the magnitude of the stimulus, so the feasible range of the fractional exponent might be between 0. 05 and 0. 2. The value of might correspond to the type, function, or location of specific neurons [38]. These regional differences are not exclusive to the cortex, for example, Purkinje cells in Lobule X and Lobules III – V show different degrees of spike adaptation [75]. Thus, different values of can be used to map the general voltage and spike time adaptation properties of neurons throughout the brain. Previous work has provided an experimental foundation to determine if a spiking neuron is a fractional differentiator [37]. Our framework provides a more general methodology to determine power-law neuronal dynamics from the sub-threshold to the spiking regime. At the sub-threshold level there are several measurements that can indicate that the membrane of the neuron follows a power-law. For example: In the same experiment a series of protocols can also be applied to determine if the spiking activity follows power-law dynamics. Some of these measurements are straightforward, others require longer recordings. For instance: Although power-laws are found at multiple scales of biological organization, their function and importance are still debated [88], [89]. In our work we propose that the membrane voltage of neurons can follow a power-law due to the emergent property of the combination of multiple active conductances. The value of the fractional derivative can be mapped to spike time adaptation dynamics taking place in multiple cell types across the brain. Computationally, a low value of results in spiking dynamics that are at the same time highly adaptable and reliable. Thus, neurons following power-law adaptation could have a large operational range while providing the reliability to filter out noisy signals and increase information capacity. The lack of a sub-threshold resonance frequency allows the neuron to filter signals homogeneously over a wide range of frequencies. In such a case, the fractional leaky integrate-and-fire model provides the basis to study the computational capacities and information processing properties of neurons showing high degree of spike time adaptation.
  The equations were coded and implemented using our recently developed fractional integration toolbox [47] and the simulation software package MATLAB [90]. The toolbox can be downloaded at www. utsa. edu/SantamariaLab. The parameters for all simulations were fixed and are described in Table 1. In order to compare to experiments we extracted the data from the referenced material using WebPlotDigitizer (http: //arohatgi. info/WebPlotDigitizer/). Then we imported the data points into Matlab. We then ran simulations varying the value of, usually between [0. 5,1. 0] at 0. 1 steps. We minimized the mean squared error between the data and the simulations, , where is the number of points. We determined the 95% confidence intervals by then varying around this minimum value of and calculated when it changed for an amount larger than 5% in either direction. In order to get only sub-threshold oscillations we used where Hz is a sigmoidal frequency function of time that varies from 0 Hz to 100 Hz in 10 seconds. The impedance as a function of frequency is defined as, where V is the membrane voltage, I is the input current, R is the resistance and X is the reactance. The absolute value of can be calculated using a fast Fourier transform and the phase as [57]. The simulations to study spike time reliability were generated by injecting a current nA, where is a Gaussian white noise with zero mean and standard deviation = 0. 03 nA. The stochastic input is filtered with an alpha function with time constant ms. The spike trains were obtained from trials, and the trail-to-trail variability of those N different responses were caused by the noise while and were fixed [68], [71]. In order to avoid initial condition effects we analyzed the spike trains of the last 5 s from 10 s simulations. The reliability measurements were computing using a correlation-based measure [68]–[71]. In brief, the spike trains obtained from N trails were smoothed with a Gaussian filter of width 3, and then pairwise correlated. The correlation-based measure reliability is defined as (10) where is the number of trials and the vectors are the filtered spike trains, filtered using ms. The values of range from 0 (lowest reliability) to 1 (highest reliability), and the reliability was computed for in the range [0. 1,1. 0]. For the quantification of the reliability in the coding of an embedded signal we injected the following current nA, where nA is the mean current which generates very low firing rate, nA is the standard deviation of the intrinsic noise, and is the standard deviation of the embedded signal also generated by a Gaussian noise and varies from 0 to 6 nA. The intrinsic and embedded signals are filtered with an alpha function where ms. We compared our fractional model to a previously developed analytical model of a fractional leaky integrate-and-fire [74]. Briefly, this analytical model is obtained with the following steps. Equation 1 can be converted to (11) where is the membrane time constant. By applying on both sides we obtained (12) Equation 12 is solved using the Fourier-Laplace transform (for details see [36], [74]) and the solution is given by (13) where is the Mittag-Leffler function [36], and for small times this function is approximated as (14) For the simulation of this model presented in our work we used the full Mittag-Leffler function instead of this approximation. In Eq. 13 the term with the Mittag-Leffler function (right side and right term) represents the memory trace. As in the classical integrate and fire model when the voltage reaches a spike is generated and V is reset to for a refractory period. The subthreshold voltage is integrated using Eq. 13 with initial voltage = and new initial time. During each integration cycle the Mittag-Leffler function restarts from 0 since the initial time reset to a new value. We call this memory reset. Thus, this model wipes out the memory trace after every spike, in contrast to our model that integrates the entire voltage trace. As a result, the inter-spike intervals of the spike train of the analytic solution are equal. If the approximation (Eq. 14) is used to simulate the voltage, the firing rate can also be approximated analytically by combining Eq. 13 and Eq. 14 (see also [74]). Let be the time when the voltage takes to increase from to and to fire. The time is given by [74] (15) Using the above the firing rate is approximated by (16) In the Results section we compare this analytical model with our model with memory reset (re-starting the memory trace after every spike) and with the full model (integrating the memory trace from the beginning of the simulation).
summary: Spike adaptation is a property of most neurons. When spike time adaptation occurs over multiple time scales, the dynamics can be described by a power-law. We study the computational properties of a leaky integrate-and-fire model with power-law adaptation. Instead of explicitly modeling the adaptation process by the contribution of slowly changing conductances, we use a fractional temporal derivative framework. The exponent of the fractional derivative represents the degree of adaptation of the membrane voltage, where 1 is the normal leaky integrator while values less than 1 produce increasing correlations in the voltage trace. The temporal correlation is interpreted as a memory trace that depends on the value of the fractional derivative. We identify the memory trace in the fractional model as the sum of the instantaneous differentiation weighted by a function that depends on the fractional exponent, and it provides non-local information to the incoming stimulus. The spiking dynamics of the fractional leaky integrate-and-fire model show memory dependence that can result in downward or upward spike adaptation. Our model provides a framework for understanding how long-range membrane voltage correlations affect spiking dynamics and information integration in neurons.
keywords: physics
  computational neuroscience
  single neuron function
  biophysics theory
  biology and life sciences
  physical sciences
  computational biology
  biophysics
  coding mechanisms
title: Neuronal Spike Timing Adaptation Described with a Fractional Leaky Integrate-and-Fire Model