Rothman and Manis (2003) developed a Hodgkin-Huxley-like neuron model that has been widely used for modeling the phasic firing of cells in the auditory brainstem. It consists of a sodium current I_Na, a high-threshold (I_KHT) and a low-threshold (I_KLT) potassium currents, a hyperpolarization-activated cation current I_h and a leak current I_lk. The current balance equation has the following expression: (1)CdVdt=−INa−IKHT−IKLT−Ih−Ilk+I(t), where C is the membrane capacitance, V is the membrane voltage and I(t) is the external input current. Each ionic current I_i is governed by an activation and/or inactivation variable, and Equation (1) can be written as

(2)CdVdt=2[−g¯Nam3h(V−ENa)−g¯KHT(0.85n2+0.15p)(V−EKHT)            −g¯KLTw4z(V−EK)−g¯hr(V−Eh)−g¯lk(V−Elk)]+I(t), where ḡ_i and E_i are respectively the maximal conductances and reversal potentials for the ionic current i = Na, KHT, KLT, h, lk. Maximal conductances and channel gating rates are multiplied by a factor of 2 and 3, respectively, as in Gai et al. (2009), to mimic the brain slices during whole cell recordings at temperature 32°C.

Of particular interest for our problem is the dynamics of activation of low-threshold potassium current I_KLT governed by the variable w and the dynamics of inactivation of the transient sodium current I_Na governed by the variable h.

To highlight the primary biophysical mechanisms for phasic firing properties and to facilitate our analysis, we proceed as in Meng et al. (2012), and develop 3 reduced versions of the 8-variable RM03 model by identifying and approximating some nonessential features for the model's excitability and spike generation mechanism. Namely, we set m = m_∞(V), we freeze the inactivation gating variable z of I_KLT and the activation gating variable r of I_h, and we remove I_KLT (see Meng et al., 2012 for more details and a justification of this simplification).

The reduced models considered are two 2-variable models, both divisive and subtractive dominant versions, and a third 3-dimensional model that combines both mechanisms.

We compute the ratio between voltage and gating variable time constants for the divisive and the subtractive models.

We define the total conductance as gtot(V,h,w)=2(gNa(V,h)+gKLT(V,w)+gl). The ratio for the divisive model is given by: r(V,h)=3Cgtot(V,h,w0)τh(V), and for the subtractive: r(V,w)=3Cgtot(V,h0,w)τw(V). Notice that values of r close to 1 indicate that the two variables evolve on a similar time scale, while values of r close to 0 indicate that the voltage evolves on a much faster time-scale than the gating variable.

We consider several types of deterministic inputs: steps, ramps and also a half-wave rectified sinusoidal input I(t) with amplitude A and frequency ω, (10)I(t)=A[sin(2πωt)]+, where [·]⁺ = max(·, 0). In some cases we inject external noisy input of the form I_noise = Cση(t), where η(t) is a “white noise” process.

We also consider conductance based synaptic-like currents. Each synaptic site generates a minimal excitatory or inhibitory postsynaptic conductance (EPSG/IPSG) that is modeled as an alpha function with time constant τ_syn = 0.3 ms, unless otherwise stated, and maximal conductance g_max: (11)gsyn(t)=gmaxtτsyne1−t/τsyn. Thus, conductance-based synaptic currents have the form: Isyn(t)=gsyn(t−ts)(V−Esyn), for t ≥ t_s, where E_syn = 0 and E_syn = −75 mV for excitatory and inhibitory synaptic inputs, respectively, and t_s is the pre-synaptic spike time.

For the case of steady (but random) synaptic inputs, the synaptic event times are Poisson distributed with a given frequency and their maximal conductances are modeled as a summation of N = 7 independent miniature EPSGs with fixed amplitude, each of which occurs with probability p = 1/2. Thus, maximal conductances for composite EPSGs are discrete and modeled with a binomial distribution, B(N, p), where N = 7 is the maximum number of mEPSGs each with probability p = 0.5 of occurence. Thus, the mean number of mEPGS is 3.5, i.e., there will be 50% probability for N ≤ 3 and 50% probability for N ≥ 4.

