Asymmetrical voltage response in resonant neurons shaped by nonlinearities

Under the right combination of ionic currents and input frequency, neurons can exhibit subthreshold resonance, i.e., the membrane voltage response peaks at a nonzero frequency of the oscillatory input current. Experimental and theoretical studies have related this phenomenon to interacting processes at the ionic level, e.g., time constants and conductances.^1–4 However, there are still open questions such as whether and, if yes, how resonance has any functionality for signal processing,^5,6 and how alterations in the biophysical properties of neurons modify the resonant properties of the voltage response.⁷ 

Resonant currents such as the hyperpolarization-activated current (⁠$Ih$⁠) contribute to the emergence of resonance. Recent studies have established the role that $Ih$ plays in determining the properties of neuronal resonance based on (i) dynamical systems analysis;^8,9 (ii) characterization of response curves, mostly measured in terms of the impedance (⁠$Z$⁠) magnitude and phase;^2,4,10 and (iii) autocorrelograms.¹¹ Another resonant current is the M-current (⁠$IM$⁠) which produces comparable effects to $Ih$⁠. In addition to the resonant currents, amplifying currents such as the potassium inward rectifier (⁠$IKir$⁠) and persistent sodium (⁠$INaP)$ do not create resonance but generate positive feedback effects, which in turn can augment the resonant effects when present.^2,12

Typically, by applying an oscillatory current with constant amplitude and linearly increasing frequency, such as the ZAP (Impedance Amplitude Profile) function,¹³ and measuring the voltage response, one can construct the input-output relationship over a range of different input frequencies. The upper and lower voltage envelope-curves can be used to demonstrate this relationship by simply connecting voltage peaks and voltage troughs as frequency increases. In the frequency domain, the impedance amplitude profile of the neuron identifies this property and is given by the ratio of the Fourier transforms of the output and the input.^2,3 The resonance frequency $fres$ is the impedance profile peak frequency (see illustration in Fig. 1). Nevertheless, the impedance profile, traditionally used to measure the response of neurons to oscillatory inputs, does not necessarily capture the properties of the neuronal response, which may require looking back into the envelope responses. When the ZAP input is applied with low amplitudes (tens of pA or less) and because the dynamics is approximately linear, in general, the voltage response oscillations are symmetric about a reference voltage line, i.e., the peak frequency of the upper voltage envelope coincides with the trough frequency of the lower voltage envelope, justifying the use of the impedance profile.

The steady state response of linear systems to oscillatory inputs is characterized by three properties: (i) the number of input and output cycles coincide and the output amplitude is uniform across cycles for a given input frequency, (ii) the output amplitude is proportional to the input amplitude rendering the impedance independent of the input amplitude, and (iii) the output is symmetric with respect to the stable equilibrium of the unforced system around which the forced system oscillates. In this work, we are implicitly assuming (i), while (ii) and (iii) are violated. In other words, in this work, we are considering nonlinear models for which (i) is satisfied, while (ii) and (iii) are not.

However, for voltage-dependent ionic currents (e.g., $Ih$⁠, $IM$⁠) it is expected that large voltage changes bias the impedance response differently between depolarized and hyperpolarized currents. An amplitude high enough to activate these currents accentuates nonlinearities leading to a neuron’s response with asymmetrical upper and lower voltage envelopes. In fact, recent work reported evidence of asymmetric voltage responses to ZAP functions, where the peak frequency of the upper voltage envelope differs from the trough frequency of the lower voltage envelope.^12,14–17

The question arises of how do the intrinsic nonlinear neuronal properties shape the voltage response (upper and lower envelopes) and the impedance profile. To address these issues, in this work, we investigate the asymmetric properties of the voltage response of neurons to oscillatory inputs and how these properties are shaped by the participating ionic currents. However, because the impedance profile does not necessarily capture asymmetries in the upper and lower envelopes, we propose a different version of the conventional impedance by looking at the upper and lower voltage envelope profiles normalized by the input amplitude, which we refer to as the upper (⁠$Z+$⁠) and lower impedances (⁠$Z−$⁠). For linear systems, the impedance is the difference between these two normalized quantities. We show how the upper and lower impedances may be used to characterize asymmetrical responses in the biophysical model. Our simulations are constrained within the biological range.¹⁸ We explain how asymmetrical responses are related to the activation curves of the modeled ionic currents, in particular, we show that the gating variables of these currents may lead to an unexpected frequency-dependent pattern. We explain our results by using dynamical systems tools (phase-plane diagrams). Finally, we use a simplified piecewise linear system subject to an oscillatory input to explain our findings in a more tractable way.

We use a conductance based neuron model with voltage-dependent currents and a leak current. The membrane potential (⁠$V$⁠) of the model obeys

where $Cm=1$$μ$F/cm$2$ is the membrane capacitance and $Ileak$ the leak current, given by $Ileak=gleak(V−Eleak),$ where $gleak$ is the leak conductance and $Eleak$ the reversal potential of the leak current. We chose the geometry of a cylinder with 70 $μ$m of diameter and 70 $μ$m of length which condenses soma and dendrite in a single compartment preserving the average capacitance of a pyramidal cell (⁠$C≈150$ pF).^4,19 In the case of other ionic currents, $Ii$ may represent the hyperpolarization-activated current (⁠$Ih$⁠), the M-current (⁠$IM$⁠), the potassium inward rectifier (⁠$IKir$⁠), or the persistent sodium (⁠$INaP)$⁠. They follow the Hodgkin-Huxley formalism²⁰ with dynamics given by

