Membrane potential resonance in non-oscillatory neurons interacts with synaptic connectivity to produce network oscillations

Several neuron types have been shown to exhibit (subthreshold) membrane potential resonance (MPR), defined as the occurrence of a peak in their voltage amplitude response to oscillatory input currents at a preferred (resonant) frequency. MPR has been investigated both experimentally and theoretically. However, whether MPR is simply an epiphenomenon or it plays a functional role for the generation of neuronal network oscillations and how the latent time scales present in individual, non-oscillatory cells affect the properties of the oscillatory networks in which they are embedded are open questions. We address these issues by investigating a minimal network model consisting of (i) a non-oscillatory linear resonator (band-pass filter) with 2D dynamics, (ii) a passive cell (low-pass filter) with 1D linear dynamics, and (iii) nonlinear graded synaptic connections (excitatory or inhibitory) with instantaneous dynamics. We demonstrate that (i) the network oscillations crucially depend on the presence of MPR in the resonator, (ii) they are amplified by the network connectivity, (iii) they develop relaxation oscillations for high enough levels of mutual inhibition/excitation, and (iv) the network frequency monotonically depends on the resonators resonant frequency. We explain these phenomena using a reduced adapted version of the classical phase-plane analysis that helps uncovering the type of effective network nonlinearities that contribute to the generation of network oscillations. We extend our results to networks having cells with 2D dynamics. Our results have direct implications for network models of firing rate type and other biological oscillatory networks (e.g, biochemical, genetic).