In order to test coincidence detection, we inject periodically modulated excitatory synaptic inputs from several fibers/sites as in Jercog et al. (2010), similar to what might occur in vivo (Joris et al., 1994). Each cycle's composite input was generated from eight small (mini) excitatory postsynaptic conductances (mEPSGs) modeled as alpha functions with fixed amplitude and event times per site drawn from a von Mises distribution independently at each cycle. The von Mises probability density function for the angle θ is given by: (12)f(θ,a,b)=ebcos(θ−a)2πI0(b), where a = 1/4 is the mean and b is the temporal coherence (b = 0 corresponds to a uniform distribution and as b increases the distribution becomes more concentrated about the angle a), and I₀(x) is the modified Bessel function of order 0. Notice that to vary b is equivalent to vary the vector strength of the periodic input train. The vector strength (VS), also known as “synchronization index” (Goldberg and Brown, 1969), is a measure of how clustered are events over a cycle. To compute the VS one associates to each event time a vector on the unit circle with a phase angle and computes the mean vector. The VS is given by the length of the mean vector. Perfect clustering is obtained when VS = 1. The relationship between VS and the temporal coherence b of the input distribution is given by: VS(b)=I1(b)/I0(b), where I₀(x) and I₁(x) are the modified Bessel function of order 0 and 1, respectively.

In our simulations we chose maximal conductances for individual mEPSGs according to the following criteria: for low input strength, we pick g_max,e so that six coincident inputs, but not less, generate a composite EPSG that exceeds threshold for spike at rest, and for high input strength we reduce this number to four coincident inputs. Notice, that since the EPSG threshold for spike is different for each model, we used different values of g_max,e for each model.

In order to test the role of timed inhibition we also inject periodic inhibitory inputs generated in the same way as the excitatory ones: composite IPSGs are obtained from mini-IPSGs (mIPSGs) modeled as alpha functions (see Equation 11) with maximal conductance g_max,i and the timings in the arrival of mIPSGs are also drawn from a von Mises distribution. The frequency of the inhibitory input will be the same as for the excitatory input, but we will vary other parameters for inhibition: the temporal coherence b_i in the von Mises distribution (Equation 12), the maximal conductance g_max,i and time scale τ_i of the mIPSGs (Equation 11). Since we want explore the role that arrival time of inhibition plays in the maximal response of the system, we also vary the phase difference φ between the inhibitory and the excitatory inputs (φ = a_i − a_e (mod 1), φ ∈ [0, 1]), where a_i and a_e are the means of the von Mises distribution for the inhibitory and the excitatory input trains, respectively. Notice that when φ is close to 0, inhibition just follows excitation, while when φ is close to 1, inhibition just precedes excitation.

RM03-like models have low input resistance. The high conductance shunts EPSCs as well as spike currents. Strong EPSGs are required to elicit a spike and distinguishing spikes from EPSPs requires care. Although spikes may vary in amplitude, we observe that there is a sharp rise in voltage response when the input exceeds a certain value (see Figure 2C), showing the presence of a regenerative process. Thus, our criterion for spike detection requires that V passes a set level (−20 mV) and that the net intrinsic current when V = −20 mV is negative, i.e., dominated by sodium current (in order to detect the regenerative process).

Equations were integrated numerically in C using an implicit 4th order Runge-Kutta method. Stochastic differential equations had a white noise term and were integrated using an Euler-Maruyama method for stochastic differential equations (Higham, 2001) with time step Δt = 0.005 ms, along with a random number generator from the GSL-GNU Scientific Library. The bifurcation diagrams were computed using the Auto feature in XPPAUT. We used Matlab and Python to analyze and plot the data.

We have developed reduced models (with 2 or 3 variables) that show type III excitability (each model fires once only at the onset of a step current and not repetitively thereafter) using two different biological mechanisms, subtractive and divisive, that we managed to isolate.

Our starting point model is a reduced version of the Rothman and Manis model (RM03) (Rothman and Manis, 2003) for phasic neurons in the auditory brainstem (see Equation 1 in Methods), that we previously showed behaves semi-quantitatively and qualitatively as the full model (Meng et al., 2012). This reduced version contains a potassium current that is partially activated at V_rest and can significantly activate for subthreshold voltages, and a transient sodium current that can significantly inactivate for subthreshold voltages. Noting that “threshold” does not correspond to a fixed voltage value for spike initiation, especially in the case of phasic models (Azouz and Gray, 2000; Platkiewicz and Brette, 2010, 2011), we use as threshold the voltage value at which sodium activation (responsible of spike generation) rises sharply (i.e., m∞3(V)), which is approximately −40 mV. This voltage value would correspond to spike initiation if the dynamic subthreshold negative feedback factors were not present (see Section Methods for our precise criterion for spike identification). For this reduced RM03 model we also remove the I_KLT current and freeze the activation variable of I_h, as well as the I_KLT inactivation variable at their resting values (V_rest = −63.6 mV).

