Frequency Preference Response to Oscillatory Inputs in Two-dimensional Neural Models: A Geometric Approach to Subthreshold Amplitude and Phase Resonance

We investigate the dynamic mechanisms of generation of subthreshold and phase resonance in two-dimensional linear and linearized biophysical (conductance-based) models, and we extend our analysis to account for the effect of simple, but not necessarily weak, types of nonlinearities. Subthreshold resonance refers to the ability of neurons to exhibit a peak in their voltage amplitude response to oscillatory input currents at a preferred non-zero (resonant) frequency. Phase-resonance refers to the ability of neurons to exhibit a zero-phase (or zero-phase-shift) response to oscillatory input currents at a non-zero (phase-resonant) frequency. We adapt the classical phase-plane analysis approach to account for the dynamic effects of oscillatory inputs and develop a tool, the envelope-plane diagrams, that captures the role that conductances and time scales play in amplifying the voltage response at the resonant frequency band as compared to smaller and larger frequencies. We use envelope-plane diagrams in our analysis. We explain why the resonance phenomena do not necessarily arise from the presence of imaginary eigenvalues at rest, but rather they emerge from the interplay of the intrinsic and input time scales. We further explain why an increase in the time-scale separation causes an amplification of the voltage response in addition to shifting the resonant and phase-resonant frequencies. This is of fundamental importance for neural models since neurons typically exhibit a strong separation of time scales. We extend this approach to explain the effects of nonlinearities on both resonance and phase-resonance. We demonstrate that nonlinearities in the voltage equation cause amplifications of the voltage response and shifts in the resonant and phase-resonant frequencies that are not predicted by the corresponding linearized model. The differences between the nonlinear response and the linear prediction increase with increasing levels of the time scale separation between the voltage and the gating variable, and they almost disappear when both equations evolve at comparable rates. In contrast, voltage responses are almost insensitive to nonlinearities located in the gating variable equation. The method we develop provides a framework for the investigation of the preferred frequency responses in three-dimensional and nonlinear neuronal models as well as simple models of coupled neurons.

the gating variable equation. The method we develop provides a framework for the investigation of the preferred frequency responses in three-dimensional and nonlinear neuronal models as well as simple models of coupled neurons.
The relation between intrinsic STOs and subthreshold resonance is still an open question. For some neuron types, STOs and subthreshold resonance have been shown to result from the same mechanism [12]. However, theoretical and experimental studies have demonstrated that they are not equivalent phenomena [11,28]: neurons may exhibit one and not the other [10,11,21]. Furthermore, standard calculations for linear models show that their natural (intrinsic) and resonant frequencies do not generally coincide except in some, rather restricted, parameter regimes [11,29] (see also our discussion in Appendix A.3).
The phase-shift (or phase) of the neuronal voltage response to subthreshold oscillatory input currents has received less attention that the corresponding amplitude response [11,14,29]. This despite the fact that phases are expected to play a major role in determining the synchronization properties of neuronal networks [30]. A zero-phase response indicates that both voltage output and current input peak at the same time, thus generating in-phase synchronized patterns. We use the term phase-resonance to refer to the ability of neurons to exhibit a zero-phase response to oscillatory inputs at a non-zero (phase-resonant) frequency. The resonant and phaseresonant frequencies do not generally coincide [29] for neuronal models (they do so for circuits in parallel but not for circuits in series as neuronal models are). In addition, resonance may occur in the absence of phase-resonance [29] (see also our discussion in Appendix A.3).
The properties of subthreshold resonance have been investigated in many systems [9,11,14,15,20,21,25,29]. Theoretical studies have focused on simulations of conductance-based models and the analysis of the impedance profiles (curves of amplitude and phase-shift as a function of the input frequency) for the corresponding linearized models around resting potential [11, 14-16, 18, 31-35]. However, the mechanisms underlying the generation of resonance and phase-resonance in neurons are not fully understood. This is partly because adequate tools are lacking.
These mechanisms can be addressed from two different, but complementary perspectives: biophysical and dynamic. The former focuses on the role of the ionic currents and their biophysical properties in shaping the neuron's voltage response, and it has been discussed in terms of the so-called resonant and amplifying ionic currents (see Sect. 2.5) [10,11]. The latter focuses on (i) the geometric properties of the neuronal models in terms of the nullclines and phase planes, and (ii) the interaction between the neuron's intrinsic time scales and the time scales associated with the input currents to produce optimal voltage responses in both amplitude and phase.
In [29] we have identified the basic biophysical mechanisms of generation of resonance and phase-resonance in two-dimensional linear and linearized conductancebased models, and we have conducted a thorough study of the properties of the voltage response in terms of the biophysical parameters. In particular, we have shown how changes in the maximal conductances affect the resonant and phase-resonant frequencies and other attributes of the impedance amplitude and phase profiles in ways that are not always intuitive.
The goal of this paper is to investigate the dynamic mechanisms that give rise to resonance and phase-resonance. Specifically, we aim to identify the basic dynamic and geometric principles that govern the generation of these phenomena in order to understand (i) how the resonant and phase-resonant frequencies are selected, (ii) how the voltage response is amplified at the resonant frequency band, (iii) how these properties are affected by changes in parameters (e.g., maximal conductances, time constants), and (iv) how resonance and phase-resonance are related to intrinsic STOs.
We use dynamical systems tools and numerical simulations, and we extend the classical phase-plane analysis approach to generate a tool (envelope-plane diagrams) that enables us to address the mechanistic issues described above. For simplicity, in this paper we focus on two-dimensional linear systems and we illustrate how the ideas we develop here can be used to investigate linearized conductance-based models and nonlinear models.
From a dynamic perspective, we view the properties of the forced oscillations as reflecting the interaction between the forced neuron's intrinsic time scales, which result from a combination of its biophysical properties and the time scales of the input. The time scales that emerge from these interactions determine the amplitude and phase of the voltage response for each value of the input frequency. The resonant and phase-resonant frequencies correspond to the emergent time scales that allow the neuron to maximize its voltage amplitude response and to peak in phase with the input, respectively. Therefore, we aim to elucidate how these emergent time scales are generated, how intrinsic and emergent time scales are related, and how they are affected by changes in the model parameter.
It is clear that for linear systems the voltage response can be computed analytically. However, due to its complexity, the dependence of the voltage response properties with the model parameters is not straightforward, not even for linear twodimensional systems [29]. Additionally, due to the same complexity, it is difficult to extract from the closed-form solutions a mechanistic understanding of the neuron's resonant properties in terms of the time scales. The use of geometric tools aids in this effort. Since the phase plane contains information about the biophysical properties of neurons through the structure of the nullclines and time scales, dynamical systems tools allow us to make statements that are valid for generic classes of biophysical models.

