Firing-rate models for neurons with a broad repertoire of spiking behaviors

Capturing the response behavior of spiking neuron models with rate-based models facilitates the investigation of neuronal networks using powerful methods for rate-based network dynamics. To this end, we investigate the responses of two widely used neuron model types, the Izhikevich and augmented multi-adapative threshold (AMAT) models, to a range of spiking inputs ranging from step responses to natural spike data. We find (i) that linear-nonlinear firing rate models fitted to test data can be used to describe the firing-rate responses of AMAT and Izhikevich spiking neuron models in many cases; (ii) that firing-rate responses are generally too complex to be captured by first-order low-pass filters but require bandpass filters instead; (iii) that linear-nonlinear models capture the response of AMAT models better than of Izhikevich models; (iv) that the wide range of response types evoked by current-injection experiments collapses to few response types when neurons are driven by stationary or sinusoidally modulated Poisson input; and (v) that AMAT and Izhikevich models show different responses to spike input despite identical responses to current injections. Together, these findings suggest that rate-based models of network dynamics may capture a wider range of neuronal response properties by incorporating second-order bandpass filters fitted to responses of spiking model neurons. These models may contribute to bringing rate-based network modeling closer to the reality of biological neuronal networks.


Introduction
The simulation of large networks of spiking neurons on the scale of cortical columns or even whole areas of the cortex has become feasible due to advances in computer  (Helias et al. 2012;Kunkel et al. 2014). In order to relate simulation results to experimental findings, it is important to employ neuron models that accurately capture actual neuron dynamics in response to realistic stimuli. Dynamical models that reproduce the responses of individual neurons to injected currents go back to the seminal work by Hodgkin and Huxley (1952). Their conductance-based model quantitatively described the action potential initiation and propagation in the squid giant axon in response to depolarizing currents and spawned many variants and simplifications that have been analyzed and used in computational neuroscience ever since. Examples are the FitzHugh (1961) and the Morris-Lecar model (Morris and Lecar 1981). On the more abstract side of neuron modeling, Lapicque's neuron model (Lapicque 1907), widely known as the leaky integrate-and-fire (IAF) neuron, models the membrane potential V (t) as a passive current integrator with leak current, emitting a spike whenever V (t) reaches a threshold value θ , followed by a membrane potential reset (Tuckwell 1988;Burkitt 2006a, b).
These simple integrate-and-fire neuron models have particular appeal to computational neuroscientists because they capture the essential function of a neuron, while still being amenable to mathematical analysis in many input and network scenarios.
Yet, the ideal model would be a neuron model that is both simple in its dynamical equations and still captures most of the actual response dynamics of a real neuron to a wide range of stimuli. To this end, Izhikevich suggested a twodimensional neuron model that is able to reproduce at least twenty different characteristic responses that are commonly used to classify neuron response types in experiments, such as tonic, phasic and rebound spiking and bursting, or adaptation (Izhikevich 2010). The response types are illustrated in Fig. 1. The stimuli used to induce these spiking behaviors are direct current injections, ramp current injections or short direct current steps or pulses as indicated at the bottom of all panels.
In a network context, however, neurons usually receive noisy input currents. Moreover, they are known to respond highly reliably to repeated injections of the same frozen noise injection, while responses vary widely across trials when neurons receive identical direct current (Mainen and Sejnowski 1996). Neurons thus respond stereotypically to certain temporal input features rather than to mere current amplitude.
Motivated by such findings, Gerstner and colleagues showed that nonlinear IAF models, including the spikeresponse model and the adaptive exponential IAF model, can succesfully be mapped to experimental spike data in a noisy input regime and even have good spiketime prediction power (Brette and Gerstner 2005;Jolivet et al. 2006). Yet, the nonlinearity and the number of parameters in general make fitting a difficult task. The International Competition on Quantitative Neuron Modeling has challenged modelers to fit their neuron models to a set of spike data recorded from neurons stimulated with noisy input currents (Jolivet et al. 2008).
The resulting neuron models were tested with a noisy input current that was not included in the training set, and the predicted spike times were compared to those of the actually emitted spikes. The multi-timescale adaptive threshold model (MAT model) introduced by Kobayashi et al. (2009), a surprisingly simple model with linear subthreshold dynamics, solved this task best. Despite its simplicity, the MAT model can generate type-I and type-II excitability, as well as burst firing. Moreover, an extended version of the MAT model, the augmented MAT (AMAT) model, which incorporates threshold dynamics that depend on the membrane-potential history, is able to reproduce all twenty spike response patterns described for the Izhikevich model (Yamauchi et al. 2011). Because of its few parameters and simple dynamics, the AMAT model has low computational cost while providing a large dynamical repertoire, and is thus highly attractive for large-scale network simulations.
In an actual neuronal network, neurons typically integrate spikes from thousands of presynaptic neurons, yet not all spikes might necessarily have a strong impact on the membrane potential. In many spiking network models, the effect of individual spikes on the membrane potential is assumed to be small, and spiking activity asynchronous and irregular. In this limit it is indeed possible to substitute the input current by, e.g., Gaussian white noise or an Ornstein-Uhlenbeck process (Johannesma 1968). However, experimental findings have repeatedly demonstrated that, even though most synapses are weak, synaptic weight distributions typically have heavy tails, with some corresponding to post-synaptic potentials of up to 10 mV (Song et al. 2005;Lefort et al. 2009;Avermann et al. 2012;Ikegaya et al. 2013). It is thus important to extend the analysis of neuronal response dynamics to input spike trains that elicit large individual post-synaptic potentials.
At an even higher level of abstraction are models that ignore specific spike times and heterogeneities in network structure, i.e., rate and field models. In contrast to high-dimensional networks of spiking neurons, such models are often easier to analyze mathematically due to their low dimensionality, and hence can offer insight into steady states of network activity and bifurcations that give rise to complex spatio-temporal phenomena, such as oscillatory dynamics, traveling waves or activity bump formation. Prominent examples are neural mass models, such as the Jansen-Rit model (Jansen and Rit 1995), and neural field models, such as the Wilson-Cowan model (Wilson and Cowan 1972), which include spatial interactions between neurons. In these models, the dynamics of large, possibly heterogeneous, populations of neurons are substituted by rate variables in a mean-field manner (Ermentrout 1998;Coombes 2005).
An important conceptual step in the derivation of these models is the substitution of the spiking activity of a neuron in response to a certain input current I (t) by an appropriate rate function 1 mapping the input history {I (s)|s ≤ t} to the response rate at time t. Common choices are abstract models such as threshold-linear or sigmoidal functions F (I (t)) depending only on the input current at time t. The threshold-linear form is often chosen because of mathematical convenience, but also because it mimics to first order the gain function of many individual neurons in experiments (Chance et al. 2002;Blomquist et al. 2009), while the sigmoidal also models the saturation at very high firing rates. Yet, parameters of the gain function such as time Response types for current input as defined by Izhikevich (2004). This illustration was created with our NEST implementation of the augmented MAT model (Yamauchi et al. 2011). Some model and stimulus parameters differ from those given by Yamauchi et al. (2011), see the Appendix. Membrane potential V m is shown in blue, the threshold V th in red, and the input current in green, while emitted spikes are shown as black bars; after Yamauchi et al. (2011). Subfigures are labeled as in Izhikevich (2004, Fig. 1) constants, activation thresholds, or slope are often chosen rather qualitatively, and it is uncertain how well they match single-neuron properties or biophysics. A first step towards a stringent comparison of spiking neuron network simulations with reduced neural mass or field models is to obtain an adequate quantitative expression for the neuronal gain function F (I (t)). It is hence of interest to understand if and how the activity of individual spiking neurons in response to arbitrary input currents can be described truthfully by a rate-model formulation. Several point neuron models are simple enough to allow for an analytical derivation of the gain function, assuming that input currents are Gaussian white noise, sinusoidally modulated input, or shot noise of a given structure (see, e.g., Gerstein and Mandelbrot 1964;Stein 1965;Brunel 2000;Brunel et al. 2001;Burkitt 2006b;Richardson 2007;Richardson and Swarbrick 2010;Roxin 2011;Ostojic and Brunel 2011). However, more complex nonlinear neuron models, such as the Izhikevich model or even the AMAT model, often render such analyses futile, especially in the presence of large-amplitude post-synaptic current events that are beyond the realm of perturbation-based theories. This holds to an even larger degree for the second step towards a stringent comparison of spiking network and neural field models, namely capturing the temporal response properties of the models. A thorough understanding of complex nonlinear models thus requires simulation studies.
We provide here an analysis of the response to spike train input of the models proposed by Izhikevich (2003b) and by Yamauchi et al. (2011), following the approach by Nordlie et al. (2010) and Heiberg et al. (2013). Both models actually represent an entire class of models that can be tuned to a wide range of reponses by adjusting model parameters. We will thus refer to the Izhikevich and AMAT model classes, respectively, when we refer to the set of equations and spike-generation rules, while we will call each of the approximately 20 different parameterization a model. Each of the two model classes comprises some 20 models.
In Section 3.1 we present how the different models respond to spike train input.
In Section 3.2, we present fits of a linear-nonlinear firingrate model to the spike responses of Izhikevich and AMAT models to stationary and temporally modulated stochastic spike trains across a range of input rates, synaptic weights, and modulation frequencies and amplitudes under different background noise regimes.
We group the different models according to the filter parameters obtained in Section 3.3, before we in Section 3.4 explore how well the linear-nonlinear rate models capture the response of their spiking counterparts to novel stimuli, such as steps in the input firing rate and more complex temporally modulated input.
Finally, in Section 3.5 we investigate whether we can generalize models fitted to a specific input regime to a broader set of stimuli, before we summarize our findings in Section 4.

