Modeling multiple time scale firing rate adaptation in a neural network of local field potentials

Abstract

In response to stimulus changes, the firing rates of many neurons adapt, such that stimulus change is emphasized. Previous work has emphasized that rate adaptation can span a wide range of time scales and produce time scale invariant power law adaptation. However, neuronal rate adaptation is typically modeled using single time scale dynamics, and constructing a conductance-based model with arbitrary adaptation dynamics is nontrivial. Here, a modeling approach is developed in which firing rate adaptation, or spike frequency adaptation, can be understood as a filtering of slow stimulus statistics. Adaptation dynamics are modeled by a stimulus filter, and quantified by measuring the phase leads of the firing rate in response to varying input frequencies. Arbitrary adaptation dynamics are approximated by a set of weighted exponentials with parameters obtained by fitting to a desired filter. With this approach it is straightforward to assess the effect of multiple time scale adaptation dynamics on neural networks. To demonstrate this, single time scale and power law adaptation were added to a network model of local field potentials. Rate adaptation enhanced the slow oscillations of the network and flattened the output power spectrum, dampening intrinsic network frequencies. Thus, rate adaptation may play an important role in network dynamics.

Appendix

1.1 Linear rate models with exponential rate adaptation

Rate models are often of the form r = f(x-a), where f is some function of the input x with a feedback current a that depends on the firing rate r via some equation a(r) describing adaptation dynamics (Benda and Herz 2003, La Camera et al. 2004). In one of the simplest such models the firing rate is a linear function of the stimulus and the negative feedback adaptation. The firing rate increases linearly with increasing stimulus, while the adaptation variable decays exponentially with time:

$$ \begin{array}{l}r=mx-ga\\ {}\frac{da}{dt}=\frac{-a}{\tau }+kr\end{array} $$

(13)

where r is the firing rate, x is the time-varying stimulus, a is the adaptation variable, τ is the relaxation time constant, and m, g, and k are gain constants. One way to see that this is a high pass system is to express the rate dependence in terms of only the stimulus, thus eliminating the adaptation variable a, and then examine the equation in the frequency domain. Expressing r(t) as only a function of x(t):

$$ \begin{array}{l}\frac{d}{dt}\left[\frac{m}{g}x-\frac{r}{g}\right]=\frac{-1}{\tau}\left[\frac{mx}{g}-\frac{r}{g}\right]+kr\\ {}\frac{dr}{dt}=-r\left(\frac{1}{\tau }+gk\right)+\frac{m}{\tau }x+m\frac{dx}{dt}\end{array} $$

(14)

With τ _eff   = (1/τ + gk) ⁻¹ this becomes:

$$ {\tau}_{eff}\frac{dr}{dt}=-r+m\left(1-gk{\tau}_{eff}\right)x+m{\tau}_{eff}\frac{dx}{dt} $$

(15)

which gives in the frequency domain:

$$ R=m\left(\frac{1-gk{\tau}_{eff}+i\omega {\tau}_{eff}}{1+i\omega {\tau}_{eff}}\right)X=m\left(1-\frac{gk{\tau}_{eff}}{1+i\omega {\tau}_{eff}}\right)X $$

(16)

and in the time domain:

$$ r(t)=m\left(x(t)-gkx(t)\ast {e}^{-t/{\tau}_{eff}}\right) $$

(17)

where 0 ≤ gkτ _eff ≤1. When gkτ _eff =1, it is clear that the system is that of adapting high pass filtering. In fact, as gkτ _eff increases, the amount of adaptation increases and the response to a step increase reflects a decay time course that is increasingly exponential-like. Notice that τ _eff   ≤ τ, as has been previously shown (Ermentrout 1998, Wang 1998, Tripp and Eliasmith 2010). The transfer function of Eq. (16) is identical to that derived in Benda and Herz (see Eq 5.8, 2003), where the gain parameters were expressed as slopes of the firing rate-input curve. The filter is:

$$ H\left(\omega \right)=m\frac{m_{ss}/m+i\omega {\tau}_{eff}}{1+i\omega {\tau}_{eff}} $$

(18)

