Weakly nonlinear response of noisy neurons

We study the leaky integrate-and-fire (LIF) model in which the subthreshold voltage $v(t)$ obeys

$$\begin{eqnarray}&&\dot{v}=-v+\mu +\varepsilon s(t)+\xi (t).\end{eqnarray} \tag{ 1 }$$

Whenever $v(t)$ hits the threshold ${v}_{T}=1$, we register a spike time (figures 1(b), (c)) and set $v(t)$ to ${v}_{R}=0$ for an absolute refractory period τ. Time is measured in units of the membrane time constant and the voltage in multiples of the distance between reset voltage and threshold voltage (for details, see [23]). The voltage dynamics are subject to a constant input current μ and to white Gaussian noise $\xi (t)$ with zero mean and autocorrelation $\langle \xi (t)\xi (t+{t}^{\prime })\rangle =2D\delta ({t}^{\prime })$. The noise represents intrinsic fluctuations (e.g., channel noise) or external background input from other neurons in the network [24]. The strength of the driving signal $\varepsilon s(t)$ (figure 1(a)) is quantified by the small parameter $\varepsilon$. We consider the instantaneous firing rate $r(t)$, i.e., the probability density of a spike at time t. The rate can be obtained numerically by averaging N trials of the spike train (figure 1(c)). For this, we divide the fraction of trials that have fired in an interval $[t,t+{\rm{\Delta }}t]$ by the size of the interval ${\rm{\Delta }}t$ (depending on parameters, we used $N={10}^{8}\mbox{--}{10}^{9}$). The rate can also be determined analytically from a perturbation approach (see below). Already for small $\varepsilon$, the rate can exhibit a strongly nonlinear response to the signal as can be seen from the higher harmonic in figure 1(d).

Zoom In Zoom Out Reset image size

The voltage $v(t)$ in equation (1) is a stochastic process whose probability density, $P(v,t)$, obeys the Fokker–Planck equation [3]

$$\begin{eqnarray}&&{\partial }_{t}P={\hat{L}}_{0}P-\varepsilon s(t){\partial }_{v}P+r(t-\tau )\delta (v-{v}_{R}),\end{eqnarray} \tag{ 2 }$$

with the operator ${\hat{L}}_{0}={\partial }_{v}(v-\mu +D{\partial }_{v})$ and appropriate boundary condition (more details in the appendix). The solution of equation (2) in principle allows us to fully determine the density $P(V,t)$ and the time-dependent firing rate $r(t)$, but in practice it can only be solved via a perturbation approach.

We also perform numerical simulations of a model of excitable membrane dynamics that was proposed by Izhikevich [25] and is similar to the two-dimensional Morris–Lecar model [26]. Following reference [27] the model is endowed with channel noise, which is due to the finite number of potassium selective channels. The voltage dynamics of the Izhikevich model are given by

$$\begin{eqnarray}&&C\dot{v}={I}_{0}+\varepsilon s(t)-{I}_{L}(v)-{I}_{{Na}}(v)-{I}_{K}(v,n),\end{eqnarray} \tag{ 3 }$$

where I₀ is a constant input current, I_L is a passive leak current, I_Na is a deterministic 'persistent sodium' current, and I_K is a potassium current gated by the number of open potassium channels $0\leqslant n\leqslant N$. The number of open channels $n(t)$ comprises a time-dependent stochastic process with voltage-dependent per capita transition rates ${\rho }_{o}(v)$ for channel opening and ${\rho }_{c}(v)$ for channel closing. In our simulations, we consider a potassium current with a high threshold and standard parameters: Passive leak current ${I}_{L}={g}_{L}(v-{E}_{L})$ with ${g}_{L}\,=8\,{\rm{mS}}/{{\rm{cm}}}^{2}$ and ${E}_{L}=-80\,{\rm{mV}}$; 'persistent sodium' current ${I}_{{Na},p}={g}_{{Na}}{m}_{\infty }(v)(v-{E}_{{Na}})$ with ${g}_{{Na}}\,=20\,{\rm{mS}}/{{\rm{cm}}}^{2}$, ${E}_{{Na}}=60\,{\rm{mV}}$, and voltage-dependent activation ${m}_{\infty }(v)=1/\{1+{\rm{\exp }}[(-20\,{\rm{mV}}-v)/15\,{\rm{mV}}]\};$ potassium current ${I}_{K}={g}_{k}(n/N)(v-{E}_{K})$ with ${g}_{K}=10\,{\rm{mS}}/{{\rm{cm}}}^{2}$, ${E}_{K}=-90\,{\rm{mV}}$, open channel number $0\leqslant n\leqslant N$, and the total number of channels N = 100. Membrane capacitance is $C=1\,\mu {\rm{F}}/{{\rm{cm}}}^{2}$. The per capita transition rate for channel opening is ${\rho }_{o}(v)=1/\{1+{\rm{\exp }}[(-25\,{\rm{mV}}-v)/5\,{\rm{mV}}]\}$ and for closing is ${\rho }_{c}(v)=1-{\rho }_{o}(v)$.