Neuron models
We study response behavior for two neuron model classes, the Izhikevich model (Izhikevich 2003b) and the augmented MAT model (Yamauchi et al. 2011). As both model classes are well described in the original publications, we just summarize them briefly in Tables 1 and 2. Both models are able to reproduce 20 of the most prominent features of biological spiking neurons in response to injected current input as illustrated in Fig. 1. These response types were first summarized in tabular form by Izhikevich (2004); see also Markram et al. (2004). Tables 3 and 4 present the parameter values required to obtain the model responses displayed for each class; note that certain models only differ in the stimulus injected, while neuron parameters are identical. All models are implemented on a fixed time grid (dt = 0.1 ms).  Table 3 Table 2 Summary of AMAT model; for parameters, see Table 4 We integrate the Izhikevich model class using the forward Euler algorithm as in the original publications on the model. Izhikevich (2003b) used a 1 ms time step, but splitting the update of the membrane potential (but not the recovery variable) into two steps of 0.5 ms "for numerical stability". Figure 1 of Izhikevich (2004), on the other hand, was generated using different time steps for different cases, ranging from 0.1 ms to 0.5 ms without substepping, as evidenced by the source code used to generate that figure (Izhikevich 2003a). We extracted model parameters as shown in Table 3 from that source code, including an external current I ext injected into the model for some variants in addition to the stimulus current.
Izhikevich's source code also revealed that model variants G, L, and R use other equations than Eqs. (1)-(2) for V (t) or U(t). We therefore excluded these variants from our study. We also excluded variants I and O, since they have the same parametes as variants A and M, respectively, and differ only in the test stimulus injected to create Fig. 1 of Izhikevich (2004). Labels refer to subfigure labels in Izhikevich (2004, Fig. 1). Models A and I, G and L, and M and O, respectively, share the same parameters and differ only in their input parameters. Instances marked with an asterisk were not included in the study due to repeated parameters sets, non-standard model equations or pathological behavior; see text for details. Common parameters: V th = 30 mV, V (t = 0) ∼ U(−70 mV, 30 mV)   Furthermore, we observed that response patterns depend on the precise time step used. In particular, the response for case T, Inhibition-induced bursting, is unstable for time steps shorter than 0.5 ms. We therefore also excluded case T from our analysis.
The Izhikevich model class is not defined with consistent units in the original publication (Izhikevich 2003b). While a time unit of milliseconds is implied and membrane potential is specified in millivolts, no units are given for the parameters or explicit constants. The model equations imply that input currents have units of mV/ms, which is rather exotic. In the spirit of Izhikevich (2003b) we therefore treat all quantities except time and membrane potential as unitless for the Izhikevich model class.
The AMAT class is implemented in NEST as model amat2 psc exp using exact integration (Rotter and Diesmann 1999). The implementation follows the NEST convention of parameterizing the membrane potential equation Eq. (4) in terms of membrane time constant τ m and membrane capacitance C m and an explicit reversal potential E L , while Yamauchi et al. (2011) parameterize their Eq. 1 in terms of τ m and membrane resistance R and define E L = 0 mV. The parameterizations are related by C m = τ m /R and a shift of the membrane potential V and the resting value of the threshold ω by E L . Some parameter values were adjusted to be able to reproduce Figs. 6 and 7 in Yamauchi et al. (2011) as discussed in the Appendix. Model variants L and R are excluded from the study as they have identical parameters to variants A and O, respectively.
In all simulations reported here, a single neuron is stimulated with spike train input. For the Izhikevich model class, this spike input results in instantaneous jumps in the membrane potential v. For the AMAT class, each incoming spike evokes an exponentially decaying synaptic current. For details, see Tables 1 and 2 and Section 2.2.
Output spikes are recorded with NEST device spike detector.

