Intrinsic parameters, g and τ_w, 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_syn, and with time in units of τ_V, the structure of the phase space of the single neuron dynamics described by Eq (1) is determined by a point in the τ_w/τ_V vs. g plane, the two parameters defining the intrinsic current, w (see Fig 1). For τ_w ≪ τ_V, 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 (2) 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_int, that varies as (3) where g_crit = (τ_w − τ_V)²/4τ_wτ_V (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 τ_w. For fast τ_w, the Ω-contour density is high and the model exhibits high parameter sensitivity, while for large τ_w 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 τ_w with slope ∝ Ω² and with a slope independent of Ω for small τ_w.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

(a) Intrinsic parameter phase diagram in (τ_w/τ_V, 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 τ_w and for small τ_w). (b) When , the phase diagram can be cast in (τ_w/τ_V, Ωτ_V)-space. Iso-g lines are shown in white. (See [29] for a similar plot).

https://doi.org/10.1371/journal.pcbi.1004636.g001

Furthermore, there is a minimum relative conductance, for which a given Ω can be achieved. The minimum shifts to increasingly short τ_w 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 Ω.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Parameter groups for each dynamics.

https://doi.org/10.1371/journal.pcbi.1004636.t001

Population firing rate dynamics.

Given a population of N neurons indexed by k, in a time window, T, each one produces a spike train, (4) with n_k spikes labeled as . The average firing rate across the population in this window is (5) For stationary input, this becomes the stationary population averaged firing rate, independent of t, in the limit T → ∞, (6) 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, (7) 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_syn, 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), (8) 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 τ_I. 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. From input to ensemble response: numerics and prediction.

Model output for the default parameter set: τ_I = 1ms, σ_I = 1, τ_V = 10ms, θ = 1, τ_w = 20ms, f_int = 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.

https://doi.org/10.1371/journal.pcbi.1004636.g002

Obtaining the response from the statistics of the voltage dynamics.

To obtain ν₁(ω) analytically, we go back to the definition of ν(t) containing s_k(t), Eq (7). s_k(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], (10) 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, (11) (12) 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, τ_s, the characteristic time, τ_c, of this high pass depends only on the stationary firing rate, ν₀.

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 , (13) where and are the respective variances. Substituting this expression into Eq (12), we obtain (14) 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, (15) Solving the integral, one obtains the linear response in the mean signal channel, (16) where ν₀ is the stationary firing rate attained in the absence of modulation around the mean input current, I₀, (17) is offset by I₀ and since I₀ ≪ θ in the fluctuation-driven regime we set I₀ 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, σ: = σ_V/θ, (18) τ_s 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. τ_s thus provides a natural time unit by which to measure the rate of output spikes, ν₀, as a function of the relative voltage fluctuations, σ. ν₀τ_s 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. 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, τ_s, is short enough to maintain ν₀τ_s ≪ 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 ν₁(ω) as the finite frequency component of its Fourier transform, (19) where we note that our definition of ν₁(ω), 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 ν₁(ω) 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 (20) We can re-express it using τ_s and σ, (21) 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/τ_c, with τ_c defined as where we have removed σ with , obtained from Eq (18).

The relative linear rate response is then (22) where the dependence on ν₀τ_s is concealed in the definition of τ_c. Thus, in units of τ_s, the high pass filter resulting from crossing the spike threshold is proportional to , with From Eq (22), we see that the characteristic frequency, 1/τ_c, 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, (23) here normalized by the stationary rate, ν₀. Taking out a factor of ν₀ 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 τ_c ≪ 1 (for a bounded range of relevant input frequencies) so the linearity assumption is better at larger values of τ_c, which means larger values of ν₀τ_s. 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.

Voltage solution.

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

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

A third convenient representation consists of effective parameters, (ω_L, Q_L, τ_w), determining the shape of the filter (27) 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_low, attained in the limit ω → 0, which gives The stability constraint, g > −1 is naturally satisfied by ω_L > 0 and keeps V_low finite. With dependence on τ_V removed in the shape representation, we must explicitly add the stability constraint, τ_V > 0, which is expressed using the definition of τ_V in this representation, (28) so that the stable regime corresponds to Q_L < ω_Lτ_w.V_low is expressed in this shape representation as (29) so that V_low > 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 ω₀ 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, (30) (31) 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 δ(ω − ω₀) 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, (τ_V, τ_w, g), by (32) We now go through the analysis of this response using this representation. Using the gain, (33) we constructed a diagram of its qualitative features in Q_L vs.ω_Lτ_w (see Fig 3).