For weak signals ($\varepsilon \ll 1$), the firing rate $r(t)$ can be approximated by the first terms of a Volterra series [28]

$$\begin{eqnarray}r(t) & = & {r}_{0}+\varepsilon \displaystyle \int {{dt}}_{1}^{\prime }\ {K}_{1}({t}_{1}^{\prime })s(t-{t}_{1}^{\prime })\\ & & +{\varepsilon }^{2}\,\,\displaystyle \int \,{{dt}}_{2}^{\prime }\ {{dt}}_{2}^{\prime\prime }\ {K}_{2}({t}_{2}^{\prime },{t}_{2}^{\prime\prime })s(t-{t}_{2}^{\prime })s(t-{t}_{2}^{\prime\prime })\,+\,{ \mathcal O }({\varepsilon }^{3}).\end{eqnarray} \tag{ 4 }$$

Here, ${r}_{0}={\left[\tau +\sqrt{\pi }{\int }_{(\mu -{v}_{R})/\sqrt{2D}}^{(\mu -{v}_{T})/\sqrt{2D}}{dx}{e}^{{x}^{2}}\mathrm{erfc}(x)\right]}^{-1}$ is the unperturbed firing rate (i.e., for $\varepsilon =0$), and K₁ and K₂ in equation (4) are the first- and second-order response kernels, respectively. In general, it is more convenient to use equation (4) in the frequency domain:

$$\begin{eqnarray}&&\tilde{r}(\omega )={r}_{0}\delta (\omega )+\epsilon {\chi }_{1}(\omega )\tilde{s}(\omega )+{\epsilon }^{2}\,\int \,d{\omega }^{\prime }\ {\chi }_{2}(\omega -{\omega }^{\prime },{\omega }^{\prime })\tilde{s}(\omega -{\omega }^{\prime })\tilde{s}({\omega }^{\prime })+{ \mathcal O }({\epsilon }^{3}),\end{eqnarray} \tag{ 5 }$$

where the response functions ${\chi }_{1}$ and ${\chi }_{2}$ are defined as the Fourier transforms of the respective kernels

$$\begin{eqnarray}&&{\chi }_{1}(\omega )=\int {dt}\ {e}^{i\omega t}{K}_{1}(t);\qquad {\chi }_{2}({\omega }_{1},{\omega }_{2})=\int \int {{dt}}_{1}{{dt}}_{2}\ {e}^{i{\omega }_{1}{t}_{1}}{e}^{i{\omega }_{2}{t}_{2}}{K}_{2}({t}_{1},{t}_{2}).\end{eqnarray} \tag{ 6 }$$

The linear response function ${\chi }_{1}(\omega )$ for the leaky integrate-and-fire model has been calculated in [8, 9]; in terms of the parabolic cylinder functions ${{ \mathcal D }}_{k}(x)$ [29], it reads [8]