where $g¯i$ is the maximal conductance and $Ei$ the reversal potential. $Ai(V,t)$ is the activation variable defined as

where $τi$ is the time constant of the activation variable and $Ai∞$ is the asymptotic value of the variable, which follows the Boltzmann formalism

where the values $V1/2$ and $k$ control the sigmoid function: the first is the voltage value for which the ionic current is half activated [i.e., $Ai(V1/2)=0.5$] and the second determines the slope of the sigmoid function. The sign of $s$ controls whether it activates with hyperpolarization or depolarization. $IZAP$ represents the ZAP current, as considered in Ref. 13, given by

where $Fstart$ (⁠$Fstop$⁠) is the initial (final) frequency of the $IZAP$⁠, and $tstart$ (⁠$tstop$⁠) is the initial (final) time limit. The amplitude of the ZAP current (⁠$IZAP$⁠) is given by $Ain$⁠, which will be shown in every figure.

A constant current was applied to keep the neuron at different $Vhold$ values (⁠$IDC$⁠). Values of $Vhold$ will be indicated in the Results section. All parameters used are displayed in Table I.

The $INaP$ time constant (⁠$τNaP$ in milliseconds) is different for values of $V$ higher and below $−40$ mV, as considered in Ref. 21, which creates a fast/slow activation depending on the voltage, and is described as

The $IKir$ time constant (⁠$τKir$ in milliseconds) is described as

where $a=6.1$/s and $b=81.8$/s.^22,23 The $Ih$ and $IM$ time constants are displayed in Table I. Note that the time constants for $Ih$ and $IM$ are usually larger (slow dynamics) than the values assumed for $INaP$ and $IKir$⁠, which are small (fast dynamics). In addition, note that in order to simplify comparison between $Ih$ and $IM$⁠, we chose to model the $IM$ current with equal parameters as $Ih$ but instead of an observable activation it deactivates with hyperpolarization, i.e., $s=−1$ instead of $s=1$ [see Fig. 7(a) for examples of activation and deactivation].

Subthreshold resonance in neurons is usually approached by identifying a peak in the impedance magnitude, which is calculated as the ratio of the Fourier transforms of the output voltage and the input current $Z(f)=FFTout/FFTin$⁠. For nonlinear systems, where the output is periodic and has frequency following the input, the impedance magnitude expression can be written as

where $Vmax(f)$ and $Vmin(f)$ are the maximum and minimum voltages obtained from the output voltage for a given frequency $f$⁠.^3,7

Here, in order to characterize the asymmetric membrane potential responses, we propose two alternative measures to the usual ratio of Fourier transforms. We take the holding potential (⁠$Vhold$⁠) as a reference voltage line so that voltages above it will be positive and voltages below it will be negative, and, for each frequency $f$⁠, measure the magnitudes of the maximum (peak $V+$⁠) and minimum (trough $V−$⁠) voltage traces normalized by the ZAP current amplitude (⁠$Ain$⁠) (Fig. 2). We refer to these two quantities, which depend on the frequency $f$ and have dimensions of impedance as upper impedance [$Z+(f)$] and lower impedance [$Z−(f)$] [Eq. (9)]. These measures stress the presence of asymmetries better than impedance profiles because they do not average upper and lower responses. They are also better than observations of the upper and lower envelopes because they allow a better comparison by highlighting differences due to the normalization

where $V+/−$ is the absolute peak/trough distance from $Vhold$ and $N+/−$ are the number of peaks/troughs in the time series, respectively. For each peak/trough $i$⁠, we check the corresponding frequency on the ZAP current [from Eq. (5)] and use it to draw $Z+/−(f)$⁠. In practice, to obtain $Z+/−(f),$ we simply interpolate the set ${Z+/−(fi)}$⁠. Notice that the impedance $Z(f)$ is an average over $Z+(f)$ and $Z−(f)$⁠.

All simulations were run in the NEURON simulator using the Python interface.²⁴ Phase-plane analysis was done using MATLAB (The Mathworks, Natick, MA).

Here, we show how asymmetrical subthreshold resonance emerges in the model when a single voltage-dependent (resonant) current is present. We start with the $Ih$ current and explore how the interplay of its time constant and the membrane potential shapes asymmetries on the impedance profile.

In Fig. 3, we show examples of $Z+(f)$ (solid lines) and $Z−(f)$ (dashed lines) computed from our simulations. For the low ZAP current amplitude [$Ain=10$ pA; Fig. 3(a)], $Z+(f)$ and $Z−(f)$ are identical. On the other hand, for the high ZAP current amplitude [$Ain=1$ nA; Fig. 3(b)], $Z+(f)$ and $Z−(f)$ display different resonance peaks. Moreover, depending on the holding potential, a resonance peak may exist in $Z−(f)$⁠, but not in $Z+(f)$⁠.

In Fig. 3(c), we present the impedance profiles [$Z(f)$] computed for the same simulations in Fig. 3(b) (colors follow the same scheme). For $Vhold=−60$ mV, $Z(f)$ does not display resonance even though in Fig. 3(b) $Z−$ has a resonance peak. On the other hand, for $Vhold=−90$ mV, both $Z(f)$ and $Z−$ have resonance peaks. Since $Z(f)$ is a combination of $Z+$ and $Z−$⁠, it may fail to capture information about hyperpolarized and depolarized voltages.