Conductance-Based Models
We consider conductance-based models of Hodgkin-Huxley type [36]. The currentbalance equation is given by where V is the membrane potential (mV), t is time measured in msec, C is the membrane capacitance (µF/cm 2 ), I app is the applied bias (DC) current (µA/cm 2 ), I in (t) is a time-dependent input current (µA/cm 2 ), is the leak current, and I j (j = 1, 2) are ionic currents of the form with maximal conductance G j (mS/cm 2 ) and reversal potentials E j (mV), respectively. The ionic currents (2) we consider here are restricted to have a single gating variable x j and to be linear in x j . This choice is motivated by the persistent sodium (I Nap ), h-(I h ) and M-type (I M ) currents found in several neurons that exhibit subthreshold resonance [15,[37][38][39] (see Appendix B). All gating variables x obey a first order differential equation of the form where x ∞ (V ) and τ x (V ) are the voltage-dependent activation/inactivation curves and time scales, respectively. For external sinusoidal inputs we use the following notation: where T = 1000 msec and [f ] = Hz.
In this paper we focus on two-dimensional models having one dynamic gating variable (x 1 ) and, possibly, an additional gating variable evolving on a fast time scale for which the adiabatic approximation x 2 = x 2,∞ (V ) is made. Additional fast currents can be including without significantly changing the formalism. The investigation of three-dimensional systems is beyond the scope of this paper.

Linearized Conductance-Based Models
We follow Richardson et al. [11] and linearize the autonomous part of system (1)-(3) The linearized equations are [11] where the effective conductances g L and g 1 , and time constantτ 1 are defined by The sign of the effective ionic conductances g j (j = 1, 2) determines whether their associated gating variables are resonant (g j > 0) or amplifying (g j < 0) [10,11]. The effective conductance g L includes information not only about the original leak conductance G L but also about the voltage-dependent ("passive") terms of the ionic currents. The ionic currents I 1 and I 2 each contribute a positive term to g L . The ionic current I 2 contributes to g L with an additional term (g 2 ) due to its instantaneous dynamics. If x 2 is amplifying and g 2 < 0 is large enough in absolute value, then g L is negative. Note that the gating variable w in (5) has units of voltage.

Rescaled Linearized Models
We rescale system (6)- (7) to aid in the analysis and to reduce the number of parameters that effectively govern its dynamics (without loss of information). The rescaling we use here focuses on the geometric properties of the system and the time-scale separation between the participating variables, and is amenable to analysis using dynamical systems tools (phase-plane analysis). A different scaling, appropriate for addressing biophysically related questions, has been used in [11,29].
We define the following dimensionless time and (dimensional) voltage variables: andÎ Substituting into (6)-(7) and dropping the "hat" signs we obtain Geometrically, the parameter α is the slope of the w-nullcline and can be thought of as representing the strength of the gain of the feedback in the linearized system. The parameter represents the time-scale separation between v and w.
Since for resonant gating variables g 1 > 0, the sign of both α and depends on whether g L is positive or negative. In the absence of fast amplifying currents (G 2 = g 2 = 0), g L > 0 and then both α > 0 and > 0. When an amplifying current is present and its contribution to g L is small enough, the sign of both α and remains positive ( Fig. 1, dashed curves). However, when stronger contributions of the fast amplifying current causes g L to be negative, the sign of both α and are also negative ( Fig. 1, solid curves). Resonance occurs in both cases [11,29]. Since resonance becomes amplified as g L decreases [29], we expect resonance to be more amplified for negative values of both and α as compared to positive ones. The cases including values of α and having different signs are excluded from this study since the underlying autonomous system is either unstable (saddle) or stable (node) but does not exhibit resonance (Fig. 3). In both cases x 1 is amplifying (g 1 < 0). The case < 0 and α > 0 requires that both g L < 0 and g 1 < 0, while the case > 0 and α < 0 requires that g L > 0 and g 1 < 0.

Impedance and Impedance-Like Functions
The voltage response of a linear system receiving sinusoidal current inputs of the form (4) is given by

Fig. 2
Schematic diagrams of the impedance (a) and phase (b) profiles (impedance and phase as a function of the input frequency f ). a1 Band-pass filter (resonance). a2 Low-pass filter (no resonance). b1 Zero-frequency phase crossing (phase-resonance). b2 Monotonically increasing and positive phase (no phase-resonance). a The resonant frequency f res is the input frequency f at which the impedance Z(f ) reaches its maximum Z max . The resonance amplitude Q Z = Z max − Z(0) measures the resonance power. The half-width frequency band Λ 1/2 is the length of the frequency interval in between f res and the input frequency value at which Z(f ) = Z max /2, and measures the system's selectivity to incoming frequencies close to f res . b The phase-resonant frequency f phas is the zero-crossing phase frequency. The minimum phase φ min measures the magnitude of the negative phase where A out is the amplitude and φ is the phase-shift (or phase), defined as the difference between the peaks of the current input I in (t; f ) and the voltage output V out (t; f ). Linear systems exhibit resonance if there is a peak in the amplitude of the impedance function Z(f ) given by (17) at some positive (resonant) frequency f res . In what follows, we will refer to impedance amplitude simply as the impedance Z(f ). In Fig. 2a we show representative graphs of the impedance function Z(f ) for a model that does (panel a1) and does not (panel a2) exhibit resonance. We characterize the impedance profiles using four parameters: (i) the resonant frequency f res , (ii) the maximum impedance Z max = Z(f res ), (iii) the resonance amplitude Q Z = Z max − Z(0), and (iv) the halfwidth frequency band Λ 1/2 , defined as the frequency interval in between f res and the input frequency value at which Z(f ) = Z max /2. Λ 1/2 is a measure of the frequency selectivity. Neurons have a higher selectivity to inputs with frequencies around f res the sharper the graph of Z(f ). In Fig. 2b we show two representative graphs of the phase φ(f ) where φ vanishes at a non-zero value of f (panel b1) and φ is always positive (panel b2). We refer to the ability of the neuron to exhibit a zero-phase frequency response at a non-zero frequency as phase-resonance and to the corresponding frequency as the phase-resonant frequency f phas . In panel b1, f phas > 0. The voltage response is "advanced" and "delayed" for lower and higher frequency inputs, respectively. Although phase advance and phase delay are ambiguous concepts to describe phase differences between inputs and outputs in oscillatory systems, we still use them since typical phase differences lie in the range (−π/2, π/2). In panel b2, f phas = 0, that is, the voltage response is delayed for all values of f . We characterize the phase profiles using two parameters: (i) the phase-resonant frequency f phas , and (ii) the minimum phase φ min . For nonlinear systems, or for linear systems with non-sinusoidal inputs, (17) does no longer provide an appropriate definition of the impedance function. Here we assume that the voltage output is periodic and has the same frequency as the input and we use the following definition: where V max (f ) and V min (f ) are the maximum and minimum of the oscillatory voltage response V out (f ) for each value of the input frequency f . For linear systems receiving sinusoidal inputs, (18) and (17) are equivalent. The resonant frequency f res is the peak frequency of Z(f ) in (18). Similarly to the linear case φ(f ), the phase is computed as the distance between the peaks of the output and input normalized by the period. Note that (18) can be thought of as a filtered version of the impedance function computed using the so-called ZAP functions [10,11,28] that sweep through a given range of frequencies continuously over time.

Resonant and Amplifying Ionic Currents
Biophysically, subthreshold resonance has been argued to result from a combination of low-and high-pass filter mechanisms that have been described in terms of neural currents [10]. RC circuits act as low-pass filters. (As the input frequency increases the voltage amplitude response of passive neurons decreases from its resistance value, for zero input frequency, to zero, for input frequencies approaching infinity.) Ionic currents, or more precisely their associated gating variables, have been classified into resonant (g 1 > 0) and amplifying (g 1 < 0) [10,11]. Resonant gating variables (e.g., hyperpolarization-activated h-currents [15][16][17] and slow potassium, Mtype currents [14]) have the ability to create resonance by opposing voltage changes (negative feedback). Amplifying gating variables (e.g., persistent sodium currents [14][15][16][17] and high-threshold calcium currents [10]) generate a positive feedback effect that enhances voltage changes but they do not create resonance [10]. Some currents such as the low-threshold calcium current I T have both resonant and amplifying gating variables [10].
The interaction between resonant and amplifying currents in shaping the neuronal voltage response (impedance and phase profiles) to oscillatory input currents is complex [29], sometimes non-intuitive, and involves not only changes in the maximal amplitude of the voltage response but also in other attributes including the resonant and phase resonant frequencies. The role that different types of resonant and amplifying currents play in shaping the impedance profile has been recently clarified in [29]. An important outcome of this study is that the standard classification "resonant vs. amplifying" does not capture in its entirety the effect that changes in their biophysical parameters have on the shapes of the impedance and phase profiles. For instance, while hyperpolarization-activated (h-or I h ) and M-type slow potassium (I Ks ) currents have qualitatively similar effects on the impedance profile as the corresponding ionic conductances (G h and G Ks ), they may have opposite effects in the presence of a persistent sodium current (amplifying). Specifically, in the latter case, an increase in G h leads to an amplification of the voltage response, while an increase in G Ks causes an attenuation of the voltage response. Additionally, an increase in the time constant associated to the gating variable, which increases the system's time-scale separation, does not only affect the resonant and phase-resonant frequencies, but it may also produce a significant amplification of the voltage response.

Stability and Resonant Properties in the α-Rescaled System
Here we review the stability and resonant properties of the rescaled equations (14)- (15) in terms of the parameters α and for later use. These properties have been discussed in [11,29] for a different rescaling based on dimensionless effective conductances.

Stability Properties for the Autonomous System and Intrinsic Oscillations
We first consider system (14)- (15) with A in = 0. The fixed point is given by (v,w) = (0, 0) and the eigenvalues are given by (Appendix A) The stability diagram in the α-parameter space is presented in Fig. 3a1. From (19), the fixed point is a focus if (1 − ) 2 − 4α < 0 and a node otherwise. Foci are stable if 1 + > 0. Nodes are stable if 1 + > 0 and (1 + α) > 0. If (1 + α) < 0 the fixed point is a saddle. The natural frequency of the damped oscillations (for stable foci) is given by

Impedance, Phase, Resonance, and Phase-Resonance
The impedance function (17) for system (14)- (15) with A in > 0 is given by (Appendix A)  (15). a Stability and resonance diagrams in the α-parameter space. a1 Stability diagram for the autonomous system. The blue curves separate between regions with different stability properties. a2 Resonance diagram. The red curves separate between regions where the system does (above and below) and does not (middle) exhibit resonance. a3 Superimposed stability (blue curves) and resonance (red curves) diagrams showing that intrinsic oscillations and resonance may occur in the absence of the other. b Natural (f nat ), resonant (f res ), and phase-resonant (f phas ) frequencies as a function of α (b1) and (b2) illustrating that these characteristic frequencies have different values (b1 and b2) and different monotonic properties (b2) as the model parameters change. c Maximum impedance (Z max ) and resonance amplitude (Q Z ) as a function of α (c1) and (c2) illustrating the two basic mechanisms of generation of resonance in 2D linear systems. c1 As α increases, resonance results from a decrease in both Z(0) and Z max with Z(0) decreasing faster than Z max . c2 As decreases (time-scale separation increases), resonance results from an increase in Z max with Z(0) fixed The resonant frequency, if it exists, is given by System (14)- (15) exhibits resonance if Ω res > 0. In order Ω res to be defined, α(α + 2 + 2) ≥ 0 and 2 α(α + 2 + 2) > 2 . The first condition is satisfied either if For the second condition to be satisfied, This yields either These regions are illustrated in the resonance diagram in Fig. 3a2. Figure 3a3 illustrates that resonance may occur in the absence of intrinsic oscillations and vice versa. The phase φ for system (14)- (15) is given by The phase-resonant frequency, if it exists, is given by The natural, resonant, and phase-resonant frequencies do not necessarily coincide (Fig. 3b). In fact, they rarely do so. Figure 3c illustrates the two basic mechanisms of generation of resonance for 2D linear systems [29] in terms of α and : (i) as α increases, resonance emerges from a combined and unbalanced decrease in both Z(0) and Z max with Z(0) decreasing faster than Z max (panel c1), and (ii) as decreases (time-scale separation between v and w 1 increases), resonance emerges from an increase in Z max with Z(0) unchanged.

The Structure of the "Oscillatory Phase Plane"
Here we consider system (14)- (15) with A in ≥ 0. For the type of analysis we present in this paper it is useful to rescale time by defininĝ t = f t (28) in order to separate the effect of the input frequency f from the input's time dependence. Substituting (28) into (14)- (15), and dropping the "hat" sign from the timet, we obtain dw where the sinusoidal input has the same frequency (1 cycle per T units of time) for all values of the input frequency f . The latter affects the speed of the trajectories in the phase plane (i.e., the speed of the system's response) without affecting the direction of the underlying vector field. Here we use T = 1000.
The v-and w-nullclines for the unforced system (A in = A t = 0) are given by N v (v) = −v and N w (v) = αv, respectively. The stability properties of this autonomous system have been discussed in Sect. 3.1.1. For A in > 0, system (29)- (30) is no longer two-dimensional. In our analysis, we will think of the projection of the zero-level curve for (29) for each value of t as the v-nullcline N v (v) for the autonomous system forced by the sinusoidal input A t : For t = 0 (or any multiple of t = 500), The moving v-nullcline reaches these two lines at t = 250 and t = 750, respectively. As N v,t moves so does its intersection with the w-nullcline N w generating a "moving fixed point" which oscillates with frequency 1 between the endpoints (v 250 ,w 250 ) (A 250 = A in ) and (v 750 ,w 750 ) (A 750 = −A in ), reaching the origin three times within a cycle at t = 0, t = 500, and t = 1000. This moving fixed point together with the moving nullclines organize the dynamics of the forced system. Specifically, the stable fixed point (either a node or a focus) acts as a moving target for trajectories, which track its motion by evolving with a f -dependent speed.
The dynamics of the forced system (29)-(30) results from the interaction between the vector field of the autonomous system (A in = A t = 0) and the forcing term (A t ) that causes the cyclic motion of both the v-nullcline and the fixed point (34) in addition to a cyclic change in the direction field. For each point in the phase plane the value of the direction field is independent of f . However, the trajectories' fdependent speed causes them to sweep different distances in a given unit of time. Since they reach different points, subject to different values (magnitude and direction) of the vector field, trajectories describe different curves for different values of f . The resonant frequency, if it exists, is the value of the input frequency f for which the corresponding limit cycle trajectory has the maximal amplitude in the v-direction, provided this amplitude is larger than the instantaneous amplitude corresponding to the effective resistance Z(0).