The model exhibits low frequency voltage gain amplification (V_low > 1) or attenuation (V_low < 1) depending on whether w is depolarizing (g < 0) or hyperpolarizing (g > 0), respectively. ω_L = 0 at g = −1 and grows with g as . Q_L = ω_Lτ_r/2 also grows with g, generating three parameter regions of qualitatively distinct low pass filter gain shapes: ω_Lτ_r < 1, 1 < ω_Lτ_r < 2 and ω_Lτ_r > 2. Indeed, in units of ω_L, the shape of the current-to-voltage filter depends only on τ_r and τ_w, 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 ω_Lτ_w in each of three distinct regimes of ω_Lτ_r. Note that relatively slow and fast intrinsic dynamics is obtained when Q_L is less than or greater than , respectively.

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

For 1 < ω_Lτ_r < 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 ω_L. For ω_Lτ_w > 1, the plateau existing for ω_Lτ_r < 1 becomes a flat-topped peak in the gain with a maximum again slightly lower than ω_L. Otherwise, the behavior is low pass. Note that ω_Lτ_r > 1 is also where the intrinsic frequency exists. However, this property does not contribute to a resonance until ω_Lτ_r > 2. Indeed, the resonance here, as in the regime ω_Lτ_r < 1, arises solely from a high pass attenuation of low frequencies sculpting a peak from a low pass, and comes alongside a region, Q_L < (2 + ω_Lτ_w)^−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_L all filters are voltage resonant.

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

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 τ_w, i.e. with the slowness of the intrinsic dynamics, for the reason that the low frequency amplification continues over a broader range the further ω_L/γ and 1/τ_w are apart. This amplification in the relative response is actually over-compensated by a broadband attenuation with τ_w, 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 τ_r (with range of elevated response fixed) because there is less dissipation.

These two resonance mechanisms contribute to the height of the response at ω_L, (34) which is resonant by definition if it is greater than V_low. The condition for voltage resonance is thus and the relative ratio of their contributions is so that at a given ω_Lτ_w the sculpting mechanism always contributes more gain than the intrinsic frequency mechanism. Indeed, this sculpting can exist in the absence of an intrinsic frequency (ω_Lτ_r < 1), so long as the intrinsic dynamics is slow enough. Conversely, even with an intrinsic frequency (ω_Lτ_r > 1), the response can lack a resonance if in addition ω_Lτ_w < 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 ω_Lτ_w ≫ 1, which itself converges to Ω for Ωτ_r ≫ 1. For smaller values of ω_Lτ_w, the location converges to a value slightly larger than ω_L.

We will make use of the representation of the current-to-voltage filter in terms of (ω_L, τ_w, Q_L) 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), (35) The auto-correlation thus requires computing an inverse Fourier transform integral of the form (36) The result is (37) where λ_± are the eigenvalues of the voltage dynamics, Eq (57), and the units are [Time³]. The correlation has two components, one decaying with τ_I and the other with τ_r. The first component is strongly suppressed for τ_r/τ_I ≪ 1. The second component exhibits damped oscillations within the exponential envelope with frequency Ω. Examples are shown in Fig 6. for increasing g and τ_r. Note that variation in g affects the width of the function around 0-delay while it is fixed over a variation in τ_r. 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 (38) and τ_eff = τ_V/(1 + g) (Eq (2)) is the maximum speed over τ_w of the voltage kinetics, approached when τ_w ≪ τ_V by the tonic conductance change induced by w. For the LIF (g = 0), τ_eff = τ_V, α_I = α_w = 1 and the variances simplify to from which the differential correlation time τ_s 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 |λ_±|² = r² + Ω², 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/|λ_±|². In all representations, the influence of intrinsic kinetics set by τ_w is negligible when τ_w is near the input timescale, τ_I. 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 (39) where the limiting value of τ_s for τ_w smaller (larger) than all other timescales is, respectively The ratio of slow and fast limiting values is , so that τ_s increases over the full range of τ_w/τ_I. 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 τ_s/τ_I vs. τ_w/τ_I over a range of τ_eff < τ_V. Referring to that figure, for τ_eff < τ_I, the curves monotonically interpolate between the limiting values, with the abscissa value at half-maximum increasing linearly with τ_V/τ_I. With τ_w/τ_I increasing from 0, τ_s first drops from to a minimum (whose depth grows with g) and then rises into a -scaling regime around τ_w/τ_I = 1, where it passes through the same value as that attained in the limitτ_s,fast, and then eventually saturates for τ_w/τ_I ≫ 1 at its maximum, . Thus, for τ_w/τ_I → ∞ 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 τ_s/τ_I ∼ 1, implying that . When τ_eff ∼ τ_I, the approximation is valid acrossτ_w < τ_I and for other τ_eff < τ_I only in ranges around the value of τ_w > τ_I for which τ_s ∼ τ_I. We also find for relatively slow intrinsic dynamics that , for τ_s ≤ Ω⁻¹. When Ω exists, we can write τ_s as (40) where , is the center frequency analyzed in the previous section. In the implicit representation and as a function of τ_r (see Fig 7), τ_s grows faster and slower than linear for τ_w/τ_I less than or greater than 1, respectively, and passes through when τ_w/τ_I ∼ 1, finally saturating at Ω⁻¹. Up to this saturation level, τ_s > τ_r for τ_w < τ_I, so that the condition ν₀τ_s ≪ 1 implies that ν₀τ_r ≪ 1 and the approximation to reset dynamics is valid. In the case τ_w > τ_I, the range of τ_r over which the ν₀τ_r ≪ 1 validity constraint is not already covered by the ν₀τ_s ≪ 1 built-in constraint is centered around τ_w = Ω⁻¹ and grows in size with τ_w/τ_I.

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. The model’s stationary response to increased input noise exhibits a cross-over from silence to linear growth around σ_I ∼ θ, 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₀ is similar for both τ_I and τ_w, growing from zero at vanishing time constants to a maximum located just below the membrane time constant. While the rate decays with increasing τ_I, it seems to saturate and even rise for slowly with τ_w for τ_w > τ_V. The stronger the flow of the dynamics around the resting state at I₀, the more the voltage fluctuations are dampened so that the the firing rate decreases with Ω. As for the I₀-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, (41) 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.

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, (42) (43) with the size of ν_high relative to ν_low, . We note that both ν_low and ν_high can be written without explicit dependence on the intrinsic timescale, it influences the limiting values only by setting the value of τ_s in the way demonstrated in the previous section.