$$\begin{eqnarray}&&{\chi }_{1}(\omega )=\displaystyle \frac{{r}_{0}i\omega /\sqrt{D}}{i\omega \,-\,1}\displaystyle \frac{{{ \mathcal D }}_{i\omega -1}\left(\tfrac{\mu -{v}_{T}}{\sqrt{D}}\right)-{e}^{{\rm{\Delta }}}{{ \mathcal D }}_{i\omega -1}\left(\tfrac{\mu -{v}_{R}}{\sqrt{D}}\right)}{{{ \mathcal D }}_{i\omega }\left(\tfrac{\mu -{v}_{T}}{\sqrt{D}}\right)-{e}^{{\rm{\Delta }}}{e}^{i\omega \tau }{{ \mathcal D }}_{i\omega }\left(\tfrac{\mu -{v}_{R}}{\sqrt{D}}\right)}\end{eqnarray} \tag{ 7 }$$

with ${\rm{\Delta }}=[{v}_{R}^{2}-{v}_{T}^{2}+2\mu ({v}_{T}-{v}_{R})]/(4D)$.

For an illustration of the second-order effects that are quantified by ${\chi }_{2}$, we consider the sum of two cosines as an input signal $s(t)=\alpha \cos ({\omega }_{1}t)+\beta \cos ({\omega }_{2}t+\phi )$, with a relative phase difference ϕ. Inserting this into equation (4), expressing the cosine functions by complex exponentials, identifying the integrals as the respective susceptibilities, and neglecting higher than second-order terms in $\varepsilon$ yields

$$\begin{eqnarray}&&r(t)\approx \underline{{r}_{0}\,+\,\displaystyle \frac{{\varepsilon }^{2}{\alpha }^{2}}{2}{\chi }_{2}({\omega }_{1},-{\omega }_{1})\,+\,\displaystyle \frac{{\varepsilon }^{2}{\beta }^{2}}{2}{\chi }_{2}({\omega }_{2},-{\omega }_{2})}\\ &&\,+\varepsilon {[\alpha | {\chi }_{1}({\omega }_{1})| \cos [{\omega }_{1}t-{\phi }_{1}({\omega }_{1})]+\beta | {\chi }_{1}({\omega }_{2})| \cos [{\omega }_{2}t+\phi -{\phi }_{1}({\omega }_{2})]]}_{{LR}}\\ &&\,+\displaystyle \frac{{\varepsilon }^{2}}{2}\phantom{\rule{}{3ex}}({\alpha }^{2}| {\chi }_{2}({\omega }_{1},{\omega }_{1})| \cos [2{\omega }_{1}t\,-\,{\phi }_{2}({\omega }_{1},{\omega }_{1})]\\ &&\,+{{\beta }^{2}| {\chi }_{2}({\omega }_{2},{\omega }_{2})| \cos [2{\omega }_{2}t+2\phi -{\phi }_{2}({\omega }_{2},{\omega }_{2})])}_{{HH}}\\ &&\,+{\varepsilon }^{2}\alpha \beta \phantom{\rule{}{3ex}}\{| {\chi }_{2}({\omega }_{1},{\omega }_{2})| \cos [({\omega }_{1}\,+\,{\omega }_{2})t\,+\,\phi \,-\,{\phi }_{2}({\omega }_{1},{\omega }_{2})]\\ &&\,+{| {\chi }_{2}({\omega }_{1},-{\omega }_{2})| \cos [({\omega }_{1}-{\omega }_{2})t-\phi -{\phi }_{2}({\omega }_{1},-{\omega }_{2})]\}}_{{MR}},\end{eqnarray} \tag{ 8 }$$

where ${\phi }_{\mathrm{1,2}}$ is the complex argument of the first- and second-order response function, respectively. The underlined terms represent the time-averaged firing rate that is affected by the cosine signals in the second order, quantified by ${\chi }_{2}(\omega ,-\omega )$. By time averaging, we mean integration over a sufficiently long time window and division by the time window. The time-averaged firing rate obtained in this way will contain only contributions from the time-independent terms in equation (8), while each of the time-dependent terms will average out to zero. The bracket ${[...]}_{{LR}}$ encloses the well-known linear response terms at the driving frequencies (the ground modes of the response) which are quantified by ${\chi }_{1}(\omega )$. Higher harmonics of each driving frequency (terms that oscillate twice as fast as the driving oscillations) are contained in ${(...)}_{{HH}}$ and are quantified by ${\chi }_{2}(\omega ,\omega )$. Last but not least, nonlinear interactions of the two signals result in a mixed response (terms in $\{...\}{}_{{MR}}$) at sums and differences of the two signal frequencies quantified by ${\chi }_{2}({\omega }_{1},\pm {\omega }_{2})$. Note that if the signal consists of two cosines with equal frequencies (${\omega }_{1}={\omega }_{2}$), the last term in $\{...\}{}_{{MR}}$ will be time-independent and will contribute to the time-averaged firing rate.