connectivity and involve the interplay of the nonlinearities and time scales present in the ionic and synaptic currents. In some cases, the network time scales directly reflect the time scales of the individual neurons. This class includes the synchronized activity of population of oscillators where the frequency band of both the network and the individual oscillators coincides.
There are other cases where the oscillatory time scales are latent (or hidden) at the individual neuron level and become apparent only at the network level. This class includes the oscillatory networks of non-oscillatory neurons that are the focus of this paper. More specifically, we investigate oscillatory networks where at least one of the participating (non-oscillatory) cells exhibits (subthreshold) membrane potential resonance (MPR), defined as the occurrence of a peak in the cell's voltage amplitude response to oscillatory input currents at a preferred (resonant) frequency (Hutcheon and Yarom 2000;Richardson et al. 2003;Nadim 2014b, 2015). Because the individual cells are intrinsically non-oscillatory, the resonant frequency reflects an oscillatory latent time scale that can be uncovered in the presence of oscillatory input currents, but not by direct observation of their spontaneous behavior.
The mechanisms of generation of sustained (limit cycle) oscillations in single neurons are reasonably well understood (Ermentrout and Terman 2010;Borgers 2017;Izhikevich 2006;Dayan and Abbott 2001). They require the interplay of negative and positive feedback effects mediated by the ionic current gating variables or related processes. Resonant ionic processes (e.g., hyperpolarization-activated mixed-cation I h current, M-type slow-potassium current I Ks and T-type calcium inactivation I CaT ) oppose changes in voltage, while amplifying ionic processes (e.g., persistent sodium current I Nap , T-type calcium activation) favor these changes.
From the oscillatory dynamics point of view, there is a hierarchy of phenomena that requires the presence of a resonant process and whose degree complexity increases with the levels of the amplifying current (Hutcheon and Yarom 2000;Rotstein 2017b) in system where sustained oscillations (subthreshold or spikes) are generated by Hopf bifurcation mechanisms (Ermentrout and Terman 2010;Borgers 2017;Izhikevich 2006). At the bottom of this hierarchy are the overshoot type of responses to square-pulse perturbations ( Fig. 1, green curves) in neurons that exhibit MPR (Hutcheon and Yarom 2000;Richardson et al. 2003;Nadim 2014b, 2015), but not subthreshold oscillations (STOs). We refer to them as resonators. For higher amplification levels the neuron may display damped subthreshold oscillations (Fig. 1, red curves). In these two cases the underlying systems may be quasi-linear in large enough vicinities of the resting potential (fixed-point) (Rotstein 2017b). (Damped oscillators may also exhibit resonance, but we do not refer to them as resonators.) At the top of the hierarchy are the sustained (limit cycle) oscillations ( Fig. 1, blue curves) that require high enough amplification levels for the development of the nonlinearities necessary for the existence of limit cycles (Rotstein 2017b). If these limit cycles represent STOs, additional amplification levels can produce spikes or depolarization block. Examples of models exhibiting this type of behavior are the Morris-Lecar model (Morris and Lecar 1981) and the I h + I Nap or I Ks + I Nap models studied in Rotstein (2017b) (see also Rotstein 2017c).
MPR has been investigated in many neuron types both experimentally and theoretically (Hutcheon and Yarom 2000;Richardson et al. 2003;Lampl and Yarom 1997;Llinás and Yarom 1986;Erchova et al. 2004;Schreiber et al. 2004;Hutcheon et al. 1996;Gastrein et al. 2011;Hu et al. 2002Hu et al. , 2009Johnston 2007, 2008;Marcelin et al. 2009;D'angelo et al. 2001Pike et al. 2000;Tseng and Nadim 2010;Tohidi and Nadim 2009;Solinas et al. 2007;Wu et al. 2001;Muresan and Savin 2007;Heys et al. 2010Heys et al. , 2012Zemankovics et al. 2010;Nolan et al. 2007;Engel et al. 2008;Boehlen et al. 2010Boehlen et al. , 2013Narayanan 2012, 2014;Fox et al. 2017;Chen et al. 2016;Beatty et al. 2015;Song et al. 2016;Art et al. 1986;Remme et al. 2014;Higgs and Spain 2009;Yang et al. 2009;Mikiel-Hunter et al. 2016;Rau et al. 2015;Sciamanna and Wilson 2011;Lau and Zochowski 2011;Rotstein 2014a, b, 2015 Response of I h +I Nap and I Ks +I Nap models to negative square pulses of current: representative dynamic scenarios. a I h +I Nap model. It includes three ionic currents: hyperpolarization-activated (h-), persistent sodium and leak (see Section 2.4 in Methods). b I Ks +I Nap model. It includes three ionic currents: M-type slow potassium, persistent sodium and leak (see Section 2.4 Methods). Both I h and I Ks are resonant and I Nap is amplifying. Increasing the levels of I Nap causes a transition from overshoot responses (green) to damped oscillations (red) to persistent (limit cycle) oscillations (blue) in both models. The gray curve is a caricature of the square wave input deflected from zero with amplitude 1. We used the following parameter values: C = 1, E Na = 42, E L = −75, E h = −26, G L = 0.3, G h = 1.5, I app = 0.55, V hlf,p = −54.7, V slp,p = 4.4, V hlf,q = −80.2 and V slp,q = 7.2 (I h +I Nap model) and C = 1, E Na = 42, E L = −75, E Ks = −96, G L = 0.3, G Ks = 1.5 and I app = 4, V hlf,p = −54.7, V slp,p = 4.4, V hlf,q = −28, V slp,q = 8 (I Ks +I Nap model) 2017a; Szucs et al. 2017). However, in contrast to single cell intrinsic oscillations, the consequences of cellular MPR on network oscillations are not well understood. Only a few studies have addressed these issues in networks having neurons that exhibit MPR Stark et al. 2013;Tikidji-Hamburyan et al. 2015;Tchumatchenko and Clopath 2014;Schmidt et al. 2016;Moca et al. 2014; Baroni et al. 2014;Rotstein et al. 2017d) or have resonant gating variables (Wang and Rinzel 1992;Manor et al. 1997;Torben-Nielsen et al. 2012). To our knowledge, no study to date has examined the detailed mechanisms of generation of oscillations in networks of non-oscillatory resonators and how the network oscillations reflect the latent time scale provided by the resonant frequency.
From the mechanistic point of view, we seek to understand how the resonant properties of individual nodes interact with the network connectivity to produce oscillations in reciprocally connected networks. We reasoned that if oscillations are to be generated in networks where the participating neurons only provide the resonant properties, then the amplification effects should result from the network connectivity. According to this hypothesis, oscillations should be generated in self-excited (Fig. 2a), but not selfinhibited (Fig. 2b) resonators and in two-cell networks of mutually inhibited or mutually excited cells that include one resonator (Figs. 2d, e and f), but not in mutually inhibited non-resonant cells (low-pass filters) (Fig. 2c). Moreover, the resonant frequency of the individual resonators should control, or at least have a direct effect, on the network frequency. Analogously to single cell oscillations, the mechanism of generation of network oscillations should involve a Hopf bifurcation and the dynamic hierarchy described above. Some of these patterns have been observed for similar systems (Manor et al. 1999) and for network models using the Wilson-Cowan formalism (Beer 1995;Ermentrout and Terman 2010;Wilson and Cowan 1972). However, the role that the filtering properties of the individual nodes (preferred frequency responses to oscillatory inputs) play in the generation of network oscillations has not been investigated.
We test these ideas using the simplest types of oscillatory networks of non-oscillatory neurons, consisting of a linear resonator reciprocally connected to a linear cell (either a low-pass filter showing no resonance or another resonator) with instantaneous graded synapses (see the motif diagrams in Fig. 2). We use linear (linearized conductancebased) models for the individual neurons to isolate the resonant (negative feedback) effects from the nonlinear amplifications that may lead to sustained oscillations.
These linearized models capture the quasi-linear dynamics of models having the passive currents and I h or I Ks , but no amplifying currents (e.g., I Nap ) (Rotstein 2017b). They also capture the dynamics of uncoupled components of the firing rate models of Wilson-Cowan type (Wilson and Cowan 1972) with adaptation (Curtu and Rubin 2011;Shpiro et al. 2009;Tabak et al. 2011).
Finally, we use graded synapses because of the subthreshold range of voltages in which they operate and because it is the type of nonlinearities used in firing rate models. They are assumed to be instantaneously fast and to have no dynamics (Wang and Rinzel 1992;Manor et al. 1997Manor et al. , 1999Ambrosio-Mouser et al. 2006;Brea et al. 2009;David et al. 2015;Curtu and Rubin 2011;Shpiro et al. 2009;Tabak et al. 2011 ) to strip them from any additional dynamic effects. The questions we ask in this paper aim to conceptually address the mechanisms by which neuronal frequency filters interact within a network. Our results in conjunction with the results of previous, complementary studies (Manor et al. 1997(Manor et al. , 1999Wang and Rinzel 1992;Ambrosio-Mouser et al. 2006;Chen et al. 2016) have implications not only for the understanding of neuronal oscillations, but also for the understanding how frequency-dependent information is communicated across neurons and networks and the phenomenon of network resonance (Ledoux and Brunel 2011;Stark et al. 2013).
The overview of the paper is as follows. In Section 3.1 we review the frequency response properties of individual neurons with one-and two-dimensional linear dynamics. In the subsequent sections we combine resonators (2D bandpass filters) and passive cells (1D low-pass filters) to analyze the circumstances under which network oscillations are generated for the various circuit motifs shown in Fig. 2 and their properties. In Section 3.2 we discuss the oscillatory properties of self-excited resonators (Fig. 2b) and show that the limit cycle oscillations monotonically depend on the resonator's resonant frequency. In Section 3.3 we show that mutually inhibitory networks consisting of one resonator and a passive cell (Fig. 2e) are able to generate limit cycle oscillations, discuss their properties and show that their frequency monotonically depends on the resonator's resonant frequency. In Sections 3.4 and 3.5 we extend these results to mutually excitatory networks consisting of one resonator and a passive cell (Fig. 2d) and a mutually inhibitory network consisting of two resonators (Fig. 2f). The "negative results" are discussed in the Appendix. There we show that self-inhibited resonators (Fig. 2a) and two-cell networks of passive cells (e.g., Fig. 2c), that do not include resonators, do not show sustained (limit cycle) oscillations. In addition to network oscillations, the two-cell networks we consider in this paper have non-oscillatory regimes (stable fixed-points) that may show linear and nonlinear resonance in response to oscillatory inputs, and could be functional in the generation of oscillations in larger networks. The investigation of these more general scenarios is outside the scope of this paper. Finally, in Section 4 we discuss our results and their implications for network dynamics.

Networks of linearized cells with graded synapses
We used linearized biophysical (conductance-based) models for the individual cells and (nonlinear) graded synaptic connections. The linearization process for conductancebased models for single cells has been previously described in Richardson et al. (2003) and Rotstein and Nadim (2014b). We refer the reader to these references for details.
The dynamics of a two-cell network are described by for k = 1, 2. In Equations (1)-(2) t is time, v k is the voltage (mV) relative to the voltage coordinate of the fixedpoint (equilibrium potential)V k , w k is the gating variable relative to the gating variable coordinate of the fixed-point w k and normalized by the derivative of the corresponding activation curve, C k is the capacitance, g L,k is the linearized leak maximal conductance, g k is the ionic current linearized conductance, τ k is the linearized time constant and I syn,k is the graded synaptic current from the other neuron in the network and given by where G syn,k is the maximal synaptic conductance, E syn,k is the synaptic reversal potential relative toV k and where the half-activation point v hlf is also relative toV k . We use the following units: mV for voltage, ms for time, μF/cm 2 for capacitance, μA/cm 2 for current and mS/cm 2 for the maximal conductances. Unless stated otherwise, we used the following parameter values: C = 1, V hlf = 0, V slp = 1, E in = −20, E ex = 60.
Note that the heterogeneity due to different values of the DC current I app,k and other biophysical parameters in the original conductance-based model is translated into the reversal potentials E syn,k and the functions S k,∞ (v) through the fixed-point (V 1 ,V 2 ). Specifically, if E syn and V hlf are the synaptic reversal potential and synaptic half-activation point of the original (not rescaled) model, then E syn,k = E syn −V k and v hlf = V hlf −V k .
(3). The phase-plane diagram is 2D. Because there is only one cell involved, we omit the subscript in the notation of the participating variables and parameters. The vand w-nullclines are given, respectively, by and

Graded network of two 2D cells: hyper-nullclines, fixed-points and dynamic phase-plane analysis
These networks are given by system (1)-(4). The phasespace diagram is 4D. The v 1 -and v 2 -nullsurfaces (obtained by making the current-balance equation for the corresponding nodes equal to zero) depend on different variables (the v 1 -nullsurface depends on w 1 and v 2 and the v 2 -nullsurface depends on v 1 and w 2 ). The w 1 -and w 2 -nullsurfaces are planes given by w 1 = v 1 and w 2 = v 2 , respectively. By substituting into the corresponding current-balance equations and rearranging terms we obtain the following equations describing curves in the v 1 -v 2 plane and These are extensions of the nullclines (5) and (6) for the networks of 1D passive cells. Their intersection (v 1 ,v 2 ) give the v 1 -and v 2 -coordinates of the 4D fixed-points (v 1 ,v 2 ,w 1 ,w 2 ) = (v 1 ,v 2 ,v 1 ,v 2 ). However, they are not nullclines, but projections of hyper-nullsurfaces onto the v 1 -v 2 plane. We refer to them as hyper-nullclines. For the hybrid networks having one 2D and one 1D cell we set g 2 = 0 in Eq. (11).

Bifurcation diagrams
As we mentioned in the previous section, the fixed-points are the intersections between the nullclines (for 2D systems) or the hyper-nullclines (for 3D and 4D systems). To determine the stability of the fixed-points we calculate the eigenvalues of the corresponding linearized system. For the 2D system of two 1D passive cells, the eigenvalues are easily calculated (see Appendix A). The expressions of the eigenvalues for the other considered networks (3D or 4D) are much more extensive and we will not show them in this work.
In all systems we can study the eigenvalue expressions when the parameter values vary, and we determine the existence of static bifurcations (such as, pitchfork and saddlenode) and dynamic bifurcations (for example, Hopf bifurcation) (Guckenheimer and Holmes 1983). If a Hopf bifurcation exists, we calculate the first Lyapunov coefficient with the MATLAB package MatCont (Dhooge et al. 2003), to determine the direction and stability of the emerging branch of cycles.
Considering the bifurcations of the fixed-points, we construct bifurcation diagrams in several parameter spaces determining regions with different dynamical scenarios. In particular, we can determine parameter values in which stable limit cycles exist.

Conductance-based models
Primarily for illustrative purposes, in some of our simulations we used biophysical (conductance-based) models (Skinner 2006;Hodgkin and Huxley 1952) to describe the subthreshold dynamics of neurons having one resonant and one fast amplifying currents. The current balance equation is given by where V is the membrane potential (mV), t is time (ms), 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 ), I L = G L (V − E L ) is the leak current, and I j = G j x j (V − E j ) are generic expressions for ionic currents (j = 1, 2) with maximal conductance G j (mS/cm 2 ) and reversal potentials E j (mV) respectively. The gating variables obey kinetic equations of the form where x j,∞ (V ) and τ j,x (V ) are the voltage-dependent activation/inactivation curves and time constants respectively. The former are given by where V hlf,x and V slp,x > 0 are constants and the sign of σ x indicates whether the curve describes an activation (σ x < 0) or inactivation (σ x > 0) process. In this paper we use voltage-independent time constants τ j,x . This assumption is mostly for simplicity since we are focusing on the subthreshold voltage regime where the time constants are typically slowly varying functions of V . The ionic currents I j we consider here are persistent sodium,

Numerical simulations
The numerical solutions were computed by using the modified Euler method (Runge-Kutta, order 2) (Burden and Faires 1980) with a time step t = 0.1 ms in MATLAB (The Mathworks, Natick, MA). Smaller values of t have been used to check the accuracy of the results.

Frequency preference response of individual cells revisited: resonators and passive cells / low-pass and band-pass filters
In this paper we consider individual cells with 1D and 2D linear dynamics having a stable fixed-point. The dynamics are described by Eqs. (1)-(2) with I syn,k = 0 (for cells with 1D dynamics g k = 0). The response of linear cells to oscillatory inputs at different frequencies f (e.g., sinusoidal functions of t) is captured by the impedance Z(f ), which is complex function with amplitude and phase. Following other authors we use the term impedance (and we write Z(f )) to refer to the impedance amplitude |Z(f )| and we refer to the corresponding curve as the impedance profile (Fig. 3). In this paper we focus on the impedance amplitude.
Linear 1D (passive) cells are low-pass filters (Fig. 3, red curve), while linear 2D cells can be either low-pass filters (not shown) or band-pass filters (Fig. 3, blue curve). Resonance refers to the ability of a cell to exhibit a peak in their amplitude response (impedance profile) at a preferred (resonant) frequency f res (Hutcheon and Yarom 2000;Richardson et al. 2003;Rotstein and Nadim 2014b). The corresponding unforced cells can have either a node, and exhibit an overshoot (as in Fig. 1, green curves), or a focus, band-pass filter -resonance low-pass filter -no resonance Fig. 3 Response of the individual cells to oscillatory inputs for representative parameter values. Impedance profile for a resonator (band-pass filter, blue) and a passive-cell (low-pass filter, red). The parameters Z max and f res are the maximal impedance and the resonant frequency, respectively and display damped oscillations (as in Fig. 1, red curves). We use the term resonator to refer to cells that exhibit resonance, but not damped oscillations. The resonant properties of 2D linear systems, including their relationship between the intrinsic properties of the unforced cells (e.g, eigenvalues, intrinsic oscillatory frequencies) and the dynamic mechanism of generation of resonance have been investigated extensively by us and other authors (Richardson et al. 2003;Rotstein 2014a, b). We refer the readers to these references for details.

Self-excited resonators can produce limit cycle oscillations and their frequency monotonically depends on the resonator's resonant frequency
The self-excited resonator model is given by system (1) Because there is only one cell involved, we omit the subscripts in the notation of the variables and parameters. The nullclines of the phaseplane diagrams are given by Eqs. (8) and (9). The individual resonator does not oscillate.
Self-excitation is the simplest mechanism of network oscillation amplification of a resonator. Mathematically, a self-excited resonator has the same structure as individual resonator+amplifying current models (e.g., I h + I Nap or I Ks +I Nap ), which are able to produce sustained oscillations for large enough amplification levels ( Fig. 1) (Rotstein 2017b). In both types of models the activation of the amplifying component (I syn and I Nap ) is instantaneous (or very fast), the shapes of their activation curves are similar, and the reversal potentials (E Na and E ex ) are above the resting potential. Models having I h or I Ks as the only active ionic currents are quasi-linear resonators (Rotstein 2017b). Therefore, it is not surprising that self-excited linear resonators are able to produce oscillations given that resonant+amplifying models can do so. However, since resonance and amplification belong to different levels of organization in self-excited resonators, we can dissociate these two effects and investigate the effects of the resonant  Self-excited (linear) resonators can produce sustained (limit cycle) oscillations, while self-excited damped-oscillators may fail to produce sustained oscillations. Phase-plane diagrams for representative parameter values. The vand w-nullclines are given by (8) and (9), respectively. The fixed-point for the uncoupled (linear) system is a stable node (f nat = 0) in panels a and b and a stable focus in panels c (f nat ∼ 48.9). a f res ∼ 17.6 for g L = 0.25, g = 1 and τ = 100. b f res ∼ 10.4 for g L = 0.25, g = 0.25 and τ = 100. c f res ∼ 55.2 for g L = 0.25, g = 1 and τ = 10. We used the following additional parameter values: , E ex = 60, v hlf = 0, v slp = 1 frequency of the individual neurons on the oscillation frequency, which we cannot do in individual cells.

Self-excited resonators can produce sustained (limit cycle) oscillations for appropriate balances among the resonance, amplification and excitation levels
Geometrically, increasing values of the excitatory maximal conductance G ex create nonlinearities of cubic type in the phase-plane diagram (Fig. 4). In single neurons this type of nonlinearities are typically created by amplifying gating variables (e.g., I Nap ) in the presence of resonant gating variables (e.g., I h or I Ks ) (Rotstein 2017b; c) (see also Ermentrout and Terman 2010;Izhikevich 2006). Fig. 4-a illustrates the effects of increasing values of G ex when the linearized resonant conductance g (= 1) is much larger than the linearized leak conductance g L (= 0.25) for a resonator (with no intrinsic damped oscillations when G ex = 0). For low values of G ex , the coupled cell shows damped oscillations as the cubic-like nonlinearities of the v-nullcline begin to develop (panel a1). Limit cycle oscillations emerge as G ex increases further (panel a2) and disappear when the fixed-point moves to the right branch of the cubic-like v-nullcline for larger values of G ex and regains stability (panel a3). As G ex increases within the oscillatory range the amplitude increases and a time scale separation between the participating becomes more prominent, generating, for large enough values of G ex , oscillations of relaxation type (Fig. 5a). Figure 4-b illustrates that oscillations are not generated when the g (= 0.25) to g L (= 0.25) ratios are relatively low. The cubic-like nonlinearities are still developed for high enough values of G ex (panels b2 and b3), but the fixed-point is located on the right branch of the v-nullcline where the fixed-point is stable, and moves further away from the knee as G ex increases. The amplification still happens, but it leads directly to depolarization block without oscillations. Similar behavior was observed when the fixedpoint of the isolated cell is a stable focus instead of a stable node. However, oscillations can be restored by increasing the value of v hlf , which moves the fixed-point to the middle branch where it loses stability (not shown).
The transition of a resonator to a damped oscillator can be achieved by decreasing the value of τ (Rotstein and Nadim 2014b). Contrary to intuition, the presence of damped oscillations in the cell does not necessarily generate sustained oscillations in the self-excited network (Fig. 4c). When it happens, the time scale separation is smaller than for the resonator and therefore relaxation oscillations are more difficult to obtain (Fig. 5b).

The intrinsic resonant frequency controls the network oscillations frequency
Self-excited resonators are the simplest models where we can investigate the effects that changes on the resonant frequency (f res ) of the individual non-oscillatory cells have on the network oscillation frequency (f ntw ). The resonator parameters that control f res (28) also control the values of other attributes of the impedance profile Z(f ) such as the maximal impedance Z max (29). In order to establish the effects of f res on f ntw it is necessary change the model parameters in such a way as to cause the minimal possible changes on the shape of Z(f ) ). In the ideal situation, changes in f res would be accompanied only by a translation of Z(f ). This is not possible for 2D linear models, but it is possible to change the model parameters in a balanced way so that f res changes, but Z max remains constant . In this way the impedance profiles are displaced with minimal changes in their shape (Fig. 6a). Oscillations in self-excited resonators: the intrinsic resonant frequency controls the network frequency. a Representative resonator impedance profiles with different resonant frequencies (f res ) and the same maximal impedance: Z max ∼ 9.5 (a1) and Z max = 3.9 (a2). b Network oscillation frequency (b1) and amplitude (b2) as a function of f res for representative values of g L . c Network oscillation frequency (c1) and amplitude (c2) as a function of f res for representative values of G ex . We used the following parameter values: E ex = 60, v hlf = 0, v slp = 1 (Fig. 6b1) and G ex (Fig. 6c1). As expected, the oscillations are more amplified the lower g L (Fig. 6b2) Nadim 2014b, 2017b) and the higher G ex (Fig. 6c2).