The previous examples (Fig. 3) suggest that $Z+(f)$ and $Z−(f)$ impedance curves depend on the biophysical properties of $Ih$⁠, more specifically $τh$ and $Vhold$⁠. To study the combined effect of these parameters on the impedance profile (low or bandpass filter), we characterized all four possible scenarios: (i) both $Z+$ and $Z−$ do not exhibit resonance, i.e, are low-pass filters; (ii) $Z+$ is a low-pass filter, and $Z−$ is a bandpass filter; (iii) both $Z+$ and $Z−$ are bandpass filters; (iv) $Z+$ is a bandpass filter and $Z−$ is a low-pass filter, which was not detected in our simulations.

A high ZAP current amplitude abolishes resonance in $Z+$ for depolarized $Vhold$ [Fig. 4(b)]. Thus low-pass or bandpass filtering behavior for depolarizing current is dependent on the current amplitude. This suggests that in a physiological context, the arrival of oscillatory synaptic inputs might be able to modulate the postsynaptic neuronal response such that high amplitude inputs are preferentially transmitted at low frequencies (a selected frequency band) if the postsynaptic neuron is depolarized (hyperpolarized).

Interestingly, the existence of resonance in $Z−$ is not affected by the amplitude of the input current [Figs. 4(a) and 4(b)]. This might be due to the fact that $Ih$ is activated by hyperpolarization. Since $Ih$ is not activated by depolarization, the only remaining effects are due to the leak current making the response of the system linear. In addition, it is expected that for small $τh$ the system also behaves in a linear way,³ which surely hampers the presence of nonlinearities responsible for the asymmetries.

For a further characterization of the resonance properties under high ZAP current amplitude, we quantified the differences between the depolarizing and the hyperpolarizing membrane resonance peaks and frequencies (Fig. 5). We calculated the difference between the resonance peaks as $ΔZ=Z+(fres+)−Z−(fres−)$ and the shift between the resonance frequencies as $Δf=fres+−fres−$⁠. We show these shifts for different combinations of $τh$ and $Vhold$ in Fig. 6. Gray area indicates that $Z+$ has no resonance.

Different combinations of $τh$ and $Vhold$ can modulate the neuron’s response in both amplitude and frequency (Fig. 6). For small amplitudes such as $Ain=0.1$ nA, the kinetics of $Ih$ controls the effect of $Vhold$ on $ΔZ$ [Fig. 6(a1)]. When $τh$ is high we see two effects: whereas hyperpolarized $Vhold$ amplifies $Z+$⁠, depolarized $Vhold$ amplifies $Z−$ [see the schemes for points (b1) and (b2) in Fig. 6 and selected voltage traces in Fig. S1]. On the other hand, lower $τh$ values weaken the $ΔZ$ sensitivity to $Vhold$⁠. As $Ain$ increases, the influence of $τh$ on $ΔZ$ decreases, which may be related to the lack of activation of $Ih$ at high amplitudes. In other words, this latter observation is related to higher amplitudes resulting in large displacements of the membrane potential consequently masking the $Vhold$ effect. These remarks suggest a high sensitivity of the impedance amplitude to intermediate current amplitudes (⁠$Ain=0.1$ nA) and a saturation for higher amplitude currents (⁠$Ain≥0.3$ nA).

In Fig. 6(a4), $Vhold$ has no apparent influence on $Δf$ for slow $τh$ but higher frequency shifts are displayed at fast kinetics, which increase with $Vhold$⁠. When we observe the panels with higher ZAP current amplitudes [Figs. 6(a5) and 6(a6)], $Vhold$ begin to influence $Δf$ even for low $τh$⁠. In general, an increase in the input amplitude corresponds to an increase in the absolute value of $Δf$⁠. For regions close to $Z+$ with low-pass filtering behavior (light gray area), $fres+$ is near zero and, consequently, $Δf$ is maximized. These effects can be visualized in the schematic representation for point (b3) in Fig. 6.

Neuronal resonance generated by $Ih$ is governed by the interplay between the $Ih$ kinetics and the voltage-dependent conductance. Our results here suggest that asymmetries in the depolarized and hyperpolarized resonances (⁠$ΔZ>0$⁠, $ΔZ<0$⁠, and $Δf<0$⁠) are also driven by the same two parameters, but mostly at high amplitudes of oscillatory inputs. More specifically, when the neuron is subjected to an oscillatory stimulation with amplitudes of the order of tens of pA, larger displacements of the membrane voltage in relation to $Vhold$ are generated, and, consequently, $Ih$ will be less activated for $V>Vhold$⁠. When $Ih$ is insufficiently activated, the effect of $Ileak$ dominates and the neuron behaves as a low-pass filter, which does not display resonance. This effect explains why higher $Vhold$ values abolish $Z+$ resonance, whereas $Z−$ resonance remains. Furthermore, we predict that currents with opposed monotonic behavior on the activation curves with respect to the voltage (activation by hyperpolarization) could reverse this phenomenon, i.e., currents with activation for depolarized voltage values would generate resonance in $Z+$ but not in $Z−$ (see, for example, resonance exclusively on upper voltage envelopes at Figs. 6A2 and 6A3 in Ref. 15, where $Z+$ would demonstrate such effect).