The first- and second-order response functions can be also determined from numerical simulations (here applied for the Izhikevich model) or from experimental measurements of neural spike trains in response to a sum of two cosine stimuli with $\alpha =\beta =1$, $\phi =0$, and ${\omega }_{1}\gt {\omega }_{2}$. To this end, we consider the Fourier transform of equation (8) (here only for $\omega \gt 0$)

$$\begin{eqnarray}\tilde{r}(\omega ) & \approx & +\displaystyle \frac{\varepsilon }{2}{[{\chi }_{1}({\omega }_{1})\delta (\omega -{\omega }_{1})+{\chi }_{1}({\omega }_{2})\delta (\omega -{\omega }_{2})]}_{{LR}}\\ & & +\displaystyle \frac{{\varepsilon }^{2}}{4}{({\chi }_{2}({\omega }_{1},{\omega }_{1})\delta (\omega -2{\omega }_{1})+{\chi }_{2}({\omega }_{2},{\omega }_{2})\delta (\omega -2{\omega }_{2}))}_{{HH}}\\ & & +\displaystyle \frac{{\varepsilon }^{2}}{2}{\{{\chi }_{2}({\omega }_{1},{\omega }_{2})\delta (\omega -({\omega }_{1}+{\omega }_{2}))+{\chi }_{2}({\omega }_{1},-{\omega }_{2})\delta (\omega -({\omega }_{1}-{\omega }_{2}))\}}_{{MR}}.\end{eqnarray} \tag{ 9 }$$

Now, in order to determine the response functions numerically, we first numerically estimate the time-dependent firing rate by simulating N realizations of a spike train and dividing the fraction of trials that have fired in an interval $[t,t+{\rm{\Delta }}t]$ by the size of the interval ${\rm{\Delta }}t$. Then, by performing a discrete Fourier transform of the estimated rate, we obtain $\tilde{r}$ and, finally, we determine the numerical estimates of the response functions by inverting equation (9):