ν_low scales the above filter shapes up or down and itself scales down linearly with τ_eff/τ_V 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 τ_c is below or above τ_eff respectively. The high pass behaviour for large g or Q_L is not due to an increase in ν_high (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_limit, 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, (44) 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 τ_s ≪ τ_eff, diverges as . For τ_s ≫ τ_eff, it falls off as . In Fig 7c, we plot τ_s as a function of ν₀. One can now use the plot in this figure to determine the high or low pass behavior for a given τ_w/τ_I and ν₀. For example, when τ_s < τ_eff (attained for instance with small τ_w and large τ_V/τ_eff), there is only low-pass behavior due to the divergence of 1/2πτ_s. The high pass region nevertheless grows quickly with τ_s > τ_eff.

The low and high input frequency limits become independent of τ_s when time is expressed in those units. Nevertheless, we can still write the critical condition independent of τ_s when expressing time in units of τ_I by combining Eqs (42) and (43) and eliminating τ_s 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: τ_w, τ_eff, τ_V, and 1/ν₀ (the latter value chosen by setting σ_I 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 τ_I, we have (45) For g = 0, we have α_w = 1 and τ_eff = τ_V 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/τ_c be no larger than either ω_L or 1/τ_w. Otherwise, the only effect of the spiking is to flatten the high frequency response beyond 1/τ_c. 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, ν_ωL: = |ν₁(ω_L)|/ν_low. The response then depends on the five shape features, ν_low, ν_high, ω_L, Q_L, and ν_ωL. Denoting ξ = τ_w/τ_c, so that and , we can re-express the response function as (46) with dynamic gain (47) When ω = ω_L, , which implicitly defines ξ in terms of ν_ωL, Q_L and ν_high and closes the representation. Indeed, with time in units of and gain values relative to ν_low, the shape of the filter depends only on this triplet: each of the six regions in (ν_∞, ν_ωL)-space defined by the boundaries ν_∞ = 1, ν_ωL = 1, and ν_∞ = ν_ωL provides filters of a qualitatively similar class (see Fig 9).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. The 6 distinct filter shapes in (ν_ωL, ν_∞)-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_L. Colored lines (blue to red) represent the Q_L-dependent boundary below which filter shapes are forbidden because of unstable dynamics. We note that ν_{ωL, ν∞ → 0} = Q_L. An intrinsic frequency exists in region above the Q_L = 1/2 boundary. A voltage resonance exists in the region above the Q_L = 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 ν_ωL, ν_∞ = 10^−1.5, 10⁰, 10^1.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.

https://doi.org/10.1371/journal.pcbi.1004636.g009

In particular, depending on the region there is a peak, dip or step at ω_L whose width varies with Q_L. 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_L. 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 ν_ωL = Q_L(1 + ν_∞) (which for different Q_L are shown as the colored lines in Fig 9).

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. The accessible region of filter shapes depends on Q_L and the relative speed of spiking to intrinsic dynamics ξ = τ_w/τ_c.

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 ν_ωL with Q_L. For relatively slow intrinsic spiking (a, b, c), there are regions of non-spiking resonant(ν_∞ > ν_ωL), but voltage resonant filters. Filters for relatively fast intrinsic dynamics (d, e, f) only exist as high pass resonant filters for large Q_L. (Left to right: , 1.1. Top row: ξ = 10. Bottom row: ξ = 0.1).

https://doi.org/10.1371/journal.pcbi.1004636.g010

Which root can also be checked by which of τ_w and τ_c 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 ν_ωL 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 ν_ωL = ν_∞ for large ν_∞. For increasing Q_L, the accessible region shifts right and up so that the band is eventually contained in ν_ωL > ν_∞ and ν_ωL > 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 ν_ωL. The lower bound is ν_ωL > Q_L(1 + ν_∞) when and when . The latter bound differs significantly from Q_L(1 + ν_∞) when Q_L > 1/2.

The voltage resonance condition can also be mapped to this space by replacing ω_Lτ_w 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 ν_ωL < ν_∞, 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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_L = 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.

https://doi.org/10.1371/journal.pcbi.1004636.g011

For ξ > 1, and Q_L decreasing from , the lower bound to the voltage resonance region interpolates across ν_∞ from the line , which rapidly approaches ν_ωL = 1 as Q_L is increased, to the lower bound of the region of stable filters, ν_ωL = Q_L(1 + ν_∞). Thus, stable filters exhibit a voltage resonance when ν_ωL 1, independent of Q_L. The absence of a spiking resonance, ν_ωL < ν_∞, 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_L < 1/2, the depolarization condition, γ⁻² < ξν_high < γ², 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 ω = ω_L 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.

Synaptic current.

I_syn contributes current terms of the form g_syn(t) (V − E), where E is the reversal potential for the synapse type and g_syn(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_syn(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_syn 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_syn(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_syn(t) on which rides a DC component described by the mean channel. We can thus write (49) where the zero-mean Gaussian process δI(t) is characterized by the variance, , and correlation time, τ_I, 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 σ_I for a desired output firing rate and measuring time relative to τ_I.

Subthreshold current.

In the most simple case (no longer exactly the Hodgkin-Huxley formalism), each somatic current, I_mem,i, contributes additively to I_mem with a term of the form (50) where g_i(V) is a voltage-dependent conductance, whose effect on the voltage dynamics depends on the driving force, V − E_i, the difference of the voltage and the reversal potential, E_i. g_i 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 (51) where v = V − V* and are the linearized variables for the voltage and subthreshold gating variable, x_i, respectively; and are the effective membrane conductances for the leak and for x_i, respectively; and τ_i = τ_i(V*) is the time constant of the dynamics of w_i. C_M is the capacitance of the membrane. The w variables have dimensions of voltage. Activation and inactivation gating variables have g_i > 0 and g_i < 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 τ_w the authors multiply the voltage equation by τ_w/C_M and analyze the behavior as a function of two dimensionless model parameters, α = g_Mτ_w/C_M and β = g_wτ_w/C_M, 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 τ_w. Thus, we retain both of the timescales of the neuron model, τ_V and τ_w. We then parametrize our model using the relative conductance g = β/α = g_w/g_M, the relative membrane time constant τ_V/τ_I = α⁻¹ = C_M/τ_I g_M, and the relative w timescale, τ_w/τ_I = β/τ_I 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], (52) We have absorbed the 1/g_M factor into the units of I_syn so that all dynamic quantities are in dimensions of voltage. We keep τ_V > 0 by setting g_M > 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_mem 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_i, 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], Δ_T, 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 Δ_T 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_T, 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_T, becomes good for Δ_T → 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_limit depending on Δ_T, the artefact can be safely ignored by considering the shape of the response only for f < f_limit. Conveniently, an upper bound on realistic values of Δ_T given by the surprisingly quick rise of real action potentials leads to a value of f_limit 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_T)/Δ_T], 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_limit, 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/(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], (53) The high frequency limit of the Gauss-Rice GIF is Eq (43). Equating these two expressions, we obtain (54) where g, τ_c, and τ_s 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⁻¹ 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 Δ_T = 3.5 mV(0.35 in our units), and a value an order of magnitude smaller, Δ_T = 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_limit.

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, τ_I, relative to the voltage correlation time defined here as the differential correlation time, . τ_s 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 τ_I/τ_s ≪ 1. In the other limit, τ_I/τ_s ≫ 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 τ_I, finding correspondence only in a fairly tight band around the membrane time constant, τ_V, from τ_I = 0.5τ_V to τ_I = 2τ_V [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, τ_I/τ_s ∼ 1, which reduces to a less strong version of the one Badel used for the LIF where . From the derivation of τ_s for the Gauss-Rice GIF exposed in the main text, the condition τ_I/τ_s ∼ 1 implies that the membrane time constant is no longer required to lie within an order of magnitude of τ_I 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 ν₀τ_r ≪ 1, where ν₀ is the firing rate and τ_r 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² > detB and when r² < 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 τ_r ≲ τ_s holds up to a saturation level, and this together with ν₀τ_s ≪ 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.