Stimulation
We briefly summarize here the sinusoidal stimulation protocol and response characterization based on Nordlie et al. (2010) and presented in detail in Heiberg et al. (2013). More general stimulation protocols are described in Section 2.5.
( 7 ) Mean rates a 0 , modulation depth a 1 , and modulation frequency f stim are varied systematically; modulation depth is limited to 0 ≤ a 1 ≤ a 0 to avoid rectification. We used NEST device sinusoidal poisson generator to generate the input spike trains.
The weights w > 0 of the synapses transmitting the stimulus spike train a(t) are varied from about 10% to about 75% of the synaptic weight w θ required if a single incoming excitatory spike shall evoke a threshold crossing from rest. For the AMAT model class, w θ is the same for all model variants and we use weights between 100 pA and 900 pA in our experiments.
For the Izhikevich model class, in contrast, model parameters do influence the response to isolated spikes. We therefore define a weight factor ξ for each model variant as the smallest weight for which a single excitatory input spike triggers the spike initiation process. Synaptic weights w are set to fractions of this value, ranging from 0.1 to 0.75, i.e., within the same range as for the AMAT model.
In addition to the resulting current stimulus, I stim (t), we consider stationary noisy background input currents I bg (t), representing unspecific weak network input. This allows us to study neuronal responses to I stim (t) in different input scenarios. The full input a neuron receives is thus given by I (t) = I stim (t) + I bg (t). We characterize the background current by its mean μ bg and standard deviation σ bg .
The NEST implementation of the Izhikevich neuron model is equipped with instantaneous current-based synapses. Assuming high rates and small synaptic strength, balanced spiking input can be approximated well by Gaussian white noise. We thus inject approximate Gaussian Fig. 2 A model neuron is driven by a spike train with sinusoidally modulated rate a(t) with mean a 0 , modulation depth a 1 , and frequency f stim , cf. Eq. (7). As a first-order approximation, the output spike train of the neuron is characterized by the sinusoidally modulated response firing rate r(t) with mean r 0 , amplitude r 1 , frequency f stim and phase φ, cf. Eq. (10). Adapted from Nordlie et al. (2010), Fig. 1 white noise realizations of defined mean μ bg and standard deviation σ bg using NEST's noise_generator 2 .
The AMAT model, as used here, has current-based exponential synapses with characteristic time constants τ syn,E = 1ms and τ syn,I = 3ms. We inject background current as Poisson spike trains through synapses with small, fixed weight w E,bg = 1 pA and w I,bg = −4/3 pA, respectively, using NEST model poisson generator.
The resulting noise input current has mean and standard deviation For given μ bg and σ bg , we obtain noise input rates by solving Eqs. (8)- (9) for ν E and ν I . We consider three background current regimes: first the case without additional background current I bg (t) = 0 pA, where all spiking activity is purely stimulus induced. In the second case, I bg (t) is chosen such that μ bg = 0 pA, and σ bg is large enough to elicit spiking activity with background input alone, i.e., if I stim (t) = 0 pA. In the third case, we consider a net inhibitory background current, with μ bg < 0 pA and sufficient standard deviation σ bg to again elicit baseline spiking in absence of I stim (t). While the first scenario can be considered a typical situation for neurons in slice preparations, the latter two mimic the situation in vivo, e.g., in cortical layer II/III where ongoing spiking activity is sparse (see e.g., Sakata and Harris 2012; Petersen and Crochet 2013) and input currents are balanced or even inhibition dominated (Haider et al. 2013).

Sinusoidal rate model
We characterize the response of the neurons by a sinusoidal rate model as illustrated in Fig. 2. For a purely linear response, r 0 represents the background firing rate of the neuron, r 1 the stimulus response amplitude (with phase shift φ 1 ), and we expect r m = 0 for all higher harmonics (m ≥ 2). Any nonlinearities in the system will typically be associated with power in the higher harmonics. We consider power at harmonics as significant (z-test, 99% confidence level) if where B is the estimated background power of the spike train between the harmonics and the weighted standard deviation of the spike train power spectrum across frequencies. For details, see Section 2.2.2 of Heiberg et al. (2013).

Linearity
We proceed as follows to characterize the linearity of the firing-rate curve in response to stationary input: We obtain the firing-rate curve r 0 = f (a 0 ) for a given neuron model, noise regime and synaptic weight by measuring the output rate r 0 as a function of stationary input rate a 0 in the absence of modulation (a 1 = 0). To characterize the linearity of f over an interval [α, β], we define the linearity measurē as the normalized mean square difference between f (x) and (x), the best linear fit to is perfectly linear, we haveL 1 = 0, whileL 1 = 1 means that the average squared distance between firing-rate curve and linear fit is equal to the mean value over the interval. Larger values ofL 1 are difficult to interpret, though. We therefore define as linearity measure. L 1 = 1 indicates perfect linearity, L 1 = 1/2 a deviation from linearity equal to the mean value, while L 1 approaches 0 for large deviations. For a given mean rate a 0 and modulation depth a 1 , we evaluate linearity over [α, β] = [a 0 − a 1 , a 0 + a 1 ], i.e., the range of rates spanned by the temporally modulated input.