$$\begin{eqnarray}&&{\chi }_{1}({\omega }_{1})\approx \displaystyle \frac{2{\rm{\Delta }}t}{\varepsilon }\tilde{r}({\omega }_{1}),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\chi }_{2}({\omega }_{1},{\omega }_{2})\approx \displaystyle \frac{2{\rm{\Delta }}t}{{\varepsilon }^{2}}\tilde{r}({\omega }_{1}+{\omega }_{2}),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\chi }_{2}({\omega }_{1},-{\omega }_{2})\approx \left\{\begin{array}{l}\displaystyle \frac{2{\rm{\Delta }}t}{{\varepsilon }^{2}}\tilde{r}({\omega }_{1}-{\omega }_{2})-\displaystyle \frac{{\chi }_{1}({\omega }_{2})}{\varepsilon }\quad {\rm{if}}\quad {\omega }_{1}=2{\omega }_{2},\\ \displaystyle \frac{2{\rm{\Delta }}t}{{\varepsilon }^{2}}\tilde{r}({\omega }_{1}-{\omega }_{2})\quad {\rm{else}},\end{array}\right.\end{eqnarray} \tag{ 12 }$$

where we neglected higher-order terms in $\varepsilon$. Employing the symmetry relations for ${\chi }_{2}({\omega }_{1},{\omega }_{2})$ that we explain below, we can extract the second-order response function for all frequency arguments.

Using a perturbation ansatz for $r(t)$ and $P(v,t)$ in equation (2), we obtain a hierarchy of equations (details in the appendix) that determine the second-order response of the LIF model (see [21] for the treatment of a similar problem). Solving this hierarchy with the appropriate boundary conditions, we arrive at the Fourier transform of the second-order response kernel

$$\begin{eqnarray}{\chi }_{2}({\omega }_{1},{\omega }_{2}) & = & \displaystyle \frac{{r}_{0}(1\,-\,i{\omega }_{1}\,-\,i{\omega }_{2})(i{\omega }_{1}\,+\,i{\omega }_{2})}{2D(i{\omega }_{1}\,-\,1)(i{\omega }_{2}\,-\,1)}\displaystyle \frac{{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}-2}\left(\tfrac{\mu -{v}_{T}}{\sqrt{D}}\right)-{e}^{{\rm{\Delta }}}{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}-2}\left(\tfrac{\mu -{v}_{R}}{\sqrt{D}}\right)\,}{{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}}\left(\tfrac{\mu -{v}_{T}}{\sqrt{D}}\right)-{e}^{{\rm{\Delta }}}{e}^{i({\omega }_{1}+{\omega }_{2})\tau }{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}}\left(\tfrac{\mu -{v}_{R}}{\sqrt{D}}\right)}\\ & & +\displaystyle \frac{i{\omega }_{1}+i{\omega }_{2}}{2\sqrt{D}}\displaystyle \frac{\left(\tfrac{{\chi }_{1}({\omega }_{1})}{i{\omega }_{2}-1}\,+\,\tfrac{{\chi }_{1}({\omega }_{2})}{i{\omega }_{1}-1}\right){{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}-1}\left(\tfrac{\mu -{v}_{T}}{\sqrt{D}}\right)}{{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}}\left(\tfrac{\mu -{v}_{T}}{\sqrt{D}}\right)-{e}^{{\rm{\Delta }}}{e}^{i({\omega }_{1}+{\omega }_{2})\tau }{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}}\left(\tfrac{\mu -{v}_{R}}{\sqrt{D}}\right)}\\ & & -\displaystyle \frac{i{\omega }_{1}+i{\omega }_{2}}{2\sqrt{D}}\displaystyle \frac{\left(\tfrac{{\chi }_{1}({\omega }_{1}){e}^{i{\omega }_{1}\tau }}{i{\omega }_{2}-1}\,+\,\tfrac{{\chi }_{1}({\omega }_{2}){e}^{i{\omega }_{2}\tau }}{i{\omega }_{1}-1}\right){e}^{{\rm{\Delta }}}{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}-1}\left(\tfrac{\mu -{v}_{R}}{\sqrt{D}}\right)}{{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}}\left(\tfrac{\mu -{v}_{T}}{\sqrt{D}}\right)-{e}^{{\rm{\Delta }}}{e}^{i({\omega }_{1}+{\omega }_{2})\tau }{{ \mathcal D }}_{i{\omega }_{1}+i{\omega }_{2}}\left(\tfrac{\mu -{v}_{R}}{\sqrt{D}}\right)}.\end{eqnarray} \tag{ 13 }$$