In this subsection, we explore the effect of other ionic currents with different combinations of activation curves and reversal potentials. We chose $IM$ which is again a resonant current but now activated by depolarization instead of by hyperpolarization as $Ih$ [Fig. 7(a)]. We also chose two amplifying currents, namely, $INaP$ and $IKir$⁠, which are activated by depolarization and activated by hyperpolarization, respectively [Fig. 7(a)]. Our hypothesis is that the shape of the activation curves of a given ionic current controls the asymmetrical behavior in $Z+$ and in $Z−$⁠. In Fig. 7(a), we show the activation curves of the ionic currents studied in this paper. A second hypothesis would be that the main responsible for the asymmetry is not the activation curve but its effect on the phase-plane of the system. If this is the case, one would obtain similar results by having nonlinearities stressed on the voltage-nullcline instead. We will check this second hypothesis in Subsection III D.

Figures 7(b) and 7(d) show $Z+$ and $Z−$ in the case of large ZAP current amplitudes $Ain=1$ nA for the model with $Ih$ and for the model with $IM$⁠, respectively. The $Z+$ (⁠$Z−$⁠) curves for one current are nearly mirrored by the $Z−$ (⁠$Z+$⁠) curves for the other current. Whereas for $Ih$ [Fig. 7(b)] resonance occurs in the bottom envelopes, for $IM$ [Fig. 7(d)] resonance occurs only in the upper envelopes. The addition of an amplifying current such as $IKir$ and $INaP$ in the models has little effect and slightly shifts $Z+$ and $Z−$⁠, as observed in the figures.

In Figs. 7(c) and 7(e), we show the amplification of the impedance due to $IKir$ and $INaP$ currents. While $IKir$ preferentially amplifies at voltages close to $Vhold=−120$ mV, $INaP$ amplifies at $Vhold=−50$ mV. We also note that the amplification is different in these cases, the latter mostly amplifies the voltage at zero frequency while the former amplifies both the voltage at a nonzero frequency (resonance amplification). This voltage-dependent amplification is a consequence of the activation curves [as depicted in Fig. 7(a)] and it has been demonstrated elsewhere.⁷ 

In this subsection, we characterize the observations above using phase-plane analysis. By using such an approach we can provide the reader a geometrical understanding of how the voltage dynamics are affected by the ionic currents, thus allowing to explain and predict these dynamics. We present our results in Fig. 8 (⁠$Ih$⁠) and Fig. 9 (⁠$IM$⁠). We investigate the dynamics of the neuron in both the low and high amplitude cases in a frequency-dependent manner. We plot the nullclines either in red (⁠$V$-nullcline for $dV/dt=0$⁠) or in green (⁠$Ai$-nullcline for $dAi/dt=0$⁠). For reference, we also plot (dashed-red) the location of the $V$-nullclines corresponding to the peak and trough values of the oscillatory input and in dashed gray the intersection of the nullclines. Note that in these cases, so that we could show frozen trajectories, we substitute the ZAP current by a sinusoidal input current of a single frequency, which we can control. In other words, this method is equivalent to using the ZAP current as in Eq. (5) but for $Fstart=Fstop$⁠, and, therefore, it delivers the same output response with the advantage of observing it without the time-dependent input frequency (see Fig. 13 in Ref. 12 and Fig. 10 in Ref. 7 for examples and details of this method).

We follow Ref. 16 and describe our observations as below: at low frequencies (first panel; $f=2.5$ Hz) the trajectories follow the $Ah$-nullcline. Although we do not explore the limit case when $f→0$⁠, this is a situation where the system stays in a point. As frequency increases, there is an observable rotation. The resonant frequency is observed when the space covered by the trajectory in the voltage domain is maximal (middle panel; $f=6.5$ Hz). For even higher frequencies, there is an observable shrinkage in the limit cycle trajectory up to the point that only a small line is observed at very high frequencies (third panel; $f=90$ Hz). For such high frequencies the system becomes quasi-one dimensional and, in the limit $f→∞$⁠, the trajectory coincides with the fixed-point and the system lies in a point.¹² 

In the second row, we see how the system behaves for high amplitudes. Given that the trajectory covers a larger voltage range the neuronal dynamics is more exposed to nonlinearities on the plane. Nonlinearities on the plane shape the neuron’s response, which differs from the ones described above (small input amplitudes). From low to high frequencies (panels from left to right; $f=2.5$ Hz, $f=5$ Hz, and $f=90$ Hz) there are noticeable differences in voltage peaks and troughs which are captured by the shape of the trajectory. The amplitude of the peaks only become smaller for increasing frequency. However, the amplitude of troughs increase up to a peak at the resonant frequency, and then decrease. This effect is related to the resonance being exclusively found on the bottom envelope in this setup [compare to blue curves in Fig. 7(b)].

The mechanism of resonance in $Z−$ and the absence of resonance in $Z+$ for $Ain=1$ nA can be geometrically explained as follows. Note that the trajectories for low frequencies do not develop further than the dashed-red lines and that the green solid line intersects the leftmost red-dashed line creating an enclosed region. Due to this region, for low frequencies, the troughs cannot develop with high amplitudes. As frequency increases a rotation is observed, the trajectory escapes from this region and troughs increase in amplitude until a maximal value (resonance). For higher frequencies, we observe shrinkage because the system tends to quasi-one dimensional.

