Complete Firing-Rate Response of Neurons with Complex Intrinsic Dynamics.

The response of a neuronal population over a space of inputs depends on the intrinsic properties of its constituent neurons. Two main modes of single neuron dynamics-integration and resonance-have been distinguished. While resonator cell types exist in a variety of brain areas, few models incorporate this feature and fewer have investigated its effects. To understand better how a resonator's frequency preference emerges from its intrinsic dynamics and contributes to its local area's population firing rate dynamics, we analyze the dynamic gain of an analytically solvable two-degree of freedom neuron model. In the Fokker-Planck approach, the dynamic gain is intractable. The alternative Gauss-Rice approach lifts the resetting of the voltage after a spike. This allows us to derive a complete expression for the dynamic gain of a resonator neuron model in terms of a cascade of filters on the input. We find six distinct response types and use them to fully characterize the routes to resonance across all values of the relevant timescales. We find that resonance arises primarily due to slow adaptation with an intrinsic frequency acting to sharpen and adjust the location of the resonant peak. We determine the parameter regions for the existence of an intrinsic frequency and for subthreshold and spiking resonance, finding all possible intersections of the three. The expressions and analysis presented here provide an account of how intrinsic neuron dynamics shape dynamic population response properties and can facilitate the construction of an exact theory of correlations and stability of population activity in networks containing populations of resonator neurons.

Introduction

Integration and resonance are two operational modes of the spiking dynamics of single neurons. These two modes can be distinguished from each other by observing the neuron’s signal transfer properties: how features in its input current transfer to features in its output spiking. The traditional approach to investigating neuronal transfer properties is to measure the stationary response: the time-averaged rate of firing of spikes as a function of the mean input current, or fI-curve. In Hodgkin’s classification [1], Type I membranes can fire at arbitrarily low rates, while the onset of firing in Type II membranes occurs only at a finite rate. This distinction arises naturally from the topology of the bifurcations that a neuron can undergo from resting to repetitive spiking [2]. In many central neurons, it is fluctuations rather than the mean input current that drive spiking, putting them in the so-called fluctuation-driven regime [3]. Many dynamical phenomena are nevertheless tightly linked to excitability type. For example, Type II neurons exhibit rebound spikes, subthreshold oscillations and spiking resonance (e.g. mitral cells, [4-6], respectively). The qualitative explanation for these phenomena is that the dynamical interplay of somatic conductances endow some neurons with a voltage frequency preference, i.e. a subthreshold resonance. This preference can contribute to a superthreshold resonance in the modulation of their output spiking [7]. How dynamic response properties of spiking dynamics such as resonance emerge can be directly assessed by considering the neuron’s dynamic gain.

Dynamic gain, first treated by Knight [8], quantifies the amount by which features at specific frequencies in the input current to a neuron are amplified or attenuated in its output spiking. It can accurately distinguish functional types and unveil a large diversity of phenomena shaping the response to dynamic stimuli [9-18]. Dynamic gain and response are also essential ingredients for theoretical studies of network dynamics in recurrent circuits [8, 12, 13, 18–49]. First, they determine the stability of the population firing rate dynamics [21, 25, 26]. Second, they determine how input correlations between a pair of cells are transferred to output correlations [28, 42, 44–49], from which self-consistent relations for correlations in recurrent circuits can be obtained.

Experimental studies have started over the past years to use dynamic gain measurements to investigate the encoding properties of cortical neuron populations [9-18]. Although theoretical studies have investigated many neuron models, very few models are known for which dynamical response can be explicitly calculated. One basic reason for this lies in the fact that Fokker-Planck equations for neuron models with two or more degrees of freedom are not solvable in general [50]. For Type II neuron models that require at least two degrees of freedom, no solvable model is known.

The simplest model capable of subthreshold resonance was introduced by Young [51] in the early theories of excitability. Later, Izhikevich formulated a structurally similar neuron [52]. Richardson and coworkers performed the first calculation of the linear response function of a neuron model capable of resonance, the Generalized Integrate-and-Fire (GIF) neuron [22, 29]. Only in the limit of relatively slow intrinsic current time constant can analytical expressions for the GIF response be obtained, however. The distinct transfer properties of resonant vs. non-resonant dynamics leads to different information transfer properties. While this has been demonstrated in the mean-driven regime [53, 54], no such results exist for the fluctuation-driven regime, in part due to a lack of exact analytical expressions for even the linear dynamic gain. Type II excitability and dynamic response thus are representative of the more general challenge posed by response properties of neurons with complex intrinsic dynamics.

In the current study, we derive and analyze the linear response function in the fluctuation-driven regime of a neuron model capable of resonance. We express it as a filter cascade from current to voltage to spiking. It is valid across all relevant input frequencies and over all relevant values of the intrinsic parameters. In particular, we apply to the GIF neuron model the Gauss-Rice approach in which the voltage reset after a spike is omitted. The methods generalize to additional intrinsic currents and to the full nonlinear response with spike generation. To understand how subthreshold features interact to determine a neuron’s filter characteristics, including resonance, we provide a two-dimensional representation of the response properties that completely characterizes all possible filter types. For this idealized model, we determine analytically and numerically a wide and biologically-relevant regime of validity of the derived expression.

The paper begins with the definition of the model and its numerical implementation. We then derive a general expression for the linear response in the mean channel of a Gauss-Rice neuron. In the next section, the analytical results for the response properties of the Gauss-Rice GIF neuron model are obtained. The final section then presents an analysis of the expression. For the sake of mathematical clarity, most calculations appear in the main text; the rest, including an exposition of model assumptions, are contained in the Methods.

Results

Definitions and methods for a population of Gauss-Rice GIF neurons

We consider the most simple hard-threshold, no-reset, GIF-type neuron capable of exhibiting resonator dynamics, whose response properties have been partially studied in [25]. A reset version of this model is treated in [29], where the population spiking response properties were calculated assuming large intrinsic time constant. In the Methods, we present a more detailed exposition of the model assumptions, and justify an additional simplification of the voltage reset after a spike. The feature that distinguishes the GIF model from the classical Leaky Integrate-and-Fire (LIF) model is that the dynamics of the voltage, V, is coupled to an intrinsic activity variable, w, where g is a relative conductance and τ and τ are the respective time constants of the dynamics. The notation denotes the derivative with respect to time of the variable x. Spikes are emitted at upward crossings of a threshold, θ. Synaptic current modeled by I drives the model whose dynamics are kept stable by keeping g > −1. When g < 0, w is depolarizing. When g > 0, it is hyperpolarizing and can lead to resonant voltage dynamics.

Intrinsic parameters, g and τ, shape the phase diagram of the intrinsic dynamics

Here we present analysis of the phase diagram of the intrinsic dynamics of the model, which is a reparametrization of Fig.1 from [29]. Beyond that work, here we analyze the Ω-contour density and scaling behavior. For a fixed, constant value of I , and with time in units of τ, the structure of the phase space of the single neuron dynamics described by Eq (1) is determined by a point in the τ/τ vs. g plane, the two parameters defining the intrinsic current, w (see Fig 1). For τ ≪ τ, w speeds up or slows down V depending on whether g is hyperpolarizing (g > 0) or depolarizing (g < 0) characterized by an effective time constant While the dissipative voltage term stabilizes the voltage dynamics, the dynamics can be effectively unstable for g < −1, and we exclude this case. For depolarizing intrinsic current, there is a region where the two eigenvalues of the voltage solution, λ ±, are complex and the model exhibits an intrinsic frequency, Ω = 2πf , that varies as where g = (τ − τ)2/4ττ (see Methods for details). For a fixed g > 0, a given value of Ω can be achieved at both a high and a low value of τ. For fast τ, the Ω-contour density is high and the model exhibits high parameter sensitivity, while for large τ the contour density is low and the model is relatively insensitive to local parameter variation. Taking the respective limits, the set of isofrequency curves are linear for large τ with slope ∝ Ω2 and with a slope independent of Ω for small τ.

Fig 1

The type of w-current depends on the values of the intrinsic parameters.

(a) Intrinsic parameter phase diagram in (τ/τ, g). w can be depolarizing (g < 0) or hyperpolarizing (g > 0). w contributes an intrinsic frequency to the model in the colored region. The dynamics are unstable if g < −1. Iso-Ω lines are shown in white ( for large τ and for small τ). (b) When , the phase diagram can be cast in (τ/τ, Ωτ )-space. Iso-g lines are shown in white. (See [29] for a similar plot).

Furthermore, there is a minimum relative conductance, for which a given Ω can be achieved. The minimum shifts to increasingly short τ with Ω. To emphasize the timescale of the intrinsic frequency when it exists, we reparametrize the model by replacing g with Ω using Eq (3), and arrive at the implicit representation shown in Table 1) (see Methods for more details). The statistical structure of the relative timings of the output spikes of the model will be affected by Ω.

Table 1

Parameter groups for each dynamics.

Neuron model	Intrinsic	Input	Voltage	Voltage filter	Spiking filter
τ_Vgτ_w}Explicit Rep.	τ_rΩτ_w}Implicit Rep.(λ_±)	τI	σ V	ωLQL=ωLτ_r/2τ_ω}Filter Rep.	ν 0
		σI	σ_V˙		τc
		ω	τ_s=σ_V/σ_V˙		νωL = ν 1(ωL)/ν 1(0)
θ	τ eff ≔ τ V/(g+1)	A	σ = σ V/θ	V_low:=V¯˜(0)	ν ∞ = ν 1(∞)/ν 1(0)

Table notes: There are three equivalent representations of the parameter space of the voltage dynamics based on the explicit parameters of the model, the implicit timescales of the model, and the filter parameters of the mean voltage filter, respectively.