Transient Dynamics for Instantaneously Forced Systems
As a first step in our investigation we describe the response of system (29)-(30) to instantaneous constant inputs (with no associated time scales) in order to understand how the system's intrinsic time scales depend on the parameters α and , and how these parameters contribute to the amplification or attenuation of the voltage response to these instantaneous changes. We will use this knowledge to explain the dynamics of the forced system when the inputs do have associated times scales.
For the purpose of this part we will consider constant values of A t in system (29)- (30). Without loss of generality, we assume that (i) A t instantaneously changes from A t = 0 to either A t = 1 or A t = −1, (ii) this change occurs at t = 0, and (iii) f = 1.

Transient Dynamics and Effective Time Scales
An instantaneous change in A t causes an instantaneous displacement of the vnullcline on the phase plane, and consequently an instantaneous change in the location of the fixed point. The vector field changes accordingly, and can also be thought of being instantaneously displaced.
The stability properties of the fixed point, however, remain unchanged since the system (29)-(30) is linear and, therefore, its associated Jacobian matrix is independent of A t .
Consequently, for different values of A t , trajectories at a given initial point have the same long term behavior as they approach the corresponding fixed points, but they have different initial transient behaviors (Fig. 4). For the purpose of this paper, we are particularly interested in the initial transient behavior caused by changes in A t . As we explain in more detail below, this transient behavior rather than the long term behavior is the one that plays an important role in determining the amplitudes of the limit cycle trajectories for the time-dependent forced systems.
The trajectories' transient behavior relevant to this paper corresponds to the time elapsed from the initial point (red dots in Fig. 4) until they cross the corresponding v-nullcline (dashed-red line), reaching their maximum value for v. The effect of different values of A t on the behavior of trajectories initially at the same point can be understood by looking at the horizontal and vertical components of the vector field respectively. As A t increases, D v increases and the motion in the horizontal direction strengthens relative to the motion in the vertical direction. Consequently, the larger A t the "more horizontal" the trajectories' transient direction of motion. From a different perspective, trajectories starting at the same initial point in the "standard" coordinate system whose origin is (0, 0), start at different initial points in translated coordinate systems whose origins are the fixed points (v t ,w t ) determined by A t through (34). The larger A t , the larger the distance between the trajectory's initial point and the corresponding v-nullcline, and so they are less affected by it and they can transiently move in more horizontal directions. Figure 4 also illustrates the roles of α and in determining both the effective time scales for system (29)- (30) and the amplitude of the voltage response (maximum value of v for any given trajectory) to constant perturbations (A t ). For = 0.01 (top panels) the time scales are well separated and the system is fast-slow. Trajectories move first fast in almost horizontal directions, then they slow down as they approach the v-nullcline, and finally they reverse direction and move towards the fixed point Each panel shows superimposed phase-planes diagrams for three different constant values of A t (= 0, 1, −1) generating three v-nullclines (solid-red for A t = 0, dashed-red for A t = 1 and A t = −1) and three fixed points (blue dots on the intersections between the red and green lines). The w-nullcline (green line) is common to all values of A t . Solid-red line: v-nullclines for A t = 0. Dashed-red lines: v-nullclines for A t = 1 (above) and A t = −1 (below). Red dots at (0, −1) and (0, −0.5): representative initial conditions. Solid-blue lines: trajectories initially located at these initial points. Each trajectory emerging from the red dots corresponds to a different value of A t and converges to the corresponding fixed point. The fixed points in panels a-top, a-middle and b-top are stable nodes and the fixed points in panels a-bottom, b-middle and b-bottom are stable foci in close vicinities of the displaced v-nullcline. This effect is more pronounced for smaller values of α (compare top panels a and b) because of the product α in (30), which affects the effective time scales. As α decreases so does the slope of the vnullcline, thus "pressing" trajectories that become more "compacted". As increases (from the top to the bottom panels), the separation of time scales decreases, trajectories become more "rounded", and reach lower maximal values of v.
An important observation is that, qualitatively, the initial transient behavior of trajectories in Fig. 4 is not dependent on whether the fixed point is a stable node or focus but rather on the magnitude of the differences between the corresponding values of α and . Specifically, trajectories with close enough values of both α and will

