Weakly nonlinear response of noisy neurons

We calculate the instantaneous firing rate of a stochastic integrate-and-fire neuron driven by an arbitrary time-dependent signal up to second order in the signal amplitude. For cosine signals, this weakly nonlinear theory reveals: (i) a frequency-dependent change of the time-averaged firing rate reminiscent of frequency locking in deterministic oscillators; (ii) higher harmonics in the rate that may exceed the linear response; (iii) a strong nonlinear response to two cosines even when the response to a single cosine is linear. We also measure the second-order response numerically for a neuron model with excitable voltage dynamics and channel noise, and find a strong similarity to the second-order response that we obtain analytically for the leaky integrate-and-fire model. Finally, we illustrate how the transition of neural dynamics from the linear rate response regime to a mode-locking regime is captured by the second-order response. Our results highlight the importance and robustness of the weakly nonlinear regime in neural dynamics.


Introduction
Linear response theory is a successful tool for analyzing how a nonlinear dynamical system reacts to external perturbations [1]. It was originally developed in physics but has also been applied in various other fields [2,3]. The response to time-dependent signals is of particular importance for nerve cells (neurons) because the response properties of these highly nonlinear and noisy elements determine the way that our brain perceives sensory stimuli, makes decisions, and generates behavior. The linear response of neurons has been measured in experiments [4][5][6] and calculated for simple models [7][8][9]; it has been used to explore, for instance, transmission of fast signals [10,11] and the emergence of stable activity in recurrent networks [12].
Despite its success, linear response theory inherently cannot fully characterize the behavior of nonlinear systems and, consequently, extensions of the theory to second and higher orders have been studied, e.g., in nonlinear optics [13][14][15] and magnetic particle imaging [16]. In neuroscience, the second-order response of nerve cells has been measured experimentally [17][18][19][20], used in the weakly nonlinear analysis of network stability [21], and calculated for a simple Poisson neuron model [22]. However, a complete theory for the second-order response for a neuron model with explicit voltage dynamics and an understanding of the significance of the second-order response for biological neurons is still missing.
In this paper, we analytically derive the second-order response of a noisy integrate-and-fire neuron to a general time-dependent external signal. Calculated in this way, the rate yields interesting single-neuron statistics that can be compared with experimental data; it can be also regarded as the instantaneous population rate of a group of uncoupled neurons subject to the same external drive. Hence, our results are relevant for population coding of time-dependent signals. We illustrate the significance of the weakly nonlinear response by means of signals consisting of one or two cosine functions. By means of extensive numerical simulations, we verify that our results are valid for a broad range of noise intensities.

LIF model
We study the leaky integrate-and-fire (LIF) model in which the subthreshold voltage v t Whenever v t ( ) hits the threshold v 1 T = , we register a spike time (figures 1(b), (c)) and set v t ( ) to v 0 R = 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 t x ( ) with zero mean and autocorrelation t t t D t 2 x x d á + ¢ ñ = ¢ ( ) ( ) ( ). 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 s t e ( ) (figure 1(a)) is quantified by the small parameter e. 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 t , + D [ ]by the size of the interval t D (depending on parameters, we used N 10 10 8 9 = -). The rate can also be determined analytically from a perturbation approach (see below). Already for small e, the rate can exhibit a strongly nonlinear response to the signal as can be seen from the higher harmonic in figure 1(d).
The voltage v t ( ) in equation (1) is a stochastic process whose probability density, P v t , ( ), obeys the Fokker-Planck equation [3] ) 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 timedependent firing rate r t ( ), but in practice it can only be solved via a perturbation approach.

Izhikevich model with ion channel noise
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 where I 0 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 n N 0   .  = (b) membrane potential; (c) spike times for different trials (those for (b) are shown in orange); (d) instantaneous firing rate from simulations (red), linear theory (dashed line), and nonlinear theory (solid black line) given in equation (8). Parameters: and the total number of channels N=100. Membrane capacitance is C 1 F cm 2 m = . The per capita transition rate for channel opening is v v