Sustained (limit cycle) oscillations are lost in self-excited 2D cells as the transition from resonators to low-pass filter
Resonance can be lost by various mechanisms (Rotstein and Nadim 2014b). One of them is having low enough values of the resonant conductance g (in the limit of g = 0 the coupled cell is 1D and therefore oscillations are not possible). Another one is having low-enough values of the time constant τ . Figure 7a illustrates how oscillations are lost as τ decreases. Note that the location of the fixedpoint is independent of τ . Figure 7b anc c illustrate that oscillations cannot be recovered by decreasing G ex (panel b) or increasing v hlf (panel c) for the same value of τ as in panel a3. In both cases, these changes move the fixed-point to the middle branch of the v-nullcline, but it remains stable.

Mutually inhibitory 2D/1D hybrid networks can generate sustained (limit cycle) oscillations and their frequency monotonically depends on the intrinsic resonant frequency
The hybrid networks we consider here consist of a linear resonator (2D, cell 1) and a passive cell (1D, cell 2) reciprocally inhibited through graded synapses. We use system (1)-(2) with g 1 > 0 and g 2 = 0 and the additional description of the synaptic connectivity presented in Section 2.
These networks can be thought of as two "overlapping" circuits, non of each able to produce oscillations on their own: the linear 2D resonator used in the Appendix B and the reciprocally inhibited passive cells discussed in Appendix A. The oscillations result from the combined Fig. 7 Self-excited 2D cells: phase-plane diagrams for representative parameter values. The vand w-nullclines are given by (8) and (9), respectively. The quantities f nat and f res refer to the natural and resonant frequencies of the uncoupled cells. The fixed-point for the uncoupled system is a stable focus. a G ex = 0.04 and v hlf = 0. b G ex = 0.015 and v hlf = 0. c G ex = 0.04 and v hlf = 1. We used the following parameter values: g L = 0.25, g = 1, E ex = 60, v slp = 1 activity of these two "sub-circuits" where the mutuallyinhibitory component acts as an amplifier of the resonant component.
For our analysis we represent the dynamics of these 3D networks using projections of the 3D phase-space (for v 1 , w 1 and v 2 ) onto the v 1 − v 2 plane and use the hypernullclines (10)-(11) (g 2 = 0) defined in Section 2.2 (e.g., Fig. 8, left columns). In order to relate the dynamics of the hybrid networks to these of the mutually inhibitory passive cells we include in the phase-plane diagrams the v-nullcline for cell 1 (dashed-red curve) for g 1 = 0 (no resonant gating variable). Figure 8 shows the oscillations generated in these networks for values of G in (= G in,1,2 = G in,2,1 ) that increase from top to bottom. Because the networks are mutually inhibited, these oscillations are not synchronized in-phase. They are created in a supercritical Hopf bifurcation (Fig. 11a1) and therefore they have small amplitude and are sinusoidal-like for small enough values of G in (Fig. 8a). The effect of the resonant gating variable w 1 is to bring the fixed-point of the mutually inhibitory (non-oscillatory) 1D/1D system (intersection between the dashed-red and green curves) to the oscillatory region where the v 2 hyper-nullcline (green curve) is non-linear.