Amplification of the Voltage Response to Instantaneous Constant Inputs
The voltage response is amplified by a given parameter if v is able to reach a higher value when this parameter changes.  Fig. 5). These results reflect the roles of the resonant and amplifying currents in determining the amplification of the voltage response through the dimensionless parameters α and given by (12). For example, an increase in time constant τ 1 causes a decrease in (time-scale separation), leading to an amplification of the voltage response. The explanation of the effects of other biophysical parameters (e.g., g L and g 1 ) through the dimensionless parameter α is less straightforward since α has been used as a voltage scale. A thorough study of these effects has been carried in [29]. Figure 4 also demonstrates that voltage responses are more amplified for negative values of α and than for positive ones (compare panels a and c). This is consistent with the fact that negative values of α and are obtained for negative values of the effective conductance g L , which in turn reflects the presence of a strong, fast amplifying current (I 2 ) in the biophysical model.

Voltage Amplitude Response to Sinusoidal Inputs at Different Frequencies:
Resonance and Phase-Resonance Here we build up on the ideas discussed in Sect. 3.3 to investigate the mechanism of selection of the resonant and phase-resonant frequencies in the two-dimensional system (29)- (30) in response to sinusoidal inputs of the form A t (31). Without loss of generality we consider the canonical case A in = 1. With this choice, for each value of f the impedance function (18) coincides with the maximum value of v; i.e., v max (f ) = Z(f ).
As t increases, A t changes and the v-nullcline N v,t (32) moves cyclically, parallel to itself between the lines N + v (v) and N − v (v) (33) (dashed-red lines in Fig. 4) with frequency equal to 1. The fixed point (34) moves cyclically with the same frequency between the two corresponding extreme values. The input amplitude is coded by the distance between the moving fixed point and the origin. Trajectories respond to changes in A t by evolving according to the dynamics dictated by the underlying vector field (35). The trajectories' speed ν is given by the product of the magnitude of the vector field and f −1 At any given point, the underlying vector field is independent of f . However, the f -dependent speed (36) causes the resulting limit cycle trajectory to be f -dependent as explained at the end of Sect. 3.2. This dependence is reflected primarily on the shape and orientation of the limit cycle as we illustrate this in Fig. 6a for a system that exhibits resonance (see Fig. 6b) and phase-resonance (see Fig. 6c). For small values of f limit cycles are elliptic and very narrow, with their major axis lying on a vicinity w-nullcline. As f increases, the limit cycle first widens and rotates, then it narrows as it continues to rotate until the major axis becomes horizontal (on the v-axis), and finally the limit cycle shrinks until it collapses to a point as f → ∞. The system exhibits resonance because the maximal value of v for f ∼ 65 is larger than the maximal value of v for f → 0. Figure 6d shows a graph containing the curve generated by the points with maximum values of v on the limit cycles for continuously changing values of f ∈ [0, ∞) and the parameter values in Fig. 6a. We refer to this curve as the upper envelope curve of the voltage response. (The lower envelope curve, not shown, is determined by the minimum values of v on the f -dependent limit cycles, and is symmetric to the upper envelope curve with respect to the v-nullcline and w-nullclines.) We refer to the diagrams containing the envelope curves together with the v-and w-nullclines (green, solid-red, and dashed-red), as in Fig. 6d, as the envelope-plane diagrams.

Envelope-Plane Diagrams: Resonant and Phase-Resonant Responses
Envelope-plane diagrams contain geometric and dynamic information about a system's frequency response to oscillatory inputs, and they are the frequency analogs to phase-plane diagrams. Trajectories in the envelope-plane diagrams (upper and lower envelopes) are curves parameterized by the input frequency as trajectories in the phase planes are curves parametrized by time. Neither f nor t are explicit in the corresponding diagrams. The red-dashed lines quantify the maximal input amplitude. For f → 0, the envelope curve is at the intersection between the upper dashed-red line N + v and the w-nullcline. For f → ∞, the envelope curve is at the origin (limit cycle shrinking to a point).
The cusp ("horizontal peak") in the envelope curve in Fig. 6d corresponds to the peak in the impedance profile (Fig. 6b), and hence to the resonant frequency. At this cusp the voltage response is the largest across input frequencies, and larger than the response for f → 0: v max = (1 + α) −1 .
The point in the envelope curve tangent to the upper dashed-red line corresponds to the phase-resonant frequency for reasons that we will clarify later in the paper.

Low Input Frequencies Generate Quasi-one-dimensional Dynamics Along the w-Nullcline
For values of f 1 both equations in system (29)-(30) are very fast (dv/dt → ∞ and dw/dt → ∞), and then the limit cycle trajectory tracks the motion of the fixed point almost instantaneously (in a quasi-steady-state fashion). In the limit f → 0, the  = 0), respectively. The dashed-red lines represent the v-nullclines for A t = 1 (above), A t = 0 (middle), and A t = −1 (below). As t increases, the red nullcline moves cyclically between the two dashed-red lines. The blue dot indicates the initial location of the trajectory (t = 0 = 1000) limit cycle trajectory moves cyclically along the w-nullcline in between the points (34) with A t = ±A in = ±1. From (34), v max /A in = (1 + α) −1 , which coincides with Z(0) in (21). In Fig. 7 we present snapshots of the evolution of the limit cycle trajectory for f = 1 (corresponding to Fig. 6a, top-left panel).

High Input Frequencies Generate Quasi-one-dimensional Dynamics Along the Horizontal v-Axis
For large enough values of f 1, both (29) and (30) are very slow (dv/dt → 0 and dw/dt → 0), and then the limit cycle trajectory evolves with a very low speed according to (36). Relative to the trajectory's speed, the v-nullcline N v,t moves very fast. The time it takes the trajectory to cover a very small distance while tracking the motion of the v-nullcline is enough for the latter to reach its maximum level N + v , reverse direction, and intersect the trajectory on its "way back". This causes the trajectory to reverse direction. As a result, the limit cycle trajectory is not able to cover any significant distance away from the origin, and so the amplitude of the limit cycle is small as compared to other values of f (Fig. 6a, bottom-right panel). In the limit f → ∞, the amplitude of the limit cycle trajectory is zero (the limit cycle shrinks to the origin), which coincides with lim f →∞ Z(f ) = 0.