The same mechanism described above is present in Fig. 9, but now since the ionic-nullcline is mirrored we see a resonance exclusively on the upper voltage envelope [see $Z+/−$ examples in Fig. 7(d)]. For low amplitudes (first row), the trajectories follow the $AM$-nullcline and, as the frequency increases, they rotate and shrink. For higher amplitudes (second row), the troughs decrease monotonically and the peaks resonate. The shape of the $AM$-nullcline and how it crosses with the $V$-nullcline defines this behavior.

Yet, observations of the second row in Fig. 9 shows a surprising unexpected dependency of the activation variables with frequency. The trajectory does not move around the fixed point for low frequencies, but only does so for higher frequencies. To the best of our knowledge, such an unexpected frequency-dependent pattern was not reported in other subthreshold resonance phase-plane analyses (see, e.g., Refs. 3,9,16,25, and 26). This is a point often overlooked in different studies due to the fact that they rely on small amplitude currents.

Our explanation of the above mentioned frequency-dependent pattern of $AM$ is threefold: it is first related to the proximity of the nullclines creating a region of slow velocity, second to the time scale separation between $V$ and $IM$⁠, and third related to the activation properties of $IM$⁠. As one can observe in the bottom row of Fig. 9, the $V$-nullcline and the $AM$-nullcline are vertically close to each other creating a region with slow velocity in the proximity of these nullclines. Since there is a time scale separation, for low frequencies, the trajectory follows preferentially the $V$-nullcline and, as frequency increases, it rotates and follows the $AM$-nullcline. However, in this situation, the $AM$-nullcline is not placed horizontally forcing the system to go through a transient rotation until it manages to reach the fixed point. Due to activation properties of $IM$ [see Fig. 7(a)], the fixed point is reached through increase of $AM$⁠. To investigate this behavior even further we explored in Fig. 10 the dynamics of the $AM$ activation variable while delivering the ZAP input to the neuron for a few parameters.

The frequency-dependent pattern of $AM$ can be controlled by adjusting the value of $V1/2$⁠. When we change $V1/2$ we move the $AM$-nullcline and consequently the fixed point. If $V1/2$ is close to the resting potential, then $AM$ is constant and it does not activate or deactivate with frequency; it is fixed at $AM=0.5$ such that $IM$ is half activated. For values of $V1/2$ above (below) the resting potential, we found that the ionic channel will activate (deactivate) toward the fixed point for increasing frequency.

In Sec. III E, we will further explore the asymmetries in the voltage response using a simplified piecewise-linear model, and will provide a geometrical explanation of some of these phenomena.

In this subsection, we investigate if the asymmetrical response can still arise due to nonlinearities associated to the voltage variable and not associated to the ionic currents as we observed in Secs. III A–III C. We ask if the activation curves are the main reason for the observed asymmetries or if we could obtain qualitatively similar behavior by other means. We chose conductance-based models from Ref. 7 with $Ih+INaP$ currents. More specifically, we select model 1 and model 2 of that article. These models differ with respect to the parameters which introduce nonlinearities on the voltage-nullcline: it is either quadratic or cubic (see phase-planes in Fig. 11).

In comparison with the model that we have been using [Eqs. (1–4)], these models follow the same equations as in the other sections, but we emphasize two differences between them: their parameters and the fact that the $INaP$ current has its activation variable approximated to its steady state as in Eq. (10). Model 1 (quadratic) has the following parameters: For the leak current, we use $gleak=0.5$ mS/cm$2$ and $Eleak=−65$ mV; for $INaP$⁠, we use $g¯NaP=0.5$ mS/cm$2$⁠, $ENaP=55$ mV, $V1/2=−38$ mV, $k=6.5$ mV, and $s=−1$⁠; for $Ih$⁠, we use $g¯h=1.5$ mS/cm$2$⁠, $Eh=−20$ mV, $τh=80$ ms, $V1/2=−79.2$ mV, $k=9.78$ mV, and $s=1$⁠; for the external current, we use $IDC=−2.5$ $μA$/cm$2$ and ZAP current with $Ain=0.1$ $μA$/cm$2$⁠. Model 2 (cubic) has the following parameters: For the leak current, we use $gleak=0.3$ mS/cm$2$ and $Eleak=−75$ mV; for $INaP$⁠, we use $g¯NaP=0.08$ mS/cm$2$⁠, $ENaP=42$ mV, $V1/2=−54.8$ mV, $k=4.4$ mV, and $s=−1$⁠; for $Ih$⁠, we use $g¯h=1.5$ mS/cm$2$⁠, $Eh=−26$ mV, $τh=80$ ms, $V1/2=−74.2$ mV, $k=7.2$ mV, and $s=1$⁠; for the external current, we use $IDC=0.01$ $μA$/cm$2$ and ZAP current with $Ain=0.1$ $μA$/cm$2$ (For a more detailed description of the model we refer the reader to Ref. 7.),

In Fig. 11, we present a comparison of the quadratic and cubic models. In both models, the asymmetric response emerges. In the first column, we see that $Z+>Z−$ and that the quadratic and cubic models show qualitatively the same behavior. Given that the $Ah$-nullcline adds a quasilinear influence on the trajectories, the quadratic/cubic $V$-nullcline is the only nonlinear influence on the response of the model. For $f→0$ and $f→∞$⁠, the trajectory follows the same behavior described above: at slow frequencies, it follows closely the $Ah$-nullcline and at fast frequencies, it converges to the fixed-point. Resonance occurs at intermediate frequencies when the trajectory covers a larger voltage range.