Now, we formulate 2-variable models that isolate the two negative feedback mechanisms. The Substractive (S) model (see Equation 3) isolates I_KLT using a non-inactivating sodium current by freezing the inactivation gating variable of the RM03 model at its resting state value. Without dynamic Na⁺ inactivation, the S model can develop bistability (coexistence of V_rest with a stable depolarized state). To prevent such bistability and preserve type III excitability we reduce ḡ_Na from 1000 to 177 nS (see Appendix B in Supplementary Material for more details). Although this is a strong reduction (by a factor of 1/5 or so), notice that Na⁺ strongly inactivates at subthreshold values of voltage and spikes ocurr with h approximately 0.1, thus the total Na⁺ conductance at spike onset is more comparable. Nevertheless, even with this reduction, the density of sodium channels for the S model is sufficient to amplify subthreshold synaptic potentials, thus retaining the excitability property. The divisive (D) model (see Equation 6) isolates Na⁺ inactivation by freezing the potassium activation variable at its resting value (w₀ = w_∞(V_rest)). However, freezing w to its resting value converts this reduced D model (as well as the original RM03) from phasic to tonic (Day et al., 2008; Gai et al., 2009, 2010; Meng et al., 2012). To counteract this effect and retain phasic firing in the D model we shift h_∞ by 6mV leftwards with respect to the original h_∞ in RM03 (see Figure 1). Moreover, in order to better compare the results with the S model we decrease ḡ_Na from 1000nS to 500nS. Finally, we combine these two mechanisms in a 3-variable model: the combined (C) model (see Equation 9) has I_KLT of the S model and I_Na of the D model.

Each of the models, C, D, S, fires phasically (only one spike, if any, at the onset of the input; see Figure 2A, inset) and its steady state V_ss depends on I, the applied current, but is stable for any I (Figure 2A). There is no bifurcaton to repetitive firing on the V_ss vs. I curve - a signature of type III excitability. For a slow modulation of the input current I there is no spiking; the voltage tracks the V_ss vs. I relation. Depolarization in the S and C models with I is modest because of I_KLT activation. In contrast to this rectification in the V_ss vs. I relation for S and C, the D model has a nearly linear relation without V-gating of I_KLT and negligible steady state I_Na. Thus, the set of voltages at steady state spans a wider range for the D model than for the S and C models.

In a phasic firing system, spikes are elicited only for transient inputs, fast enough and strong enough inputs. For step current input our models respond sensitively to a step's amplitude (Figures 2A(inset),C), effectively showing a threshold-like behavior for this stimulus class. The C model has higher current threshold than the D model – they have the same ḡ_Na but C has two negative feedback processes. S has the highest current threshold, reflecting the decreased ḡ_Na. Notice that the amplitude of spikes is modulated accordingly.

Phase plane analysis and relative time scales between voltage and gating variable help us understand dynamic features of spike generation (see Figure 2B and Methods Section). The trajectories corresponding to the time courses in Figure 2A (inset) reveal the characteristics of excitability: amplification of V for an adequate stimulus and then recruitment of dynamic negative feedback (Figure 2B). For these phasic systems, the steady state is always on the “left branch” for each I-value, guaranteeing its stability (Rinzel and Huguet, 2013). For steady strong input the V-nullcline's cubic or N-shape is greatly diminished corresponding to suppressed excitability (strongly activated I_KLT or strongly inactivated I_Na). Notice that this change is more dramatic in the D model than in the S model (the V-nullcline for I = 0.6 nA in the D model is hardly visible).

When the external input I increases, the steady state drifts downward along the left branch and the phase point approaches the steady state. In the S model, activation of I_KLT precludes excessive depolarization, limiting the steady state voltage to values no higher than −60 to −50 mV. Close to the resting state, the voltage and the gating variable have “comparable” time-scales (within a factor of 3), and the phase point does not escape rightward directly (red trajectories in Figure 2B). Only for large enough input I, the phase point can enter the region with a clear time-scale separation (blue region for values of V above −40 mV), and escape rightward toward the right branch of the V-nullcline, generating a spike (green trajectories in Figure 2B). Repolarization of the membrane potential after a spike is only due to the Ohmic leak current in the D model, unlike S and C models that have the V-gated potassium current, and there is no hyperpolarization when the spike terminates.