Rate model description
The response of a linear, time-invariant (LTI) system to any input can be calculated as a convolution of the input and the impulse response of the system. A wide class of non-linear systems can be described by a linear convolution of the input with a kernel h(t) followed by a non-linear activation function g(·), so that the response is given by To test how well this applies to the neuron models studied here, we fit linear-nonlinear firing-rate models to the responses of the spiking neuron models and compare firingrate predictions from the linear-nonlinear models to those of the fitted spiking models. We summarize the derivation of the firing-rate model below, based on Heiberg et al. (2013) and Nordlie et al. (2010).
For each neuron, we find the activation function g(·) and the kernel h(t). For constant input, a(t) = a 0 , the convolution becomes the identity operation, provided the kernel is normalized ( h(t)dt = 1). We determine g(·) by measuring the response to stationary input, r 0 = g(a 0 ) for a range of a 0 and obtain a continuous representation of g(·) by interpolation (linear B-spline).
To obtain the kernel h(t), we linearize the activation function around a given working point (a 0 , r 0 ). The response to a(t) = a 0 + a 1 s(t) can then be expressed as where the linear impulse response function combines the normalized kernel with the linear gain γ . In general, this approximation is only valid for small deviations from the working point. However, the limits are not known a priori. For brevity of notation, we will usually drop the explicit reference to stimulus parameters a 0 and a 1 below. We obtain the transfer function, i.e., the Fourier transform of the linear impulse response h 0 (t), from the model responses to sinusoidally modulated input (s(t) = sin 2πf stim t, cf. Eq. (7)) where r 1 (f stim ) and φ(f stim ) are the Fourier amplitude and phase of the response, respectively. In Nordlie et al. (2010) and Heiberg et al. (2013), firstorder low-pass filters with delay provided adequate fits to the empirical frequency responses. Here, more complex filter models are needed to fit additional response features. In particular, we expect a second filter time constant τ c = 1/(2πf c ) to be needed to model some of the response types illustrated in Fig. 1. We choose to combine low-and high-pass components of the filter as a sum: This form allows for a representation of the filter through a system of linear differential equations, see Section 2.4.1. 3 3 We also explored combining the terms in product form but did not observe significantly different results.
The filter kernels were fitted to the empirical transfer function to capture it with as few parameters as possible. For each set of stimulus parameters (a 0 , a 1 , w, μ, σ ), we obtained fits for the parameters γ 1 , f c,1 , γ 2 , f c,2 and . Fitting was performed using basin-hopping optimization provided by the SciPy Optimize toolbox with L-BFGS-B minimization (Jones et al. 2001). To avoid pathological solutions, we imposed the following constraints If a fit resulted in f c,1 > f c,2 , we swapped frequencies and gain coefficients so that f c,1 and f c,2 , respectively, are always the lower and upper characteristic frequencies of the filter. For each parameter combination, we performed 60 independent fits from different starting points and retained the best fit. We also performed 15 independent fits for a pure lowpass filter, but these never yielded better results than fits to the bandpass filter defined by Eq. (18). We now define the linear-nonlinear rate model as with the normalized kernel and correspondingly in Fourier space. We take the maximum solely to avoid negative rates that may result in rare cases from extrapolation of the activation function g(·). This kernel depends on the stimulus parameters a 0 and a 1 used to construct it, but we drop this dependence for clarity of notation; the actual range over which a kernel is useful is explored in Section 3.

Differential-equation representation
The filterH 0,SUM (f ) corresponds to a sum of low-pass filters in the time domain. For this model, the linearnonlinear model of Eq. (14) can be mapped to a set of delay differential equations using the linear chain trick (Nordbø et al. 2007). In particular, in the time domain the filter is given by with Heaviside step function (t). We introduce with u 1 (t) = (a * h 0,1 )(t) and u 2 (t) = (a * h 0,2 )(t) . (28) Straightforward differentiation of the two temporal kernels yieldṡ From this we can solve for u(t) = u 1 (t) + u 2 (t) and the full rate dynamics then follows by application of the nonlinearity g The advantage of the differential representation Eq. (29) lies in the fact that it is local in time, whereas the representation by convolution in general relies on knowledge of the complete history of the dynamics.

Tests against spike trains
We compare the response properties of our rate-based model against spiking neuron models as follows. We use synthetic (Section 2.5.1) or experimentally recorded (Section 2.5.2) spike trains S(t) as test input. Spiking neuron models are driven by these trains directly and their output spike trains R(t) are recorded as described in Section 2.6. We then use the fixed-kernel density estimation method by Shimazaki and Shinomoto (2010) with 0.05 ms bin width to estimate a continuous output firing rate r spike (t). This is the reference against which we test the rate-based model. To obtain the response of the rate-based model, we either use the known rate of the synthetic input spike trains or obtain a continuous input rate function a(t) from the input spike trains S(t) using the fixed-kernel density estimation method. Applying Eq. (14) to this rate yields the reponse of the rate model r rate (t).
We repeat each simulation experiment with five different random seeds and retain only results for which the optimal kernel width obtained by the densitiy estimation methods is 15 ms or less, as wider kernels would lead to an undue smoothing over time.
The difference between responses obtained from ratebased and spiking models is then defined as the mean squared error normalized by the variance of the response of the spiking model (Pillow et al. 2005) Rates change instantaneously at interval boundaries wherer spike is the average response rate of the spiking model. Corresponding to the linearity measure L 1 (Section 2.3.2), we define as quality measure. E r = 1 indicates perfect agreement, E r = 1/2 an error equal to the variance, while E r approaches 0 for large deviations between spiking and rate model response.

Tests with synthetic spike trains
We first test models using a Poisson spike train input with step changes in rates. Stimulus parameters are given in Table 5. Spike responses are obtained by simulating a population of 4,096 independent model neurons driven by one Poisson process each. The resulting 4,096 output spike trains are pooled to estimate the output rate r spike (t).

Tests with realistic spike trains
To test the performance of the rate models in response to realistic spike trains, we drive model neurons by spike trains recorded from retinal ganglion cells (RGCs) in cat (Casti et al. 2008). Their data set contained 128 spike trains of 8 s duration recorded during different trials, characterized by low baseline firing and fast transients as illustrated in Fig. 3. Because trains from the first trials in the dataset have noticeably lower average firing rates than those from later trials, we only use the last 96 spike trains with an average firing rate of 18.3 ± 1.3 spikes per second. The Izhikevich models in particular responds weakly to these spike trains in many cases. We therefore increase the rate of the input spike trains by merging pairs of spike trains, resulting in a total of 48 input spike trains with average rates of 36.6 spikes per second. We then drive 48 model neurons independently with one spike train each for 8000 ms and pool the resulting output spike trains for output rate estimation.

Simulation
Simulations for all model configurations are performed with the NEST Simulator (Gewaltig and Diesmann 2007;Plesser et al. 2013).
In practice, we simulate N trials by creating N mutually independent Poisson-generator-neuron pairs in a single NEST simulation. Membrane potentials are randomized upon network initialization and data collection is started Simulations underlying model fitting are performed using NEST 2.3.r10450, while some scoring of model responses according to Eq. (33) was performed using NEST 2.8.0. Trials are configured using the Neuro-Tools.parameters package (Muller et al. 2009). Data analysis is performed using NumPy 1.7.1-1.11.1, SciPy 0.18.1, Pandas 0.11.0-0.18.1 and Matplotlib 1.2.1-1.5.3 under Python 2.7.