Weakly nonlinear response theory
For weak signals ( 1 e  ), the firing rate r t ( ) can be approximated by the first terms of a Volterra series [28] r t r dt K t s t t dt dt K t t s t t s t t , . 4 is the unperturbed firing rate (i.e., for 0 e = ), and K 1 and K 2 in equation (4) are the firstand second-order response kernels, respectively. In general, it is more convenient to use equation (4) in the frequency domain: where the response functions 1 c and 2 c are defined as the Fourier transforms of the respective kernels dt e K t dt dt e e K t t ; , The linear response function 1 c w ( )for the leaky integrate-and-fire model has been calculated in [8,9]; in terms of the parabolic cylinder functions x k  ( ) [29], it reads [8] r i D i e e e 1 7 For an illustration of the second-order effects that are quantified by 2 c , we consider the sum of two cosines as an input signal s t t t cos cos with a relative phase difference f. 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 e yields where 1,2 f is the complex argument of the firstand 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 , 2 c w w -( ). 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 timeindependent 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 1 c w ( ). Higher harmonics of each driving frequency (terms that oscillate twice as fast as the driving oscillations) are contained in ... HH ( ) and are quantified by , 2 c w w ( ). 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 , ). Note that if the signal consists of two cosines with equal frequencies ( 1 2 w w = ), the last term in ... MR { } will be time-independent and will contribute to the time-averaged firing rate.

Numerical measurement of the firstand second-order response functions
The firstand 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 Now, in order to determine the response functions numerically, we first numerically estimate the timedependent firing rate by simulating N realizations of a spike train and dividing the fraction of trials that have fired in an interval t t t , + D [ ]by the size of the interval t D . Then, by performing a discrete Fourier transform of the estimated rate, we obtain r and, finally, we determine the numerical estimates of the response functions by inverting equation (9): where we neglected higher-order terms in e. Employing the symmetry relations for , 2 1 2 c w w ( ) that we explain below, we can extract the second-order response function for all frequency arguments.

Weakly nonlinear response of the LIF model
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 The function 2 c obeys the symmetries , , where the asterisk denotes the complex conjugate), which implies in particular that , 2 c w w -( )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 secondorder 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) Because the magnitude of 2 c 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.

Excitation of harmonic oscillations
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 ( 1 a = , 0 b = in equation (8)). The amplitude of the second harmonic exhibits one pronounced peak at f 0.42 1 » (at the unperturbed firing rate) and another at f 0.21 1 » (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).

Change of the time-averaged firing rate
The time-averaged firing rate does not change in linear response. Second-order effects for a single cosine signal (8)) are illustrated in figure 3 for the two different firing regimes of the model.
In the subthreshold regime ( 0.9 m = ), 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 ( 1.1 m = ), 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 0.9 m = , a lowfrequency signal (f 0  ) decreases the time-averaged rate in marked contrast to the subthreshold case. The 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 m ( ). 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 r 0 > will speed up the firing, while a signal with f 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.

Interaction between two periodic signals
In order to study the mixed response, we now consider two cosines ( 1 a b = = in equation (8)). For the subthreshold regime ( 0.9 m = ), 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 ( 1.1 m = ), 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 f r 1 2 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 0 , 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.

Weakly nonlinear response of the Izhikevich model with ion channel noise
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-andreset rule. As for the LIF model, we find that for a suprathreshold mean current (I 4.51 0 > ), the firing rate of the Izhikevich model exhibits harmonic oscillations ( figure 4(d)), i.e., the response r at 2 1 w w = is equally strong or even stronger than the response at the driving frequency itself ( 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 e 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 e much more than the values we have used here.
3.6. Transition from linear response to stochastic mode locking via the weakly nonlinear regime 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)). ), 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 0.15 e < < ) 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 0.15 e > , 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%.
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 secondorder 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].