Parameter notation: τ , membrane time constant; g, relative conductance; τ , w time constant; θ, voltage threshold; τ , relaxation time of voltage dynamics; Ω, intrinsic frequency; λ ±, eigenvalues, when ω τ > 1, λ ± = r ± iΩ, with real part, ; τ , effective membrane time constant; τ , input noise correlation time; σ , input noise standard deviation; ω, input signal frequency; A, input signal amplitude; σ , standard deviation of voltage; , standard deviation of voltage time derivative; τ , differential correlation time; σ, relative standard deviation of voltage; ω and Q , center frequency and Q-value, respectively, of low pass component of current-to-voltage filter; V , low frequency voltage response; ν 0, stationary firing rate; τ , characteristic time of voltage-to-spiking filter; ν ≔ ν 1(ω )/ν 1(0), spiking response at ω ; ν ∞ ≔ ν 1(∞)/ν 1(0), high frequency spiking response.

Population firing rate dynamics

Given a population of N neurons indexed by k, in a time window, T, each one produces a spike train, with n spikes labeled as . The average firing rate across the population in this window is For stationary input, this becomes the stationary population averaged firing rate, independent of t, in the limit T → ∞, In the other limit, taking T → 0 while keeping NT constant, such that there is a statistically invariant number of spikes in the time window, the integrand of Eq (5) is a well-defined time-dependent ensemble average, the instantaneous population firing rate, where denotes the population average of a single neuron quantity, x. Note that this population firing rate can exhibit time dependence on arbitrarily fast timescales.

Populations of fluctuation-driven neurons

The input to the neuron, I , from Eq (1) arrives from many, weak synapses. The total drive will thus resemble a continuous stochastic process. The system can then be solved under this assumption by directly simulating the corresponding stochastic differential system of Eq (1), where is the time-dependent mean input and δI(t) is a zero-mean noise process. Solutions give the output spike times, which averaged over an ensemble give the population firing rate, ν(t). Under the diffusion approximation, discussed in more detail in the Methods, the stochastic drive, δI(t), can be taken as an zero-mean Ornstein-Uhlenbeck process with variance and correlation time τ. The resulting stochastic dynamics were simulated by numerical integration via a Runge-Kutta scheme (see ref. [55] for details).

To illustrate the dynamic ensemble response, we show in Fig 2 an example of input, intrinsic, and output variable time series produced by the model for two choices of signal in the mean channel, : a weak oscillation of amplitude A and frequency ω and, separately, a step of height Δ. In addition, we show the corresponding population firing rate dynamics obtained from a histogram of the spike times of the sample ensemble produced by the two inputs. Code to produce this plot can be found in the supplemental material. The input modulation structures the spike times produced by the ensemble relative to the stationary response in a way that only becomes salient at this population level. We motivate consideration of the analytical expressions for the linear response function (Eq (41)) and the step response function (Eq (63)), obtained later in this paper, by plotting their curves, which accurately overlie the profile of the two respective measured histograms. While the input oscillation produces modulation in the output spiking at only one frequency, the step input produces a response that has power across a broad band of frequencies.

Fig 2

From input to ensemble response: numerics and prediction.

Model output for the default parameter set: τ = 1ms, σ = 1, τ = 10ms, θ = 1, τ = 20ms, f = 20Hz (g = 3.15). Left: in in the case of an oscillation of amplitude A = 0.05 and input frequency ω = (2π)20 rad/s. Right: in the case of a step of height A = 0.1. The example realization shown is the one with the maximum number of spikes from the sample ensemble. The red line is the response calculated using the analytical expressions for the oscillation and step response, Eqs (41) and (63), respectively.

Approaches to obtaining the population response

Response theory captures the population response to input signals with arbitrary frequency content and so we now turn to it, and linear response theory in particular, in the pursuit of understanding the population firing rate dynamics of the GIF neuron model.

The formal, implicit definition of the linear response function, ν 1(ω), arises from a weak oscillatory modulation of amplitude A and frequency ω in the mean input, and an expansion of the response, ν(t), in powers of A, where ν 0 is the stationary response, Eq (6). In the Methods, we restate how the linear response can be obtained (Eq (59)) directly from the spike times using the complex response vector, . Below, we show the classic formulation that shows that it can also be obtained from the voltage dynamics.

Obtaining the response from the statistics of the voltage dynamics

To obtain ν 1(ω) analytically, we go back to the definition of ν(t) containing s (t), Eq (7). s (t) can be rewritten as where Θ is the Heaviside theta function defined as Θ(x) = 0 for x < 1 and Θ(x) = 1 for x > 0. appears since spikes are only generated at upward threshold crossings of the voltage. The factor results from the coordinate change in the argument of the δ-function. When combined with , the absolute value can be omitted. For a population of such neurons, we can then obtain the population-averaged firing rate as the rate of upward threshold crossings known as Rice’s formula [56], The underlying ensemble of the population is captured by the distribution of voltages and voltage time derivatives at a given time, . When each neuron’s state is identically and independently distributed, the average over k neurons is an average over this distribution at fixed t, This time-varying expectation value over the statistics of the voltage dynamics in the population is the central time-domain quantity in the response theory for neuronal populations. It is in general analytically intractable.

Subthreshold dynamics can be approximately linear and the many, weak inputs to each neuron can permit a diffusion approximation to a Gaussian process input. In this situation, a model of voltage dynamics that omits the nonlinear voltage reset gives a voltage statistics that is also Gaussian and can be treated analytically. This is the Gauss-Rice approach, to our knowledge first published by Jung [57], where it was used to calculate correlation functions. The first application of the approach to dynamic gain appears in Supplementary Note 3 of ref. [37]. Note that the lack of a voltage reset required for this approach restricts its range of applicability (see Methods and Discussion).

Dynamic gain for the mean channel of a population of Gauss-Rice neurons

Here, we derive the dynamic gain using the Gauss-Rice approach, similar to [19, 37, 42]. The emphasis of the derivation here differs in that, first, we go directly to the linear response by linearizing around the mean voltage, and second we express the results in terms of the voltage transfer function to emphasize the additional filtering of the voltage dynamics by the spike threshold. The resulting expression, Eq (22), applying to a generic population of neurons specified only by the Gaussian statistics and frequency response of their mean voltage dynamics, simply adds a first-order highpass filter to the frequency response of the voltage. In units of the differential correlation time, τ, the characteristic time, τ, of this high pass depends only on the stationary firing rate, ν 0.

Because at zero-lag the voltage and its time derivative are uncorrelated for a stationary variance channel, , the Gaussian probability density function of the voltage dynamics factorizes over V and , where and are the respective variances. Substituting this expression into Eq (12), we obtain This expression can be computed in terms of error functions to obtain the full nonlinear dynamic response, e.g. for the Gauss-Rice LIF neuron model [19, 37, 42].

For a transparent analytical treatment of the mean channel in the fluctuation-driven regime we consider the linear response. That is, for case of weak mean input we expand, for each time t, this expression in terms of the resulting weak deviations to the ensemble mean voltage and to its derivative . To linear order, Solving the integral, one obtains the linear response in the mean signal channel, where ν 0 is the stationary firing rate attained in the absence of modulation around the mean input current, I 0, is offset by I 0 and since I 0 ≪ θ in the fluctuation-driven regime we set I 0 to 0 without loss of generality, so that . This expression can then be rewritten using only two quantities: the differential correlation time , and the size of voltage fluctuations relative to threshold, σ: = σ /θ, τ is the width of the quadratic approximation to the correlation function around zero delay. It serves as the summary timescale determined by the joint effects of all intrinsic timescales and we study in detail this dependence for the GIF model in a later section. τ thus provides a natural time unit by which to measure the rate of output spikes, ν 0, as a function of the relative voltage fluctuations, σ. ν 0τ is then interpreted as the number of spikes in a correlated window of voltage trajectory, and according to Eq (18) rises with σ, saturating for large σ at (2π)−1 < 1. Fluctuation strength is less than the voltage difference between resting and threshold for most physiological conditions, , in which case the useful bound, , holds. (Large output firing rates can nonetheless be achieved so long as the voltage correlation window, τ, is short enough to maintain ν 0τ ≪ 1.) Spike-generating voltage excursions are thus on average well-separated in time so that the produced spiking exhibits low temporal correlations.

According to Eq (16), we can then identify ν 1(ω) as the finite frequency component of its Fourier transform, where we note that our definition of ν 1(ω), Eq (16), that has the amplitude of the input modulation, A, factored out implies that A has been factored out of the voltage response. All response quantities are implicitly defined as these A-independent versions. This expression can be simplified further by pulling out the time-derivative operator. In the Fourier domain, this is just multiplication by iω so that the factors out and calculation of ν 1(ω) requires only the first two voltage moments, as any statistic derived from a stationary Gaussian process should. is the mean voltage response and the variances, and , are computed from the correlation function of the stationary unperturbed voltage correlation function, , obtained from the voltage noise spectrum δV(ω). The latter provides only the variances, and so in the space of correlation functions, only directions along which these quantities change affect the rate response [43]. The relative response can then be written We can re-express it using τ and σ, The ensemble response of a population of Gauss-Rice neurons to a small modulation in the mean input is thus simply a first-order high pass filter of the ensemble mean voltage response with characteristic frequency 1/τ, with τ defined as where we have removed σ with , obtained from Eq (18).

The relative linear rate response is then where the dependence on ν 0τ is concealed in the definition of τ. Thus, in units of τ, the high pass filter resulting from crossing the spike threshold is proportional to , with From Eq (22), we see that the characteristic frequency, 1/τ, shifts to lower values for larger output firing rates, as the prefactor, , further attenuates the low frequency response. One consequence is that the effect of the low pass voltage characteristics are made negligible by the differentiating action of the spike at high firing rate.

The dynamic gain of this complex-valued linear rate response function is its modulus, here normalized by the stationary rate, ν 0. Taking out a factor of ν 0 in Eq (9), we see that the strength of the linear term and thus the quality of the linear approximation of the response is then controlled by the size of the right hand side of Eq (23) relative to 1. The effect of this spiking filter contributes a factor that scales as when τ ≪ 1 (for a bounded range of relevant input frequencies) so the linearity assumption is better at larger values of τ, which means larger values of ν 0τ. The quality of the approximation will also depend on the size of . We also note that focusing on the linear response neglects boundedness features of the population firing rate such as its non-negativity. Nevertheless, once a voltage dynamics is specified, Eq (23) gives the explicit dependence of the dynamic gain on the underlying parameters of the single neuron model.