Voltage Response Amplification at the Resonant Frequency Band
For small and large enough values of f the corresponding limit cycle trajectories are constrained to move along quasi-one-dimensional directions (w-nullcline and vaxis, respectively). For intermediate values of f there is a transition in the shapes of the limit cycle trajectories between these two limit cases. Specifically, limit cycle trajectories are neither too fast nor too slow, and then, while they are "left behind" by the moving v-nullcline, they can take advantage of the two-dimensional vector field without being constrained to move in quasi-one-dimensional directions. For the appropriate values of α and this degree of freedom allows limit cycle trajectories to reach values of v max (f ) larger than v max (0), and so to exhibit resonance.
In Fig. 8 we show a sequence of snapshots of the evolution of the limit cycle trajectory for f = f res = 65 (Fig. 6a, bottom-left). The initial snapshot (t = 0) corresponds to a point (blue dot) on the limit cycle trajectory and A t = 0. As t increases the vnullcline moves and the limit cycle trajectory tracks its motion. The limit cycle' shape and amplitude result from the combined effect of A t , f , and the model parameters (α and ).
In Fig. 8 we are using a relatively small value of ( = 0.1). While, small values of are enough to account for the almost horizontal transient trajectories in autonomous systems (Fig. 4), they are not enough to account for the almost horizontal directions of motion for limit cycles trajectories in forced systems such as in Fig. 8. In fact, for the same value of but a different value of f (f = 1 instead of f = 65) the limit cycle trajectory moves in the (oblique) direction of the w-nullcline (Fig. 7).
To better understand how A t , f , and α affect the direction of motion of the limit cycle trajectory in Fig. 8, it is useful to go back and look at the initial transient segment of the evolution of trajectories for the autonomous systems presented in Fig. 4 (from the initial, red point until they cross the v-nullcline). Figure 4a (middle) corresponds to the same parameter values (α = 1 and = 0.1) as in Fig. 8. The increase in the constant input from A t = 0 to A t = 1 causes both an increase in the distance between the tip of the trajectory and the v-nullcline (from the solid line to the dashed-red line) and a strengthening of the horizontal component of the vector field D v (35) relative to its vertical component D w . As a result, during the initial transient interval, the larger A t > 0, the more horizontal is the trajectory's direction of motion.
A similar, but dynamic effect is partially responsible for the determination of the direction of motion in Fig. 8. Trajectories move faster in more horizontal directions as A t increases during the ascending phase (t = 0 to t = 250) because the continuous increase in A t causes a continuous increase in D v while D w remains almost unchanged. However, differently from the static case, in the dynamic case there is an opposite, dynamic effect caused by the decreasing distance between the tip of the limit cycle trajectory and the moving v-nullcline N v,t since D v decreases as the limit cycle trajectory approaches the v-nullcline. In the autonomous case, the nullclines are fixed, so the distance between the trajectory and N v,t depends only on the trajectory's speed (36). In the forced case, on the other hand, the trajectory and the v-nullcline N v,t are simultaneously moving, with different speeds. The resonance frequency f res is the input frequency for which their relative speed is optimal in the sense that it allows the limit cycle trajectory to move in a direction that maximizes (over the range of input frequencies f ) the distance (in the v direction) the limit cycle trajectory can cover before intersecting the moving v-nullcline, at which time they are forced to reverse direction. This causes a maximization (over the range of input frequencies f ) of the maximum value v max on the limit cycle trajectories. In Fig. 8 (f ∼ f res ), this intersection occurs for v max (∼ 1). The remainder of the limit cycle can be explained using similar ideas.
For values of f < f res , the limit cycle trajectory moves faster than for f = f res , and so the direction of motion is less horizontal causing this cycle trajectory to intersect N v,t at a higher point; i.e., for a lower value of v max (e.g., Fig. 6a, f = 48). For values of f > f res , the limit cycle trajectory moves slower than for f = f res . The rela-tive speed between this trajectory and N v,t is large enough to allow the limit trajectory to move in an almost horizontal direction, but, due to the limit cycle trajectory's lower speed, N v,t "returns" before this trajectory covers a large enough distance. Thus, the intersection between the trajectory and N v,t occurs for lower values of v max (e.g., Fig. 6a, f = 200).

Resonance and Intrinsic Oscillations Are Generated by Related, but not Identical Mechanisms
Intrinsic STOs and subthreshold resonance have been proposed to result from the same underlying mechanism [12]. However, in linear systems, resonance may occur in the absence of intrinsic oscillations and vice versa [11] (see Fig. 3a). Even for parameter values for which linear systems exhibit both resonance and intrinsic oscillations, the resonant and natural frequencies do not necessarily coincide (Fig. 3b).
The phase-plane analysis discussed above demonstrates that the resonant properties are not expected to be directly linked to the stability properties of the underlying autonomous system. The time scales corresponding to intrinsic oscillations, if they exist, are determined by the stability properties of the stable foci, which reflect the long term behavior of trajectories. In contrast, the time scales corresponding to the resonant frequency do not involve the stability properties of the fixed points of the underlying autonomous system (nodes or foci) but rather the initial transient behavior of trajectories as described above. This initial transient behavior is the link between the autonomous and forced systems. This transient dynamics is associated with either the onset of intrinsic oscillations (foci) or a "sag" (node) typically observed in the response of h-currents to constant inputs [11]. Importantly, as we demonstrated in Sect. 3.3, autonomous systems may display qualitatively similar transient behavior for nodes and foci, even though the long time behavior of the corresponding trajectories is qualitatively different (see Fig. 4b).

Phase Advance, Phase Delay and Phase-Resonance
In Fig. 3b we showed that the resonant (f res ) and phase-resonant (f phas ) frequencies do not necessarily coincide, and that resonance may occur in the absence of phaseresonance. For the parameters values in Fig. 6, f res = 65 and f phas = 48. Typical phase profiles show "phase advance" for f < f phas and "phase delay" for f > f phas (see Fig. 6c).
Geometrically, for the output and input to be synchronized in phase, the intersection between the limit cycle trajectory and the moving v-nullcline must occur when the v-nullcline reaches its highest level A t = A in (when the solid-red line in the moving phase-plane diagrams reaches the dashed-red line at t = 250). The phase is advanced when the limit cycle trajectory "peaks" before the input does, which requires the trajectory to evolve fast enough as compared to the speed of the v-nullcline. Conversely, the phase is delayed when the limit cycle trajectory intersects the v-nullcline "on its way back" (descending phase), after the v-nullcline has reached its highest level. For f = f res = 65 (Fig. 8), the limit cycle trajectory is still behind the v-nullcline at the time this v-nullcline reaches its highest level (t = 250), and the  = 0), respectively. The dashed-red lines represent the v-nullclines for A t = 1 (above), A t = 0 (middle), and A t = −1 (below). As t increases, the red nullcline moves cyclically between the two dashed-red lines. The blue dot indicates the initial location of the trajectory (t = 0 = 1000) phase is therefore delayed (the intersection occurs at a time in between the top-right and middle-left panels, t = 250 and t = 375, respectively).
For f = f phas = 48 (Fig. 9), on the other hand, the limit cycle trajectory intersects the v-nullcline when the latter is at its highest level (t = 250). Since the point of intersection corresponds to v = v max (on the limit cycle trajectory), the corresponding point on the envelope curve "touches" the top dash-red line. For all other frequencies, the intersection between limit cycle trajectories and the moving v-nullcline occurs away from the top dashed-red line, therefore the phase-resonant frequency corresponds to the point in the envelope curve that is tangent to the top dashed-red line. This tangent point is not necessarily the cusp that corresponds to the resonance frequency f res .
Specifically, the phase-resonance phenomenon depends on the ability of the limit cycle trajectory for f = f phas to track the v-nullcline fast enough so to reach v max before the v-nullcline reaches its highest level. But this does not preclude limit cycle trajectories for frequencies f = f phas to reach higher values of v due to the different directions of motion generated by these values of f . In fact, values of f > f phas cause limit cycle trajectories to evolve slower than for f = f phas , and so (within some range of values of f ) they can move along more horizontal directions, and thus reach higher values of v before intersecting the "returning" v-nullcline. This is the case for f = f res in Fig. 8. 3.5 Effects of α and on Resonance, Phase-Resonance, and Other Attributes of the Impedance and Phase Profiles Geometrically, the value of α determines one of the boundaries of the region (triangle) in the envelope-plane diagram where the envelope curves live (Fig. 6d). The other two boundaries are the v-nullcline for A t = 1 (dashed-red line) and the v-axis. The value of α also determines Z(0) = v max (0)/A in = (1 + α) −1 , which is the v-coordinate of the intersection between the w-nullcline and the v-nullcline for A t = 1. As we mentioned earlier, for resonance to occur, v max (f ) > v max (0) (on the envelope curve) for some range of values of f around f res . The value of plays a key role in determining the direction of motion of limit cycle trajectories analogous to its effect on the initial transient behavior of trajectories for the autonomous systems discussed in the context of Fig. 4. For fixed values of α, the smaller the more horizontal (and less rounded) is the direction of motion of the initial portion of the trajectories, and the larger the value of v these trajectories reach; i.e., the voltage response is amplified.
From the biophysical point of view, α and convey information about the effective conductances g 1 and g L , the capacitance C, and the time constantτ 1 (gating variable x 1 ) through (12). In [29] we use numerical simulations to conduct a through analysis of the effects of resonant and amplifying biophysical conductances on the determination of f res , f phas and other attributes of the impedance and phase profiles.
Here we use the envelope-plain diagrams developed in Sect. 3.4 to explain how the parameters α and affect the resonant and phase-resonant properties of (14)- (15), and how they contribute to the amplification of the voltage response. In Figs. 10 and 11 we examine the effects of changes in and α, respectively. In Figs. 12 and 13 we investigate the additional amplifying effects for negative values of both α and .