As observed in the models from Secs. III A–III C, the amplitude has a strong effect on the asymmetrical response. In the case of the quadratic model, however, if the ZAP current amplitude is enhanced action potentials develop. This shows a limitation to the increase of $ΔZ$⁠. Nonetheless, in the case of the cubic model a further increase of the ZAP current amplitude is possible because the trajectory will develop around the cubic $V$-nullcline. Note this effect in Fig. 12 in the first column.

A second effect from the cubic $V$-nullcline concerns the position of the fixed-point which can be controlled by the $Ah$-nullcline. We perform this study which is presented in the second column of Fig. 12. Clearly, there is an inverted behavior: now $Z−>Z+$⁠.

These experiments confirm that asymmetrical response may emerge either from ionic currents or from the voltage nullcline. But, most importantly, it is a phenomenon which emerges from nonlinear systems.

So far, we have shown two different, but related phenomena: (i) asymmetric responses that are a property of nonlinear systems, and that may or may not be pronounced. The asymmetries become pronounced mostly as a function of the amplitude of the oscillatory input, the time scale separation or holding voltage; (ii) the frequency-dependent pattern of the gating variable, which is related to a combination of the voltage and ionic current properties.

In this subsection, we will use a simplified model to geometrically test and explain (i), but we are not able to explore (ii) because this model does not contain any activation dynamics. First, we identify what are the most important parameters that control the asymmetric responses and how the latter depends on the input amplitude. For simplicity, we chose a piecewise-linear model (PWL) to explain the asymmetry through a geometrical approach on the phase-plane.¹⁶ Nonlinearities are added by a sudden change in the slope of one of the nullclines. In this respect, we believe that the PWL system is simple enough to be addressed as a reduced system. The PWL system reads

where the functions $hv(v)$ and $hw(v)$ are described by

where $α=0.4$⁠, $η=−1$⁠, and $ηr=−0.2$⁠. The values of $ε$ (time scale separation) and $β$ (capturing the bias DC current) we use in our simulations will be indicated as needed.

As we have shown in Sec. III C, $Vhold$ plays a relevant role in determining the asymmetric response. Thus, here we search for parameters in the PWL system that represent changes in $Vhold$⁠. As one can observe from the schematic $Ih$ activation curve in Fig. 13(a), which mimics the $Ah$ curve, as $Vhold$ increases the $V$-nullcline gets closer to a region on the plane where nonlinearities are more prominent. We make a distinction between “linear region” and “nonlinear region.” The areas where asymmetries are not present are denoted as “linear regions,” and are found in the middle of the activation curve. The areas where asymmetries are present are denoted as “nonlinear regions,” and are found where nonlinearities are more prominent on the plane. These regions can be mimicked in the PWL system by the breaking point that separates the two linear pieces. In this regard, in the phase planes of the PWL system in Figs. 13(b) and 13(e), we have observed that as the nullclines are shifted by changing the values of $β$⁠, the breaking point (⁠$v=0.2$⁠) becomes more distant or closer, the latter one emphasizes asymmetries because trajectories get closer to the “nonlinear region” [Figs. 13(e) and 13(f)]. In fact, we show that this is a general phenomenon not only observed in a biophysical neuron model but also found in systems where nullclines are placed in such a way that they allow a certain competition between linear and nonlinear effects, i.e., between a “linear region” and a “nonlinear region.”

Nonetheless, if asymmetries are created by the fact that oscillatory trajectories lie around the nonlinearities in the $Ih$ activation curve (or equivalently, close to the breaking point in the PWL system), other mechanisms could achieve the same effect. Here, we were able to retrieve asymmetries by increasing the ZAP current amplitude which forces the trajectories of PWL system to the nonlinear region again [see dotted trajectory in the schematic representation in Fig. 13(a)]. This simulation uses the same system of Figs. 13(b) and 13(c) for higher amplitude and is presented in Fig. 14(a) where we retrieved the asymmetries.

Further we analyzed the effect of varying the time-scale separation $ε$ in the PWL system (see Sec. III A where we discuss time-scale separation controlled by $τ$⁠). In addition, as previously reported,¹⁶ an increase in the time-scale separation plays an important role in amplifying the voltage response. We tested the time-scale separation with $ε$ because it corresponds to the inverse of the time-scale of the system. Our simulations confirm that resonance emerges for $0<ε<1$ and is more pronounced the smaller the value of $ε$ within some range (similar results can be found in Ref. 3). At the same point, the asymmetries observed show maximum of $ΔZ$ [see Fig. 14(b)].

The results using the PWL system show an asymmetrical response emerging close to the breaking point, but we remind the reader that the $Ih$ activation curve posses two nonlinear regions: one for hyperpolarized $Vhold$ and the other for depolarized $Vhold$⁠. At the same time Fig. 6 shows that $ΔZ$ exhibits an inversion of signals at $Vhold≈82$ mV, which could be attributed to the trajectory of the system being closer to another nonlinear region of the curve. To understand this behavior we used the PWL system where the asymmetrical response was evaluated under $ηr$ variation in order to simulate both extremes of the $Ih$ activation curve. Note that for $ηr<−1.0$ (⁠$ηr>−1.0$⁠) the PWL system is comparable with the already mentioned hyperpolarized (depolarized) nonlinear portion of the $Ih$ activation curve. This is demonstrated in Fig. 15 where we show how different angles determined by $ηr$ can influence $ΔZ$⁠.