Oscillations can be generated in 2D/1D hybrid networks and are amplified by increasing levels of mutual inhibition
The oscillations amplitude increases with increasing values of G in as the limit cycle trajectories evolve in small vicinities of the v 2 hyper-nullcline ( Fig. 8b and c). This amplification is accompanied by the development of a separation of time scales. For large enough values of G in the oscillations are of relaxation-type (Fig. 8c). This partially reflects the time constant of the resonators, which needs to be slow enough, but it is a network effect since linear models do not display sustained oscillations.
For low values of G in within the oscillatory regime, the network has only one fixed-point. As G in increases, additional fixed-points are created (Figs. 8c1 and 11c1) in a Pitchfork bifurcation (Fig. 11a1), but they are not stable and they do not obstruct the presence of oscillations. However, as G in increase further, these new fixed-points become stable by subcritical Hopf bifurcations and coexist as attractors with the limit cycle (Fig. 11a1). The oscillations are abruptly terminated when the stable limit cycle collides with an unstable limit cycle generated in one of the mentioned subcritical Hopf bifurcations (Fig. 11c1). Without oscillations the attractors that remain in the network are the fixed-points corresponding to one of the cells inhibited ( Fig. 11a1 and c1).
Similarly to the mutually inhibited passive cells discussed above, for other parameter regimes the pitchfork bifurcation can be transformed into a saddle-node bifurcation (Fig. 11a2) without causing significant qualitative changes to the network dynamics (Fig. 11c2). Figure 11b illustrates that existence of network oscillations requires balanced combinations of G in and g 1 .
The generation of oscillations requires certain heterogeneity in the underlying mutually inhibitory 1D/1D system. For the oscillations in Fig. 8, g L,1 = 0.25, g L,2 = 0.5. This has been observed also for the related system studied in Manor et al. (1999). Oscillations are not possible for the hybrid 2D/1D network when g L,1 = g L,2 and C 1 = C 2 unless there is heterogeneity in the synaptic connectivity (G in,2,1 > G in,1,2 ) (not shown).