Response to spike train input
To gain a first impression of the basic response properties of the models, we show the spike responses to stationary and sinusoidally modulated Poisson input in Figs. 4 and 5 for Izhikevich and AMAT models, respectively. Each raster plot shows the response of 30 unconnected neurons, half driven by stationary and half by sinusoidally modulated Poisson spike trains after an equilibration phase of 1000 ms. Each of the 30 neurons receives different realizations of input spike trains and noise, but the same trains and noise are used for all models.
As spiking and bursting variations are included as separate response types in the model classification (Fig. 1), we illustrate the burstiness of the responses by marking spikes fired within dT = 5 ms of each other as belonging to a burst, corresponding to the upper limit of intra-burst intervals in LGN (Funke and Wörgötter 1997, p. 71).
The models that exhibit their characteristic behaviour ( Fig. 1) based on "simple" excitatory input current shapes (e.g., steps, ramps, pulses) generally behave as expected when driven by Poisson spike trains; spiking neurons primarily spike and bursting neurons burst, but the nuances of individual models are less visible in the spiking patterns (e.g. tonic vs phasic) due to the input variability. Models that are based on more specific input current pattens or consistent inhibitory input (i.e., bottom rows) do to a lesser extent receive the required input and respond in a less characteristic manner, some even seem erratic (e.g., Fig. 5Q). Note, however, that the figures illustrate responses at a single input rate and noise regime combination and that the models to varying degree are sensitive to these conditions.
In contrast to the 20 markedly different responses to current injections (Fig. 1), responses to spiking input show more similar patterns across models, differing in the overall response rate and the proportion of spikes belonging to bursts.
While some Izhikevich and AMAT models that show identical responses to current injections also respond similarly when driven by spiking input (e.g., top two rows), we observe some with very different response patterns (e.g., depolarizing after-potential (Q) and inhibition-induced spiking (S)) across the two model classes.

Linear-nonlinear models
We now obtain the linear-nonlinear models as defined by Eq. (24).

Activation functions
We first obtain the activation function g(·) by fitting a linear B-spline to the response to stationary input, r 0 = g(a 0 ), varying a 0 from 0 s −1 to 1000 s −1 in steps of 10 s −1 . This yields one activation function fit for We thus obtain a total of 210 activation functions for the Izhikevich model class and 270 for the AMAT model class.
The top row of Fig. 6 shows the activation function for the Tonic Spiking and the Phasic Bursting models for the Izhikevich and AMAT model classes, respectively, for all three noise regimes. The spike rates obtained from the simulations are fitted very well by the B-splines. This holds for all models except model S in the Izhikevich class (inhibition-induced spiking), which has rather noisy activation curves with extremely high rates under certain conditions (above 1000 s −1 ); data not shown.
Before we investigate the response of the models to temporally modulated stimuli, we briefly explore the lin- input rates a 0 , and small modulation amplitudes a 1 give the most linear responses. Larger weights and mean rates not only increase the mean input of the Poisson input current, but also its variance. This leads to a linearization of the activation function and moves the activation threshold towards smaller rates (see also Chance et al. 2002). fitted filter H 0 (f ), thin solid lines: second harmonic r 2 , dotted lines: significance level r crit . Mean input rate a 0 and modulation range a 0 ±a 1 are marked gray in the top row. Fit parameters are given in Table 6. Third row: Response of spiking model (thin solid lines) and rate-model prediction (light thick lines) to Poisson spike trains with rate 100 s −1 for t < 700 ms and 300 s −1 for t ≥ 700 ms. Fit quality E r shown as inset. Bottom row: Response to realistic spike trains, 400 ms section starting at 3000 ms (cf. Fig. 3) with the same line types as for step responses. Connection weight w from left to right: 0.75, 700, 0.75, 500 Furthermore, firing-rate modulation amplitudes are more likely to stay within a single region of the sigmoidal firing rate curve for small a 1 , and are thus more likely to adhere to a linear fit. The stationary linearity metric L 1 for the augmented MAT model indicates overall more linear behavior, but the same general pattern of parameter dependence can be seen (Fig. 8). One notable difference is the saturation of the AMAT model at output rates of 500 s −1 -due to the absolute refractory time of 2 ms-that adds another source of non-linearity in the firing rate curves for some neurons.

Transfer function and linear filters
We obtain empirical transfer functions according to Eq. (17) for 20 combinations of working point and modulation depth (a 0 , a 1 ) for each model, noise regime and synaptic weight using the approach described in detail in Heiberg et al. (2013, Section 2.2.2), measuring the model response at 28 different stimulation frequencies f stim , logarithmically spaced from 1 Hz to 1000 Hz. We then fit the linear filterH 0,SUM (f ) according to Eq. (18) as described in Section 2.4, obtaining fit parameters (f c,1 , f c,2 , γ 1 , γ 2 , ) for each model and stimulation parameter combination. Note that γ 1 is fully captured by the activation function, and therefore does not explicitly enter the linear-nonlinear model we construct here, cf. Eq. (25).
The second row of Fig. 6 shows the resulting transfer functions and fitted kernels for the same models and conditions as the activation functions discussed above; see Table 6 for fit parameters. The examples reveal bandpass behaviour for three out of four cases, which also show significant power in the second harmonic r 2 . The exception is Tonic spiking for the AMAT class, which shows lowpass behavior and no significant power in r 2 . Phasic bursting shows a second peak in the spectrum around 200 Hz (Izhikevich) and 500 Hz (AMAT), which our fitted bandpass filter models (thick light lines) cannot capture by construction. These peaks occur as refractory effects regularize firing patterns at high rates. The fitted bandpass filters capture the frequency response of the model neurons well, except for Izhikevich phasic bursting case, where the amplitude of the fitted filter is significantly larger than the power in the first harmonic.
We found that not all model variants responded sufficiently to periodic stimulation under all stimulation conditions to provide sufficient spike data to fit a kernel. Therefore, we only obtained kernel fits for approximately three-quarters of all conditions for the Izhikevich class (3262 out of 4200 possible) and about 90% of all conditions for the AMAT class (4843 out of 5400).

Grouping of models
To systematize model responses, we cluster the kernel fit parameter sets 4 with k-means clustering using Scikit-Learn (Pedregosa et al. 2011). We cluster Izhikevich and AMAT filters independent of each other. To find a suitable clustering, we ran the clustering algorithm searching for four, five, six, and seven clusters. For each run, we clustered starting from 100 initial conditions to avoid local minima. Since we are clustering fit parameters obtained for a wide range of simulation conditions, while we are interested in grouping the 14 and 18 model variants, respectively, we assign each model variant to a model group as follows: We count how often each model occurs in each k-means cluster and assign each model to the cluster to which it is assigned most often. Each cluster to which at least one model is thus assigned forms a model group. Low number of clusters led to model groups containing models with clearly different spiking behavior, e.g., for the AMAT model, the strongly bursting case J would be grouped with cases A, H and K, which show almost no bursts at all, jf Fig. 5 when using six clusters. We thus used seven clusters in all cases which yielded six model groups. Figure 9 shows the counts for each model in each cluster, with a maximum possible count of 300 if fit parameters are available for all combinations of μ, σ, w, a 0 , a 1 and are assigned to the same cluster. The resulting grouping into six groups per model class is shown in Table 7, with median kernel parameter values for each group in Table 8. Grouping is clearly different for the Izhikevich and AMAT classes, supporting our previous observation that these models respond differently to spike input even though they show identical responses to the current injection protocol of Fig. 1.
Comparing the grouping of models to the spike responses shown in Figs. 4 and 5, we can roughly identify the groups found by k-means clustering of filter parameters to firing patterns, as indicated in the right column of Table 7. This classification is far from perfect, as several models show firing patterns different from the groups into which they have been placed, especially for the AMAT class. It should also be noted that the firing patterns are for a single stimulus configuration only and that models may behave differently under other conditions; the k-means clustering, on the other hand, is based on a wide range of stimulus conditions.