An Increase in the Time-Scale Separation (Decreasing Values of ) Causes an Amplification of the Voltage Amplitude Response
In Fig. 10 we compare the voltage responses for various representative values of and α = 1. As decreases, both Z max and Q Z increase, f res and f phas decrease, and the neuron becomes more selective (Λ 1/2 decreases) (Fig. 10b). These differences are captured by the shapes of the envelope curves, which are sharper the smaller the . For = 0.01 (panel a, left) the system is fast-slow causing trajectories to move fast along horizontal fibers as compared to larger values of (middle and right panels), and thus the limit cycle trajectory corresponding to f res reaches a higher value of v max (the voltage response is more amplified). The envelope-plane diagrams predict that f res and f phas are very close since the envelope curve is tangent to the dashed-red line almost at the peak (panel a).    = 0), respectively. The dashed-red lines represent the v-nullclines for A t = 1 (above-right) and A t = −1 (below-left). The solid-blue lines represent the trajectories of the forced system for a single period (T = 1000). b Impedance profile. c Phase profile. d Envelope-plane diagram. The solid-blue line represents the envelope curve: Each point on this curve is the maximum point on the limit cycle response to sinusoidal inputs parametrized by the input frequency f which increases from f = 0 (blue-square at the intersection between upper dashed-red and green curves) to f → ∞ (blue dot at the origin). (The v-coordinates of the envelope curve are the impedance function Z, since A in = 1.) Solid-red and -green lines: v-and w-nullclines for the unforced system (A t = 0). Dashed-red lines: v-nullclines for A t = 1 (above) and A t = −1 (below) coincides with the intersection between the green and dashed-red lines, and hence the system exhibits no phase-resonance (f phas = 0), while it still exhibits resonance.
In the limit of → 0, the system becomes quasi-one-dimensional and the dynamics occurs almost exclusively along the v-axis for almost all values of f , except a small range close to f = 0. For = 0 the dynamics is fully one-dimensional, thus generating a low-pass filter response (no resonance). For f → 0 (and = 0) in (29) and (30) the maximum value of v approaches the intersection between the dottedred line and the v-axis in the envelope-plane diagram. The maximum value of v (along the v-axis) decreases as f increases. As soon as > 0, the dynamics becomes two-dimensional. The envelope curves unfold in the two-dimensional space and ac- quire "triangular-like" shapes similar to these shown in Fig. 10a (left), but peakier the smaller . Regardless of how small is , there is always a range of values of f < such that as f → 0, the limit cycle trajectory evolves fast, "almost along" the wnullcline (green line) following the motion of the fixed point. Therefore, the point on the envelope curve for f = 0 lies on the intersection of the dotted-red line and the wnullcline. The remaining of the envelope curve will be as in Fig. 10a, but the distance between the envelope curve and the dotted-red line will be larger the smaller . This distance is determined by the minimum phase, which increases in absolute value as decreases.

Increasing Values of α Generate Resonance and Phase-Resonance and Cause an Amplification of the Voltage Response
In Fig. 11 we compare the voltage responses for various representative values of α and = 1. The value of α determines the slope of the w-nullcline, which is steeper for larger values of α, and v max (0) = (1 + α) −1 . The smaller slope of the w-nullcline for lower values of α (right panel) reduces the trajectories' freedom of motion by constraining them to evolve in close vicinities of both nullclines (the w-nullcline and the moving v-nullcline). The consequent decrease in speed prevents the limit cycle trajectory from reaching large enough values of v max . This together with the fact that v max (0) increases with α prevents the system from exhibiting resonance.
For larger values of α (left and middle panels), v max (0) is smaller. Although an increase in α increases the horizontal component of the vector field (through the product α ) and hence reduces the effective time-scale separation, this is not enough to prevent limit cycle trajectories from moving beyond v max (0), and the system is able to exhibit resonance.
The envelope-plane diagrams predict the generation of phase-resonance as α increases above some critical value (in between these for α = 2 and α = 0.5). For α = 0.2 and α = 0.5 (middle and right panels) the envelope curve is tangent to the dashed-red at its intersection with the w-nullcline, and hence the system exhibits no phase-resonance. In contrast, for α = 2 (left panel), the tangent point occurs in between this intersection and the peak of the envelope curve, and thus the system does exhibit phase-resonance.

Negative Values of Both α and Amplify the Voltage Response
From (12), negative values of both α and reflect the presence of a strong enough amplifying current that makes the effective conductance g L < 0 and a resonant current (g 1 > 0). For the underlying autonomous system to be stable the values of are constrained to be in the region > −1 (see Fig. 3a). We present a representative example in Fig. 12 (α = −2 and = −0.5). Comparison with the examples presented in Figs. 10 and 11 demonstrates that the voltage response is more amplified for negative than for positive values of both α and . This is true even for negative values of that are not very small in absolute value (time-scale separation no necessarily small).
Geometrically, negative values of α cause the slow of the w-nullcline to be negative (Fig. 12a)  There are three main differences between the two cases (positive and negative values of both α and ). First, the phase-resonant frequency f phas = 138 is larger than the resonant frequency f res = 108 for negative parameter values (Fig. 12b), while f phas < f res for positive parameter values. This difference is also captured by the envelope-plane diagram in Fig. 12d. Second, the voltage response is amplified by decreasing, rather than increasing levels of the time-scale separation (given by | |), for negative parameter values (Figs. 13a and 13b). Finally, the phase profiles are initially at φ = −π and increase towards φ = π/2 (Fig. 13c) for negative values of both α and (compare with Figs. 10 and 11), while φ never decreases below −π/2 for positive values of these two parameters.

Extension of the Envelope-Plane Diagram Approach to Simple Nonlinear Systems
Linearization around resting potential produces a good approximation to the impedance and phase profiles for weakly forced linear systems. However, signifi-cant departures from the linear prediction are expected for larger input amplitudes in the presence of strong nonlinearities, in particular when they are close to the voltage threshold for spike generation. The question arises of how nonlinearities affect the voltage response and, in particular, how the resonant and phase-resonant frequencies depend on the type of these nonlinearities and the time-scale separation between the participating variables. Linear responses to sinusoidal inputs are characterized by (i) the coincidence of the input and output frequencies, (ii) the proportionality between the output and the input signals that renders the impedance function independent of the input amplitude A in , and (iii) the symmetry between positive and negative values of the voltage response. In terms of the envelope-plane diagrams, the information about the voltage response is captured by the upper envelope curve for A in = 1. Here we illustrate how the ideas developed in the previous sections can be extended to the investigation of the mechanisms underlying the generation of resonance and phase-resonance in nonlinear systems.
We use simple piecewise-linear (PWL) systems of the form where the functions h v and h w are continuous PWL functions with two linear pieces. We consider two representative cases: (i) h v is PWL and h w is linear, and (ii) h v is linear and h w is PWL. In both cases, the nullclines intersect at the origin. By construction, in a vicinity of the fixed point (v,w) = (0, 0) both h v and h w are linear, and system (37)- (38) has the form (14)-(15).