Development of relaxation oscillations for large mutual inhibition levels
In order to understand the mechanisms of generation of network oscillations and their properties in terms of the model parameters it is useful to consider the v 1 -nullsurface parametrized by constant values of w 1 , N 1c (v 1 ) = N 1c (v 1 , c), and track the motion of the trajectory as time progresses and the values of w 1 change. This will cause the v 1 -hyper-nullclines in Fig. 8 to move as the trajectory evolves following the dynamics of w 1 . For the second cell the curve (10) for g 2 = 0 is time-independent and therefore remains fixed. Note that Eq. (15) is Eq. (10) before w 1 is substituted by v 1 . In order to uncover the presence of nonlinearities of cubic type and to further capture the effect of the model's geometric properties that give rise to the different types of oscillations we use an adapted version (Fig. 9) of the phaseplane diagram discussed in Fig. 8 where the hyper-nullclines and trajectories are plotted relative to the v 2 -hyper-nullcline. In the adapted phase-plane diagram, the v 2 hyper-nullcline (green curve) is the zero-level line and the v 1 -hypernullclines (red curves) are cubic-like. The dashed-red curves represent the maximum (lower curve) and minimum (upper curve) levels of w 1 during the oscillations. The red curve moves in between these two dashed-red curves as the oscillation progresses. The intersections between the green line and the moving red curve generated "transient fixedpoints", which are not fixed-points of the 3D system, but they serve as targets for the evolution of the trajectories. The v 1 -and v 2 -hyper-nullclines are given by Eqs. (10) and (11), respectively. Black dots indicate stable nodes and gray dots indicate unstable foci. The dashed red curve represents the v 1 nullcline for cell 1 for g 1 = 0 (no resonant gating variable). v 1 and v 2 as a function of t). a G in,1,2 = G in,2,1 = 0.112. The network frequency is f ntw = 6.1. b G in,1,2 = G in,2,1 = 0.14. The network frequency is f ntw = 5.4. c G in,1,2 = G in,2,1 = 0.22. The network frequency is f ntw = 2.4. We used the following parameter values: Specifically, in the v 1 -v 2 plane presented in Fig. 9 (left panels), trajectories move towards the transient fixed-points with negative slope and their speed depends on the distance between the moving red curve and the green line. The existence of oscillations imply that the local extrema of the dashed-red curves do not intersect the zero-level reference green line. We emphasize that this is not a standard phaseplane diagram and it captures only specific aspects of the dynamics.