Performance of rate models
We evaluate the perfomance of the linear-nonlinear firing rate models by testing them against the corresponding spiking model as described in Section 2.5, using the fit quality E r as criterium, with E r = 1 indicating a perfect fit.
The third row of Fig. 6 shows the response to a Poisson spike train with a step in rate from 100 s −1 to 300 s −1 . We Table 7 Models grouped by k-means clustering of linear filter parameters as illustrated in Fig. 9 Group Firing pattern Izhikevich AMAT Regular isolated spikes H/Class 2 P/Bistability P/Bistability Q/Depolarizing after-potential* Group numbers are arbitrary and do not correspond to cluster numbers in Fig. 9. Models with a different firing pattern in Figs. 4 and 5 than their group are marked with an asterisk. Filter kernel parameters for the groups are shown in Table 8 use the filters fitted for the same noise regime and synaptic weight and a 0 = 200 s −1 and a 1 = 100 s −1 , corresponding to the step height. For the Tonic spiking case, the firing rate models capture the spiking neuron response very well, with E r > 0.9 in all cases (see legend). For the Phasic bursting models, we find that the rate models overshoot massively for the Izhikevich variant with no or balanced noise, while the rate models "undershoot" somewhat for Table 8 Median values of parameters for filter kernels fitted to the six groups described in Table 7 Table 5, where E r = 1 indicates a perfect fit. Each major column corresponds to one model, each major row to one of 15 input conditions; thus, each major cell corresponds to one firing-rate estimate r spike (t). Each entry inside a major cell represents the fit against one of 20 filter models obtained for the same μ, σ, w and different combinations of a 0 , a 1 as in Fig. 7. Missing results (grey) are either due to models excluded (grey labels) or insufficient responses the AMAT variant. The stationary rate attained after the step is captured well in all cases. These examples also provide an illustration of how to interpret the fit quality measure E r . Figures 10 and 11 show the fit quality observed for responses to Poisson input with piecewise constant rates as described in Section 2.5.1. For each model, we simulate responses under 15 input conditions (three noise regimes Fig. 11 Fit quality E r for AMAT model class responses to a Poisson process with firing rate changing in steps; all else as in Fig. 10 (μ, σ ) and five different synaptic weights w), yielding 15 firing-rate estimates r spike (t). We then test each firingrate estimate against the 20 linear-nonlinear rate models obtained for the same μ, σ, w and all a 0 , a 1 combinations, yielding 20 fit quality values E r .
For the Izhikevich class, the Inhibition-induced spiking and bursting models, as well as most models with the lowest weight, w = 0.1, produce too few spikes to confidently estimate firing rates from the spiking model. We also observe very poor responses for the Bistability model. Class 2 excitable stands out with poor scores, E r < 0.5, while the remaining models provide reasonable fits, E r > 0.7 at least for most cases with sufficiently strong weights (w ≥ 0.4).
The AMAT model class performs significantly better: Results are available for almost all stimulus conditions except for w = 100 pA in the absence of noise and all models except the Depolarizing after-potential model yield excellent fits (E r > 0.9) for almost all conditions. Figures 12 and 13 show the fit quality E r for responses to real spike trains from cat retinal ganglion cells as described in Section 2.5.2. We again stimulate under 15 different conditions and obtain the fit quality for each of 20 different linear-nonlinear model fits.
For Izhikevich-class models we find noticeably worse fit quality, mostly E r < 0.7, with the worst results mostly for the same conditions that also yielded low fit quality in response to Poisson input with piecewise constant rate. The main differences are that the Depolarizing after-potential model, which fitted stepped Poissonian input very well does not perform better than other models for the real spike trains, and that we obtain quality of fit values, albeit very poor ones, for the Inihibition-induced bursting model.
AMAT class responses to real spike trains show all over better fit quality than the Izhikevich class, but also for the AMAT class fit quality is lower in response to real spike trains than to stepped Poisson input. The distribution of good and bad fits is similar to the one observed for stepped Poisson input, with the worst performance for the Depolarizing after-potential model. Furthermore, more models require w ≥ 300 pA to yield a fit quality result for real spike trains.
We summarize these observations in Fig. 14, which shows the cumulative distribution P (E r ) of all individual results from Figs. 10-13 (thin lines). This clearly shows that our linear-nonlinear rate models are much more faithful for the AMAT class than for the Izhikevich class, and that within each class, responses to stepped Poisson input are rendered more faithfully than to real spike trains.
These observations pertain to all (up to 20) linearnonlinear models obtained for each input configuration μ, σ, w. In practice, we would only be interested in the optimal linear-nonlinear model for each input configuration, Fig. 12 Fit quality E r for Izhikevich model class responses to the real spike trains from Fig. 3; all else is as in Fig. 10  Fig. 13 Fit quality E r for AMAT model class responses to the real spike trains from Fig. 3; all else as in Fig. 10 i.e., for given μ, σ, w we would choose the model with the highest E r across all a 0 , a 1 combinations. The cumulative distribution of fit quality for these optimal models is shown as thick lines in Fig. 14  for each input configuration μ, σ, w. If all fits were perfect, P (E r ) would hug the x axis until jumping to 1 for E r = 1 noticeably better fit quality for stepped Poisson and real train responses for both model classes. Table 9 shows the proportion of cases for which we reach high fit quality (E opt r ≥ 0.8) for the optimal models. We find that the linear-nonlinear rate models perform well for the majority of conditions for the AMAT model class, but mostly poorly for the Izhikevich model class.