Derivation of the dynamic gain of a population of Gauss-Rice GIFs

In this section, we take the general result of the previous section, Eq (23), and go through its explicit calculation for a population of Gauss-Rice GIF neurons to obtain the result Eq (41). A work taking a similar approach, partly inspired by this work, though with with less intermediate analysis has recently appeared [25]. Our novel findings arise from an exhaustive characterization of the parameter dependence across the phase diagram of the voltage response, Fig 3. We calculate the current-to-voltage filter, expressing it in each of the three representations listed in Table 1, Eqs (25, 26 and 27) respectively. We show (Fig 4) how the low pass component of the filter undergoes a qualitative change from second-order low pass to first-order low pass to resonant as Q is increased. We find the voltage resonance condition, , where the resonance has a contribution from slow adaptation and from the frequency, Ω. Either can exist without the other (Fig 5). We then compute the voltage correlation function, Eq (37)), whose envelope depends on the relaxation time, τ (Fig 6). From this, the variances are calculated and an expression for the differential correlation time, τ, Eq (39) is obtained. We show a characteristic dependence on the ratio τ/τ (Fig 7). Finally, we show in Fig 8 how the stationary firing rate has unimodal dependence on the time constants, τ and τ, monotonic rise with input variance, , and monotonic decay with intrinsic frequency, Ω.

Fig 3

Regimes of the current-to-voltage transfer function.

(a) Phase diagram of the transfer function. The region of depolarizing w (low frequency amplifying, V > 1) is shown in purple and voltage resonance in green. The filter is unstable in the blue region. An intrinsic frequency exists above the dotted line, Q = 1/2. Note that there is a region with Q > 1/2 and no voltage resonance, and vice versa. The star and circle denote the example values of (ω τ, Q ) used in (b) and (c), respectively. (b) An example of the current-to-voltage filter in the case of resonance with no intrinsic frequency (τ = 10, τ = 100, g = 1.2). (c) An example of the current-to-voltage filter in the case of no voltage resonance despite the existence of an intrinsic frequency (τ = 10, τ = 5, g = 0.5). The rising and falling dashed lines in (b) and (c) denote the contributions of the high pass, , and the low pass, , respectively. Their combination forms the current-to-voltage filter, which are shown as solid lines.

Fig 4

The qualitative shape of voltage response depends on Q .

Here we classify the current-to-voltage filter shapes shown as colored solid lines in (a), (b), and (c), which show the three Q -regimes with respective examples for Q = 0.1, 0.75, 10. In each plot, the high pass component of the voltage response is shown as the colored dashed lines, one for each of three representative values of its characteristic frequency, ω τ = 102 > γ(blue), ω τ = 1(green), and ω τ = 10−2 < γ −1(red). The solid black line is the low pass component of the voltage response. For the regime shown in (a), the green case can not be achieved when w is hyperpolarizing (g > 0) and the example red case cannot be achieved because it violates the stability condition Q < ω τ.

Fig 5

A resonance frequency emerges in the voltage response in one of two ways depending on the intrinsic timescale.

For slow intrinsic current (), a response exhibiting a maximum at ω already exists at Ω = 0. For fast intrinsic current (), a resonance emerges at finite Ω, whose value converges for vanishing ω τ.

Fig 6

Emergence of oscillatory behavior in the voltage dynamics and a consequent ringing appears in the voltage correlation functions.

(a) The frequency of the ringing increases with the strength of the intrinsic current (g = 0, 1, 5 shown; τ in units of τ). (b) The envelope of the ringing widens with τ (dashed lines are ;ω τ = 1, 2, 5 shown and τ in units of ). Lines are Eq (37); dots are numerics.

Fig 7

Differential correlation time depends on intrinsic parameters.

(a)τ increases (not always monotonically) with τ. For the sake of comparison, we show τ normalized by its small-τ limiting value, , vs.τ/τ across τ/τ = 100, 101, 102, 103 104 (from blue to red) with g adjusted so their large-τ limiting value, τ/τ = 1 + g = 102. Shapes are sigmoidal for τ/τ > 1 (e.g. green to red) and include an initial dip for τ/τ < 1 (blue to green). The dot-dashed line denotes . (b) τ follows the relaxation time, τ, (the dashed line is τ = τ) and saturates at Ω−1. Colors indicate the value of τ/τ on a logarithmic scale from 10−1(red) to 101(purple). (c) General shape of τ vs. ν 0. Values in the blue region are forbidden due the maximum rate achievable in a Gauss neuron. The thick black line denotes the boundary between high and low pass.

Fig 8

Effect of model parameters on the fluctuation-driven stationary response.

The stationary firing rate, Eq (18) for I 0 ∼ 0 (a) increases monotonically with the strength of input fluctuations and (d) decreases monotonically with the intrinsic frequency. Across each of τ and τ ((b) and (c) respectively), the rate exhibits a maximum. Insets are the mean input dependent expression for the stationary response, Eq (17), valid in the regime I 0 ≪ 1. Inset color refers to the value of the parameter (σ , τ, τ and Ω) at the location of the colored dots in the main plots. Parameters were otherwise set to their default values.

Voltage solution

For arbitrary input, I(t), the system in the Fourier domain is Multiplying the first equation by (1 + iωτ ) and eliminating w(ω) one obtains so that the solution for any g > −1 is, with respect to the representation of the model by its explicit parameters (g, τ, τ),

When the neuron exhibits an intrinsic frequency, Ω, we can use and the definition of the complex eigenvalues by their real and imaginary parts, |λ ±|2 ≡ r 2 + Ω2 (see Methods), to substitute Ω into the denominator of Eq (25) after expanding: with and . r defines the relaxation time of the dynamics, τ = −r −1. Thus, in the representation of the model based on the implicit time scales, (Ω, τ, τ), the solution is expressed as