The Voltage Response Is Amplified by Nonlinearities in the Voltage Equation
Here we consider system (37)- (38) with with η = −1, η r = −0.4 and α = 1. Our results are presented in Fig. 14. Figure 14a shows the envelope-plane diagrams (including both the upper and lower envelopes) for representative values of α and . The v-nullcline for the unperturbed system (solid-red lines) breaks at v = 0.8. For the forced system, the location of the vnullcline changes with t and so does the breaking point. The impedance and phase profiles for the parameters in Figs. 15a ( = 0.01) and 15b ( = 0.1) are presented in Figs. 15d1 and 15d2, respectively.
For small values of A in (Figs. 15a1 and 15b1) the dynamics of the PWL system is governed by the underlying linear system and the voltage response is linear. The envelope-plane diagrams are symmetric and similar to the ones discussed in previous sections. For larger values of A in the voltage response exhibits an amplification  (Figs. 15d1 and 15d2). As A in increases, the values of both Z max and Q Z increase, the resonant frequency f res decreases, and the selectivity increases. Notably, this nonlinear amplification of the voltage response is more pronounced the smaller the value of (larger levels of time-scale separation).
The nonlinear amplification of the voltage response can be understood by comparing the envelope-plane diagrams in Fig. 14a. Geometrically, the nonlinearities are reflected by a sudden change in the slope of the v-nullcline. Similarly to the linear case, as t progresses, the nonlinear v-nullcline moves cyclically between the two dashed-red curves corresponding to ±A in . The envelope curves live in the region bounded by the w-nullcline, the moving w-nullcline (dashed-red line) and the v-axis (horizontal w = 0 line). As A in increases, the area of this region increases as the triangle in panel a1 becomes a (non-triangular) polygon in Figs. 15a2 and 15a3. These changes are accompanied by an increase in the voltage amplitude of the regions that allow trajectories to reach higher values of v. Fig. 15 Voltage response of the PWL system (37)- (38) with h v and h w given by (40) to sinusoidal inputs. a A in = 0.8. b A in = 1. c A in = 1.5. We used the following parameters: The nonlinearities do not affect the dynamics of the limit cycle trajectory for low and large values of f , which behaves as explained above for linear systems, but they do affect the dynamics of the limit cycle trajectories for intermediate values of f . For A in = 0.8 (Fig. 15a1), the voltage response is quasi-linear since the nonlinearity is not present in the triangular region where the envelope curve is located. The limit cycle trajectories cannot move beyond the piece of the v-nullcline with slope η = −1.
For larger values of A in the nonlinearity moves above the horizontal axis (Figs. 15a2 and 15a3) during the ascending phase of the oscillatory input, thus changing the shape of the envelope-curve region (to a non-triangular polygon).

The Voltage Response Is Amplified by the Time-Scale Separation in the Presence of Nonlinearities in the Voltage Equation
The nonlinear amplification of the voltage response results from the ability of the limit cycle trajectories to move beyond the dashed-red line with slope η = −1.
Clearly, this phenomenon is more pronounced for larger values of A in . As expected from our previous discussion for linear systems, this phenomenon is also more pronounced for smaller values of (compare Figs. 14a and 14b) since limit cycle trajectories move in more horizontal directions for smaller values of , and they are able to reach larger voltage values. As increases, limit cycle trajectories move in less horizontal directions and, consequently, they do not take advantage of the extra available region of motion in the envelope-plane created by the nonlinear v-nullcline. We emphasize that the nonlinearities are present, and they are the same, for all values of , but the underlying vector field causes the limit cycle trajectories to "ignore" these nonlinearities and move inside the linear region (Fig. 14d).
The asymmetry in the system's nonlinearity is reflected in the asymmetry in the voltage responses, which is captured by the differences between the upper and lower envelope curves in the envelope-plane diagrams. Clearly, the difference between these two curves increases as A in increases.
We note that information about this asymmetric voltage response is not captured by the impedance function.

The Voltage Response Is Almost Insensitive to Nonlinearities in the Gating Variable Equation
Here we consider system (37)- (38) with with η = −1, α = 1 and α r = 0.4. Envelope-plane diagrams for representative values of are shown in Fig. 15. The nonlinearities are reflected in the sudden change in the slope of the w-nullcline. This change slightly affects the dynamics of the limit cycle trajectories for low values of f , since they evolve in close vicinities of the wnullcline but they do not significantly affect the dynamics of the limit cycle trajectory for larger values of f . Even for large values of A in (not shown), the nonlinearities in the w-nullcline do not cause any significant change in the envelope-curve region in the envelope-plane, and thus the impedance function does not significantly change with increasing values of A in . The resonant frequency does decrease with increasing values of A in but significantly less than for the case considered above (not shown).