The PWL system could be used to describe other asymmetrical properties. An example is found in Ref. 15, where resonance was found only in the upper envelope and not in the bottom envelope. As we have shown in Sec. III C, this is related to activation curves such as $IM$ which activates with respect to the voltage instead of the deactivation observed for example, in $Ih$ [see Fig. 7(a)]. By working with the nullclines extracted from Eqs. (11) and (12) and producing phase-planes with qualitatively mirrored images we were able to recover situations where resonance is only found in the upper envelope confirming that the simple PWL system is a good geometrical approximation to generate asymmetrical responses.¹⁶ 

Upon receiving an oscillatory input with a varying frequency, a neuron is said to exhibit resonance when the voltage response peaks at a nonzero input frequency. In general, the oscillatory input is injected with relatively low amplitude therefore producing a quasilinear voltage response, which is almost symmetric with respect to the holding potential. In addition, a quasilinear response has its voltage number of cycles equals the number of input cycles. However, when the amplitude is increased, asymmetries in the voltage response envelope arise. Since the impedance profile is essentially an average of the upper and lower voltage response envelope profiles, in the presence of asymmetries, this impedance fails to capture important properties of these envelopes. Of particular importance is the case where the maximum of the upper envelope and the minimum of the lower envelope occur at different frequencies. In these cases, the impedance and the upper envelope predict different characteristic frequencies at which the subthreshold responses are communicated optimally to the spiking regime. Although asymmetries in the voltage response to oscillatory inputs have been observed both theoretically and experimentally,^14,17,35 their properties and the mechanisms that give rise to these asymmetric envelope profiles have not received much attention. To our knowledge, this is the first systematic study that addresses these issues.

Understanding how the neuronal intrinsic properties control a given response pattern can have implications for neuronal information processing in addition to an accurate prediction of the response spiking patterns to similar inputs (communication of oscillatory inputs to the spiking regime). One of them is for the postinhibitory rebound (PIR) phenomenon by which inhibitory inputs give rise to spiking activity in an otherwise silent neuron. PIR and resonance have been associated in the literature³⁶ since the currents that produce resonance (e.g., $Ih$⁠) are involved in PIR. Experimental observations agree that network resonance can emerge by inhibitory communication by means of PIR.³⁷ While the precise link between these two phenomena is not known, we speculate that the characteristic time scale associated to the trough of the lower response envelope plays a role in determining the PIR time scale, and therefore sheds light on the cellular mechanisms controlling the communication of inhibitory inputs to the spiking regime. More research is necessary to develop these ideas in a more precise fashion.

We unfolded the impedance profile into two quantities, the upper and lower impedance profiles, by taking the absolute values of the differences between the upper/lower voltage envelopes and the holding potential and normalizing them by the input amplitude. These metrics allowed a further characterization of asymmetry properties in terms of the model parameters. More specifically, we manipulated the nonlinearities from three different perspectives: (i) biophysical, by looking at the different representative ionic currents (⁠$Ih$⁠, $IM$⁠, $INaP$⁠, and $IKir$⁠); (ii) nonlinear voltage dependencies, by looking at the type of nonlinearities present in the model as captured by the voltage-nullclines (quadratic, cubic); and (iii) nonlinear mechanisms, by looking at the role that the nonlinearities play in different voltage ranges by using a piecewise-linear model that allows for changes of the nonlinearity type in one region but not in others, therefore allowing a more detailed analysis. This combined approach produced rules determining the effects that different types of ionic currents have on the upper and lower impedance profiles. These predictions are amenable for experimental testing using the dynamic clamp technique.^38,39

From the biophysical point of view we focus on how different resonant (⁠$Ih$ and $IM$⁠) and amplifying (⁠$INaP$ and $IKir$⁠) currents shape the upper/lower impedance profiles. We chose these four currents because they are representative of the four different type of scenarios. $INaP$ and $IM$ are depolarization-activated, but $INaP$ is inward, while $IM$ is outward.⁴⁰ Similarly, $Ih$ and $IKir$ are hyperpolarization-activated, but $Ih$ is inward, while $IKir$ is outward.⁴⁰ In particular, we characterized the effect that the holding potential has in determining the upper/lower impedance properties since the holding potential controls how much the neuron dynamics is affected by the activation curve of the currents. Depending on where the neuron’s voltage is located, the ionic current has a higher effect either depolarizing or hyperpolarizing the neuron and eventually causing resonance or amplification on the upper or lower response envelope, respectively. Notably, the gating variables showed an unexpected frequency-dependent pattern which is related with the half-activation value of the sigmoid function: if the half-activation value is above or below the resting potential, the ionic current will be either activated or deactivated with increasing frequency, respectively. This effect is not completely understood and is probably related to strong nonlinearities associated to the gating variables. However, to the best of our knowledge, this is the first time such a phenomenon is discussed and it could potentially elucidate new functions for the ionic currents. This calls for additional research.