For a time-varying input the V-nullcline and target steady state move with the stimulus. Thus, for a slow stimulus ramp, the phase point tracks slowly the drifting steady state and there is no spike. For a fast enough ramp the phase point can escape rightward and lead to a single spike and then settle to a stable depolarized steady state. For this reason, phasic neurons are often called differentiators because they respond to fast-rising, but not to slow-rising, inputs (Ferragamo and Oertel, 2002; Svirskis et al., 2002; McGinley and Oertel, 2006; Izhikevich, 2007; Ratte et al., 2014). Our D model is a less sensitive differentiator; its critical value for the ramp slope to elicit an action potential is lower (see Figure S1A) and, in response to a half-wave rectified sinusoidal input (Equation 10), it can fire at lower frequencies than C and S (see Figure S1B).

Type III excitability occurs in each of our representative models but more or less robustly for mixtures of I_Na and I_KLT. We identified regions in the parameter space (ḡ_Na, ḡ_KLT) where the models behave phasically, i.e. the fixed point is stable for any value of I (Figure 3A, gray region). For other ḡ_Na−ḡ_KLT combinations the models show destabilization of the fixed point through a Hopf bifurcation and repetitive firing (Figure 3A, red region) and the possiblity of multiple steady states, including a stable steady state of high voltage, “lockup,” (Figure 3A, white region).

The divisive mechanism is key to guarantee the robustness of the phasic properties of the system with respect to changes in channel density. Removal of Na⁺ inactivation, as in the S model, strongly constrains the parameter range for type III excitability. To avoid “lockup” in the S model (large white region in Figure 3A) we chose its value of ḡ_Na smaller than for D and C models; consequently, a stronger input was required for the S model to be excited (see Figure 2C). By comparing the parameter spaces of the D and C models, we observe that when ḡ_Na in the D model is increased, the system transitions to a regime with oscillatory behavior, while the presence of I_KLT narrows this region (compare Figure 3A for D and C).

The threshold for eliciting a spike with a brief synaptic excitatory conductance input (EPSG) depends on ḡ_Na and ḡ_KLT (see Figure 3B). If ḡ_Na is small, then the system is not excitable (white region in Figure 3B). Of course, as ḡ_Na increases the system becomes more excitable and the EPSG threshold decreases. Moreover, if ḡ_KLT is large, more inward current (I_Na plus the excitatory current) is needed to overcome the opposing potassium current and EPSG threshold increases. Thus, in Figure 3B, along contours, ḡ_Na increases as ḡ_KLT increases.

Neurons and models with type III excitability respond to input transients but not to slowly varying inputs. In contrast, type I and type II excitable systems fire repetitively for steady input current. Their input-output relations, firing frequency vs. current (f-I curve), typically show monotonic increase over much of the I-range and then at high I-values loss of oscillation with either gradually decreasing amplitude or with a sudden drop in frequency and amplitude (Rinzel and Ermentrout, 1998; Borisyuk and Rinzel, 2005; Rinzel and Huguet, 2013). For type III excitability, an f-I relation quantifies the response to a fluctuating input: the mean firing frequency vs. the mean of a noisy current. The f-I curve is non-monotonic, as seen in some illustrative cases (Higgs et al., 2006; Lundstrom et al., 2008, 2009; Gai et al., 2009), and particularly distinctive from typical (noise free) type I and II cases. The firing frequency (firing probability per unit time) decreases smoothly toward zero for high mean I and the f-I relation shows strong sensitivity to input variance.

Our three phasic models have non-monotonic f-I relations and show strong sensitivity to noise (Figure 4A). The f-I curves are similarly shaped for the S and C models but less symmetric than D's f-I curves. There is no repetitive firing without noise; the increase of firing probability with noise level evidences the models' sensitivity to input transients (different colors in Figure 4A). The S model fires less frequently than D or C to a given noise level, partly because ḡ_Na is smaller for the S model requiring stronger fluctuations to generate spikes (see Figure 3B). The C model shows intermediate sensitivity to noise.