A third convenient representation consists of effective parameters, (ω , Q , τ), determining the shape of the filter where the second order low pass filter has been re-expressed using its center frequency, at which its contribution to the gain is its quality factor, (with when , and we have pulled out the broadband voltage response, V , attained in the limit ω → 0, which gives The stability constraint, g > −1 is naturally satisfied by ω > 0 and keeps V finite. With dependence on τ removed in the shape representation, we must explicitly add the stability constraint, τ > 0, which is expressed using the definition of τ in this representation, so that the stable regime corresponds to Q < ω τ.V is expressed in this shape representation as so that V > 0 is satisfied by the stability constraint.

Each of the three expressions for the voltage response filter, Eqs (25, 26 and 27), is instructive in understanding the dependence on the contained parameters. To motivate this analysis in the context of population response, we first specify the input, I(ω). An input oscillation of frequency ω 0 will produce an oscillation in the mean input expressed as . In the frequency domain, the spectrum of the mean input, , and power spectral density of the noise, δI(t), is, respectively, with noise strength, , in the latter. Because of the linearity of the dynamics, we can solve the system for mean and fluctuating input separately. In the next paragraph, we employ Eq (30) to obtain the mean voltage response, and in the following paragraph we employ Eq (31) to obtain the voltage correlation function.

Mean voltage response function

The population mean voltage response, , required for Eq (20) is obtained by inserting the expression for the mean input, Eq (30), into the voltage solution. For the remainder of the paper, we omit the factor δ(ω − ω 0) and denote the frequency of the mean input by ω. The mean response in the three representations is then Eqs (25, 26 and 27), respectively, with the I(ω) factor dropped and with an additional a factor of .

The mean voltage response is given in the filter shape representation, (τ, τ, g), by We now go through the analysis of this response using this representation. Using the gain, we constructed a diagram of its qualitative features in Q vs.ω τ (see Fig 3).

The model exhibits low frequency voltage gain amplification (V > 1) or attenuation (V < 1) depending on whether w is depolarizing (g < 0) or hyperpolarizing (g > 0), respectively. ω = 0 at g = −1 and grows with g as . Q = ω τ/2 also grows with g, generating three parameter regions of qualitatively distinct low pass filter gain shapes: ω τ < 1, 1 < ω τ < 2 and ω τ > 2. Indeed, in units of ω , the shape of the current-to-voltage filter depends only on τ and τ, and so in the next paragraphs and with reference to Fig 4, we describe this 2D parameter space completely by considering qualitative differences in the full filter shape across ω τ in each of three distinct regimes of ω τ. Note that relatively slow and fast intrinsic dynamics is obtained when Q is less than or greater than , respectively.

For ω τ < 1 (see Fig 4a), the low-pass gain contribution can be factored into a contribution arising from two first order low pass filters, where γ = γ(Q )≥1 is the solution to Q = γ/(γ 2 + 1). The low pass gain thus begins falling as ω −2 after ω /γ and then as ω −4 after γω . The intermediate region, ω/ω ∈ (γ −1, γ), is given by the inequality Q < γ/(γ 2 + 1) and disappears as Q approaches 1/2 where γ and ω τ approach 1. The region of depolarizing w (g < 0) shown in Fig 3 satisfies in this representation, whose solution in ω τ is also the range (γ −1, γ). Thus, response shapes in this intermediate region (see Fig 4b) are only achievable by depolarizing w, and w must be depolarizing for any response exhibiting such shapes. Consequently, the three qualitatively distinct shapes of the current-to-voltage filter for ω τ < 1 are determined by the location of ω τ relative to 1/γ and γ, with the middle regime, (γ −1, γ), only achievable for depolarizing w. For ω τ > γ, the filter first rises with ω after 1/τ, is flattened at ω /γ, and then falls after γω . The result is an intermediate, raised plateau of width (γ − γ −1)ω . The condition for this voltage resonance is or in terms of Q , Q > (2 + ω τ)−1/2. For 1/γ < ω τ < γ, the response attenuates first and so the plateau is now an intermediate, downward step of width (γ − 1)ω . For ω τ < 1/γ, there is only low pass behavior and the high pass only acts to pull up the ω −4-falloff up to a ω −2-falloff. As ω τ approaches 1 from below, γ also approaches 1, and the qualitatively distinct region between ω /γ and γω shrinks as the two roots coalesce into one and the low pass expression forms a perfect square. In the case that ω τ > γ, this leave a well-defined maximum located just before ω . The slight offset arises simply because the second order low pass begins falling significantly before ω at ω τ = 1.

For 1 < ω τ < 2 (see Fig 4b), the impact of the high-pass on the shape of the filter is determined simply by whether its characteristic frequency is above or below ω . For ω τ > 1, the plateau existing for ω τ < 1 becomes a flat-topped peak in the gain with a maximum again slightly lower than ω . Otherwise, the behavior is low pass. Note that ω τ > 1 is also where the intrinsic frequency exists. However, this property does not contribute to a resonance until ω τ > 2. Indeed, the resonance here, as in the regime ω τ < 1, arises solely from a high pass attenuation of low frequencies sculpting a peak from a low pass, and comes alongside a region, Q < (2 + ω τ)−1/2, that lacks resonance. This latter region is upperbounded in general by , and specifically for stable filters by so that above these values of Q all filters are voltage resonant.

For ω τ > 2 (see Fig 4c), by definition a resonant peak emerges in the low pass filter. If ω τ < 1, this contributes a de novo resonance in the current-to-voltage filter located near ω . Otherwise, it simply acts to sharpen the existing resonance that appears progressively over 0 < ω τ < 1, and again with a peak slightly to the left of ω .

Of the two mechanisms for resonance just described, the contribution of the first ‘sculpting’ mechanism leads to a linear increase in the response height and input frequency range of elevated response with τ, i.e. with the slowness of the intrinsic dynamics, for the reason that the low frequency amplification continues over a broader range the further ω /γ and 1/τ are apart. This amplification in the relative response is actually over-compensated by a broadband attenuation with τ, so that the actual effect is the carving out of a resonant peak using adaptation, i.e. a low frequency attenuation of an otherwise low pass filter.

The second low-pass resonance mechanism emerges in the expression when the low pass filter exhibits a maximum, which itself emerges when the two low pass characteristic times of the low pass coalesce. From the point of the view of the voltage dynamics, this occurs from a sufficiently strong and negative feedback interaction between v and w, whose timescales are sufficiently similar so that the delayed feedback is constructive. In the time domain voltage solution, this occurs when the two eigenvectors align. The height of the resonant response grows linearly with τ (with range of elevated response fixed) because there is less dissipation.

These two resonance mechanisms contribute to the height of the response at ω , which is resonant by definition if it is greater than V . The condition for voltage resonance is thus and the relative ratio of their contributions is so that at a given ω τ the sculpting mechanism always contributes more gain than the intrinsic frequency mechanism. Indeed, this sculpting can exist in the absence of an intrinsic frequency (ω τ < 1), so long as the intrinsic dynamics is slow enough. Conversely, even with an intrinsic frequency (ω τ > 1), the response can lack a resonance if in addition ω τ < 2, demonstrating that an intrinsic frequency is not a sufficient condition for resonance. These two cases become apparent in a plot of the resonance frequency as a function of the intrinsic frequency (Fig 5), where we observe x- and y-intercepts because of the preexisting or absent resonance, respectively. The location of the maximum converges to for ω τ ≫ 1, which itself converges to Ω for Ωτ ≫ 1. For smaller values of ω τ, the location converges to a value slightly larger than ω .

We will make use of the representation of the current-to-voltage filter in terms of (ω , τ, Q ) to understand the full response. What is left to calculate, however, is the voltage correlation function.

Voltage correlation function and the variances, and

For the correlation function of V, we perform the calculation in the implicit representation and add Eq (31) to the modulus squared of Eq (26), The auto-correlation thus requires computing an inverse Fourier transform integral of the form The result is where λ ± are the eigenvalues of the voltage dynamics, Eq (57), and the units are [Time 3]. The correlation has two components, one decaying with τ and the other with τ. The first component is strongly suppressed for τ/τ ≪ 1. The second component exhibits damped oscillations within the exponential envelope with frequency Ω. Examples are shown in Fig 6. for increasing g and τ. Note that variation in g affects the width of the function around 0-delay while it is fixed over a variation in τ. These results were checked against numerical autocorrelation functions computed from the voltage time series output of the numerically implemented model. The correspondence is excellent.

In the model representation, the variance of the voltage and that of the time derivative of the voltage are given by where, for notational convenience, we have defined and τ = τ/(1 + g) (Eq (2)) is the maximum speed over τ of the voltage kinetics, approached when τ ≪ τ by the tonic conductance change induced by w. For the LIF (g = 0), τ = τ, α = α = 1 and the variances simplify to from which the differential correlation time τ for the LIF can be read off as . We consider this quantity more generally in the next paragraph. In the intrinsic representation, the variances can be written as where |λ ±|2 = r 2 + Ω2, and r < 0 ensures that the values are positive. Note that the only difference between the expression for the two variances is a factor of and a factor of 1/|λ ±|2. In all representations, the influence of intrinsic kinetics set by τ is negligible when τ is near the input timescale, τ. Even then, w affects the variances via g or Ω.

The differential correlation time and the stationary response

From the correlation function providing the variances, the differential correlation time is calculated with . The Gauss-Rice GIF differential correlation time for the model representation is where the limiting value of τ for τ smaller (larger) than all other timescales is, respectively The ratio of slow and fast limiting values is , so that τ increases over the full range of τ/τ. In particular, the curves of Eq (39) have a characteristic shape for the non-trivial (g ≠ 0) cases. We focus on the hyperpolarizing case. In the left panels of Fig 7, we plot some example shapes of τ/τ vs. τ/τ over a range of τ < τ. Referring to that figure, for τ < τ, the curves monotonically interpolate between the limiting values, with the abscissa value at half-maximum increasing linearly with τ/τ. With τ/τ increasing from 0, τ first drops from to a minimum (whose depth grows with g) and then rises into a -scaling regime around τ/τ = 1, where it passes through the same value as that attained in the limitτ, and then eventually saturates for τ/τ ≫ 1 at its maximum, . Thus, for τ/τ → ∞ the g ≠ 0 case is equivalent to the g = 0 case and we conclude that any novel features attributable to the extra degree of freedom are washed out in this limit by the relatively slow intrinsic dynamics. As discussed in the Methods, the validity of the no-reset approximation lies around τ/τ ∼ 1, implying that . When τ ∼ τ, the approximation is valid acrossτ < τ and for other τ < τ only in ranges around the value of τ > τ for which τ ∼ τ. We also find for relatively slow intrinsic dynamics that , for τ ≤ Ω−1. When Ω exists, we can write τ as where , is the center frequency analyzed in the previous section. In the implicit representation and as a function of τ (see Fig 7), τ grows faster and slower than linear for τ/τ less than or greater than 1, respectively, and passes through when τ/τ ∼ 1, finally saturating at Ω−1. Up to this saturation level, τ > τ for τ < τ, so that the condition ν 0τ ≪ 1 implies that ν 0τ ≪ 1 and the approximation to reset dynamics is valid. In the case τ > τ, the range of τ over which the ν 0τ ≪ 1 validity constraint is not already covered by the ν 0τ ≪ 1 built-in constraint is centered around τ = Ω−1 and grows in size with τ/τ.

Next, we compute the stationary firing rate of the neuron model Eq (1) as a function of the two input parameters and the two intrinsic parameters. It is shown in Fig 8. We focus on the parameter dependence at I 0 = 0. The model’s stationary response to increased input noise exhibits a cross-over from silence to linear growth around σ ∼ θ, simply due to the higher propensity of threshold crossings. In subsequent analyses in this paper, we explore the parameter dependence at fixed stationary output firing rate by adjusting the input variance accordingly (see Methods for this mapping). The rate dependence at I 0 is similar for both τ and τ, growing from zero at vanishing time constants to a maximum located just below the membrane time constant. While the rate decays with increasing τ, it seems to saturate and even rise for slowly with τ for τ > τ. The stronger the flow of the dynamics around the resting state at I 0, the more the voltage fluctuations are dampened so that the the firing rate decreases with Ω. As for the I 0-dependence, we see that all curves rise monotonically simply because the average voltage moves closer to the threshold.

Expression for the complex response function

With the variances and the mean voltage response in hand, we can write down the complex linear frequency response, This biquad filter is composed of a two-step cascade of a combined 1st-order high-pass and 2nd-order low-pass current-to-voltage filter followed by a first-order high-pass voltage-to-population firing rate filter. In the remaining part of the paper, we analyze the properties of this filter.

Analysis of the dynamic gain function of a GIF ensemble

In this section, we characterize the qualitative features of the response function, Eq (41), again with a focus on completeness. We first show that the high and low input frequency limits of the response constrain the parameter sets that can achieve high and low pass behavior and we give an expression, Eq (45), of the critical stationary rate separating these two regions in terms of the other parameters. We then reparametrize the expression for the response, Eq (47), using the height of the response at its center frequency, ν and high frequency limit, ν ∞, both relative to its low frequency limit. The two-dimensional shape parameter space give responses with a peak, dip or step at ω whose width varies with Q . The additional high or low pass nature of the filter give six classes of filter shape. The constraint of stable voltage dynamics restricts the area accessible to the model to ν ≥ Q (1 + ν ∞).

The ω → 0 and ω → ∞ limits simply determine a high/low pass criterion

The matched order between the high and low pass filter components of Eq (41) implies that there are finite limiting values of the dynamic gain at low and high input frequencies, with the size of ν relative to ν , . We note that both ν and ν can be written without explicit dependence on the intrinsic timescale, it influences the limiting values only by setting the value of τ in the way demonstrated in the previous section.

ν scales the above filter shapes up or down and itself scales down linearly with τ/τ and thus with g. The boundary in the parameter space between low and high pass is defined implicitly by ν ∞ = 1 providing the simple criterion for low or high pass behavior as whether τ is below or above τ respectively. The high pass behaviour for large g or Q is not due to an increase in ν (in which g does not appear) but in fact a consequence of the low frequency attenuation. Recalling that the approximation to a hard threshold keeps the response flat to arbitrarily high frequencies, while in fact it eventually decays (beyond f , as discussed in the Methods section), the high pass case here implies a large elevated high frequency band up to this cut-off, while the low pass condition implies a large intermediate downward step. The low/high pass criterion implies a critical relative variance , and in turn the critical output firing rate, at which the response changes from low to high pass. Both of these values are intrinsic properties of this model whose dependence on the input relies only on the units of time taken. For τ ≪ τ, diverges as . For τ ≫ τ, it falls off as . In Fig 7c, we plot τ as a function of ν 0. One can now use the plot in this figure to determine the high or low pass behavior for a given τ/τ and ν 0. For example, when τ < τ (attained for instance with small τ and large τ/τ), there is only low-pass behavior due to the divergence of 1/2πτ . The high pass region nevertheless grows quickly with τ > τ.

The low and high input frequency limits become independent of τ when time is expressed in those units. Nevertheless, we can still write the critical condition independent of τ when expressing time in units of τ by combining Eqs (42) and (43) and eliminating τ altogether by substituting in the expression for (cf. Eq (44)) to get the high-low pass condition explicitly and solely in terms of the four timescales: τ, τ, τ, and 1/ν 0 (the latter value chosen by setting σ appropriately using Eq (60)). Setting any three of these determines the critical value of the remaining one above, across which the model changes from high to low pass behavior. For example, when time is measured in units of τ, we have For g = 0, we have α = 1 and τ = τ and this reduces simply to , the expression presumably underlying results for the LIF in [42].

As for the limiting behavior of the phase response, Φ(ω), the model gives zero delay for both high and low frequencies. At low frequencies, this is because the input changes slowly so that the model dynamics can directly follow the oscillation. At high frequencies, the return of the lag to 0, just like the flat high-frequency gain, is an artifact associated with the hard threshold.

There are six qualitatively distinct filter shapes

When g = 0 (LIF), the filter, Eq (41), simply reduces to single order. The intermediate behavior is then only the respective monotonic decay or rise beginning and ending around the smaller and larger of the two characteristic frequencies.

For g ≠ 0, the voltage modulation by the current, w, comes into play. To analyze the effect of the high pass voltage-to-spiking filter on the current-to-voltage filter we employ a similar exhaustive characterization as was done above in the analysis of the current-to-voltage filter, i.e. by going through all the cases arising from distinct orderings of the characteristic times of the components of the combined filter. The ordering can give simple information about the filter shape. For instance, any contribution of the voltage-to-spiking filter to the qualitative behavior of the complete filter beyond just low or high pass requires that 1/τ be no larger than either ω or 1/τ. Otherwise, the only effect of the spiking is to flatten the high frequency response beyond 1/τ. In general, however, there are many possible shapes. To further facilitate the classification of these shapes, we present a single parameter space representation in which they are all simply mapped.

For this general case, we can introduce the relative quality factor for the full filter, ν : = |ν 1(ω )|/ν . The response then depends on the five shape features, ν , ν , ω , Q , and ν . Denoting ξ = τ/τ, so that and , we can re-express the response function as with dynamic gain When ω = ω , , which implicitly defines ξ in terms of ν , Q and ν and closes the representation. Indeed, with time in units of and gain values relative to ν , the shape of the filter depends only on this triplet: each of the six regions in (ν ∞, ν )-space defined by the boundaries ν ∞ = 1, ν = 1, and ν ∞ = ν provides filters of a qualitatively similar class (see Fig 9).

Fig 9

The 6 distinct filter shapes in (ν , ν ∞)-space.

(c, f) Region 1–6 denote the regions exhibiting qualitatively similar filter shapes. E.g. spiking resonance is by definition region 1 and 2. Not all of these six regions are accessible for a given Q . Colored lines (blue to red) represent the Q -dependent boundary below which filter shapes are forbidden because of unstable dynamics. We note that ν = Q . An intrinsic frequency exists in region above the Q = 1/2 boundary. A voltage resonance exists in the region above the Q = 1 boundary. We show the accessible subset of corresponding filter shapes at representative positions within the regions (located at and ) and at the border between regions (located at ν , ν ∞ = 10−1.5, 100, 101.5). (f) Same type of plot as (c), but for the phase response. π/2 and −π/2 are shown as top and bottom bounding dashed lines for the set of phase responses at each location. The gain and phase for the position denoted by the circle are shown in (a) and (c), and for the star in (b) and (e), respectively.

In particular, depending on the region there is a peak, dip or step at ω whose width varies with Q . The additional high or low pass nature of the filter gives the six classes of filter shape.

While the possible shapes are simply represented in this space, the constraints are no longer represented in a plane since they depend additionally on Q . We now dissect the effects of the stability and voltage resonance constraint on determining which filter shapes are allowed where. A main conclusion that can be drawn is that a lower bound for accessible filters is ν = Q (1 + ν ∞) (which for different Q are shown as the colored lines in Fig 9).

With reference to Fig 10, the stability constraint, Q < ω τ, translates into with the correct root of ξ given by the values of Q , ν ∞, ν , and ω .

Fig 10

The accessible region of filter shapes depends on Q and the relative speed of spiking to intrinsic dynamics ξ = τ/τ.

The purple region marks the region of voltage resonant filters. This region is contained in the red region of stable filters, whose lower bound moves to larger ν with Q . For relatively slow intrinsic spiking (a, b, c), there are regions of non-spiking resonant(ν ∞ > ν ), but voltage resonant filters. Filters for relatively fast intrinsic dynamics (d, e, f) only exist as high pass resonant filters for large Q . (Left to right: , 1.1. Top row: ξ = 10. Bottom row: ξ = 0.1).

Which root can also be checked by which of τ and τ is larger. This constraint breaks into branches when combined with the other constraints.

For ξ < 1 so that the intrinsic dynamics is faster than the spiking dynamics, the region exhibiting stable filters is constrained to a sliver, , with an additional constraint on the lower bound, , so that stable filters only exist for and . For values of ν ∞ and ν increasing from this lower bound point, the accessible region forms a band whose vertical thickness grows with ν ∞ and it extends out parallel with the line ν = ν ∞ for large ν ∞. For increasing Q , the accessible region shifts right and up so that the band is eventually contained in ν > ν ∞ and ν > 1 region, i.e. only high pass, resonating filters are allowed.

For ξ > 1 so that the spiking is faster than the intrinsic dynamics, the region exhibiting stable filters has no upper bound in ν . The lower bound is ν > Q (1 + ν ∞) when and when . The latter bound differs significantly from Q (1 + ν ∞) when Q > 1/2.

The voltage resonance condition can also be mapped to this space by replacing ω τ by giving . For both roots of ξ, all stable filters are voltage resonant when .

For ξ < 1, and the voltage resonant filters exist at large ν ∞ only for ν < ν ∞, i.e. only for non-spiking resonant filters, possible because the high pass limit is brought up by the additional high pass filter above the peak of the resonance, e.g. Fig 11. Conversely, the spiking resonant filters here lack a voltage resonance because the spiking resonance arises not from the voltage resonance but from the lower frequency amplification due to the high pass spiking filter.

Fig 11

An example of filter shaping: attenuation at high frequencies uncovers an amplified band of intermediate frequencies.

(a) The shape space representation showing the region of accessible filters (white) for Q = 0.1. The blue regions exhibit unstable filters. Filters obtained from points above the thick black line are spiking resonant. Filters obtained from points above the black dashed line are voltage resonant. The arrow illustrates a path in shape space along which ν ∞ is decreased. (b) and (c) show the beginning and end filters along the path in (a). For (b) and (c), blue dashed lines are the high and low pass components of the current-to-voltage filter, which itself is shown in solid blue. Shown in red is the voltage-to-spiking filter which combined with the current-to-voltage filter gives the full filter, shown in black.

For ξ > 1, and Q decreasing from , the lower bound to the voltage resonance region interpolates across ν ∞ from the line , which rapidly approaches ν = 1 as Q is increased, to the lower bound of the region of stable filters, ν = Q (1 + ν ∞). Thus, stable filters exhibit a voltage resonance when ν 1, independent of Q . The absence of a spiking resonance, ν < ν ∞, however holds over a large sub region of these stable, ξ > 1, and voltage-resonant filters, for same reason as in ξ < 1 that the high pass limit is brought up by the additional high pass filter above the peak of the resonance, thus covering it.

For Q < 1/2, the depolarization condition, γ −2 < ξν < γ 2, also excludes some regions for hyperpolarizing w (see Fig 9).

The phase response across this representation is shown in Fig 9f. We find 0 lag when ω = ω so that the input and the response are synchronous. For the spiking resonance region, we always find a delay for slower and an advance for faster input frequencies. For non-resonant cases, it is possible to observe delays or advances for both faster and slower input frequencies.

Discussion

A neuron’s dynamic gain constrains its signal processing capabilities. Our analysis provides the first complete analysis of an expression for dynamic gain of a resonator neuron model. The level-crossing approach used here has been previously applied to 1D models to study correlation gain [23, 42, 44, 57], dynamic response [19, 37, 42], and Spike-Triggered Averaged stimulus and variance [23, 42]. Consistent with conditions for the validity of the approach [19], experiments have directly demonstrated that Gauss-Rice neurons can provide a surprisingly accurate description of cortical neurons [44, 58]. We find that the space of gain functions contains six types, two of which are resonant. The height of a resonant response is strictly dominated by intrinsic adaptation, while its sharpness is controlled by the strength of the subthreshold resonance. In particular, sharper peaks arise for higher intrinsic frequencies. We determined the parameter region where an intrinsic frequency exists and where subthreshold and spiking resonance are exhibited. We find that all possible combinations of the presence or absence of these three features have finite volume in parameter space. We expect profitable applications of our results to the study of the connection between intrinsic properties and population oscillations.

Model limitations

Neuron models with hard-thresholds, such as the LIF and GIF, have been unexpectedly successful in modeling cortical neurons [58]. They are obtained from more complex models by a series of reductions.

In Methods, we gave a rationale for the reduction to a no-reset, hard-threshold model, where we state the additional limitations imposed by lifting the voltage reset. First, these models do not apply to mean-driven situations and so do not cover phenomena such as the masking of a subthreshold resonance by a resonance at the firing rate [29]. Second, to avoid extremely bursty spike patterns, we extend previous work [19] and argue that the correlation time of the input, τ, and the correlation time of the voltage statistics, τ, can not be too different. This precludes analysis involving white current noise but implies that satisfaction depends additionally on intrinsic parameters through their dependence on τ. For example, since τ ≤ τ, the rough validity condition 1 ∼ τ/τ ≲ τ/τ so that the timescale of the input fluctuations, τ, should not be much slower than the membrane time constant, τ. Third, for correspondence with threshold models the voltage relaxation time, τ, should fall within the average inter-spike interval, ν 0τ ≪ 1. Last, these models should only be considered in the irregular firing regime, ν 0τ ≪ 1. We found that τ ≤ τ for τ > τ, so that this last constraint is in fact implied by the combination of the second and third.

To verify the validity of the no-reset model within the prescribed range, we made a direct, quantitative comparison to a canonical model with an active-spike generating mechanism. The dynamic gain of the two models coincides up to the high frequency limit, f , beyond which the low pass effect of the finite action potential rapidness dominates. Thus, all of the 6 distinct types of response shapes are altered by additional low pass behavior at high frequencies. For a previously used value of the rapidness, the intermediate frequency behavior is affected, while for a higher, and perhaps more accurate value it is not, and the artificially flat high frequency response is brought down by the realistic finite onset rapidness. In summary, these results show that the simplification to a no-reset, hard threshold is an adequate approximation when response features are slower than the speed of action potential onset.

A topic of related future work regards the possibility of accelerated kinetics of auxiliary currents during a spike [59]. To study such a scenario, one could numerically compute the gain for a model where the auxiliary current, w, undergoes a jump at spike times.

In this study of the Gauss-Rice GIF neuron and a previous on the Gauss-Rice LIF [42], exponentially-correlated Gaussian noise was used as an example of a Gaussian input statistics with non-trivial temporal correlations. These input statistics will not in general produce self-consistent firing statistics. It is therefore important to note that the approach to the linear response taken here admits arbitrary temporal correlations in the input, so long as their effect on the short-delay features of the temporal correlation of the voltage can be calculated, since that is what determines τ and thus the effect of temporal correlations on the response properties. We also note that since the voltage correlation affects the response properties only through τ, there is an equivalence class structure over the space of input correlation functions based on how they influence τ.

Relation to previous work on Type II membrane excitability

Excitable membranes are classified by the type of bifurcation that they undergo from resting to spiking, with Type I and II referring to super and sub critical Hopf bifurcation, respectively. The respective set of eigenvalues around the resting state are real and complex, with the imaginary part of the latter providing an intrinsic frequency. In this case, the voltage impulse response exhibits decaying oscillations and the voltage response function can exhibit a resonant peak near the intrinsic frequency. The mean-driven stationary spiking response rises continuously from 0 for Type I while firing in Type II neurons begins only at a finite frequency. The dynamic gain of the spiking response of Type II neurons can exhibit a superthreshold resonance arising from such subthrehsold resonance.

Frequency-sweeping ZAP input currents have revealed resonant responses from neurons in the inferior olive [60, 61], thalamus [62], hippocampus [63], and cortex [64]. Consistent with the type classification, these cells often display Type II membrane excitability properties such as subthreshold oscillations with power at similar frequencies as the spiking resonance (for a review, see ref. [7]). Type II stationary spiking responses have been measured in cortical interneurons [65]. Direct measurements of the dynamic gain of resonator neurons are lacking, however. Moreover, these existing measurements used the mean input to drive the neurons to spike. Resonator response properties in the in vivo fluctuation-driven regime remain unmeasured.

Numerical simulations of resonator models containing the minimally required currents can nevertheless reproduce the peaked voltage and ZAP response and bimodal ISI distributions in both mean and fluctuation-driven regimes [66-68]. Inspired in part by the research presented here, Tchumatchenko and Clopath [25] used similar methods as those used here on excitatory and inhibitory GIF networks where they investigated the role of subthreshold resonance and electrical synapses on the emergence of network oscillations for a particular choice of model parameters, in which they also confirm the correspondence between the response properties with and without voltage reset. The remaining few analytical results for the stationary and linear response have so far been restricted to the long intrinsic time constant limit, τ ≫ 1 [22, 29]. In this paper, we are able to obtain exact results for the stationary and linear response for all values of τ, something not possible in ref. [22] due to the difficulty of the analytics of the Fokker-Planck approach used there. For large τ and the fluctuation-driven regime, our results qualitatively match their high noise results, where σ ∼ 0.1 − 1. Since Gauss-Rice models apply only to the fluctuation-driven regime, there is no meaningful mean-driven, deterministic limit attained in the limit of vanishing noise strength with which to compare to the mean-driven results of [22], such as the shift in the resonant frequency with increasing, small amounts of noise. Their phase diagram of subthreshold behavior is essentially the same as ours, up to reparametrization. We also note that the low frequency limit will differ slightly between the models due to the slightly differing slopes of their fI-curves. These small quantitative discrepancies between idealized models should not, however, be emphasized over their ability to provide a qualitative explanation of the phenomena. Finally, we note the GIF Spike-Triggered Averaged can be obtained from our expression for dynamic gain. It has also been computed through other methods [69].

Uses of the dynamical response in the theory of recurrent networks

Explicit expressions for the linear response, such as Eq (41) obtained above, are essential ingredients for the analysis of the collective states in recurrent networks. First, they are the key quantity in the evaluation of population stability [21]. The dynamics of the population firing rate linearized around one of its fixed points is defined by the linear response function. Second, knowledge of the response function additionally reveals the correlation gain in the mapping of input current correlations to output spiking correlations. Recurrent networks exhibiting such gain will generate self-consistent patterns of inter-neuron correlations [47, 49, 70]. In the Gauss-Rice approach used here, the linear response providing the population stability and correlation gain is tractable for arbitrary Gaussian input current. Many networks generate such input statistics, most prominently balanced networks [71, 72]. We expect that the correlation gain and population firing rate stability of these networks can be theoretically investigated using the expressions for the linear response derived here.

One target application area is in understanding the connection between circuit oscillations and single cell excitability. Subthreshold resonance is often neglected in modeling studies of the PING and ING mechanisms for population oscillations [73]. This is despite the ample suggestive evidence of phase locking between subthreshold oscillations and gamma band population oscillations [7]. This connection has been studied in the olfactory bulb where mitral cells display a host of resonator properties such as subthreshold oscillations [5, 6], rebound spikes [74], and Type II phase resetting curves [75]. The role of this resonance in sustaining the population oscillation has not been directly assessed in detailed network models of resonating mitral cells [76], though it should play a role in either of two existing hypotheses for the origin of the oscillations [77]. Combining subthreshold and PING mechanisms has been studied in other contexts [78].

The demonstrated subthreshold resonance in inhibitory interneurons in cortex likely also contributes to the population oscillation observed there (as suggested by the numerical results of [79] and [78]) and could be investigated using the expression for dynamic gain that we provide. A first of such studies inspired by an unpublished version of the work presented here has already appeared [25], where the Gauss-Rice GIF response gain was also derived.

Finally, an ad hoc dynamic response filter of the same form as the one derived here [80] has been successful in modeling responses of cortical neurons (personal communication O. Shriki). The explicit dependence in our derived expression on the parameters of an underlying neuron model can be used to extend those studies, in particular, by inferring from the fitted values the properties of the intrinsic dynamics of the measured cells.

Response properties depend on the differential correlation time

The differential correlation time, τ, was used in a variety of ways throughout this paper.

First, it appeared in expressions for other important quantities in the theory. It appears most prominently in our expression for the fluctuation-driven voltage autocorrelation function for exponentially-correlated Gaussian input current. The result for a Type II GIF, Eq (37), gives exponentially enveloped, oscillatory decay, with a decay constant equal to the relaxation time of the model and oscillation frequency given by the intrinsic frequency, Ω. Despite these oscillations, we find that the dynamic gain depends only on the initial falloff behavior away from 0-delay, a feature that can be shown to define, τ. From the perspective of the response then, voltage correlation functions differ only insofar as they exhibit different τ. The characteristic time, τ, and thus also the attenuation of the spiking filter scales linearly with τ, influencing the high or low pass nature of the filter accordingly.

Second, τ appears in the validity conditions for the model. Namely, the range of valid firing rates for all Gauss-Rice neurons must lie below .

Third, model parameters such as the intrinsic time scale, τ, have an effect on dynamic response features, such as the high and low frequency limits, only through τ. The analysis of their effect on τ provides insight as to their role in sculpting the response properties. In Fig 7c for example, τ grows with τ, and for large τ saturates at τ → τ = τ, so that τ can only be made shorter, not longer, than the membrane time constant, τ, by intrinsic and synaptic current parameters.

The central role of τ could be tested by applying a variety of input correlation functions with significant differences only away from the fall-off at 0-delay so that they provide the same τ. Our model predicts no significant change in the response properties. Such a large number of experiments could be performed by methods of high-throughput electrophysiology currently under development.

The six filter types of the Gauss-Rice GIF

We re-expressed the response expression, Eq (20), using the center and high frequency response relative to the low frequency response, ν and ν ∞ respectively. We find six qualitatively distinct filter shapes distributed around (1, 1) in the (ν ∞, ν ) plane, with the value of Q determining which of the six are accessible. Depending on the region there is a peak, dip or step at ω whose width is determined by Q . We summarize below the constraints on the accessible shapes set by Q . For Q < 1/2, all six filters shapes are possible for fast relative spiking (τ < τ). There are no high pass resonating shapes in the limit of vanishing Q for slow relative spiking (τ > τ). For all accessible shapes have elevated response at the center frequency, ν > 1. For Q > 1, all allowed filter shapes are resonating, that is ν > ν ∞. There are no low pass resonating filters for slow relative spiking and so a sharp resonance, i.e. a high Q , is only possible when the overall filter is high pass.

Neither voltage nor spiking resonance strictly imply the other in this model. First, there can be voltage resonance with no spiking resonance because the spiking high pass pulls up the response in the high input frequency range above the elevated response around the intermediate-range resonant input frequency. The high frequency limitation of the approach (e.g. Fig 12) implies that the elevated response extends up to the speed of the action potential, leaving a broad resonant band at high input frequencies. Second, there can be spiking resonance with no voltage resonance because of a low frequency attenuation by the spiking high pass filter of a low pass current-to-voltage filter.

Fig 12

Correspondence of response between analytical result of no-reset model (blue line) and the numerical result of its EIF version (black circles).

The correspondence holds up to a high frequency cut-off, f (Eq (54)), due to finite rise time of action potential controlled by Δ = 0.35, 0.035. The EIF-version was simulated with V = 1.15, 3, and V = 0.8, −1 (the latter was adjusted to keep ν 0 = 2Hz fixed). The black dashed lines correspond to the high frequency limit of the response of the EIF-type model (Eq (53)). The no reset model had the default parameters.

In addition, neither voltage nor intrinsic resonance strictly imply the other. First, the existence of an intrinsic frequency does not imply voltage resonance in general because the response at ω where Ω becomes finite is Q = 1/2 and is thus still attenuated relative to the response at low input frequencies. This response only becomes resonant at Q = 1. Second, there can be a voltage resonance with no intrinsic resonance for the same reason that a high pass with low characteristic frequency (this time from relatively slow intrinsic dynamics) can sculpt a peak from the low pass component of the full filter.

Finally, we found that the strength of the spiking resonance (∼ν ) is composed of a contribution from the intrinsic timescale, τ and from the intrinsic frequency, Ω. Nevertheless, ν is dominated by the attenuation at low input frequencies associated with the high pass effect of large τ, while the unique effect of Ω is to sharpen this resonance.

The cascade representation of the dynamic response

The effect of spiking in the Gauss-Rice formulation of the response is as an explicit first-order high pass filter of the voltage dynamics (see Eq (20)). We note that this high pass behavior associated with spiking is distinct from that discussed in the literature as arising from sodium channel inactivation [81]. This has nothing to do with the Gauss-Rice high pass arising in this paper. In this work, we always consider the threshold fixed. Closed form expressions are thus obtained for the low frequency limit and characteristic time of this filter in terms of the parameters of the model. When the characteristic frequency is high, the filter has the effect of flattening an otherwise decaying voltage response. The flattening effect is physiologically meaningful up to frequencies at which the spike-generator cut-off appears. It thus sculpts a plateau of constant response at high frequencies that can be elevated or depressed relative to the low frequency response. On the other hand, when the characteristic frequency is low, the resulting effect is a low frequency attenuation that carves out a resonant peak. The high pass characteristics are then also dependent on the intrinsic timescales.

Methods

Reduction from conductance-based models

Here, we detail how one arrives at a model like the one used in this paper from simplifications made to the synaptic, subthreshold, spiking, and spiking reset currents of a Hodgkin-Huxley type neuron model for the dynamics of the somatic transmembrane voltage potential, V (here measured in mV), where C is the membrane capacitance, I is the sum of all membrane currents and I is the total synaptic current arriving at the soma. Our exposition of the reductions to synaptic and subthreshold currents is standard. To the exposition of the reductions of spiking currents we add analysis determining the high frequency limit, f , below which the approximation to a hard threshold is valid. To the exposition of the reduction of reset currents, we add more detailed consideration of the mechanisms through which the no reset approximation breaks down.

Synaptic current

I contributes current terms of the form g (t) (V − E), where E is the reversal potential for the synapse type and g (t) is the time-varying, synaptic input conductance for that class of synapse whose time course is determined by presynaptic activity. For a neuron embedded in a large, recurrently-connected population, this presynaptic activity arises from both the recurrent presynaptic pool of units (numbering K ≫ 1 on average) and any external drive. In networks with sufficient dissipation, the external drive acts to maintain ongoing activity. The measured activity of networks in this regime is asynchronous and irregular and can be achieved robustly in models by an approximate -scaling of the recurrent coupling strength, J. This scaling choice has the effect of balancing in the temporal average the net excitatory and inhibitory input to a cell, leaving fluctuations to drive spiking. In this fluctuation-driven regime, the mean-field input to a single neuron resembles a continuous stochastic process. In the limits of (1) many, (2) weak, and (3) at most weakly correlated inputs, a diffusion approximation of I (t) can be made such that it obeys a Langevin equation [82, 83]. While not yet developed for the Gauss-Rice neuron approach, analytical tools for computing the response in the case of the shot noise resulting when (1) fails are appearing [84]. Strong inputs do exist in the cortex where synaptic strengths can be logarithmically distributed. Nevertheless, many strengths are weak, and we treat only (2) here. Finally, an active decorrelation in balanced networks justifies (3). Expanding I to leading order in the conductance fluctuations reduces the input to additive noise yielding the Gaussian approximation to the voltage dynamics, also known as the effective time constant approximation [84, 85]. The quality of this approximation depends on the relative difference between the reversal potential and the voltage. Somas receive input from two broad classes of synapse: excitatory ones for which the difference is large, and inhibitory ones for which the difference is smaller so that they are less well-approximated. The two types can also differ in their kinetics. While both are generally low pass, their characteristic times can be different. Their combination can thus have qualitative effects on the response [19]. We retain only a single synapse type so as to concentrate on the shaping of the filter properties by the intrinsic currents of the neuron model.

In this approximation to additive Gaussian noise, the time-dependent ensemble from which the input signal, I (t), is sampled is completely described by a variance channel carrying the dynamics of the fluctuations of the network activity, and a mean channel carrying the dynamics of the mean network activity. More complicated compound input processes described by higher order statistics offer more channels but they are negated by the diffusion approximation to a Gaussian process. The variance channel determines the fluctuations of I (t) on which rides a DC component described by the mean channel. We can thus write where the zero-mean Gaussian process δI(t) is characterized by the variance, , and correlation time, τ, of the fluctuations, both of which can in general vary in time, and is the time-dependent population mean. The population mean of a quantity, x, will be denoted by a bar so that , where k indexes the neuron. For stationary input, the time average of is due to the balance. In this paper, we consider deterministic changes in the mean channel, , produced for example by a global and time-dependent external drive. We compute the resulting frequency and phase response, and leave the analysis of the variance channel to a forthcoming work. For much of the paper, we will also remove explicit dependence of the model’s behavior on the input by setting σ for a desired output firing rate and measuring time relative to τ.

Subthreshold current

In the most simple case (no longer exactly the Hodgkin-Huxley formalism), each somatic current, I , contributes additively to I with a term of the form where g (V) is a voltage-dependent conductance, whose effect on the voltage dynamics depends on the driving force, V − E , the difference of the voltage and the reversal potential, E . g obeys kinetic equations based on channel activation whose specification is often made ad hoc to fit the data since the details of the conformational states and transitions of a neuron’s ion channels is often unknown or at least not yet well understood. Nevertheless, for voltages below the threshold for action potential initiation the voltage dynamics can be well-approximated by neglecting spike-generating currents and linearizing the dynamics of the subthreshold gating variables around the resting potential, V*. Following ref. [29], the resulting subthreshold dynamics is then given by where v = V − V* and are the linearized variables for the voltage and subthreshold gating variable, x , respectively; and are the effective membrane conductances for the leak and for x , respectively; and τ = τ(V*) is the time constant of the dynamics of w . C is the capacitance of the membrane. The w variables have dimensions of voltage. Activation and inactivation gating variables have g > 0 and g < 0, respectively. We denote the linearized voltage by V instead of v throughout the paper to better distinguish it from the firing rate, ν.

With the addition of a hard (i.e. sharp and fixed) voltage threshold and a reset rule to define the spiking dynamics, this defines the GIF class of models [29]. Among the models considered in ref. [29], the simplest has only one additional degree of freedom, with spikes occurring at upward crossings of the threshold, θ. With time in units of τ the authors multiply the voltage equation by τ/C and analyze the behavior as a function of two dimensionless model parameters, α = g τ/C and β = g τ/C , upon which the qualitative shape of the current-to-voltage filter for white noise input depends.

We consider correlated noise input that introduces an additional time scale which serves as a more natural time unit. We are also interested in the explicit dependence on τ. Thus, we retain both of the timescales of the neuron model, τ and τ. We then parametrize our model using the relative conductance g = β/α = g /g , the relative membrane time constant τ/τ = α −1 = C /τ g , and the relative w timescale, τ/τ = β/τ g. Input variance is independently fixed in order to achieve a desired firing rate. We thus make a slight alteration to the model in ref. [29], We have absorbed the 1/g factor into the units of I so that all dynamic quantities are in dimensions of voltage. We keep τ > 0 by setting g > 0, that together with g > −1, this gives stable voltage dynamics. This model is the same as the one stated at the beginning of the Results section, Eq (1).

The approximation to a hard threshold from a set of spike-generating currents that are in principle contained in I but are not considered explicitly in [29] involves some assumptions and approximations that have since been nicely formalized in [24] and so we include them in the following section.

Spike-activating current

The formulation of spike-activating currents can be simplified using the fact that all the information that the neuron provides to downstream neurons is contained in the times of its action potentials and not their shape. Only the voltage dynamics contributing to this time is retained in the model; namely, the summed rise of voltage-gated activation of the spike-generating x , summed into a single function, ψ(V), dependent only on the voltage when its dynamics is relatively fast [24]. ψ(V) then appears as a term in the voltage dynamics and, when supralinear in V, acts as the spike-generating instability that, in the absence of superthreshold, hyperpolarizing currents, causes the voltage to diverge in finite time. These latter currents are simply omitted and the time at which the voltage has diverged is used in these models as the spike time. The socalled spike slope factor [24], Δ, is the inverse curvature of the I-V curve near threshold and sets the slope of the rise of the action potential, with smaller values giving steeper rise. Its value should be measured at the site of action potential initiation, the precise location of which is not yet known in general. An upper bound on the realistic range of Δ is, however, likely smaller than that achievable by conventional Hodgkin-Huxley-like models, even with multiple compartments [37, 39], and this speed has motivated neuron models with fast action potential onset rapidness [86].

The time between the crossing of a fixed threshold voltage, V , defined implicitly by , and the spike time vanishes quickly with , so that the further approximation to a hard threshold, i.e. for omitting ψ(V) altogether by setting the spike time at V , becomes good for Δ → 0. However, the instantaneous rise in voltage in this limiting approximation introduces artefactually fast population responses at high input frequencies, denoted by f, raising the scaling behavior to and constant for white and colored noise, respectively [20]. Nevertheless, since the discrepancy begins above some f depending on Δ, the artefact can be safely ignored by considering the shape of the response only for f < f . Conveniently, an upper bound on realistic values of Δ given by the surprisingly quick rise of real action potentials leads to a value of f well beyond the range of input frequencies over which realistic filtering timescales act. As a result, the approximation to a hard threshold does not alter the sub-spiking timescale response properties of the full model.

For concreteness, a popular choice for ψ(V) is ψ(V) = exp[(V − V )/Δ], the family of so-called exponential integrate-and-fire (EIF) models [87] for which the difference between the threshold crossing and the spike time vanishes very fast as . Its high frequency response falls off as 1/f, with a high frequency cut-off . We consider an EIF version of our model defined having an additional, superlinear term in the -equation, . We note that a similar comparison is made in [25]. The approximate upper limit of input frequencies, f , below which the no-reset approximation is valid is given implicitly by the intersection of the response of the simplified model computed in this paper and the analytical high frequency response of the EIF version of the full model, computed from an expansion of the corresponding Fokker-Planck equation in ω −1 = 1/(2πf). We choose examples where the intrinsic dynamics are slow relative to the cut-off so we use the high frequency limit result of the EIF with no additional degree of freedom calculated in [24], The high frequency limit of the Gauss-Rice GIF is Eq (43). Equating these two expressions, we obtain where g, τ, and τ are parameters defined later. We check this condition through numerical simulations of the EIF-version of the model. Instead of the heuristic constraints for choosing the integration time step dt as specified in [24], we more simply obtain the f −1 fall-off by raising the numerical voltage threshold for spiking, allowing the speed of the action potential to play a role at higher frequencies. While this gives an artifact in the phase response (not shown), the high frequency limit of the gain is correct. Two example gain functions are shown in Fig 12 for a widely used value of Δ = 3.5 mV(0.35 in our units), and a value an order of magnitude smaller, Δ = 0.35 mV(0.035 in our units). The former value gives a cut-off slow enough that it affects the resonant feature, while the latter value gives a cut-off high enough that it does not. The features of the filter in this case are thus well below f .

Notably, the LIF FP methods have been used to obtain the linear response to a piecewise linear models [33, 34]. In these works, the high frequency artifacts induced by the hard threshold are treated explicitly and removed.

Resetting current

Models that neglect the downward part of the action potential require the addition of, or have already built-in a reset voltage to which the voltage is reset after a spike. The reset makes the dynamics discontinuous and a closed form expression for the frequency response for more-than-1D models appear intractable. We forgo this reset rule in order to open up the problem for deeper analysis. With this simplification, however, come three issues that we avoid by narrowing the scope of the analysis.

First, without the reset and for the case of mean-driven activity, the mean voltage is taken into an unrealistic, super-threshold range. Thus, only fluctuation-driven activity with low, subthreshold mean input is covered by this approximation, leaving out mean-driven phenomena such as the masking of a subthreshold resonance by a resonance at the firing rate as shown, e.g. in ref. [29]. This is nevertheless the operating regime of cortical networks that we wish to study. We thus set the mean input to 0.

Second, the lack of reset produces periods of artefactually high and low firing rates for respectively small and large values of the input correlation time, τ, relative to the voltage correlation time defined here as the differential correlation time, . τ is the quadratic approximation to the voltage correlation function around 0-delay (discussed in detail in the main text). This definition precludes the use of white noise input whose correlation function is non-differentiable around 0-delay. Indeed, the fractal nature of the voltage traces when the no-reset model is driven by white noise endows the model with the problematic feature that every threshold crossing has in its neighborhood infinitely many such crossings [57]. A version of this effect explains the discrepancy between reset and non-reset dynamics even in the finite realm where τ/τ ≪ 1. In the other limit, τ/τ ≫ 1 means that the voltage stays super threshold for long spans of time and so must also be excluded. Badel compares the stationary response of the LIF with and without reset across τ, finding correspondence only in a fairly tight band around the membrane time constant, τ, from τ = 0.5τ to τ = 2τ [19]. Given that the stationary response of the LIF also deviates from more realistic models, in this paper we do not aim for exact correspondence with the LIF but rather analyze the more general and less strong condition, τ/τ ∼ 1, which reduces to a less strong version of the one Badel used for the LIF where . From the derivation of τ for the Gauss-Rice GIF exposed in the main text, the condition τ/τ ∼ 1 implies that the membrane time constant is no longer required to lie within an order of magnitude of τ but that the validity now holds around a manifold in the space of intrinsic parameters of the model.

Third, for those neurons that do exhibit reset-like dynamics, this approach can nevertheless provide a good approximation so long as the model dynamics allow for the sample paths of the voltage trajectory after a spike with and without reset to converge onto one another before the next spike occurs. The formal condition for this is ν 0τ ≪ 1, where ν 0 is the firing rate and τ is the relaxation time of the deterministic dynamics of the voltage, i.e. the negative of the largest real part of the eigenvalues of the solution to the linearized voltage dynamics. For the case of 2D linear dynamics considered in this paper, with differential matrix operator B, when r 2 > detB and when r 2 < detB, where r is the real part of the complex eigenvalue (see next paragraph for details). For relatively fast intrinsic kinetics, this constraint limits the range of parameters and output firing rates over which the no-reset model approximates reset dynamics to within some tolerance. However, we will show that, for relatively slow intrinsic kinetics, the condition τ ≲ τ holds up to a saturation level, and this together with ν 0τ ≪ 1 (a condition that all healthy Gauss-Rice neurons must satisfy) guarantees the near equivalence of reset and no-reset dynamics, independent of the other parameters. In other words, the approximation is valid in this regime if the relaxation time falls within a correlated window of voltage trajectory as this is a lower bound to the time between spikes. Indeed, for any temporally correlated dynamics, it always takes some time for the state to move some fixed amount. In this context, that effect induces an relative refractory period in reset dynamics as the state must move from reset to threshold again in order to spike. It is not absolute because this time depends on the firing rate. The same type of refractoriness emerges in non-reset dynamics as the voltage must fall back below threshold in order to pass it from below again.

Parametrization of the model

When the eigenvalues of the solution to the voltage dynamics are complex, we can re-express the denominator of Eq (25) using the intrinsic frequency, i.e. the imaginary part of the eigenvalues of the voltage solution. We first obtain the eigenvalues. For the linear matrix evolution operator the eigenvalues are obtained via the identity where as the negative reciprocal of the harmonic mean of the two time constants, τ, and where ω is the center frequency of the voltage filter. When ω τ < 1, the magnitude . When ω τ > 1, the eigenvalues are complex with as the real part. We define the imaginary part that plays the role of the intrinsic frequency, Ω > 0, via λ ± = r ± iΩ, so and now the magnitude is . We can substitute the expression for ω , obtaining the relation between g and Ω, Eq (3), where is the condition for complex eigenvalues (see Fig 1).

Obtaining the response function directly from spike times

Here we rederive the linear relationship between the vector strength and the linear response. ν 1(ω) from Eq (9) can be expressed using the complex response vector, where in the last step we use n ≈ ν 0 T, good when T is made much larger than . Taking the ensemble average, Using the decomposition of the response into its gain and phase, ν 1(ω) = |ν 1(ω)|e , the dynamic gain is thus obtained from the norm of the ensemble-averaged response vector, called the vector strength, where here we have simplified the notation by having m run over all the spikes from the full ensemble. We computed this expression using the spike times obtained directly from numerical simulations of the stochastic dynamics generated by the neuron model. We use the result to confirm the validity of the analytical gain function derived below, whose utility goes far beyond the numerical result because it provides the explicit dependence on the model parameters.

Computing the input variance for given firing rate

Rearranging the expression for ν 0 and then substituting in the σ -dependent expression for the voltage fluctuations, σ, we have When we study the model’s behavior we will use this relation to set the input variance for a chosen output firing rate so that the dimensions of the parameter space to be explored are the four time scales in the problem, (τ, τ, τ, 1/ν 0) and when Ω exists, (τ, 1/Ω, τ, 1/ν 0).

Step response

The firing rate response derived in this paper allows us to compute the response to any weak signal and we demonstrate that in this section where we derive the response to step-like input. The time-domain version of linear frequency response, ν 1(t), is the impulse response function, which when convolved with any input times series gives the corresponding response time series, where has units of [Time]−2[Current]−1. If there is an accessible frequency representation of the input, the interaction can be made in the frequency domain and then the result transformed back to the time domain, We used this definition to study the response to step-like input, I(t) = AΘ(t), with step height, A, and with frequency domain expression for the Heaviside theta function, Applying the inverse Fourier transform to the product of this with the linear frequency response gives the expression for the response. The relative response is then, with . We can express Eq (63) in terms of r and Ω, where and depend on the parameters. Taking the limit t → 0+, the relative instantaneous jump height is , consistent with the notion that higher characteristic cutoff frequencies, i.e. , imply stronger instantaneous transmission. The exponent of the subsequent decay is , providing an envelope that funnels into the relative asymptotic response, , attained in the limit t → ∞. Since the oscillation amplitude scales as 1/Ω while the asymptotic response scales with 1/Ω2, there will be a tapering envelope for Ω > 1. Within this envelope the response oscillates at the intrinsic frequency and with a phase that is explicitly dependent on the neuron parameters, as well as implicitly though τ. This function was used to calculate the step response shown in Fig 2.

A zipped code package is supplied with which the data shown in Fig 1 can be generated, outputted, and plotted.

Details are given in the readme.txt contained in the package.

(ZIP)

Click here for additional data file.