Discussion
In this study, we derived an analytical expression for the second-order response function of the LIF model and illustrated the significance of the weakly nonlinear regime. Combined with the linear theory, our result (equation (13)) predicts the firing rate of a neuron (or, equivalently, the population rate of uncoupled neurons) for an arbitrary signal up to the second order in the signal amplitude. Since the temporal modulation of the rate is thought to be an essential channel in neural information transmission [28], understanding and predicting it theoretically for a broad range of signal frequencies is crucial.
The response of noisy systems to weak periodic stimuli can be successfully described by linear response theory. In particular in computational neuroscience, this theory has been used to investigate, for example, the influence of noise on signal transmission [8,9,[34][35][36], the stability of the autonomous activity in recurrent neural networks [12], and the magnitude and speed of information transmission in neural populations [10,11]. However, in the case of strong periodic stimulation, linear response theory fails and tools from nonlinear dynamics are predominantly employed; in particular, measures of mode locking and synchronization [37]. For general nonlinear systems, some features of the transition between these extremes, linear response and mode locking, can be described by the theory of stochastic synchronization [38,39]. However, in the context of neural oscillators, the complete transition is analytically poorly understood.
Here, we described in detail the transition between the linear-response regime and the stochastic mode-locking regime for a noisy neuron. This transition emerges in two ways: either with increasing signal amplitude or by  (8), where we updated the parameter e in every time step; for the theoretical prediction in figures 5(e)-(h), we also employed equation (8) and updated the parameter D in equations (7) and (13) in every time step.
decreasing the intrinsic noise of the neuron (both illustrated in figure 5). Both manifestations of this transition are well captured by the second-order response that we derived in this paper. We also demonstrated a number of striking features of the second-order response such as a higher harmonic that can exceed the ground mode (figure 1), or the highly nonlinear response to a sum of cosine signals ( figure 3(c)). Beyond the analytically tractable LIF model with white Gaussian current noise, we also measured numerically the second-order response function in a biophysically more realistic spiking neuron model with channel noise. We found the features of this response function to be very similar to those calculated for the LIF model. Given the differences in the noise model (discrete multiplicative temporally correlated noise instead of additive uncorrelated Gaussian noise) and the differences in the neural dynamics, this similarity is surprising. It suggests that the discussed features do not hinge on the particularities of the model and are, thus, general.
Many papers [4,[40][41][42][43][44][45] have studied the response of biological neurons to periodic stimulation in different firing regimes, and it is crucial to find a theoretical description that can capture the regimes' response properties. We expect the weakly nonlinear regime to be relevant in experimental preparations where neurons are subject to weak noise, e.g., in vitro where a small-amplitude periodic current stimulation is applied and where channel noise is the main source of fluctuations. However, the weakly nonlinear response may also be relevant in situations where network fluctuations are weak, e.g., in networks with very low firing rates [46] and where network mechanisms lead to the emergence of either one or multiple global oscillations [47].
Possible extensions of the framework presented here are the application of the Richardson method [48] for an efficient numerical estimation of higher-order response functions, the computation of the second-order response for neural models with synaptic filtering [10,36], and the inclusion of finite synaptic potentials [49].

Acknowledgments
This work was supported by the BMBF (FKZ: 01GQ1001A) and the DFG (research training group GRK 1589/2).

Appendix A. Calculation of the second-order response function
Here, we present the calculation of the second-order response function for an LIF neuron with white Gaussian noise. For completeness, we also sketch the first-order response (see [50] for a detailed presentation or [36] for an alternative approach to the perturbation calculation). The nonlinear response for one cosine signal (an important but not the complete second-order response) was already calculated by Brunel and Hakim in [21].
The voltage v t ( ) in equation (1) is a stochastic process whose probability density, P v t , ( ), obeys the Fokker-Planck equation [3] ) and an absorbing boundary at the threshold which leads to the boundary conditions P v t , 0 where the prime denotes a partial derivative with respect to v (see [50]). Additionally, the density is continuous everywhere, integrable and normalized to 1. For the computation of the firstand second-order response functions in equations (4) and (8), we choose a signal that consists of a sum of cosines. Note that in general, the firstand second-order response functions describe the response to signals with arbitrary time dependence. The signal here (a sum of two cosines with different frequencies), however, is chosen because it is amenable to an analytic calculation and still completely probes the second-order response properties of the neuronal dynamics (i.e., it allows computation of the response function 2 c for all possible combinations of frequency arguments). This is also the main difference to a similar calculation performed for the second-order response function [21] where, due to a different choice of probing signal, the mixed response (see equation (8)) is not present and hence the second-order theory is incomplete. With s t t t cos cos where c c . . denotes the complex conjugate. In analogy to equation (A.2), we make the following ansatz for the density P v t , ( ): Inserting equations (A.2) and (A.3) into equation (A.1) and sorting out the terms with the corresponding time dependence, we obtain the following hierarchy of equations:  where we again dropped the frequency arguments for the sake of notation and consider the homogeneous equation A.13 are parabolic cylinder functions [29] and linear independent solutions of  First, for the sake of simplicity, we employ the notation f v e ,