We highlight the different effects of subthreshold V-dependent conductance gating by replotting the input-output relations of Figure 4A in terms of mean voltage, 〈V〉 (Figure 4B). For S and C the maximum firing rate occurs for 〈V〉 around V_rest, slightly depolarized for D. The firing rate of S and C decreases for increased 〈V〉 and falls abruptly to zero for 〈V〉 just below −50 mV, less than the activation voltage for I_Na. This abrupt falloff reflects the shunting effect of strongly activated g_KLT and reduced amplitude of V fluctuations. On the other hand, firing rate rises gradually with 〈V〉 below V_rest with substantial firing probability 〈V〉 well below V_rest. At these hyperpolarized levels the input resistance is high (g_KLT is deactivated) and V-fluctuations are substantial. In contrast to the asymmetry in the f-V relations for S and C the relation for D is more symmetric, like that for f-I relation of D. Ohmic leak is the only conductance at the foot and tail of the V-range for D. Therefore, the input resistance and voltage fluctuations are comparable.

The noisy current input of Section 3.4 may be viewed as an idealization of random synaptic input delivered to a cell. Here, we consider the response to trains of excitatory synaptic conductance inputs (EPSGs). The event times are Poisson distributed and the EPSG amplitudes are binomially-distributed (see Section Methods) with mean amplitude that will just elicit a spike from the resting state. The three models show non-monotonic input-output curves, firing rate vs. EPSG input rate (Figure 5A, black curves), qualitatively as seen in Figure 4A for noisy current injection but some features differ. On the low input side the models fire for current injection (Figure 4A) even for non-positive 〈I〉, since by chance some positive current fluctuations are strong enough to cause a spike. But for synaptic input, firing only occurs for positive EPSG rate (Figure 5A). This feature is also seen in Figure 5B (black curves) where the firing rate for stochastic EPSG input drops abruptly to zero as 〈V〉 decreases below V_rest; this corresponds to zero EPSG rate, i.e. no synaptic input. Substantial mean input, either noisy current injection or stochastic ESPGs, leads to depolarization and reduced firing probability. But for synaptic input, high EPSG rate leads to increased membrane conductance, thus EPSGs that might be superthreshold for lower input rates are not for higher rates; the EPSG threshold has effectively increased thereby attenuating the effect of input variance and smoothing the drop in firing rate (see more details in Figure S2). For the D model, the conductance nature of synaptic input accounts for the loss of symmetry in its input-output relation.

Synaptic inhibition differently affects our three models while they are undergoing stochastic EPSG input. For the S and C models, not surprisingly (although, see below), inhibition decreases the firing rate for a given EPSG rate (Figure 5A, red and green curves). In contrast, for the D model the relationship between firing rate and inhibition can be non-monotonic (Figure 5A), depending upon the EPSG rate. Specifically, the firing rate decreases with g_inh for low EPSG rate, whereas for high EPSG rate firing rate increases. Consequently, at moderate input rates the firing rate behaves non-monotonically, increasing then decreasing with g_inh (see more details in Figure S3).

There are some notable features of the input-output relations, firing rate vs. 〈V〉 (Figure 5B). Overall, the relations appear qualitatively similar for the 3 models, non-monotonic with a sharp rise and gradual fall for increasing 〈V〉. The “left branch” extends to lower 〈V〉 with g_inh; this is because g_inh lowers the effective resting voltage, where firing rate descends to zero as 〈EPSG〉 rate decreases to zero. For the S and C models, in case of either stochastic EPSGs, here (Figure 5B), or for noisy current injection (Figure 4B), firing rate decreases from its maximum to near zero for the same 10 mV range: (−60, −50 mV). Notably, the maximum firing rate occurs at the same mean voltage independently of the inhibition level; this “sweet spot“ for firing is again located about the same voltages as for the noisy input (Figure 4B): ~−60mV for S and C and at a higher depolarized state (~−55mV) for D. From the perspective of a given 〈V〉 we see (except for 〈V〉 less than, say, −60 mV) that the firing rate decreases with g_inh. The effect is monotonic decreasing in each model. The reason being that in order to maintain a fixed 〈V〉 as g_inh increases, the excitatory input, EPSG rate, must also increase. That is, to keep a balanced state of fixed 〈V〉 both synaptic excitation and inhibition increase or decrease together. The amount of increased g_inh that compensates for increased excitatory input to maintain a balanced 〈V〉 is smaller for the D model than for S and C due to D's larger input resistance (see Figure S3).