The function ${\chi }_{2}$ obeys the symmetries ${\chi }_{2}({\omega }_{1},{\omega }_{2})={\chi }_{2}({\omega }_{2},{\omega }_{1})$ and ${\chi }_{2}({\omega }_{1},{\omega }_{2})={\chi }_{2}^{* }(-{\omega }_{1},-{\omega }_{2})$ (where the asterisk denotes the complex conjugate), which implies in particular that ${\chi }_{2}(\omega ,-\omega )$ is real. These symmetries become apparent in figure 2, where we plot the absolute value and the complex argument for two distinct firing regimes of the LIF model [23]. Values along the white solid diagonal (in figures 2(a), (b), (c), (d)) quantify the second-order contribution of a cosine signal to the mean firing rate (see underlined terms in equation (8)), values along the white dashed diagonal quantify the second-order contributions of a single cosine to higher harmonics in the output (see terms contained in ${(...)}_{{HH}}$ in equation (8)), and the area off the white diagonals quantifies the nonlinear interactions of two cosine signals (see terms contained in ${(...)}_{{MR}}$ in equation (8)). Values of the second-order response functions along the black lines in figures 2(a) and (c) ($\ {\chi }_{2}(\omega ,0)={\chi }_{2}(0,\omega )\$) are proportional to the derivative of the linear response function w.r.t. the base current, $d{\chi }_{1}(\omega )/d\mu$ (see Appendix B for details). For a subthreshold base current $\mu \lt 1$, noise required to reach the threshold and spike trains are typically irregular [interspike intervals have a high coefficient of variation (CV)] whereas, for a suprathreshold base current $\mu \gt 1$, the model exhibits sustained firing even in the absence of noise and spike trains are more regular (the CV is low). In the subthreshold regime (${r}_{0}\approx 0.14$, ${CV}\approx 0.4$), the dependence of ${\chi }_{2}$ on the two frequencies is limited to two broad peaks that for our choice of parameters are located at ${f}_{1}={f}_{2}\approx \pm 0.2$ (figure 2(a), (b)). In the suprathreshold regime (${r}_{0}\approx 0.42$, ${CV}\approx 0.02$) the second-order response displays a rather different shape in the form of triangular structures (figures 2(c), (d)).

Zoom In Zoom Out Reset image size

Because the magnitude of ${\chi }_{2}$ is much larger in the suprathreshold regime than in the subthreshold regime (see figures 2(a), (c)), we mainly focus our subsequent analysis on the suprathreshold regime and proceed to discuss quantitatively the nonlinear effects of the sum of cosine signals on higher harmonics, the time-averaged firing rate, and the mixed response.

We compare the amplitude of the harmonic oscillation to the amplitude of the ground mode (orange circles and blue diamonds, respectively, in figure 2(e)) that is excited by a single cosine ($\alpha =1$, $\beta =0$ in equation (8)). The amplitude of the second harmonic exhibits one pronounced peak at ${f}_{1}\approx 0.42$ (at the unperturbed firing rate) and another at ${f}_{1}\approx 0.21$ (half the unperturbed firing rate), where it even slightly exceeds the amplitude of the ground mode, implying that for these frequencies the excitation of the harmonic oscillation is most pronounced. This is exactly the case that we illustrated in figure 1. The linear response theory is in this case particularly misleading because the higher harmonic is stronger than the ground mode. Note that the theory put forward by Ostojic and Brunel [30], which is based on the incorporation of a static nonlinearity (linear-nonlinear model), cannot capture the dominance of the higher harmonics for certain values of the driving frequency. In fact, for the parameter set used in figure 1, the Ostojic/Brunel theory yields practically the same curve as the linear response (not shown) because the stationary firing rate (i.e., the static nonlinearity) is essentially linear over the relevant range of signal values. The observed nonlinear effect in our theory and in the simulations arises from the true dynamic second-order response to signals with a time scale that is comparable to the time scale of the system (inverse firing rate).

The time-averaged firing rate does not change in linear response. Second-order effects for a single cosine signal ($\alpha =1$, $\beta =0$ in equation (8)) are illustrated in figure 3 for the two different firing regimes of the model.

Zoom In Zoom Out Reset image size

In the subthreshold regime ($\mu =0.9$), a signal with a low or intermediate frequency can only increase the time-averaged rate (figure 3(a)). High-frequency signals have practically no influence on the rate, as expected from the low-pass nature of the LIF model. A maximum is attained at an intermediate frequency that is larger than the firing rate and depends on the specific choice of μ and D.

In the suprathreshold regime ($\mu =1.1$), the signal can both increase and decrease the rate depending on the signal frequency (figure 3(b)). While fast signals do not change the time-averaged rate, as for $\mu =0.9$, a low-frequency signal ($f\to 0$) decreases the time-averaged rate in marked contrast to the subthreshold case. The different effects of slow periodic signals on the time-averaged rate can be understood from a Taylor expansion of the firing rate. For slow signals, the time-dependent firing rate can be written as the unperturbed firing rate with a modulated base current, ${r}_{0}(\mu +\varepsilon s(t))$. A Taylor expansion then yields