where m and m _SS are the slopes of the firing rate-input curve at an initial unadapted time and an adapted steady state, respectively. Here, m _SS   = m(1-gkτ _eff ) and since τ _eff   = (1/τ + gk) ⁻¹ one can see that τ _eff   = τ m _SS /m. The degree of negative feedback provided by adaptation is expressed in the reduction of τ _eff from τ, the more adaptation the smaller τ _eff . Thus, the model of Eq. (16) is functionally equivalent to the universal exponential adaptation model derived in Benda and Herz (2003) from a standard spike-generating model of membrane potential and currents. Here, linear dynamics were assumed initially rather than as an approximation during the derivation. By claiming the adaptation variable a can be described as a function of the firing rate r, rather than of individual spikes, one is assuming that τ is much greater than 1/r. In order to separate the fast spike generator from slow adaptation, one assumes that fluctuations in a do not markedly affect the time course of the spike generator (Benda and Herz 2003).

1.2 Linear rate models with power law rate adaptation

One notable function that can be approximated by exponentials is the power law. This unique function has the special property of being scale invariant such that its shape is unchanged despite scaling the x-axis. In other words, power laws do not have a characteristic time scale. Adaptation in neocortical neurons has been found to perform a function that can be approximated by power law filters (Lundstrom et al. 2008a, Lundstrom et al. 2010). Power laws can be approximated by an infinite sum of exponentials as can be seen by using the definition of the gamma function (Thorson and Biederman-Thorson 1974, Fairhall and Bialek 2002):

$$ \varGamma (k)\equiv {\displaystyle {\int}_0^{\infty }{\lambda}^{k-1}{e}^{-\lambda }d\lambda }. $$

Further, the Fourier transform of a power law is obtained by using the gamma function and setting λ = iωt and dλ = iω dt:

$$ {\left(i\omega \right)}^{-k}=\frac{1}{\varGamma (k)}{\displaystyle {\int}_0^{\infty }{t}^{k-1}{e}^{-i\omega t}dt} $$

Thus, by using the definition of the Fourier transform, the Fourier transform of a power law is:

$$ {\left(i\omega \right)}^{\alpha}\overset{F}{\leftrightarrow }k{t}^{-\left(\alpha +1\right)} $$

(19)

where in the time domain t ≤0 = 0. Experimentally, (iω) ^α is a reasonable filter that describes the effect of adaptation on the stimulus (Lundstrom et al. 2008a, Lundstrom et al. 2010) with a model of adaptation as follows, shown below in the frequency domain:

$$ R\left(\omega \right)=k{\left(i\omega \right)}^{\alpha }X\left(\omega \right)+b $$

(20)

where k and b are constants and α describes degree of adaptation. In the time domain, the equivalent filter for H(ω) = iω ^α with a finite number of frequencies can be expressed as

$$ h(t)=\frac{1}{\varGamma \left(1-\alpha \right)}\frac{d}{dt}{\displaystyle {\int}_0^t{\left(t-\tau \right)}^{-\alpha }}\frac{ \sin nt}{\pi t}d\tau, \kern0.75em \left(0\le \alpha <1\right) $$

where the following approximation for the delta function with finite range of frequencies is used (Arfken and Weber 1995):

$$ {\delta}_n(t)=\frac{ \sin nt}{\pi t}=\frac{1}{2\pi }{\displaystyle {\int}_{-n}^n \exp (iwt)dw}. $$

However, this form of the filter is difficult to use. An easier-to-use filter can be found by using the approach described above, where resulting dynamics are such that either the adaptation variable or the response to a step impulse has the form of a power law. For example, assuming power law dynamics of the adaptation variable with h(t) = kt ^-α, the appropriate model would be:

$$ r(t)=mx(t)-mk{\displaystyle \int x\left(\tau \right){\left(t-\tau \right)}^{-\alpha }d\tau +b} $$

(21)

with rate reponse r to a time-varying stimulus x with constants m, b, and k, with α controlling the degree of adaptation. Alternately, if a model that has precisely a power law decay of power α to a step increase is desired, the derivative of the power law yields the result:

$$ r(t)=-m{\displaystyle \int x\left(\tau \right){\left(t-\tau \right)}^{-\alpha -1}d\tau }+b $$

(22)

with constants m and b and α controlling the degree of adaptation, as displayed in Fig. 7. Practically, discrete time power law filters as in Eqs. (21) and (22) are limited approximations for the frequency domain filter of (iω) ^α. They are not defined at time zero, and thus the initial response to a stimulus, which should not be dependent on any adaptation dynamics, may be discontinuous with the subsequent decay. In addition, a power law has an infinitely long tail, which can only be approximated by a finite length filter.

Fig. 7

a A power law filter h with α =0.2, according to Eq (22). The magnitude of the filter was normalized such that its first point was equal to 1, which is not displayed for convenience. b The response of h convolved with a step function. Note that the initial peak response is determined by the arbitrary first point of the filter. The power law does not have a defined steady state, and thus eventually decays to zero given an infinite filter length, in contrast to an exponential filter which decays to some steady state

1.3 The relationship between magnitude and phase in minimum phase systems

Typically there is no precise relationship between the frequency-response magnitude, or gain, for a linear time invariant system and its phase. However, for systems characterized by a rational response function there is a relationship and for a subset of rational systems termed minimum phase systems, specifying the phase determines the magnitude to within a single scale factor, and vice versa (Oppenheim et al. 1999). A rational LTI system that is causal, stable, and that in the Laplace domain has all zeros on the left side of the s-plane (or inside the unit circle of the z-plane for the z-transform) is a minimum phase system (Ulrych and Lasserre 1966, Oppenheim et al. 1999). Causality assumes that the filter is right-sided, that is, it is equal to zero for negative time points, meaning that the output cannot precede the input. Stability implies a bounded output sequence for every bounded input sequence, which suggests that the discrete time impulse response is absolutely summable:

$$ {\displaystyle \sum_{n=-\infty}^{\infty}\left|h\left[n\right]\right|}<\infty $$

(23)

The last requirement for a minimum phase system, which amounts to assuming that the inverse system H is also causal and stable, is equivalent to requiring all τ _n in the adaptation filter models to be nonnegative. Specifically, in the Laplace domain Eq. (10) has a numerator of the form:

$$ \left(1/{\tau}_1+s\right)\left(1/{\tau}_2+s\right)\ldots \left(1/{\tau}_n+s\right) $$

(24)

where zeros can be seen to be negative, that is, on the left side of the s-plane, as long as all τ_n are nonnegative. This implies that the magnitude and phase of the system are related through the Hilbert transfrom (Ulrych and Lasserre 1966, Oppenheim et al. 1999):

$$ \angle H\left(\omega \right)=-\frac{1}{\pi }PV{\displaystyle {\int}_{-\infty}^{\infty}\frac{ \log \left|H(x)\right|}{\omega -x}dx} $$

(25)

where ∠ H(ω) is the phase and PV signfies the principal value of the Cauchy integral of the Hilbert transform. The Hilbert transform can be understood as convolving a function with the filter 1/(π x), where the function in this case is the logarithm of the magnitude of H(ω). Approaching negative and positive infinity as values approach zero from the negative and positive side, respectively, 1/x is effectively a differentiating filter. Thus, from Eq. (25) the phase is related to the derivative of the magnitude. The phase is positive for high pass filters, meaning that rate adaptation gives rise to phase leads.

1.4 Implementing the Jansen and Rit model of EEG oscillations

The standard Jansen and Rit model (Jansen and Rit 1995) specified as a system differential equations was implemented with a fourth order Runge Kutta solver as:

• $$ {y}_1^{\prime }={y}_4 $$
• $$ {y}_2^{\prime }={y}_5 $$
• $$ {y}_3^{\prime }={y}_6 $$
• $$ {y}_4^{\prime }={k}_e/{\tau}_e\left(S+{c}_2Sgm\left({c}_1{y}_3\right)\right)-2/{\tau}_e{y}_4-1/{\tau}_e^2{y}_1 $$
• $$ {y}_5^{\prime }={k}_i/{\tau}_i{c}_4Sgm\left({c}_3{y}_3\right)-2/{\tau}_i{y}_5-1/{\tau}_i^2{y}_2 $$
• $$ {y}_6^{\prime }={k}_e/{\tau}_eSgm\left({y}_1-{y}_2\right)-2/{\tau}_e{y}_6-1/{\tau}_e^2{y}_3 $$