Neurons that mainly fire due to rapid depolarization to threshold caused by time coincident inputs (like the models considered herein) are often called differentiatiors (Izhikevich, 2007; Prescott et al., 2008a; Lundstrom et al., 2009; Ratte et al., 2014). Although some neurons may display features of integrators and differentiators depending on the temporal pattern of the inputs (Gutkin et al., 2003; Rudolph and Destexhe, 2003), here we have focused on pure neural differentiatiors which, because of their intrinsic properties, are incapable of responding to temporally integrated, slowly varying input. This property is strongly linked to the concept of type III excitability. That is, the system approximately tracks a slowly varing input as a pseudo steady state -there is no bifurcation to oscillatory behavior with the input treated as a parameter. Input variability is the essential cause of spiking (Lundstrom et al., 2009, 2008). This excitability type complements the classical type I and type II excitability classes, where a stable fixed point destabilizes through a SNIC or Hopf bifurcation, respectively (Rinzel and Ermentrout, 1998; Borisyuk and Rinzel, 2005; Izhikevich, 2007; Prescott et al., 2008a; Meng et al., 2012; Rinzel and Huguet, 2013).

Neurons with type II excitability (emergence of repetitive firing via a Hopf bifurcation) are typically classified as resonators since the emergence of repetitive firing is via a Hopf bifurcation and thereby associated with subthreshold damped oscillations (some type III models can also show damped oscillations and thus resonator properties Prescott et al., 2008a,b; Mikiel-Hunter et al., 2016). However, such a system may display differentiator-like features when operating near but outside its repetitive firing regime (Prescott et al., 2008a; Ratte et al., 2014). The sensitivity to timing of multiple inputs, coincidence detection, is naturally inherited from the time windows of opportunity provided by subthreshold resonance. Some differentiator properties that we studied here are found exclusively in type III excitable systems. For instance, the fact that, in the presence of noise, pure differentiators have highly non-monotonic input-output curves and its firing rate is sensitive to variance and insensitive to mean input in the noise free case (see Figure 4). Moreover, in the present study, we added an extra ingredient: when the non-monotonic input-output curves are plotted as a function of the mean voltage (see Figures 4B, 5B), the firing rate decay occurs for subthreshold values of voltage (below ~−40 mV; when m is half activated), irrespective of the negative feedback mechanism considered. This property does not occur for type II neurons, for which the firing rate increases with mean input and the decay in firing rate may occur but at a much higher activated state.

Previous studies on specialized coincidence detectors in the auditory brainstem (Manis and Marx, 1991; Rathouz and Trussell, 1998; Svirskis et al., 2002) have highlighted the relevance of I_KLT in conferring phasic properties. Our early finding that Na⁺ current in MSO neurons inactivated at low-voltages suggested that this feature contributes to enhancing coincidence detection (Svirskis et al., 2004). Here, we underscore that Na⁺ inactivation plays a strong role in guaranteering the functional properties of the neuron model in terms of both spike preservation and phasicness. Thus, if Na⁺ inactivation is disabled, as in the S model, sodium conductance must be tightly constrained in order to preserve the spike and avoid the “lockup” state (see Figure 3). We tried but were not able to counteract “lockup” and retain phasic behavior with an additional potassium current that activated at high voltage (see Appendix B in Supplementary Material). With the divisive mechanism alone, spiking is preserved if h_∞ is left-shifted to lower values of voltage – I_Na must be strongly inactivated at rest – to retain type III excitability and avoid repetitive firing in response to steady current. However, having both types of negative feedback mechanisms in the neuron model (C model) provides a larger parameter space for conductances over which the system preserves phasicness, without compromising the spike or the excitability of the system. Thus, phasic models equipped with two mechanisms of different nature are more robust to changes in channel density (see Figure 3A).

In our analysis we considered extreme cases for our models. Indeed, it is not necessary that a C model, being phasic and performing properly, remain phasic when the subtractive mechanism is disabled. For example, the original RM03 model is phasic but behaves tonically when I_KLT conductance is frozen (Meng et al., 2012). Actually, there is a range of admissible positionings of Na⁺ inactivation for a C model that guarantee similar phasic features as our C model, even when the D model (frozen w) is not phasic. However, as h_∞ is right-shifted, the phasic region for the C model is reduced and when the shift is too dramatric, it resembles our S model. In this sense, our S model is a limiting case of the C model. Thus, our parameter choice for models allow us to, on one hand, clearly label them as S or C and on the other hand, keep the parameter values as similar as possible for a fair comparison. See more details in Figure S6.