$$\begin{eqnarray}&&{r}_{0}(\mu +\varepsilon s(t))\approx {r}_{0}+\displaystyle \frac{{{dr}}_{0}}{d\mu }\varepsilon s(t)+\displaystyle \frac{1}{2}\displaystyle \frac{{d}^{2}{r}_{0}}{d{\mu }^{2}}{\varepsilon }^{2}s{(t)}^{2}.\end{eqnarray} \tag{ 14 }$$

For the time-averaged rate of an LIF driven by a slow cosine signal, this leads to

$$\begin{eqnarray}&&{\langle r(t)\rangle }_{T}\approx {r}_{0}+\displaystyle \frac{{\varepsilon }^{2}}{4}\displaystyle \frac{{d}^{2}{r}_{0}}{d{\mu }^{2}}.\end{eqnarray} \tag{ 15 }$$

Hence, the change of the time-averaged rate due to periodic low-frequency signals is related to the curvature of the sigmoidal rate function ${r}_{0}(\mu )$. The curvature of this rate function is positive for the subthreshold regime (at the foot of the sigmoid) and negative for the suprathreshold regime (close to the plateau of the sigmoid).

Another feature in the suprathreshold regime is a precursor of frequency locking, which has already been studied in the deterministic [31] and stochastic [32] versions of the LIF model. A signal with a frequency close to the model's unperturbed firing rate can entrain the model within a limited frequency range (inset in figure 3(b)). A signal with $f\gt {r}_{0}$ will speed up the firing, while a signal with $f\lt {r}_{0}$ will slow down the model. Hence, the second-order response kernel provides a link between the deterministic nonlinear behavior of the LIF model and its stochastic firing rate.

In order to study the mixed response, we now consider two cosines ($\alpha =\beta =1$ in equation (8)). For the subthreshold regime ($\mu =0.9$), the second-order response function has a pronounced peak on the white dashed diagonal (quantifying the amplitude of the higher harmonics) but only weak contributions off the white diagonals, indicating a weak mixed response (figure 2(a)). However, in the suprathreshold regime ($\mu =1.1$), the second-order response function exhibits pronounced peaks and stripes off the white diagonals, indicating a strong mixed response (figure 2(c)).

In the latter regime, we first stimulate the neuron by the sum of two cosines (figures 3(c), (d)) and record the rate response (red noisy curve in figure 3(e)). The two signal frequencies are chosen such that they lie on the off-diagonal stripe in figure 2(c), satisfying ${f}_{1}+{f}_{2}\approx {r}_{0}$. In order to extract the significance of the mixed response (terms indicated by $\{...\}{}_{{MR}}$ in equation (8)), we then stimulate the neuron by the two cosines separately. The grey noisy curve in figure 3(e), obtained by summing up the respective responses and subtracting the mean firing rate, r₀, is well represented by the linear theory (keeping only zero and first-order terms in equation (8)), indicating that the harmonics induced by the cosines are negligible. The striking difference between grey and red lines shows that the response to the sum of two cosines can be strongly nonlinear even in a regime where the responses to the single cosines are linear. Linear response theory in this case (dashed line) fails to grasp the modulation of the firing rate (red noisy line) both with respect to its overall amplitude as well as to the timing of local minima and maxima of the rate. The weakly nonlinear theory (equation (8), solid line) in contrast captures both these aspects quite well. The mixed response is the strongest second-order effect we observe in the weakly nonlinear regime of the LIF model.

In order to test the robustness of the results obtained for the LIF neuron model, we performed numerical simulations for a biophysically more complex neuron model with ion channel noise, equation (3), that is driven by a cosine stimulus (figure 4(a)). In contrast to the integrate-and-fire neuron, the membrane voltage (figure 4(b)) of the Izhikevich model exhibits excitability (i.e., it generates spikes) and does not rely on a fire-and-reset rule. As for the LIF model, we find that for a suprathreshold mean current (${I}_{0}\gt 4.51$), the firing rate of the Izhikevich model exhibits harmonic oscillations (figure 4(d)), i.e., the response $\tilde{r}$ at $\omega =2{\omega }_{1}$ is equally strong or even stronger than the response at the driving frequency itself ($\omega ={\omega }_{1}$).