Right. Voltage traces (curves of
Similarly the self-excited resonator discussed above, increasing amplification levels are characterized by more pronounced cubic-like nonlinearities. Here the amplification levels are provided by the levels of mutual inhibition that are measured in terms of the values of G in (compare panels a1 and b1 in Fig. 9).
When w 1 = w 1,min the red curve is at its highest level and the trajectory moves to the right (F 1 ), towards the only transient fixed-point with relatively high speed (jump up).
As this happens w 1 increases, causing the red curve to shift down with the consequent motion of the transient fixedpoint to the left. The variable v 1 reaches its maximum when the trajectory crosses the transient fixed point and is forced to reverse direction (S 1 ). As the red curve continues to shifts down, the stable and unstable transient fixed-points collide and disappear leaving only one transient fixed-point (on the leftmost side, for lower values of v 1 ), which becomes the new target for the trajectory. The trajectory moves towards this target fixed-point, but it does so on a very slow time scale (S 1 ) due to the ghost effect of the "defunct" fixedpoints until it reaches the (jump down) region of fast motion (F 2 ). The process repeats to complete the cycle through.
Relaxation oscillations are created when difference between the two local extrema on each dashed-red curve is large enough (well separated). This occus in Fig. 9a, but not in Fig. 9b, where the oscillations do not show any separation of time scales.

The resonator's intrinsic resonant frequency controls the network oscillations frequency
Similarly to the self-excited resonator networks discussed above, the functional role of cellular resonance is to determine the frequency of the network oscillations. This is illustrated in Fig. 10 for various representative parameter set values. We followed the same protocol as in Section 3.2.2 ( Fig. 6): for each value of f res , the values of g 1 and τ 1 are balanced so to maintain Z max constant. In all cases, f ntw increases with increasing values of f res (left panels). The oscillation amplitude increases with increasing values of G in (= G in,1,2 = G in,2,1 ) and is more variable than for the self-excited resonator (Fig. 11).
The oscillatory active f res band (the range of values of f res for which network oscillations are possible) is relatively small as compared to the self-excited resonator network and it depends on the value of Z max and g L,1 . All other parameters fixed, decreasing values of Z max (from Fig. 10a-b) causes the oscillatory active resonant frequency band to slide to the right. Figure 10c shows that the size active frequency band can be increased by decreasing g L,1 and increasing Z max . A proper comparison would involve changing one parameter at the time, but decreasing values of g L,1 require increasing values of Z max for the oscillations to be present.

Mutually excitatory 2D/1D hybrid networks can generate limit cycle oscillations and their frequency monotonically depends on the intrinsic resonant frequency
In Section 3.2 we showed that self-excited resonators can produce limit cycle oscillations, their frequency monotonically depends on the resonator's resonant frequency, and relaxation oscillations develop for high enough levels of self excitation. Here we extend our results to include two-cell networks. Because self-excited resonators may be thought of as representing a population of synchronized in phase cells, we expect our results from Section 3.2 to hold of these networks. However, the presence of nonlinearities of cubic type are not apparent from either the model equations or the phase-space diagrams and need to be uncovered using the method developed in Section 3.3. Figure 12a1 show the small amplitude oscillations generated in a Hopf bifurcation (Fig. 14) for low enough values of G ex (G ex,1,2 = G ex,2,1 ). This oscillations are not identical, because the cells are not identical, but they are synchronized in phase. Figure 12b1 shows that increasing values of G ex lead to oscillations of relaxation type. The dynamic mechanisms of oscillation amplification ( Fig. 12a2 and b2) as well as the cubic-based mechanisms of generation of relaxation oscillations (Fig. 12a3 and b3) are analogous to the mutually inhibitory networks discussed in Section 3.3.

The resonator's intrinsic frequency controls the network oscillations frequency
Our results are presented in Fig. 13 for (i) values of Z 1,max that increase from panel a to b (for a fixed-value of g L,1 ), and (ii) values of g L,1 that decrease from panel b to c (for a fixed values of Z 1,max ). In contrast to the mutually inhibitory networks, by increasing Z 1,max , the oscillatory active resonant frequency band is increased, while the onset of oscillations occurs for lower values of f res , similarly to mutually inhibitory networks. The opposite behavior is observed for decreasing g L,1 with all the other parameters fixed. The behavior of the oscillations amplitude is similar to the mutually inhibitory networks (Fig. 14).

Graded mutually inhibitory or excitatory 2D/2D resonator networks generate sustained (limit cycle) oscillations and their frequency interact to control the network frequency
Here we extend our results from Sections 3.3 and 3.4 to networks having two mutually connected 2D resonators (that are not damped oscillators). We consider heterogeneous networks of non-identical resonator in order to test the effects of the interaction band-pass filters with different frequency bands. Because the mechanisms of generation of oscillations are similar to those discussed in Sections 3.3 and 3.4, we focus on the effects of the resonant frequencies of the participating resonators on the network oscillation frequency. Our results are presented in Fig. 15. The gray curves corresponds to networks of resonators with the same frequency band. The network model is given by system (1)-(4) with g 1 , g 2 > 0. Figure 15 shows the dependence of the network frequency on f 1,res for representative values of f 2,res and other model parameters for mutually inhibitory (panels a and b) and mutually excitatory (panels c) networks. As expected, in all cases the network frequency monotonically depends on the resonant frequency of both oscillators.
The range of values of f 1,res for which network oscillations are possible increase with increasing values of f 2,res as does the network frequency. However, the network frequency and the resonant frequency of the oscillators is no longer one-to-one, as was the case for the 2D/1D networks investigated in Sections 3.3 and 3.4, but it depends on the a g L,1 = 0.25, Z 1,max = 3.9. b g L,1 = 0.25, Z 1,max = 3.7. c g L,1 = 0.1, Z 1,max = 6. We used the following parameter values: g L,2 = 0.5, E in = −20, v hlf = 0, v slp = 1 Fig. 11 Bifurcation diagrams for mutually inhibitory resonatorpassive cell networks (2D/1D) for representative parameter values. The shadowed region corresponds the existence of sustained (limit cycle) oscillations. The green-lined region corresponds to multistability (limit cycle and/or fixed-points). The inset trajectory diagrams indicated the dynamics within the regions bounded by the solid and dashed curves (except the solid green curve): stable nodes, stable foci, unstable foci and unstable nodes (from left to right). The inset diagrams correspond to the 3D linearized system for the fixed-point before the static bifurcation. H 0 , H 1 and H 2 note the Hopf bifurcation branches, P F notes the pitchfork bifurcation branch and SN notes the saddle-node branch. a Bifurcation diagram in G in -τ 1 parameter space. Cell 1 is a resonator for values of τ 1 > τ 1,res (dashed-black horizontal line). a1 g L,1 = 0.25 and g 1 = 0.25. a2 g L,1 = 0.25 and g 1 = 0.3. b Bifurcation diagram in G in -g 1 parameter space for g L,1 = 0.25 and τ 1 = 100. c Bifurcation diagram with G in as bifurcation parameter. The solid-and dashed-blue curves represent stable and unstable fixed-points, respectively. The solid-and dashed-black curves represent the stable and unstable limit cycle branches created at the Hopf bifurcations (red dots). c1 g L,1 = 0.25, g 1 = 0.25 and τ 1 = 40. c2 g L,1 = 0.25, g 1 = 0.3 and τ 1 = 40. We used the following parameter values: g L,2 = 0.5, E in = −20, v hlf = 0, v slp = 1 complex interaction between the two resonators. The oneto-one dependence between the network frequency and the resonant frequencies of the individual oscillators occurs when the two resonators have the same resonant frequency (black dots). The slopes of the network frequency curves for non-identical resonators are smaller than for identical resonators indicating that the network frequency is larger (smaller) than the resonant frequency to the left (right) of the black dot. This is independent of the mechanisms of amplification (mutual inhibition or mutual excitation).
Increasing values of Z 2,max for fixed values of Z 1,max and f 2,res causes the range of values of f 1,res for which oscillations exist increase, but the slope of the network frequency curves remains almost unchanged, independently of the value of g L,1 used and whether the mechanisms of amplification is based on mutual excitation or mutual inhibition (not shown).

Discussion
Network oscillations emerge from the cooperative activity of the intrinsic properties of the participating neurons (e.g., ionic currents) and the synaptic connectivity. Their G ex,1,2 =G ex,2,1 =0.04 Fig. 12 Oscillations generated in mutually excited hybrid 2D-1D networks. Cell 1 is a resonator with f res = 8 (f nat = 0) and cell 2 is a passive cell. Left. Voltage traces (curves of v 1 and v 2 as a function of t ). Middle. Phase-plane diagrams. The v 1 -and v 2 -hypernullclines are given by Eqs. (10) and (11), respectively. The dashed red curve represents the v 1 nullcline for cell 1 for g 1 = 0 (no resonant gating variable). Right. Adapted phase-plane diagrams relative to the v 2 -hyper-nullcline N 2c (v 1 ) (green curve in the phase-plane diagrams in the left panels). The red lines are the differences between the v 1 -and v 2 -hyper-nullclines in the left panels parametrized by constant values of w 1 . The solid-red curve corresponds to an intermediate value of w 1 . The dashed-red curves correspond to the maximal w 1,max (lower) and minimal w 1,min (upper) values of w 1 . The trajectories (blue curves) are also plotted relative to the v 2 -hyper-nullcline N 2c (v 1 ). a G ex,1,2 = G ex,2,1 = 0.032. The network frequency is f ntw ∼ 7.1. b G ex,1,2 = G ex,2,1 = 0.04. The network frequency is f ntw ∼ 4.3. We used the following parameter values: g 1 = 1.8, g L,1 = 0.1, g L,2 = 1, τ = 750, E ex = 60, v hlf = 0, v slp = 1 generation and dynamics involve the nonlinearities and time scales present in the circuit components and the time scales that emerge from their interplay. Several neuron types have the intrinsic ability to generated membrane potential oscillations under blockade of all synaptic connectivity. Other neuron types do not show intrinsic membrane potential oscillations, but exhibit membrane potential resonance (MPR). MPR is a property of the interaction between oscillatory inputs and the intrinsic neuronal properties (intrinsic resonant and amplifying processes) that uncovers a circuit latent time scale associated to the resonant frequency (MPR can be observed in the absence of intrinsic damped oscillations Richardson et al. 2003;Rotstein and Nadim 2014b). This hidden time scale (provided by the resonant frequency) is encoded in this impedance profile, and therefore, the impedance profile is the object of study for a non-oscillatory resonant neurons.
MPR has been investigated both experimentally and theoretically in many neuron types (Hutcheon et al. 1996(Hutcheon et al. , 2000Richardson et al. 2003;Lampl and Yarom 1997;Llinás and Yarom 1986;Gutfreund et al. 1995;Erchova et al. 2004;Schreiber et al. 2004;Gastrein et al. 2011;Hu et al. 2002Hu et al. , 2009Johnston 2007, 2008;Fig. 13 Oscillations in mutually excited hybrid 2D-1D networks: the intrinsic resonant frequency controls the network frequency. Left columns. Network oscillation frequency as a function of f res . Right columns. Network oscillation amplitude (oscillator 1) as a function of f res . The synaptic conductances G ex,1,2 = G ex,2,1 are equal to the values reported in the figure. a g L,1 = 0.1, Z 1,max = 9.2. b g L,1 = 0.1, Z 1,max = 9.87. c g L,1 = 0.08, Z 1,max = 9.87. We used the parameter value g L,2 = 1, Fig. 14 Bifurcation diagrams for mutually excitatory resonatorpassive cell networks (2D/1D) for representative parameter values. The shadowed region corresponds the existence of sustained (limit cycle) oscillations. The green-lined region corresponds to multistability (limit cycle and/or fixed-points). The inset trajectory diagrams indicated the dynamics within the regions bounded by the solid and dashed curves (except the solid green curve): stable and unstable foci. H 0 and H 1 note the Hopf bifurcation branches. a Bifurcation diagram in G ex -τ 1 parameter space for g L,1 = 0.1 and g 1 = 2. Cell 1 is a resonator for values of τ 1 > τ 1,res ∼ 0.205: b Bifurcation diagram in G ex -g 1 parameter space for g L,1 = 0.1 and τ 1 = 100. We used the following parameter values: g L,2 = 1.2, E ex = 60, v hlf = 0, v slp = 1 any functional role for network oscillations or is simply an epiphenomenon is largely an open question. A few studies have investigated the oscillatory properties of networks including neurons that exhibit MPR Stark et al. 2013;Tikidji-Hamburyan et al. 2015;Tchumatchenko and Clopath 2014;Schmidt et al. 2016;Moca et al. 2014;Baroni et al. 2014;Rotstein et al. 2017d) or have resonant gating variables (Wang and Rinzel 1992;Manor et al. 1997;Torben-Nielsen et al. 2012). But the role that MPR plays in the generation of network oscillations and how the latent time scales affect the properties of the oscillatory networks in which they are embedded remained to be understood.
Addressing these questions is not straightforward. The concept of a hidden (latent) time scale is somehow abstract in the sense that it is the response the neuron would have in the presence of an externally imposed oscillatory input and not a property directly measurable in the individual neurons as, say, intrinsic membrane potential oscillations. It is also difficult to manipulate, because, as we discuss in the paper (see also Chen et al. 2016), a proper comparison between the effects of different resonant frequencies would require sliding the impedance profile along the resonant frequency line while keeping the impedance profile shape unchanged (or as unchanged as possible), and this requires changing more than one model parameter, contrary to the standard mechanistic approach of changing one parameter at the time, while keep the remaining ones unchanged. Far from being simply a theoretical issue disconnected from the biological reality, the same approach and method should be used for the experimental determination of the role of resonance for network oscillations , particularly to experimentally test the predictions of our work (e.g., using the dynamic clamp technique).
In this paper we set out to investigate these issues using minimal network models consisting of non-oscillatory resonators mutually coupled to either a low-pass filter neuron or another band-pass filter (resonator). In this way we could separate the different effects that give rise to network oscillations in two different levels of organization that can be manipulated separately. The resonator provides the negative feedback and the network connectivity provides the amplification. Because we leave out resonators that can be also damped-oscillators, the network oscillations are not inherited from the individual cell level, but are created by the combination of the individual cell and connectivity properties.
We showed that oscillations can be generated in networks of increasing complexity: (i) self-excited band-pass filters, (ii) mutually inhibited band-and low-pass filters, (iii) mutually excited band-and low-pass filters, (iv) mutually inhibited band-pass filters, and (v) mutually excited band-pass filters. The presence of a resonator is necessary to generate oscillations in these networks; if the resonators are substituted by low-pass filters, network oscillations are not possible. However, what characterizes the oscillatory activity of a resonator is the resonant frequency, which cannot be assessed in the absence of oscillatory inputs. By showing that the network frequency monotonically depends on the resonant frequency of the individual band-pass filters, we provide a direct link between MPR and the generation of network oscillations. To our knowledge, this is the first time such a link is provided. A similar result was obtained in Fig. 15 Oscillations in mutually inhibitory or excitatory resonator cell networks (2D-2D): the intrinsic resonant frequencies interact to control the network frequency. Left columns. Network oscillation frequency as a function of f res . Right columns. Network oscillation amplitude (oscillator 1) as a function of f res . The gray curves correspond to a network of identical cells for fixed values of g L,1 = g L,2 and Z 1,max = Z 2,max . The colored curves correspond to a fixed cell 2 with the resonance frequency f 2,res indicated in the figure. a g L = 0.1, Z max = 6, G in = 0.1. b g L = 0.25, Z max = 3.7, G in = 0.1. c g L = 0.1, Z max = 9.2, G ex = 0.03 electrically coupled networks ), but in these cases, the network oscillations were driven by one of the nodes that was a sustained oscillator. Network oscillations have been shown to emerge as the result of the interaction of damped oscillators (Torben-Nielsen et al. 2012;Loewenstein et al. 2001;Manor et al. 1997), but in these cases, the network oscillations are inherited from the oscillatory activity of the individual intrinsically oscillatory nodes.
The existence of sustained oscillations in networks of non-oscillatory neurons is not without precedent. The inferior olive oscillatory network studied in Manor et al. (1997) is composed of electrically coupled neurons that, when isolated, are damped oscillators. In this case, the individual neurons are nonlinear and include both resonant and amplifying effects, but the connectivity is linear. The model investigated in Wang and Rinzel (1992) involves nonlinear neurons reciprocally inhibited with graded synapses. For baseline values of the DC input current, the neurons are quasi-linear and are at most damped oscillators. The nonlinearities developed for negative input current values combined with the dynamics resulting from the mutual synaptic inhibition result in the postinhibitory rebound mechanism underlying the observed network oscillations. Post-inhibitory rebound (PIR) and subthreshold resonance are closely related phenomena since both require the presence of a negative feedback effect, but they are different in nature. The mechanisms investigated in Wang and Rinzel (1992) depend crucially on the effective pulsatile nature resulting from the dynamic interaction between cells and synaptic connectivity. Models having an h-current also show PIR. Even in the presence of an (additive) amplifying current, such as the I h + I Nap model, the functional connectivity in these models is not PIR-based, but rather resonance-based as we show in this paper (unpublished observation). In Manor et al. (1999) and Ambrosio-Mouser et al. (2006) oscillations emerge in two reciprocally inhibited passive cells where one of them is self-excited, thus providing additional dynamics to the network. The model studied in Chen et al. (2016) consists of an oscillator electrically coupled to a follower resonator whose intrinsic resonant frequency directly affects the network frequency while the shape of the impedance profile remains almost unchanged.
The minimal models we used in this paper serve the purpose of establishing the role of MPR for the generation of network oscillations. Other types of models could include resonant properties at the network level. Moreover, there are alternative possible scenarios where, for example, amplification occurs at the single cell level and the negative feedback effect occurs at the network level. These types of networks are beyond the scope of this paper. The understanding of the oscillatory properties of such networks requires more research.
The types of models we used could be argued to be too simplistic and not realistic. We used these models precisely because of their simplicity in order to understand some conceptual points that can be generalized and applied to more realistic networks. However, one should note that the type of models we used are very close to the firing rate models of Wilson-Cowan type (Wilson and Cowan 1972) with adaptation (Curtu and Rubin 2011;Shpiro et al. 2009;Tabak et al. 2011), which are essentially resonators (unpublished observation). In this models, the nonlinearity is similar to the one we used (sigmoid type, instantaneously fast). Therefore, our results can be easily generalized to these models. One example are the networks of OLM cells and fast spiking (PV + ) interneurons (INT) that have been shown to be able to produce network oscillations (Gillies et al. 2002;Rotstein et al. 2005). OLM cells show MPR , while the presence of MPR in INT is debated (Pike et al. 2000;Zemankovics et al. 2010). Although our models are simplistic, they make predictions that can be tested using the dynamic clamp technique (Sharp et al. 1993;Prinz et al. 2004).
The results presented in this paper advance the conceptual understanding of the oscillatory interaction among nodes in a network, particularly when there are hidden time scales, and proposes ideas to understand the dynamics of these networks. We open several questions regarding the ability of networks of band-pass filters to generate oscillatory patterns and how the properties of these patterns depend on the properties of these filters. More research is required to address these issues.

A.1 Linearization and eigenvalues
The linearization of system (1) with g k = 0 (k = 1, 2) and I syn,k given by Eqs. (3) and (4) reads where The eigenvalues (r 1 and r 2 ) are given by The first two terms in Eq. (22) are always negative (provided g L,1 > 0 and g L,2 > 0). The second term in the radicand is positive if F v 2 and G v 1 have the same sign and negative if F v 2 and G v 1 have different signs. Therefore, the fixedpoint for networks with the same type of connections (both excitatory or both inhibitory) can be either stable nodes or saddles, while the fixed-points for excitatory-inhibitory networks can be either stable nodes or stable foci (e.g., Fig. 16).

A.2 Absence of limit cycles
We compute The v 1 -and v 2 -nullclines are given by Eqs. (5) and (6), respectively. Black dots indicate stable fixed-points (nodes) and gray dots indicate unstable fixed-points (saddles). Cells and connectivity are identical. The parameter G in represents G in,1,2 = G in,2,1 . As G in increases (from a1 to a3), the v 1 -and v 2 -nullclines transition from quasi-linear to nonlinear and the system undergoes a pitchfork bifurcation as a stable fixed-point (a1) looses stability and two additional stable fixed-points are created (a3). Heterogeneous networks have non-symmetric phase-plane diagrams and show a qualitatively similar behavior, but bistability results from saddle-node bifurcations. We used the following parameter values: g L,1 = g L,2 = 0.25, E in,1 = E in,2 = −20, v hlf = 0, v slp = 1. Excitatory-inhibitory networks of passive cells.
The v 1 -and v 2 -nullclines are given by Eqs. (5) and (6), respectively. Black dots indicate stable nodes and gray dots indicate stable foci. The parameter G ex,1,2 = 0.01 is fixed. As G in,2,1 increases (from b1 to b3), the fixed-point transitions from stable nodes to stable foci and back to stable nodes. We used the following parameter values: g L,1 = g L,2 = 0.25, E in,1 = −20, E ex,2 = 60, v hlf = 0, v slp = 1 Since U < 0 for all v 1 and v 2 (provided g L,1 > 0 and g L,2 > 0), then by the Bendixson-Dulac theorem (Guckenheimer and Holmes 1983), there are no limit cycles in the (v 1 , v 2 )-plane. This argument breaks down when either g L,1 < 0 or g L,2 < 0 and small enough, indicating a strong positive feedback effect generated by an ionic process.

Appendix B: Dynamics of autonomous and forced 2D cells
We consider the following system where the parameters g L , g, C and τ are as in system (1)-(2) by omitting the subindex (k), A in is the input amplitude, and f is the input frequency. We assume here that all intrinsic parameters (g L , g, C and τ ) are positive. The constraint g > 0 indicates that the ionic current that the term g w linearizes is a resonant process (negative feedback) (Richardson et al. 2003;Rotstein and Nadim 2014b). The linearized parameter g L captures the effects of the biophysical leak current and possibly another ionic amplifying process (positive feedback) provided by an additional current. Strong enough amplifying processes may cause g L to be negative.
System (1)-(2) (A in = 0) has a uniquefixed-point (v,w) = (0, 0). This fixed-point is stable provided g L τ + C > 0. It is a stable node if the radicand in Eq. (26) is nonnegative and a stable focus otherwise. We refer the reader for details on the dependence of the fixed-point type (node or focus) with the model parameters to Rotstein and Nadim (2014b).
An important aspect to note, relevant for this paper, is that resonance can occur in the absence of damped oscillations; i.e., when the fixed-point is a stable node. In this paper we focus on resonators that do not show damped oscillations.

B.3 Quasi-displacement of impedance profiles: fixed peak values and changing resonant frequencies
From Eq. (29) we can compute the value of g as a function of Z max and the other model parameters g = (Z 2 max + τ 2 − Z 2 max g 2 L τ 2 ) 2 4 Z 2 max τ (−τ 2 + Z 2 max (1 + g L τ ) 2 ) .
Equation (30) relates the model parameters of a forced 2D linear system of the form (24)-(25) for which the impedance peak Z max is constant. In order to calculate the balanced values of g and τ , if they exist, for given values of Z max and g L (fixed) we proceed as follows. First we take values of τ within certain range and compute the corresponding values of g using Eq. (30). For these values of g L , g and τ we compute ω res using (28) and C = 1. In this way we have g = g(τ ) and ω res = ω res (τ ) for a given value of Z max . values provided S ∞ (v = E in ) is small enough (the sigmoid function S ∞ changes fast enough around v hlf and is negligible at v = E in ).
From our discussion above (Section B.1), the uncoupled system (G in = 0) has a stable fixed-point. We expect this to persist for small enough values of G in .

C.1 Fixed point
The fixed-points of the self-inhibited 2D system are the zeros of whose derivative is given by The first two terms in Eq. (32) are negative, while the third one is negative provided v > E in . However, for v < E in this third term is negligible. Therefore, H (v) is a decreasing function for all v. Because H (v) < 0 for large enough values of v, a fixed-point exists if H (v) > 0 for some v. The first term in Eq. (31) is positive for negative values of v and so is the second term provided v < E in . Therefore, the self-inhibited cell has a unique fixed-point (v * ,v * ). Since H (0) = G in S ∞ (0)E in < 0 and H (E in ) = −(g L + g) E in > 0, then E in <v * < 0. The stability properties of the fixed-point (v * ,v * ) are determined by looking at the equation for the eigenvalues (26) with g L substituted by g * L = g L + G in S ∞ (v * )(v * − E in ) + G in S ∞ (v * ) > g L . Therefore, the stability of the fixed-point is preserved. If (v,v) is a stable node, then (v * ,v * ) is a node for all values of G in (Fig. 17a). In contrast, if (v,v) is a stable focus, then (v * ,v * ) remains a stable focus for small enough values of G in , but it transitions to a stable node for large enough values of G in (Fig. 17b). The vand w-nullclines are given by (8) and (9), respectively. Black dots indicate stable nodes and gray dots indicate stable foci. a g L = 0.25. The fixed-point for the uncoupled system is a stable node. b g L = 0.01. The fixed-point for the uncoupled system is a stable focus. We used the following parameter values: g = 0.25, τ = 100, E in = −20, v hlf = 0, v slp = 1

C.2 Absence of limit cycles
Because the self-inhibited 2D systems is relatively simple we do not expect the existence of limit cycles. We address this in the region R = {(v, w) ∈ R 2 : v > E in }. We compute By substituting S ∞ (v) = S ∞ (v)(1 − S ∞ (v)) > 0 where S ∞ , given by Eq. (4), we obtain If v > E in (assuming g L > 0), then U < 0. Thus, by the Bendixson-Dulac theorem (Guckenheimer and Holmes 1983), there are no limit cycles lying entirely in the region R.