We have considered the idealized situation of a point neuron model, without regard for potential effects of differential spatial distributions of ion current mechanisms for phasic behavior. MSO neurons provide a case for such considerations. In MSO neurons somatically recorded spikes are typically quite weak in vitro (Scott et al., 2007) and distinguishing what might be spikes or synaptic transients is challenging from in vivo extracellular recordings. Spikes that are generated in the axon initial segment are presumably shunted by the considerable conductance of I_KLT and I_h in the soma-dendritic membrane (Khurana et al., 2011), even though Na⁺ channels (inactivated near rest) are present in the soma (Scott et al., 2010). A thorough modeling study of how MSO phasic firing may depend on both g_KLT and Na⁺ inactivation could include the spatial distribution of the various V-gated currents (Ko et al., 2016); the influence of spatial distribution should likewise be considered in other phasic firing neurons. MSO neurons are superb at computing time differences. In gerbil, they receive fast excitation and fast, soma-targeted, inhibition (Myoga et al., 2014). Perhaps the sodium current while inactivated at rest may be recruited (de-inactivated) by transient inhibition to respond selectively to well-timed inhibitory-then-excitatory event pairings as in post-inhibitory facilitation (Dodla et al., 2006).

Our findings can be extended to analyze the role of negative feedback processes in a wide range of biological rhythms. Indeed, the presence of two mechanisms of different nature opposing to depolarization is not exclusive of phasic systems. Thus, the original Hodgkin-Huxley model has also Na⁺ inactivation and K⁺ activation, which in this case, act mainly superthreshold. To counteract a strong non-inactivating inward current with a strong outward current constrains the allowable channel density. However, controlling the effect of autocatalysis by a divisive mechanism, makes the system less sensitive and more flexible to changes in channel density (Sengül et al., 2014). Spiking, crucial for neuronal communication and computation, needs to be preserved with respect to variations in channel density. A similar statement about “two better than one” can be applied to a network of excitatory and inhibitory neurons, say as in a firing rate framework like the Wilson-Cowan equations. In the classical formulation, inhibition acts as a subtractive mechanism to counter recurrent excitation (equivalent to sodium activation). Dynamic synaptic depression on the excitatory-excitatory interaction behaves divisively like sodium inactivation by directly controlling and reducing autocatalysis. Having both mechanisms enables a network model to behave phasically (Tabak et al., 2006, 2011).

In further regard to robustness and parameter choices in our study, we acknowledge some compromises. Since we needed to choose a low value of ḡ_Na for the S model, the excitability for this model was reduced. Thus, in order to keep a fair comparison with the other models, we also reduced ḡ_Na for D and C models with respect to the original RM03 model, although this was not necessary since D and C still are type III for a larger range of ḡ_Na (see Figure 3A). Lowering Na⁺ conductance reduces, of course, the excitability properties of the system, and leads to lower amplitude spikes. Lowering ḡ_Na was also suggested in Lundstrom et al. (2008) as a mechanism to turn an integrator into a differentiator.

Transient inhibition can reduce or, surprisingly, enhance firing probability for peri-threshold random synaptic input events, depending on the timing and precision of the inhibition (Dodla et al., 2006). This property of post-inhibitory facilitation (Dodla et al., 2006) is a transient analog for post-inhibitory rebound after prolonged hyperpolarization. It can be expected in any system with a dynamic excitability-suppressing factor that is partially activated at rest and reducible by transient hyperpolarization or in spontaneously firing conditions for randomly arriving excitatory and inhibitory events (Dodla and Rinzel, 2006). As a model, we considered the response sensitivity for periodically delivered, compound excitatory input events (say multi-synaptic inputs) that are near to threshold but by chance are sometimes subthreshold/superthreshold and more or less time dispersed. Our models (S and C) with I_KLT express the behavior. Inhibition reduces the excitability-suppressing I_KLT current and although simultaneously positioning membrane potential further from threshold it opens a time window of enhanced excitability if inhibitory conductance decays fast enough (see Figures 5, 10). However, when excitability suppresion is primarily due to Na⁺ inactivation as in the D model, hyperpolarization from a transient inhibitory input can remove some Na⁺ inactivation and thereby effectively lower the spike threshold in spite of hyperpolarization. A normally subthreshold excitatory input may then trigger a spike even if inhibition has not completely decayed away (see Figures 5, 10). The C model should behave more as the S model or the D model depending on the contribution of each mechanism. Thus, the PIF phenomenon for the D model is less sensitive to the arrival time of the inhibitory inputs, compared to S and C, and this is especially noticeable at high frequencies.