Zoom In Zoom Out Reset image size

The numerically measured second-order response function (figures 4(e)–(g)) shows a striking resemblance with the second-order response function that we calculated analytically for the LIF model (figure 2), indicating a high robustness of the weakly nonlinear theory for the response of spiking neurons. In particular, we observe the following similarities: the triangular structure of the second-order response amplitude (with peaks located at half the firing rate) in the suprathreshold case, the broadly peaked structure in the subthreshold case, and the overall dependence of the phase on the two frequency arguments in both regimes. However, there are also differences: (1) for the Izhikevich model, the amplitude is larger in the subthreshold than in the suprathreshold regime (the other way around for the LIF model); (2) the detailed behavior of the phase along the diagonal in the suprathreshold regime, which is almost constant for the Izhikevich model but decreases by π in the same range for the LIF model; (3) in the suprathreshold regime, the peak on the diagonal is larger than the peaks on the frequency axes. Points (1) and (2) could be due to the specific choice of parameters or could reflect more substantial differences in the two models. Note in particular that we used the same number of channels N for the two regimes for the Izhikevich model in figure 4 but different values of the noise intensity D for the LIF model in figure 2. The difference (3) could be caused by the finite-amplitude simulation from which the second-order response is estimated for the Izhikevich model. As we have verified also for the LIF model, a too-large signal amplitude $\varepsilon$ leads to an underestimation of the second-order response at the maxima on the frequency axes (not shown). For the computationally costly Izhikevich model, it is not possible to reduce $\varepsilon$ much more than the values we have used here.

Let us return to the LIF model and illustrate under which circumstances the weakly nonlinear response becomes manifest in the activity of spiking neurons.

First, we consider an LIF neuron that is driven by a cosine stimulus of increasing amplitude (figure 5(a)) and is subject to a background noise of constant strength (figure 5(b)). (Because we cannot directly plot the white noise $\xi (t)$ (it possesses an infinite variance), we use in figures 5(b) and (f) a filtered version of the noise, $\eta (t)={\int }_{t}^{t+{\rm{\Delta }}t}{{dt}}^{\prime }\sqrt{2D}\xi ({t}^{\prime })$.) For weak signal amplitudes ($\varepsilon \lt 0.08$), the firing rate of the neuron (red noisy curve in figure 5(e)) is well described by the linear response theory (dashed black line). For large signal amplitudes ($\varepsilon \gt 0.15$), the firing of the neuron shows pronounced periods of silence which become apparent from the raster plot (figure 5(c)). In this regime, which we call stochastic mode locking, the neuron fires twice within each period of the driving signal (2:1 mode locking) and noise only shifts the spike times but does not induce any skipping of firing. The transition ($0.08\lt \varepsilon \lt 0.15$) between the regimes of linear response and stochastic mode locking is well described by the second-order theory (black solid line in figure 5(d)). For $\varepsilon \gt 0.15$, the second-order theory still predicts the position of the peaks in the firing rate quite well but attains unphysiological negative values between the peaks (results not shown). The regions of the different regimes cannot be sharply defined but here are chosen to be determined by computing the relative mean squared error between theory and prediction (computed over one signal period), and terminating the theory where the error exceeds 10%.

Zoom In Zoom Out Reset image size

Another transition in which the weakly nonlinear response becomes manifest is depicted in figures 5(e)–(h). Here, we consider an LIF neuron that is driven by a cosine stimulus of constant amplitude (figure 5(e)) but where the strength of the intrinsic noise decreases slowly with time (figure 5(f)). While the rate response is well described by the linear theory for large noise strengths (dashed line in figure 5(h)) and the activity exhibits 2:1 mode locking for low noise strengths, the transition between the two regimes is well captured by the second-order theory (black solid line in figure 5(h)). The transition from the linear to the nonlinear regime controlled by noise is consistent with the known linearizing effect of fluctuations [33].