with sigmoid Sgm(v) = e ₀ / [1 + exp(0.56(6 - v), external stimulus (S), constants c _1–4  = [135 108 33.75 33.75], and parameters k _e =3.25 mV, k _i =22 mV, τ _e =10 ms, τ _i =20 ms, and e ₀   = 5 Hz, unless otherwise noted. Primes indicate temporal first-order derivatives.ODEs were solved by a fourth-order fixed step Runge Kutta solver with dt =1–5 ms, with identical results obtained regardless. The stimulus was either constant or a sine wave. To stimulate rate adaptation with three time scales, the following system was implemented:

• $$ {y}_1^{\prime }={y}_4 $$
• $$ {y}_2^{\prime }={y}_5 $$
• $$ {y}_3^{\prime }={y}_6 $$
• $$ {y}_4^{\prime }={y}_{10} $$
• $$ {y}_5^{\prime }={y}_{11} $$
• $$ {y}_6^{\prime }={k}_e/{\tau}_eSgm\left({y}_7\right)-2/{\tau}_e{y}_6-1/{\tau}_e^2{y}_3 $$
• $$ {y}_7^{\prime }={y}_8 $$
• $$ {y}_8^{\prime }={y}_9 $$
• $$ {y}_9^{\prime }=\left({y}_{10}^{\prime }-{y}_{11}^{\prime}\right)+{k}_1\left({y}_{10}-{y}_{11}\right)+{k}_2\left({y}_4-{y}_5\right)+{k}_3\left({y}_1-{y}_2\right)-{k}_4{y}_9-{k}_5{y}_8-{k}_6{y}_7 $$
• $$ {y}_{10}^{\prime }={k}_e/{\tau}_e\left(S^{\prime }+{c}_1{c}_2{y}_6Sgm^{\prime}\left({c}_1{y}_3\right)-2/{\tau}_e{y}_{10}-/{\tau}_e^2{y}_4\right) $$
• $$ {y}_{11}^{\prime }={k}_i/{\tau}_i{c}_3{c}_4{y}_6Sgm^{\prime}\left({c}_3{y}_3\right)-2/{\tau}_i{y}_{11}-1/{\tau}_i^2{y}_5 $$

with the following additional parameters:

• $$ {k}_1=1/{\tau}_1+1/{\tau}_2+1/{\tau}_3 $$
• $$ {k}_2=1/\left({\tau}_1+{\tau}_2\right)+1/\left({\tau}_1+{\tau}_3\right)+1/\left({\tau}_2+{\tau}_3\right) $$
• $$ {k}_3=1/\left({\tau}_1{\tau}_2{\tau}_3\right) $$
• $$ {k}_4={k}_1=k{g}_1+k{g}_2+k{g}_3 $$
• $$ {k}_5={k}_2+k{g}_1/\left({\tau}_2+{\tau}_3\right)+k{g}_2\left({\tau}_1+{\tau}_3\right)+k{g}_3/\left({\tau}_1+{\tau}_2\right) $$
• $$ {k}_6={k}_3+k{g}_1/\left({\tau}_2{\tau}_3\right)+k{g}_2\left({\tau}_1{\tau}_3\right)+k{g}_3/\left({\tau}_1{\tau}_2\right) $$

where kg _n and τ_n as parameters governing the gain and time constant of each exponential filter. Similar smaller systems of ODEs were simulated for the cases with one or two exponential filters with identical results as those obtained by setting the appropriate gain parameters to zero above. This model can also be simulated with the corresponding filter or integral equations rather than the above differential equations. The overall output of the model is the difference between the post-synaptic excitatory and inhibitory membrane potentials to the pyramidal neurons (y ₁ – y ₂ ), consistent with what is thought to primarily underlie the signal of electroencephalography (Nunez and Srinivasan 2006).