While timed inhibition can facilitate signal detection, it can also enhance coincidence detection and precision. Several studies have shown that timed inhibition close to excitation, either preceeding or following excitation, can improve temporal precision of neurons (not only neurons with type III) by narrowing the window for coincidence detection, sometimes at the price of reducing the firing rate (Brand et al., 2002; Grothe, 2003; Ingham and McAlpine, 2005; Kuenzel et al., 2011). Indeed, we observed that timed inhibition close to excitation can deselect some spikes that do not correspond to coincident inputs and improve precision, especially for the S and C models, but also modulate the rate encoding and dynamic range of temporal coherence for coincidence detection (interval of b-values for gradation of firing probability), creating windows over different b ranges (see Figure 11).

Many of our results are demonstrated for a specific conductance-based model (adapted from Rothman and Manis, 2003). However, the geometrical analysis and underlying mathematical structure of our reduced versions of the model suggest that qualitative aspects of phasicness will be found in a more general class of phasic models. Take for instance, Clay's model for a healthy squid giant axon (Clay et al., 2008), that was obtained from the original Hodgkin-Huxley (HH) model by steepening and left-shifting the activation of I_K. The modified model can be seen as our C model, since HH has sodium inactivation. Another possibility is to turn the standard HH model into a type III excitable model by lowering the conductance of I_Na from 120 to 83 mS/cm2 (Lundstrom et al., 2009). In this case, the model can be seen as our D model but with a potassium current that activates superthreshold or close to threshold.

Indeed, we may have kept the high threshold potassium current I_KLT (present in the original RM03 model) in the D model. The role of this current is to contribute to the repolarization of the membrane potential once the spike has occurred – recall that in the D model the repolarization is only due to the leak current – but it has no effect on the spike generation, since this current is not active at subthreshold values of voltage. For this reason, the results described for the D model will also follow if I_KLT is added, since they mainly involve the subthreshold dynamics. Some differences will be observed if the activation of I_KLT is shifted to lower values, closer to threshold (around −40 mV for instance). In this case, I_KLT will behave as an I_KLT current, so the model will become a C model, with proper scaling of conductances.

Similarly, the I_KLT current can also be included in the C and S models. For the C model, the inclusion of I_KLT does not change the behavior of the model significantly (Meng et al., 2012). The S model is different since the inclusion of I_KLT might compensate the persistent I_Na, and allow for a larger value of ḡ_Na. Unfortunately, when I_KLT is incorporated in the S model, the system then switches to type II excitability (see Appendix B in Supplementary Material for a detailed discussion on this topic).

In many phasic models one can find a hyperpolarization-activated cation current I_h (Rothman and Manis, 2003; Khurana et al., 2011). The current I_h is mostly activated below V_rest with an activation time constant of 200−300 ms and deactivation by membrane depolarization on the order of tens of milliseconds. In our reduced models we have frozen g_h to its resting value, which is small. Allowing dynamic I_h may have an influence only on those results that involve dynamics below V_rest. This will constraint the range of voltages for negative currents, because I_h prevents excessive hyperpolarization, in the same way as I_KLT prevents excessive depolarization. Moreover, when I_h is activated the input resistance is lower at values of voltage below V_rest, thus the responsiveness to noise is smaller. The effects of a dynamic I_h are more noticeable in the C and S model than in the D model. Indeed, g_KLT is frozen at the resting value for the D model, therefore, when voltage hyperpolarizes, dynamic I_KLT deactivates for S and C, while for the D model it remains partially activated causing a similar effect on input resistance as if I_h were present: input resistance is smaller and there is less probability to fire due to noise effects. Thus, the presence of I_h (with dynamic rather than frozen-at-rest conductance) will have little effect on the essential results regarding coincidence detection properties. Even in the presence of inhibition, when I_h may reduce the IPSP size, the results will remain valid but for a stronger inhibitory input.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.