Model generalizations
As we have shown above, for the AMAT model class our linear-nonlinear rate models can capture the responses of spiking neuron models to real spike trains quite accurately. For the Izhikevich class, on the other hand, fits were poorer. Unfortunately, to find the optimal linear-nonlinear model for each input configuration μ, σ, w, we had to test a set of 20 different linear-nonlinear models to then pick the best one. This is impractical. We will now consider how to generalize our linear-nonlinear rate models, so that we can select an optimal model a priori.
We consider four different types of generalization:  For the M and MN generalizations, we exploit that the activation functions g(a) for many models and conditions scale roughly linear in the synaptic weight. We thus pool the scaled activation function data g(a)/w for a given model across all input conditions (M) or just all synaptic weights for given noise (MN) and fit a single splineḡ(a) to the pooled data. We then use wḡ(a) as activation function in the linear-nonlinear model. For MNW and MNWS we use the original splines fitted directly against measurements.
To generalize the linear kernels, we take the median value for each of the kernel fit parameters f c,1 , f c,2 , γ 1 , γ 2 , d and use these median parameters as parameters of our generalized kernel; using the median instead of the mean avoids problems with outliers. For M generalization, we take the median across all μ, σ, w, a 0 , a 1 combinations, for MN across all w, a 0 , a 1 for given μ, σ and for MNW across all a 0 , a 1 for given μ, σ, w.
For MNWS generalization, we proceed differently: For each combination of μ, σ, w we select the filter parameters f c,1 , f c,2 , γ 1 , γ 2 , which yielded the highest fit quality E r = E opt r in response to the stepped Poisson input, our test stimulus.
While these generalizations, especially of the M and MN type, may seem rather crude, they perform reasonably, as . This is consistent with the overly large amplitudes of the fitted filters. The corresponding generalized model filters fit the experimental data better, avoid the overshoot and thus track the response of the spiking model better.
To systematically quantify the quality of our generalizations, we compute the fit quality in response to stepped Poisson and real spike train input for all input configurations μ, σ, w for each generalization variant. Then measures how close the model generalization X is to the optimal linear-nonlinear model, where ρ X = 1 is best.
Results are shown in Figs. 16 and 17. The coarser the generalization, the more frequently do we observe low generalization quality ρ X . Interestingly, differences between the various generalizations are larger for the AMAT class than for the Izhikevich class, and generalization seems to fail for the AMAT class mostly for biased noise and Inhibitioninduced spiking. For the Izhikevich class, on the other hand, generalization mostly fails for low synaptic weights.
The most important observation, though, is that MNWS generaliztion works well, with ρ MNWS > 0.9 in almost all cases for both model classes. This means that by selecting a filter model based on a fixed stepped Poisson protocol, we will obtain a linear-nonlinear rate model that is close to the optimal model for given noise regime and synaptic weight when applied to real neuronal dynamics.
Combined with the observation from Table 9 that the optimal model will provide a good approximation to actual neuronal firing rates in roughly two thirds of all conditions, we can thus use our fitting approach together with the stepped Poisson test to select a reasonably reliable linearnonlinear neuron model.

Discussion
In this paper we numerically investigated the response properties of two neuron model classes, the Izhikevich model and the AMAT model, to noisy spiking input. Both neuron models can reproduce a wide range of experimentally observed spike response patterns when stimulated with current injections. However, how these neurons behave with more natural synaptic inputs has so far not been studied systematically. We considered three different background noise regimes, one with no background noise at all, one balanced and one biased with enough background noise to put neurons in a spontaneously active state at low output rates. The first scenario can be considered to represent the situation in slice preparations, the other two correspond to neurons embedded in a network with ongoing excitatory and inhibitory activity. The stimulus spikes were modeled as stationary and sinusoidally modulated excitatory Poisson input spike trains, mimicking afferent inputs from sensory pathways with different synaptic connection strengths w.

Responses to spike input
We found that the response complexity observed under current injection collapses to only a few response types when the neurons are driven by stationary or sinusiodally modulated Poisson input. This is not entirely surprising, since some of the models are parametrically quite similar, and variations in response behavior to current stimulation depend on very specific current injection patterns that are not realizable in terms of Poisson spike inputs. Still, actual neurons receive inputs that often are well-described by Poissonian statistics and this can thus be considered the functionally more relevant input scenario. It is hence of interest to see which, possibly quite different, neuron models behave approximately equivalent.
The respective groupings for Izhikevich and AMAT are all-in-all very different. In particular, direct comparison of the individual corresponding neuron models reveals completely different response properties for most neuron models. This is in part explained by the differences in subthreshold dynamics which are linear for the AMAT model but nonlinear for the Izhikevich model. Individual spikes thus have quite different effects in the two models: Any input spike to an AMAT neuron will always evoke the same postsynaptic membrane-potential response and these responses simply superimpose due to the subthreshold linearity. For the Izhikevich models, on the other hand, the postsynaptic response depends intricately on the value of the dynamic variables, such as the membrane potential, and the effect of several incoming excitatory spikes of same weight at one moment might be smaller than that of just one such spike at another moment. These differences hence make it hard, or even impossible, to set up synaptic weights w for the two neuron model classes that are directly comparable.
We therefore chose to gauge synaptic strengths in terms of the minimal weight w θ needed to evoke a spike from rest, cf. Section 2.2, and to use weights spanning roughly from 10% to 75% of w θ . This allowed us to quantify and compare input coupling strength within and between model classes. In general we observed that output rates for Izhikevich neurons were much lower than for AMAT neurons for the same input frequency and relative synaptic strength. It is therefore possible that model classes might become more similar if the Izhikevich neurons were driven at higher input rates or at other background noise levels, although we did not observe such a trend.
We observed here that neuron models can show very similar responses to spike input, even though they show very different responses to current injections, and in particular that models of different mathematical nature, showing identical current responses can respond very differently to spiking input. Given that neurons are mainly driven by spike input in vivo, this raises the intriguing question of how valuable a classification of neuronal response types based purely on current injection experiments is. While in vitro characterization using carefully crafted current injections is an important tool to classify neuronal cell types, it appears that a systematic classification based on a neuron's response to spiking input may be required to select suitable neuron models for spiking and rate-based network models.