Discussion
Subthreshold (membrane potential) resonance has been investigated in a number of neuron types both experimentally and theoretically [10, 12-22, 29, 40], and has been implicated in the generation of preferred neuronal firing rates [11] and network oscillations [23,25,27,[41][42][43]. Phase-resonance has also been investigated in neuronal systems [11,14,29], although to a lesser extent than subthreshold resonance, and it has been proposed to play an important role in neuronal synchronization [30]. The voltage response to oscillatory inputs in neurons and other electric circuits is typically characterized in terms of the impedance and phase profiles [10]. For linear systems receiving sinusoidal inputs, these curves can be computed analytically. For nonlinear systems or linear systems receiving non-sinusoidal oscillatory inputs, analytical calculations are not possible, and these curves are computed numerically.
Even when analytical expressions for the impedance and phase are available, the information they provide about the dynamic mechanisms leading to resonance is limited. For instance, the impedance and phase profiles fail to provide the necessary insight into the roles played by the interaction between the neuron's intrinsic time scales and the time scales associated with the input currents in the generation of resonance. Furthermore, linearized models fail to capture important nonlinear properties of the voltage response such as its nonlinear amplification for significant levels of the time-scale separation (small values of ).
In this paper we have developed and used dynamical system tools to investigate the dynamic mechanisms of generation of subthreshold and phase resonance in twodimensional linear and linearized (conductance-based) models, and we have extended these tools to include the effect of simple, but not necessarily weak, types of nonlinearities. The mechanistic analysis and the envelope-plane diagrams we developed in this paper provide a framework for the investigation of the preferred frequency responses in three-dimensional and nonlinear models. They can also be used as tools complementary to both the numerically computed and the experimentally measured impedance and phase profiles. More research is needed to adapt these tools to small neuronal networks.
Mechanistic studies of subthreshold resonance have mainly focused on the role that resonant and amplifying currents (and their associated gating variables) play in shaping the neuronal voltage response [10,11,29]. As we have shown in [29], the determination of the attributes and phase profiles as the result of the interaction between these two current types is more complex and not as straightforward as previously thought [10] (see our discussion of resonant and amplifying currents in Sect. 2.5). This is further emphasized by the fact that resonance may occur in the absence of intrinsic oscillations and vice versa [11,29] (see also Fig. 3).
In this study, we set out to understand the dynamic mechanisms underlying the generation of resonance and phase-resonance for a generic class of linearized biophysical models. We were particularly interested in issues that cannot be addressed solely by considering the impedance and phase profiles. These include (i) the identification of the mechanisms of amplification of the voltage response, (ii) the mechanisms of selection of the resonant and phase-resonant frequencies, (iii) the identification of the roles played by the intrinsic time scales and the time scales associated to the current inputs, (iv) the mechanisms that govern their interaction, (v) how all this is affected by changes in the model parameters, (vi) the relation between intrinsic STOs and subthreshold resonance, and (vii) the relation between subthreshold resonance and phase-resonance.
Perhaps our intuition on the resonant properties of neuronal models is based on the dynamics of the so-called λ-ω systems (55)-(56) (in Appendix A.3), which were used to build the resonate-and-fire models introduced in [31]. These systems display intrinsic oscillations with natural frequency Ω nat = ω for all values of λ (58). They also exhibit resonance and phase-resonance. The natural, resonant and phaseresonant frequencies (59) coincide for λ = 0 and they are different for other values of λ. Although λ-ω systems have been used to investigate the dynamics of resonant neurons [31], they are not representative of linearized neuronal models [11]. They rather correspond to cases where the voltage and gating variables evolve with comparable rates ( = 1) (see (61)-(62) in Appendix A.3), and leave out the large class of neuron types that have a strong time-scale separation ( is small) in the subthreshold regime.
The fact that resonance may occur in the absence of intrinsic STOs for nonnegligible parameter regimes [11,29] implies that resonance is not necessarily reflecting the amplification of an existing intrinsic oscillation by an oscillatory input current. Instead, resonance is uncovering the ability of the neuron to operate at time scales that are neither intrinsic nor imposed by the oscillatory input, but they emerge as the result of the interaction between a neuron's intrinsic time scale (determined by the neuron's intrinsic properties) and the time scales of the oscillatory input. These emergent time scales are likely to be the ones that play a significant role in network interactions.
To address these mechanistic issues from a geometric/dynamic perspective, we have extended the classical phase-plane analysis approach to include the effects of oscillatory inputs. This approach consists of projecting the three-dimensional space (for v, w and t) onto the two-dimensional plane (for v and w) and viewing the projection of the two-dimensional v-nullsurface (spanned by the v-nullcline for the autonomous system and t) as a moving one-dimensional v-nullcline as t progresses. Trajectories track the motion of the v-nullcline and the fixed point with a speed that depends on the underlying vector field and the input frequency f . The shapes of the limit cycle trajectories depend on f . A system exhibits resonance if for some value of f the amplitude of the limit cycle trajectory in the v-direction is larger than for f = 0. Qualitative predictions on the effect of changes in parameters on both resonance and phase-resonance can be made by looking at the phase-plane and the associated envelope-plane diagrams. Approaches including "moving nullclines" have been used before to investigate the mechanisms of synchronization of neuronal [44][45][46][47] and other systems [48,49] but they have not been used before in the context of the analysis of subthreshold and phase resonance.
From a geometric perspective, we view both resonance and phase-resonance as the result of the interaction between the time scales of the neuron (determined by the geometry of the unperturbed phase plane and the intrinsic time-scale separation) and the oscillatory input (captured by the moving v-nullcline). The resonant frequency, if it exists, is the input frequency at which the voltage responds optimally to the oscillatory driving current, which is captured by the cyclic movement of the v-nullcline. The trajectory's response is neither too fast nor too slow so that the voltage is able to reach a higher value than for other input frequencies. The time scale associated to the zero-phase frequency is also an emergent time scale, and is optimal in the sense that the trajectory is neither too fast nor too slow so that the trajectory intersects the moving v-nullcline at the exact time at which the latter reaches its maximum; i.e., both the voltage response and the input current peak at the same time. We showed that these two phenomena are captured by the shape of the envelope-plane diagrams not only for linear models, but also for nonlinear ones.
The concept of time scale for neural models is not always easy to quantify. In some limiting cases the time constants provide reliable information about the rate of change of the participating variables. This is the case for the one-dimensional passive membrane equation and for the so-called slow-fast systems [47]. However, time constants do not always capture the effective time scales. For the purpose of our study, we have qualitatively characterized the effective time scales of the isolated neurons by looking at their response to instantaneous DC (tonic) inputs and the shapes of the corresponding trajectories. We have identified the effective time scales of the isolated neurons that are relevant for their interaction with oscillatory inputs with the time scales that govern the behavior of the initial portion of the trajectory from its initial point until it reaches its highest voltage value. Geometrically, this is determined by the crossing point with the v-nullcline displaced by A in units (dotted-red v-nullclines). The effective time scales defined in this way may or may not represent the dynamic behavior of the autonomous trajectories for all subsequent times t, depending on the parameter values, but they are appropriate to describe the interaction between neurons and oscillatory inputs. As we showed, trajectories reverse direction after crossing the moving v-nullclines.
This notion of effective time scales for the underlying autonomous system is connected with its eigenvalues. However, our discussion of resonance is independent of the behavior of the autonomous system for large values of t. Therefore, the mechanisms of generation of resonance are in general independent from those responsible for the generation of intrinsic STOs.
The mechanism of generation of intrinsic STOs in linear systems is well understood [50]. The natural frequency (Ω nat ) corresponds to the eigenvalues' imaginary part. The eigenvalues' real part affects the amplitude of solutions, and consequently the evolution of trajectories in the phase plane, without affecting the oscillation frequency. The initial time interval of transient behavior that determines the effective time scales depends on both the eigenvalues' real part and Ω nat , and not only on Ω nat . Only when the effect of the eigenvalues' real part is small enough, both the generation of intrinsic STOs and resonance are governed by a similar underlying mechanism as occurs for the λ-ω systems discussed above.
To our knowledge, no previous study has addressed questions on subthreshold resonance in nonlinear systems from an analytical perspective. In this paper, we have extended our analysis for linear systems to include simple, piecewise-linear types of nonlinearities in both the equations for v and w. Our main question was: to what extent does the linear prediction capture the nonlinear effects? We found that when the v-nullcline is nonlinear the differences between the nonlinear response and the linear prediction increase not only with increasing values of the input amplitude A in but also with increasing levels of the time-scale separation between the voltage and the gating variable (decreasing values of ). These differences almost disappear when both equations evolve at comparable rates. In contrast, voltage responses are almost insensitive to nonlinearities located in the gating variable equation. In these two latter cases, the nonlinearities are present in the system but the voltage response does not detect them. More research is needed to understand whether these findings play out in more realistic nonlinear models.
Resonance has been first studied in the damped harmonic oscillator subject to sinusoidal forcing (see Appendix A.4). Similar to the λ-ω system, the natural and resonant frequencies coincide when the system is undamped and they are different in the damped case. Unlike the neural models discussed in this paper, the over-damped system (real eigenvalues) does not exhibit resonance. The damped harmonic oscillator can be rewritten as a system of two first order equations (64)-(65). However, differently from the systems considered in this paper, the sinusoidal forcing is located in the second equation rather than the first. In other words, the forcing term does not directly affect the dynamics of the variable for which we compute the response to the oscillatory input (analogous to v), but its derivative (analogous to w). Therefore, one should not necessarily expect resonance to occur in the absence of intrinsic oscillations. The geometric ideas developed in this paper can adapted to these systems. The moving nullcline will be the analog to the w-nullcline (oblique with negative slope) instead of the v-nullcline. The latter will be horizontal and fixed.

Competing Interests
The author declares that he has no competing interests.  (49) and The solution (48) can be written as The impedance amplitude and phase are given by (see Appendix A.5) and respectively. The impedance Z peaks at the resonant frequency provided the quantity inside the radical is positive.

A.4 The Harmonic Oscillator with Sinusoidal Forcing
The equation for the harmonic oscillator with sinusoidal forcing reads where β ≥ 0 and γ > 0. By defining y = −dx/dt, (63) can be rewritten as the following two-dimensional linear system of ODEs: dy dt = γ x − βy − C in sin(Ωt).
Differently from the systems considered in this paper, and the general form (41)- (42), the sinusoidal forcing is located in the second equation rather than the first. We note that the two equations are not interchangeable since we are following the convention that the first equation describes the dynamics of the variable (x) for which we compute the response to the oscillatory input. The eigenvalues of the autonomous system are given by The natural frequency of the autonomous system is given by provided 4γ − β 2 > 0. The impedance is given by Note that the formula (54) is not applicable in this case and the impedance has to be calculated separately. The impedance peaks at provided the quantity inside the radical is positive. For β = 0, Ω nat = Ω res = √ γ . For other values of β for which both Ω nat and Ω res are defined, Ω nat = Ω res . If the eigenvalues are real, the over-damped harmonic oscillator does not exhibit resonance since the inequalities 4γ − β 2 < 0 and 2γ − β 2 > 0 cannot be simultaneously satisfied.
Solving for A out and φ we obtain This model has been proposed in [37]. It has a persistent sodium current and a twocomponent (fast and slow) h-current given by This model has been proposed in [37]. It has a persistent sodium current and a slow potassium (M-type) current given by I p = G p p(V − E Na ) = G p p ∞ (V )(V − E Na ) and I Ks = G q q(V − E x k ) with E Na = 55 mV and E k = −90 mV. The voltagedependent activation/inactivation and time constants are given by