From the dynamics point of view, we provided a geometrical interpretation of the role of the asymmetries in shaping the voltage responses using a dynamical systems analysis of both the biophysical model and a simplified piecewise-linear model. We show how the asymmetric responses are linked to the model nonlinearities, which are captured by the nullclines in the phase-plane diagram, and the time scale separation between the participating variables. We found that higher current amplitudes are more likely to generate asymmetric responses because the systems’ trajectories are more prone to reflect nonlinearities from the plane since they cover larger areas of the phase-plane where the nullcline nonlinearities on both sides of the holding potential differ more than in a close vicinity of it. Moreover, by taking advantage of the ability to modify one branch of the voltage nullcline, while keeping the other intact in piecewise linear models, we learned that asymmetries exist even in a system without currents with kinetics (activation curves). Markedly, this result demonstrates that the gating variables are not essential to obtain asymmetrical responses as they naturally emerge in simplified systems as long as the input amplitude, the nullclines position, and the time-scale separation allow such effect.

We have carried out simulations with longer ZAP and the same frequency range for representative parameter values and we have observed no difference with the results discussed in the paper (not shown). We have also carried out simulations in the reverse frequency order (starting with the higher frequencies and ending with the lowest ones) and, again, we have observed no differences (not shown). This indicates that, as expected, there are no intrinsic cellular processes that are too slow as to fail to be captured by the ZAP currents we have used in this work.

Our results are in line with recent experimental observations.^14,17 In Ref. 14, for example, the authors extensively studied how different currents shape the upper and lower envelopes of the voltage profile in neurons in the crab pyloric CPG. Recently, a clear asymmetric subthreshold voltage response was experimentally recorded in Ref. 35 (see Figs. 3 and 4 of that article). There, by applying ZAP currents to frog vestibular neurons, the authors observed asymmetrical voltage responses where the resonant frequency of the upper envelope increases with depolarization, i.e., $Δf$ changes with $Vhold$ as observed in our simulations. That work allowed a better classification of neurons in the frog vestibular system by means of their response properties. Another experimental example where asymmetries were recorded is found in Ref. 41, where the disruption of M-currents with the selective blocker XE991 in medial entorhinal cortex cells mediated such pattern (see Fig. 7 of that article). The latter example indicates a clear channel dependency on the asymmetrical response of a neuron as we have reported here by the use of computational models. Another key point is that in all the models we have used, the number of output cycles coincide with the input cycles meaning that this linearity principle is maintained. Further research should be conducted to clarify impedance profiles in which this linearity principle is broken.

Our observations suggest important implications for signal processing in neurons. First, since the resonance peak difference $ΔZ$ can be either negative or positive, we speculate that an oscillatory input to a neuron could lead it to preferably respond with a hyperpolarizing or a depolarizing displacement in the membrane potential. Second, the existence of the resonance frequency shift $Δf$ suggests that the resonance properties of a neuron may be related to additional functions than increasing neuronal excitability. For $Δf≠0$⁠, the system has two characteristic frequencies, one depolarizes, the other hyperpolarizes. However, we found the $Δf$ values to be not very large, which raises the question of whether there are other channel combinations, possibly including calcium concentration dependent channels, that could further increase $Δf$ or amplify its effect when the neuron is embedded in a network. In this regard, we propose a role for neuronal response modulation under oscillatory inputs for different values of the $Ih$ kinetics parameters and its voltage regulated activation. $Ih$ displays a high variability with the time constant of its kinetics spanning from tens of milliseconds to several seconds.^27–31 This could be a potential mechanism for output modulation.

The characterization of $Z+(f)$ and $Z−(f)$ has implications for model parameter estimation of neuronal systems. When measuring the properties of ionic currents, the information about each one is usually obtained separately by means of experimental procedures (e.g., patch clamp and voltage clamp). However, these experiments do not provide information about the interaction between currents and it has been shown that the effects of these interactions cannot be understood in terms of the sum of the effects of the participating currents. In addition, because the number of model parameters are typically larger than the information obtained by applying constant pulses of current/voltage, the problem is unidentifiable (parameter values cannot be uniquely determined). Resonance experiments provide additional, independent dynamic information that can be used to help characterizing the ionic currents and their interactions, more so when the responses exhibit asymmetries.

As it has been shown, the impact of subthreshold resonance ranges from suprathreshold neuronal properties (e.g., firing rate) to network behavior and learning.^26,32–34 Aligned with this, our results suggest mechanisms by which a network could switch between different states not only according to its activity and input frequency but also to hyperpolarization and depolarization. This can have implications on neuronal spiking properties in a wide range of circumstances. These phenomena deserve further exploration.

See the supplementary material for selected voltage traces from the same simulation in Fig. 6.

This paper was developed within the scope of the IRTG 1740/TRP 2015/50122-0, funded by DFG/FAPESP. This work was partially supported by the Research, Innovation and Dissemination Center for Neuromathematics (FAPESP Grant No. 2013/07699-0). This work was partially supported by the National Science Foundation Grant No. DMS-1608077 (HGR, RFOP). R.F.O.P. was also supported by a FAPESP Ph.D. scholarship (Grant No. 2013/25667-8), C.C.C. is supported by a CAPES Ph.D. scholarship, V.L. is supported by a FAPESP M.Sc. scholarship (Grant No. 2017/05874-0), R.O.S. is supported by a FAPESP Ph.D. scholarship (Grant No. 2017/07688-9), and A.C.R. is partially supported by a CNPq fellowship (Grant No. 306251/2014-0). This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES)—Finance Code 001.