Firing-rate models
In the second part of the paper, we made use of the measured stationary and frequency responses to fit linear-nonlinear firing-rate models to the data. It has previously been shown that the firing-rate dynamics in response to complex spiking input can be well described by such models (Paninski et al. 2004;Ostojic and Brunel 2011;Weber and Pillow 2017;Østergaard et al. 2018). In particular, Nordlie et al. (2010) studied simple leaky integrate-and-fire (LIF) models with strong current-based synapses. They showed that a lowpass fit to the frequency response together with the nonlinear activation function yielded linear-nonlinear rate models that predicted responses to arbitrary inputs with high accuracy. Heiberg et al. (2013) adapted this approach and studied two LIF-like models, one with current-based, the other with conductance-based synapses, that were fit to actual data recorded from cat and macaque LGN in response to retinal stimulation. They found the performance of linear-nonlinear rate models to be good as well.
Here, we presented results of the same basic approach for the Izhikevich and AMAT neuron model classes. Frequency responses were in most cases more complex than simple lowpass behavior and we employed fits to bandpass filters that better capture the non-monotonous passband structure observed in simulations. We then used novel test stimuli, i.e., step responses and more structured, highly variable spike input sampled from actual recordings of retinal ganglion cells (Casti et al. 2008) to study rate-model performance. The main finding is that the AMAT neuron model class is approximated much better than the Izhikevich class by our linear-nonlinear rate models: for the former, good rate model responses (E r ≥ 0.8) were obtained in 64% of all cases tested, while the latter provided such good results in only 15% of cases tested. This difference might again be explained by the fact that the AMAT model class has subthreshold linear dynamics. However, the AMAT class is not completely linear either, because its firing threshold depends on the history of the membrane-potential dynamics. Therefore, neuronal transfer is not expected to be linear in any model class. Some of the model variants gave consistently poor results, typically those that show very nonlinear behavior in response to direct current stimulation, e.g., the Bistability, Inhibitioninduced spiking and Inhibition-induced bursting models for the Izhikevich class, and the Depolarizing afterpotential and Inhibition induced spiking models for the AMAT class.
To estimate the effects of linearity on rate-model performance, we measured the linearity of the stationary response function r 0 (a 0 ) in terms of L 1 , cf. Eq. (12). If the stationary response function is linear, the activation function g(·) is also linear and only its slope is relevant, independent of the working point, cf. Section 2.4. We computed L 1 for all background noise regimes as a function of synaptic strength w and working point a 0 , and find that the AMAT model generally is more linear than the Izhikevich model with respect to L 1 . We further find that the linearity measure L 1 does not predict rate-model performance (data not shown). Thus, a nonlinear activation function does not imply poor rate-model performance, nor does linearity in terms of L 1 necessarily predict a good rate-model performance.
Furthermore, despite exploration of many potential performance predictors, we were unable to identify any single quantity or group of quantities that reliably predicted whether the response of a neuron model in a given input regime could be captured well by a linear-nonlinear rate model. We found, though, that a relative simple protocol, testing the rate model's performance in response to a Poisson spike train input with piecewise constant rate (stepped Poisson), allowed us to reliably identify rate models that render spiking neuron model responses to realistic spike input with reasonable accuracy.

Application to network modeling
We have shown that the firing-rate responses of the widelyused Izhikevich model (Izhikevich 2003b) and in particular of the award-winning augmented multi-timescale adaptive threshold (AMAT) model (Yamauchi et al. 2011;Jolivet et al. 2008) can be captured by linear-nonlinear firingrate models with bandpass filters fitted to spiking neuron responses through a systematic, automated process. We have further shown that these models can be generalized across a wide range of input conditions without excessive loss of fidelity, and that optimal parameter sets can be chosen using simple test stimuli. Furthermore, since we use a bandpass filter in sum form, it can be represented by a system of two first-order differential equations, which is straightforward to integrate into standard formalisms for rate-based network models (see, e.g., Nordbø et al. 2007).
This suggests the following approach to improve ratebased neuronal network models based on our findings. Assuming a model with two neuronal populations, start by selecting for each population the AMAT (or Izhikevich) model variant that best matches the response of individual neurons of each population to spiking input. Then apply the fitting procedure described in this paper (Section 2.4) to obtain the parameters of the nonlinear activation function and the linear bandpass filter using test stimuli covering the expected dynamic range of your network model. From the large set of fits obtained, either select individual fits based on a simple test protocol (Section 3.4) or generalized at a suitable level (Section 3.5), and apply the model parameters thus obtained to the differential-equation representation of the bandpass filter (Section 2.4.1).
The approach presented here may thus contribute to bringing rate-based network modeling closer to the reality of biological neuronal networks. While a systematic delineation of the range of validity of the linear-nonlinear models described here for network modeling is beyond the scope of this paper, we consider the generalization results in Figs. 16 and 17 an indicator: Linear-nonlinear models will not be useful in cases where responses are insufficient (gray areas in Figs. 10-13 and 16-17); we also observed poorer performance for spiking patterns that deviate strongly from tonic behavior (e.g., phasic bursting, bistability, depolarizing afterpotential). But good generalization scores, in particular where coinciding with good scores for individual conditions (Figs. 10-13) suggest that in these cases our linear-nonlinear model provides a faithful representation for the rate dynamics of the underlying spike responses, and hence a strong potential for good performance also on the network level. the Heaviside step function. We will commonly refer to I E as excitatory and I I as inhibitory input. For efficient integration of the synaptic currents (Plesser and Diesmann 2009), we describe them by differential equationṡ I X + τ X I X = j w X,j δ(t −t X,j ) .
We can summarize the resulting system of eight firstorder linear differential equations aṡ with state and input vectors and system matrix I ext (t) represents the jumps in the piecewise constant external input current.
We can solve this system exactly (to the limits of machine precision) on a fixed time grid t j = jh for h > 0 using exact integration (Rotter and Diesmann 1999), provided that I ext (t) only changes at grid points t j , i.e., I ext (t) = j δI j δ(t − t j ). Starting from the initial state y 0 = y(t = 0), exact integration updates the state according to y j+1 = Ay j + x j+1 (51) where A is the propagator matrix The δ(t j −t k )-functions in the input vector x(t) are replaced by Kronecker symbol δ j,k upon discretization, restricting spike times to the time grid; see Morrison et al. (2007) for an extension to off-grid spikes. The propagator matrix P can be obtained numerically, e.g., using the expm functions provided by SciPy or Matlab, using an algorithm due to Higham (2005). These methods can, though, fail under certain circumstances (Moler 2012; Al-Mohy and Higham 2009) and we have not performed any systematic tests regarding their reliability with respect to model neuron dynamics. We used the Mathematica symbolic algebra system (Wolfram 1999) to obtain an explicit expression for the propagator matrix P and generate C++-code for the matrix elements. Note that the expression obtained in this way requires that all time constants except τ 1 and τ 2 differ, i.e., τ m = τ V = τ E = τ I . Other expressions for P pertain if any two time constants are equal. The resulting model is implemented in NEST as amat2 psc exp.
While our implementation of the AMAT2 model is based on the equations given by Yamauchi et al. (2011), we found that we needed to modify some model parameters to reproduce the responses in Figs. 6 and 7 of that paper: -all models: membrane time constant τ m = 10 ms instead of τ m = 5 ms subthreshold oscillations: α 1 = 1, β = 0.2 instead of α 1 = 10, β = 0.1 resonator: β = 0.5 instead of β = 0.1 threshold variability: β = −0.5 instead of β = −0.1 Parameters used in our model are given in Table 4.
We further needed to make some adjustments to stimulus parameters to reproduce Figs. 6 and 7 in Yamauchi et al. The changes in stimulation currents are only relevant for the illustrative Fig. 1 and do not directly affect the remainder of the results presented here.
We can only speculate about why these parameters changes were necessary. We inferred the time axis of Figs. 6 and 7 of Yamauchi et al. (2011) from the scale bars given in panel A of each figure. This clearly indicates that pulse/ramp durations given in Table 1 of that paper are inconsistent with the figures for, e.g., bistability (Fig. 6F) and accomodation (Fig. 7I). Concerning the membrane time constant τ m , the authors tested various values of this parameter (Yamauchi et al. 2011, p. 2, right column) and this may have led to a mix-up.