## Abstract

Several types of intrinsic dynamics have been identified in brain neurons. Type 1 excitability is characterized by a continuous frequency-stimulus relationship and, thus, an arbitrarily low frequency at threshold current. Conversely, Type 2 excitability is characterized by a discontinuous frequency-stimulus relationship and a nonzero threshold frequency. In previous theoretical work we showed that the density of K_{v} channels is a bifurcation parameter, such that increasing the K_{v} channel density in a neuron model transforms Type 1 excitability into Type 2 excitability. Here we test this finding experimentally, using the dynamic clamp technique on Type 1 pyramidal cells in rat cortex. We found that increasing the density of slow K_{v} channels leads to a shift from Type 1 to Type 2 threshold dynamics, i.e., a distinct onset frequency, subthreshold oscillations, and reduced latency to first spike. In addition, the action potential was resculptured, with a narrower spike width and more pronounced afterhyperpolarization. All changes could be captured with a two-dimensional model. It may seem paradoxical that an increase in slow K channel density can lead to a higher threshold firing frequency; however, this can be explained in terms of bifurcation theory. In contrast to previous work, we argue that an increased outward current leads to a change in dynamics in these neurons without a rectification of the current-voltage curve. These results demonstrate that the behavior of neurons is determined by the global interactions of their dynamical elements and not necessarily simply by individual types of ion channels.

- bifurcation
- dynamic clamp
- ion channel density
- pyramidal neuron
- threshold dynamics

the threshold dynamics of a neuron is a key determinant of its spiking characteristics and, mainly via influencing synchronization, of its interaction with other neurons (Gouwens et al. 2010; Izhikevich 2010; Ratte et al. 2013; Zeberg et al. 2010, 2012). Consequently, the threshold dynamics of single neurons is an essential factor in determining the spiking activity of brain networks in general.

One way to classify the spiking of a neuron is according to its frequency-stimulation relation. Already in 1948, Hodgkin did this in a study of crab axon activity. Spiking characterized by a continuous relationship between frequency and stimulation current (*f-I*_{stim}), and thus by an arbitrarily low frequency at threshold current, he named Class 1. Spiking characterized by a discontinuous *f-I*_{stim} relationship and a distinct minimum threshold frequency he named Class 2. This classification is now generally used, although often referred to as Type 1 and Type 2 (Tateno and Robinson 2006). A remarkable phenomenon that Hodgkin reflected upon in this study was that a dynamical system composed of processes with time constants in the millisecond spectrum can operate on the scale of seconds. The introduction of dynamical systems theory into the study of cellular neurophysiology has revealed how this “slowness” can emerge (FitzHugh 1961; Izhikevich 2010). Essentially, the spiking properties of neurons can be explained in terms of four generic bifurcations. Type 1 behavior can be described with a dynamical system with three fixed points or equilibria, two of which merge into a quasi-stable equilibrium at threshold (leading to a saddle-node bifurcation on an invariant circle, SNIC) (Arhem et al. 2006; Arhem and Blomberg 2007; Ermentrout and Terman 2010; Izhikevich 2010; Zeberg et al. 2010). The existence of the quasi-stable equilibrium point in the trajectory leads to spiking at frequencies that are arbitrarily close to zero at threshold. This type of spiking behavior can be found in cortical pyramidal neurons (Izhikevich 2010; Tateno and Robinson 2006). More complex spiking patterns, comprising bursting and adaptive processes, have also been described in some pyramidal cells (Izhikevich 2010), and it seems that an interplay between the axon initial segment and the first node of Ranvier is essential for neuronal burst firing in layer 5 “intrinsic bursting” cells (Kole 2011).

Fast-spiking (basket morphology, parvalbumin expressing) interneurons often show another basic spiking behavior that is generically described by a dynamical system that has a single equilibrium point that becomes unstable at increased stimulation through a local subcritical Andronov-Hopf bifurcation. Already at lower stimulation this leads to repetitive spiking via a global fold limit cycle. This spiking pattern is characterized by a discontinuous (Type 2) stimulus-frequency curve, often with a minimum frequency of ∼20–50 Hz (Tateno et al. 2004). It should be noted, however, that saddle-node (but not SNIC) and supercritical Andronov-Hopf bifurcations also lead to Type 2 behavior, and that even if fast-spiking interneurons mainly spike with a minimum frequency of 20–50 Hz, Type 2 excitability can generate low firing rates in the presence of noisy membrane voltage fluctuations (Higgs et al. 2006) and near the transition to Type 1 dynamics (Zeberg et al. 2010).

Different types of interneurons show diverse spiking patterns, such as irregular spiking, delayed onset, and bursting. Nevertheless, regular and fast spiking can be considered the predominant patterns of periodic spiking in the mammalian neocortical network. Table 1 summarizes some electrophysiological and dynamical system properties of these two basic spiking patterns.

In previous theoretical studies (Arhem et al. 2006; Arhem and Blomberg 2007; Zeberg et al. 2010), we have shown that for a dynamical neuron model that is characterized by a SNIC bifurcation at threshold, an increase in the density of the K_{v} channels could lead to a reduction in the number of near-threshold equilibrium points to 1, leading to a subcritical Andronov-Hopf bifurcation and change in the threshold firing rate from essentially 0 to a distinct minimum value, typically 20–50 Hz.

The excitability and bifurcation patterns of this neuron model were shown to form two main regions in the Na and K channel density plane, one being associated with a nonmonotonic *I*_{ss}-*V* curve (where *I*_{ss} denotes the steady-state current and *V* denotes the membrane potential) and with one or three equilibrium points and the other being associated with a monotonic *I*_{ss}-*V* curve and with one equilibrium point. In the monotonic region the oscillations always had a distinct onset, and the stability was lost via a fold limit cycle. The fold limit cycle bifurcation could be associated with a subcritical Andronov-Hopf bifurcation of the equilibrium point, or the bifurcation could occur because of a purely global bifurcation without the change of local stability of any associated equilibrium point. In the region with the possibility of three equilibria, we found regions of excitability and repetitive spiking. The oscillation could arise through a SNIC bifurcation or through a more complex geometry (see Zeberg et al. 2010). At the border between the saddle-node bifurcation area and the Andronov-Hopf bifurcation area, the bifurcations are of a codimension 2 type [Bogdanov-Takens bifurcations (Takens 1974)]. In this region the threshold frequency is discontinuous but typically low (see Zeberg et al. 2010). It transpires that this, somewhat exotic, region in our previous work will be of great importance for pyramidal neurons, as we demonstrate below.

The present study is an experimental test of these theoretical findings. We chose to test cortical pyramidal cells showing a Type 1 spiking pattern (Tateno et al. 2004). Using the dynamic clamp technique (Destexhe and Bal 2009; Economo et al. 2010; Robinson and Kawai 1993; Sharp et al. 1993) on pyramidal cells in rat neocortical brain slices, we analyzed the effects of injecting currents corresponding to those of a naturally occurring delayed-rectifier channel, referred to here as the slow K_{v} channel (Korngreen and Sakmann 2000; see methods). As mentioned above, a Type 1 spiking pattern is due to a SNIC bifurcation and requires three equilibrium points and an N-shaped *I*_{ss}-*V* curve. An N-shaped curve can be straightened by an increased delayed-rectifier K_{v} current, reducing the number of equilibrium points from three to one and altering the spiking pattern to Type 2. Injecting delayed-rectifier current is therefore predicted to switch the Type 1 pyramidal neurons to Type 2 neurons. Such a mechanism has been used to explain the switch in spiking type observed when leak current (shunt) is injected, or when an increase of M-type K_{v} current (K_{v}7) or Ca-activated K current are induced, in CA1 hippocampal pyramidal neurons (Prescott et al. 2008).

In agreement with the theoretical predictions, we found that increasing the slow K_{v} current in Type 1 pyramidal neurons gave rise to a robust switch to Type 2 behavior. However, a low-dimensional model revealed that the shift from Type 1 to Type 2 spiking does not occur exclusively because of a rectification of an N-shaped *I–V* curve but, at less K_{v} current injection, because of a surprisingly complex mechanism. The mechanism suggested here also holds for a switch induced by an increased leak current, and might possibly explain, at least partially, the observations by Prescott et al. (2008).

## METHODS

#### Ethics statement.

All animals were handled in strict accordance with good animal practice as defined by UK Home Office regulations and killed according to UK Home Office-approved Schedule 1 procedures. All animal work was approved by the University of Cambridge.

#### Cell and animal characteristics.

The study comprises electrophysiological recordings from 33 cells, 18 of which were chosen for a detailed analysis. Cells were characterized as pyramidal cells on the basis of their morphology. The 15 cells that were not analyzed either were of low quality or were native resonators (*n* = 5). The age of the analyzed cells was 13–17 days (except 1 cell from a 8-day-old rat), and the mean age was 14.5 days. We found no substantial difference in firing frequency over the temperature range studied (22–34°C). The kinetic description of the K channels is based on measurements at room temperature (Korngreen and Sakmann 2000). The mean membrane capacitance among the 18 cells was 111 pF with a standard deviation of 63 pF.

#### Slice preparation and electrophysiological recordings.

Wistar rats were killed by cervical dislocation. Sagittal slices of the somatosensory cortex (300 μm) were prepared with a vibratome (DSK Microslicer Zero 1, Dosaka EM, Kyoto, Japan) in a chilled solution composed of (in mM) 125 NaCl, 25 NaHCO_{3}, 2.5 KCl, 1.25 Na_{2}PO_{4}, 2 CaCl_{2}, 1 MgCl_{2}, and 25 glucose and oxygenated with 95% O_{2}-5% CO_{2} gas. Slices were held at room temperature for at least 30 min before recording. The tissue was visualized with an Olympus BX50WI upright microscope (Olympus UK, London) using infrared differential interference contrast videomicroscopy. During recording, slices were perfused with oxygenated solution identical to the slicing solution. Whole cell recordings were made from the somas of pyramidal neurons in cortical layers 2/3 and 5. Regular-spiking pyramidal cells were identified by their morphology and threshold firing frequency. Patch pipettes of 4- to 8-MΩ resistance were pulled from borosilicate capillary glass and filled with an intracellular solution containing (in mM) 105 K-gluconate, 30 KCl, 10 HEPES, 10 phosphocreatine, 4 ATP, 4 MgClB_{2}B, and 0.3 GTP adjusted to pH 7.3 with KOH. Current-clamp recordings were performed with an Axon Multiclamp 700B amplifier (Axon Instruments, Foster City, CA). Membrane potentials were corrected for the prenulled liquid junction potential. Signals were filtered with a four-pole low-pass Bessel filter, with a cutoff frequency at −3 dB of 5 kHz, sampled at 20 kHz or faster, and recorded with custom software written in MATLAB (The MathWorks, Natick, MA).

#### Conductance injection.

Virtual ion channels were introduced with conductance injection (see Robinson and Kawai 1993; Sharp et al. 1993). Based on Ohm's law *I* = *g*(*V* − *E*_{rev}) where *g* is the conductance, *V* is the membrane potential of the cell, and *E*_{rev} is the reversal potential of the conductance, an effective conductance was inserted in the recorded cell by injecting the corresponding current. Conductance injection software running on a dedicated digital signal processing board (SM-2, Cambridge Conductance, Cambridge, UK) implemented this equation with a response time of <50 μs, producing the current command signal for the current-clamp amplifier. The density of the injected conductance and the density of the injected stimulation current were calculated from estimates of the membrane area based on measurements of time constants and steady-state values of subthreshold current step responses(*C* = τ × *I*/*V*) and the assumption of a capacitance density of 0.01 F/m^{2}.

#### Ion channel kinetics.

The injected K channel conductance was modeled to mimic a tetraethylammonium-sensitive slow K current, characterized in voltage-clamp experiments in the somata of neocortical pyramidal cells in young rats (postnatal days 13–15) (Korngreen and Sakmann 2000). The K current was implemented as *I*_{Kv} = *ḡ*_{Kv}*n*^{2}*h*(*V* + 85.5 mV), where the time evolution of *n* and *h* are given by first-order kinetics with the following rate functions:

The equations given in Korngreen and Sakmann (2000) could reproduce neither the fitted functions nor their voltage clamp data; therefore, the rate functions were changed slightly (see Fig. 1). Moreover, the slow inactivation had two time constants that were reduced to a single component. However, as can be seen in Fig. 1, the above functions satisfactorily reproduce the voltage-clamp data of Korngreen and Sakmann (2000).

#### The model.

The reduced pyramidal cell model is described by the following equations:

where α and β are given above and the gating of the sodium channels is given by the following two functions:

where

where the superscript −1 denotes the inverse function. The parameters describing the Na current (*V*_{1/2} and *s* for *m* and *h*, respectively) were taken directly from measurements of Na_{v}1.6 kinetics (Mercer et al. 2007), which represents the most abundantly expressed Na channel subtype in the CNS (Goldin 2001). The magnitude of the conductances are within the physiological range (Tables 2 and 3). The reversal potential for the leak current was set to match the resting potential (−80 mV). The reversal potential for the K current was chosen to be the same as that used for the dynamic clamp equations (−85.5 mV), and the reversal potential for the Na current was calculated with the Nernst equation, given the Na^{+} concentrations (+66.3 mV). The value chosen for the capacitance (150 pF) was within the range estimated in this study (111 pF ± 63 pF, *n* = 18). All parameter values used are given in Table 2. For Fig. 2, uniform noise was added at a relatively low level that was sufficient to explain the experimental results. The amplitude of the noise was set to 0.03% of the *n-V* plane in both directions, and a waveform of noise was added to the differential equations. Noise values were generated from a uniform distribution between these limits, updated every 5 ms.

#### Stability analysis and path-following algorithm.

The stability of the equilibrium points was analyzed by linearizing the system, i.e.,

where **J** is the Jacobian matrix evaluated at the equilibrium points whose eigenvalues determine the stability of the system (that the system's behavior is locally qualitatively the same as the linearized version is ensured by the Hartman-Grobman theorem). For a more detailed description of the stability analysis, see Zeberg et al. (2010).

A path-following algorithm for the bifurcation diagrams (see Fig. 8) was devised as follows: a starting point *p*_{0} of a border of interest was located, and custom-written software then searched in a circle with a sufficiently small radius around *p*_{0} in steps of 0.05° to find *p*_{1}. To locate *p*_{n+1}, a similar search was conducted around *p*_{n}. To make the algorithm more efficient, the search was limited to a ±2° deviation from the tangent.

#### Fourier analysis.

Fourier analysis was performed on interspike intervals (ISIs) ≥ 100 ms. A 10-ms segment was removed from the beginning of the interval [to account for the afterhyperpolarization (AHP) of the previous spike], and 1 ms was removed from the end of the ISI (to remove the upstroke of the following action potential). The linear trend in the resulting waveform was removed, and the mean was subtracted.

#### Cluster analysis.

For our analysis of the firing frequency, we first averaged the data for each *I*_{stim} and *g*_{K} in a pooled data set consisting of 21,567 ISIs from 18 cortical pyramidal neurons. In theoretical models (Izhikevich 2010; Zeberg et al. 2010) the *f-I*_{stim} relationship is discontinuous, but at applied currents above threshold the *f-I*_{stim} relationship exhibits a continuum. However, our experimental data demonstrated a higher variability due to spike failure and subthreshold oscillations. To obtain robust estimates, we used DBSCAN (density-based spatial clustering of applications with noise; Ester 1996) to search for clusters. DBSCAN was implemented as a function call to ELKI 0.5 (Achtert 2012) from MATLAB. The definition of a cluster is established on the concept of density reachability. A data point is density-reachable if it is not beyond a given distance ε and is surrounded by sufficiently many points *k* in its ε-neighborhood. In Fig. 6*A*, ε was set to 5% of the space and *k* was set to 3.

## RESULTS

#### K_{v} channel effects on spiking patterns.

Figure 2 shows the effect of injecting a dynamic conductance, which mimics the electrical effect of K_{v}2.1 channels, in rat cortical pyramidal cells during current step responses. The description of the kinetics is based on voltage-clamp recordings in cortical pyramidal cells (Korngreen and Sakmann 2000). Injection of this dynamic conductance switched the spiking from one behavior, comprising *1*) relatively long time delays to the first spike and relatively ISIs at near-threshold stimulation, to another behavior, comprising *2*) short delays to the first spike, subthreshold spikelets, shorter ISIs, and faster AHPs, suggesting that the neuron had switched from Type 1 to Type 2 excitability. According to the phenomenological definition of Type 1 dynamics, it should be possible to demonstrate essentially zero frequency spiking in Type 1 neurons. In practice, it is difficult to stay at the metastable potential, the so-called attractor ruin (Izhikevich 2010), for longer than a second, which is why it is difficult to obtain firing frequencies below 1 Hz. Similarly, the phenomenological definition of Type 2 excitability requires tonic high-frequency firing above threshold. However, the injected neurons were rarely able to sustain firing for the entire duration of the stimulus, leading to a pattern called stuttering (Izhikevich 2010). Nevertheless, the consistent existence of a subthreshold spikelet suggests Type 2 excitability.

#### A computational model.

To understand the dynamics of the K_{v} channel-modified neurons, we found it helpful to analyze the phase space geometry of realistic models that capture the features of the neurons. We began by constructing a four-dimensional model comprising Na_{v} channels of the Na_{v}1.6 type (Mercer et al. 2007), slow K_{v} channels of the injected type, and ohmic leak channels. This “full” model was described in classical Hodgkin-Huxley terms, comprising the potential *V* and the activation/inactivation parameters *m*, *h* (Na current), and *n* (K current):

Important elements in a geometrical analysis include orbits and equilibrium points, stable or unstable, and nullclines. This approach is most direct in two dimensions, partly for the obvious reason that two dimensions easily project onto paper, partly because two dimensions exclude chaotic solutions (cf. Poincaré-Bendixson theorem), and partly because the eigenvalues of the stability matrix are given directly by a quadratic equation. By the method described below we obtained a two-dimensional model of the following form:

To reduce the dimensionality, we took advantage of the fact that Na channel activation is faster than inactivation, which is faster than K channel activation, i.e., that

where τ_{i} denotes the time constant for Na channel activation and inactivation (*m*,*h*) and for K channel activation (*n*). To obtain *m*, we used the conventional assumption that *m* can be taken to be instantaneous, i.e., *m* = *m*_{∞}(*V*). To obtain *h*, we constructed an algorithm based on the time evolution of *m* and *n* and the steady-state curves of the *m*, *n*, and *h* parameters, i.e., *m*_{∞}(*V*), *n*_{∞}(*V*), and *h*_{∞}(*V*). The often-used assumption that *n* + *h* is constant did not apply to the present case, in contrast to the squid axon model; see FitzHugh (1961). Instead, we explored a way of estimating an effective “inactivation voltage” *V*_{h} by combining the voltages that are predicted by assuming that *m* and *n* are at their steady-state values (*V*_{m} and *V*_{n}). Using inverse functions we map *m* and *n* to potentials (*V*_{m} and *V*_{n}) that are then weighted using a function *f* to obtain a potential (*V*_{h}), which is used as an argument in *h*_{∞}(*V*). Thus we let

where the superscript −1 denotes the inverse function. Figure 3 depicts a function diagram of the reduction procedure. It is preferable if *f* is the identity function when the two arguments are the same, because then the fixed points are preserved. This preservation property can be shown as follows (Fig. 3). Let *V*_{s} be a stationary potential, and we have the following:

which shows that the fixed points are preserved. More attractive features of this reduced model are its experimentally interpretable parameters and the preservation of the steady-state curves from the full four-dimensional system. We are left to specify *f*. We chose the following function:

which, as we demonstrate, is surprisingly effective at reproducing the experimental spiking patterns (and setting *V*_{m} and *V*_{n} to *V*_{s} yields *V*_{h} = *V*_{s}). ρ is a constant between 0 and 1 and gives the relative speed of the kinetics (ρ = 1 represents the case of a sodium channel inactivation as fast as the activation, and ρ = 0 represents the case of an inactivation as slow as the potassium channel activation). We found ρ = 0.2 to be a reasonable value. To mimic a stuttering spike pattern, we added noise (see methods).

Figure 2, *C* and *D*, and Figure 4, *C* and *D*, show computed spiking patterns for models emulating noninjected and injected neurons. Details are given in the figure legends, and the full model is presented in methods. In conclusion, the reduced model presented here captures the features of the experimental recordings surprisingly well.

#### K_{v} channel effects on action potential.

The characteristics of the two spiking patterns at threshold are illustrated in Fig. 4, which compares the emergence of single impulses in a pyramidal cell during control conditions (Fig. 4*A*) and during a slow K_{v} conductance injection (Fig. 4*B*). This highlights the slow AHP and long latency to first spike in the control case (Fig. 4*A*) and the faster and more pronounced AHP and shorter latency time in the test case (Fig. 4*B*). In the latter case, the characteristic stimulation-dependent graded spikelet is observed to emerge with an increasing stimulation current, suggesting that it does not have a well-defined spiking threshold. This feature has been associated with Type 2 excitability and is the consequence of the absence of all-or-none responses. From a mathematical point of view, it has been related to the absence of a saddle equilibrium (FitzHugh 1961). The apparent illusion of threshold dynamics and all-or-none responses in Type 2 models has led the phenomenon to be called a “quasithreshold.” Our dynamical model captures the experimental findings, as shown in Fig. 4, *C* and *D*.

#### Latency to first spike.

Figure 5 shows that the latency to first spike at threshold decreases with increased *ḡ*_{K}, a finding predicted when a SNIC bifurcation is replaced by a fold limit cycle/subcritical Andronov-Hopf bifurcation (Izhikevich 2010; Wang et al. 2013). The same figure shows the latency in our model, which correspond well with the experimental data. As can be inferred from Fig. 5, the latency as a function of *ḡ*_{K} is a continuum, although strictly mathematically, the bifurcation between dynamical types occurs at ∼1 μS. Because we have avoided making assumptions about the native *ḡ*_{K}, there is a difference between the model and the experiments. The gray line in Fig. 5 is a function of the absolute *ḡ*_{K}, unlike the data, which are a function of the added K conductance, Δ*ḡ*_{K}. However, because we most likely injected more than the native *ḡ*_{K} (estimated as 0.5–1 μS, see below), the agreement between the model and the data is reasonable and the model qualitatively captures the response to increased *ḡ*_{K}.

#### Relation between spiking frequency and stimulation current.

As repeatedly mentioned above, the dynamical properties of a neuron determine its excitability, i.e., its stimulus current-frequency relation. Type 1 excitability is defined by a continuous *f-I*_{stim} curve; a discontinuous *f-I*_{stim} curve implies Type 2 excitability. As noted above (Fig. 2), most injected neurons did not produce sustained spiking but a stuttering pattern with regular interruptions, as if produced by a dynamical system close to a stable fixed point.

To obtain a *f-I*_{stim} curve for the spiking frequency during the bursts of activity, we based our analysis on measurements of the ISIs. We applied a clustering filtering technique in which ISIs that did not belong to a continuum in the *f-I*_{stim} plane were filtered out, thereby circumventing the problem with the irregular spikeless intervals.

Figure 6 shows the *f-I*_{stim} relations at different K channel densities, given as normalized *ḡ*_{K}, using estimated areas based on capacitance measurements described in methods. The figure is based on >20,000 measured ISIs. The data were obtained from measurements at two temperature ranges (<30°C and >30°C; see legend of Fig. 6). No systematic difference between the two data sets was found (see *Influence of temperature* below). The resulting data demonstrate that the onset frequency increases with higher K conductance (*ḡ*_{K}) values, with a switch in the dynamical type at ∼1 μS. A shift to higher onset frequency (15–25 Hz) can be observed at *ḡ*_{K} values of ∼2–3 μS (Fig. 6).

#### Influence of temperature.

The frequency-stimulation relationship was largely temperature independent. This can be inferred from Fig. 6*B*, which depicts data obtained from experiments conducted at high temperature (>30°C, red circles) and low temperature (<30°C, blue circles). This is not unexpected. The various components of the dynamics of a neuron are known to have differential temperature dependencies (Frankenhaeuser and Moore 1963; Jonas 1989). The capacitance, the leak current, and the reversal potential are only weakly dependent on temperature (e.g., linear as in the Nernst equations), whereas the ion channel kinetics are highly temperature dependent (exponential). The system follows Kirchhoff's first law and can be described as a sum of currents (see the model in methods).

When more *ḡ*_{K} is injected, the K current as a whole (native plus injected current) becomes less temperature dependent because the injected current is computed based on experimentally determined kinetics at room temperature. The Na current has a Q_{10} of ∼2.5 (Frankenhaeuser and Moore 1963; Jonas 1989). However, given the separation of timescales, this has a minor effect on the system; the activation of the Na current is so fast that it can be taken as instantaneous at 22°C.

#### Subthreshold oscillations.

As shown in Figs. 2 and 4, the conductance-injected Type 2 pyramidal cells showed a relatively smooth time course between spikes, suggesting that the system is not governed by a strong focus. As the bifurcation analysis reveals, the equilibrium point, although it is a focus, is in the proximity of a stable node. However, at higher *ḡ*_{K} injected neurons can be classified as true resonators with subthreshold oscillations, suggesting that the equilibrium point has switched to a stable focus. We used Fourier analysis to quantify such subthreshold oscillations (Prescott et al. 2008) (see methods). Figure 7*A* shows such subthreshold oscillations in a neuron injected with a relatively high *ḡ*_{K}, and Fig. 7*B* shows its frequency spectrum. Figure 7*C* shows that the power in the ISI increases with injected *ḡ*_{K}, a phenomenon that is consistent with the hypothesis of an underlying Andronov-Hopf bifurcation.

#### Threshold and bifurcation diagrams.

The results can be summarized in a diagram relating the dynamics to the *ḡ*_{K}-*I*_{stim} plane. Such a diagram maps the regions with different forms of dynamics based on the configuration and number of equilibrium points, defining the excitability or type of resting to spiking bifurcations. Figure 8 demonstrates this by detailing how the dynamics of the studied pyramidal neurons depends on the *ḡ*_{K} level. At zero injection and up to ∼1 μS, the excitability depends on a classical SNIC bifurcation. The Bogdanov-Takens bifurcation occurs in the middle of the three-solution region (∼1 μS), and subsequently there is a large span between 1 and 2 μS where the loss of stability is due to subcritical Andronov-Hopf bifurcation around three equilibria (a rarity in our previous modeling). From 2 to 4 μS of injection, the excitability depends on a classical fold limit cycle/subcritical Andronov-Hopf bifurcation in the one equilibrium point region.

#### Phase-plane analysis.

To more fully understand the dynamics of the studied pyramidal neurons, we analyzed the geometry of the computational model (described in methods) in the phase plane. The model demonstrates that it is possible to shift from a SNIC to a subcritical Andronov-Hopf bifurcation with the addition of K_{v} channels. Figure 9*A* shows the phase portrait of Fig. 2*C*, and Fig. 9*B* shows the phase portrait of Fig. 2*D*. The analysis is in good agreement with our previous theoretical study (Zeberg et al. 2010), and the detailed border region that we described proved to be of great relevance in explaining the present experimental results. The geometry of excitability was found to be complex, with global bifurcations around three equilibrium points (Fig. 8*B*). Two instructive examples (with associated eigenvalues) are presented in Fig. 9.

#### Effects of leak conductance.

Are the above investigated mechanisms for a shift in threshold dynamics specific for voltage-gated channels? Or can similar changes be achieved by leak channels? Figure 10*A* investigates this in the *g*_{leak}-*ḡ*_{Kv} plane. As can be seen, there are no qualitative differences: shifts can occur both through an increase in leakage current and through an increase of voltage-gated channel current. However, there are some quantitative differences. The magnitude differs, and the boundaries between the different regions are nonlinear. *Region I* represents a nonfiring regime, *region II* represents an area where a saddle-node bifurcation takes place, *region III* represents where oscillations starts with limit cycles around three fixed points, and *region IV* represents where oscillations arise via a subcritical Andronov-Hopf bifurcation.

Given that a rectification of the *I–V* curve prevents a saddle-node bifurcation, it is natural to assume that rectification is the mechanism behind the conductance-induced alterations. However, our model suggests otherwise: a change occurs before rectification is achieved. As can be inferred from Fig. 10*A*, *region III* that separates the saddle-node region (*II*) from the Andronov-Hopf region (*IV*) is substantial. A common scenario in the proposed model is thus oscillation around three fixed points. In this scenario the oscillations start with a certain starting frequency. To obtain further support for this explanation we examined *I–V* curves both from our model (which has been fitted to experimental data in the present investigation) and from dynamic *I–V* recordings on layer 5 pyramidal neurons conducted by Badel et al. (2008) (see Fig. 10, *B* and *C*). The derivate of the curve in Fig. 10*B* reaches around −100 nS, whereas in Fig. 10*C* it reaches around −1 μS. Therefore an addition of a leak current on the order of 10 nS is 1–2 orders of magnitude too small to straighten the *I–V* relationship. The dashed lines in Fig. 10, *B* and *C*, show the change induced by adding 1-μS slow voltage-gated K current, and the dotted lines show the effect of an extra 10-nS leak current with a reversal potential of −70 mV (as used by Prescott et al. 2008). The addition of 10-nS leak current causes a 45% reduction in input resistance around threshold, very much in line with the experimental value of 56% obtained by others (Prescott et al. 2008). In summary, neither the addition of a 1-μS voltage-gated K current nor the addition of 10-nS leak current is sufficient to straighten the *I–V* relationship.

Loss of stability around three fixed points has earlier been studied by us (Zeberg et al. 2010) in the context of a neuronal model based on voltage-clamp data on cultured hippocampal cells and axonal data. In that analysis this bifurcation has a minor role, whereas it plays a central role in the present study. The mathematics behind this complex transition requires an analysis of both the local properties (stability of the fixed points) and global properties (orientation of vector field and its trajectories). The type of bifurcation that occurs has been referred to as a (subcritical) big-saddle homo-clinic bifurcation (Izhikevich 2010). The appendix addresses this bifurcation in more detail for the interested reader.

## DISCUSSION

In previous theoretical work (Arhem et al. 2006; Arhem and Blomberg 2007; Zeberg et al. 2010), we showed that the membrane density of K_{v} channels in the membrane is a bifurcation parameter, i.e., that increasing the K_{v} channel conductance constants in a dynamical neuron model transforms Type 1 excitability into Type 2 excitability. We could also, using dynamical systems theory as a theoretical framework, analyze the underlying geometry (i.e., the vector field).

The present investigation is an experimental extension of the previous study. We confirmed that the number of naturally occurring slow K_{v} channels in the membrane soma of pyramidal cells was a bifurcation parameter and that we could switch the excitability from integrator to resonator dynamics by injecting synthetic slow K_{v} channel currents with the dynamic clamp technique. By constructing a two-dimensional dynamical model, we could also analyze the underlying geometry. Basically, we recognized three *ḡ*_{K}-dependent stages; uninjected neurons showed Type 1 excitability and SNIC bifurcations. Neurons injected with slow K_{v} channel currents at levels higher than ∼2 μS showed Type 2, and a phase plane with a single equilibrium point underwent typical fold/subcritical Andronov-Hopf bifurcations. At intermediate levels of slow K_{v} channel current injection (*ḡ*_{K} below 1–1.5 μS), however, the neurons showed more complex dynamics, caused by subcritical fold/Andronov-Hopf bifurcations in a phase plane with three equilibrium points.

Our results suggest that the dynamics of cortical pyramidal neurons are different from other previously investigated neurons. In a previous study (Zeberg et al. 2010), we showed that hippocampal neurons display relatively clear *ḡ*_{K}-dependent Type 1 and Type 2 excitability but very little intermediate dynamics. The Hodgkin-Huxley and Frankenhaeuser-Huxley axons show exclusively Type 2 excitability. We speculate that this difference has functional and adaptive consequences, restricting axonal dynamics and resulting in quantitative and not qualitative changes under stress (Fahlstrom et al. 2012; Shi et al. 2013). A previous study (Prescott et al. 2008) studied the effects of adding a leak conductance of 10 nS to CA1 pyramidal neurons. They also studied adaptation, using prolonged depolarizing current injections to activate M-type K currents and calcium-activated K currents. Both shunting and adaptation could cause a shift from Type 1 to Type 2. Their explanation is that of a rectified *I–V* curve in both cases. However, we believe (see Fig. 10) that the presently offered explanation also explains, at least to some extent, their results.

#### Channel-based description of the dynamics switch.

An intuitive understanding of the dynamical switch in terms of activity of the channel populations is as follows. The added K_{v} channel conductance is partially activated at subthreshold voltages, so that it partially counteracts the effect of the current injection that stimulates the neuron. This means that the current threshold for spiking increases substantially. It also means that for just subthreshold current injections, the voltage increases to a peak (graded with stimulus current) and then relaxes back to more negative voltages even with continuing current injection. Spike threshold is only reached for substantially larger currents than in noninjected neurons, and once threshold is reached the combination of stimulating current and Na current (which may recover from inactivation faster than in control due to more negative AHPs resulting from spike-evoked K_{v} channel conductance) leads to shorter ISIs and faster firing frequency at threshold, the frequency increase depending on how much K_{v} channel conductance was injected (Fig. 1 in Zeberg et al. 2010). A key element in the shift of spiking behavior is at which voltage the added channel starts to open. The activation curve of the slow delayed-rectifier channel, the current of which was injected in the present experiments, is shown in Fig. 1. Somewhat surprisingly, the midpoint of the activation curve is approximately −10 mV, although the dominant delayed-rectifier channel in somata and proximal dendrites of pyramidal neurons in the cortex and hippocampus has been reported to be of the K_{v}2.1 type (Misonou et al. 2005b), the activation curve midpoint of which being more positive than −10 mV (e.g., see Nilsson et al. 2003). This suggests that other K_{v} channels contribute significantly to the K current of the soma (possibly K_{v}1 and K_{v}3 channels). The fact that rather large conductances were needed in the study most likely reflects the low opening probability around threshold (Du et al. 2000; Guan et al. 2013; Johnston et al. 2008). Both data and modeling, however, suggest that this opening probability is enough to provide a significant current for a shift in dynamics.

#### Principles of conductance-based shifts of threshold dynamics.

As shown above, dynamical system theory provides us with a biophysical and theoretical framework for understanding altered dynamics in a wider sense and is not limited to a specific ion channel. Certain channels have gained particular interest, e.g., K_{v}3 channels in fast-spiking interneurons (Erisir et al. 1999; Lien and Jonas 2003) and K_{v}7 channels in pyramidal cells (Prescott et al. 2008; Stiefel 2003). The current injected in this work is not as slowly activating as the K_{v}7 channel (which was injected as a shunt in Prescott et al. 2008) and not as fast as the K_{v}3 ion channel.

The following formal argument clarifies why an increased K_{v} channel density eventually leads to Type 2 dynamics. Let *f*(*V*) be the sum of the steady-state currents. Any equilibrium point will involve a potential that solves the equation

where *I*_{stim} is the applied current and *I*_{i} is the ionic currents. A saddle-node bifurcation requires three equilibria and therefore three roots to *f*(*V*) = 0. If the function is strictly monotone, we know that the zero is unique (*f* is said to be injective) and that a saddle-node bifurcation is not possible. We can now see why the magnitude of the conductance is so important. For *f*(*V*) to be strictly increasing, we have

where *g*_{out} and *g*_{in} are slope conductances associated with outward and inward currents (i.e., d*I*/d*V*), respectively. The outward conductances are positive (if they open with potential and we define the current as positive for an outward cation flux) and inward conductances are consequently negative, and each conductance is proportional to its ion channel density. Therefore, increasing any noninactivating K current (usually an outward current) relative to the inward currents eventually prevents a saddle-node bifurcation (or in our case, more precisely a SNIC) and regular spiking (Type 1 excitability). Moreover, the steady-state currents are independent of the time constants. Therefore, a sufficient criterion for a K current to inhibit regular spiking at a sufficiently large conductance is that it opens with potential. Consistent with this argument, A-type K current (K_{v}A) can have the reverse effect (Ermentrout and Terman 2010).

However, as demonstrated by our model (Fig. 8) and our previous work (Zeberg et al. 2010), a shift from Type 1 to Type 2 can occur, even with an N-shaped *I*_{ss}-*V* curve and three equilibrium points. A look at the eigenvalues provides an insight into this mechanism, which is different from that of Prescott et al. (2006, 2008). For an Andronov-Hopf bifurcation to occur, the eigenvalues need to be complex, which is the case in our two-dimensional model if

where *g* is the slope conductance, τ_{n} the time constant of the voltage-dependent K current, *C* the membrane capacitance, *I* the sum of the ion currents, *V* the membrane potential, and *n* the opening probability of the voltage-dependent K current. Adding a leak term will affect the slope conductance *g*, whereas an increase in an active K current will also affect ∂*I*/∂*n*. A more thorough mathematical analysis of the different effects on the stability matrix of leak and active conductances is beyond the scope of the present study.

Thus it can be observed that this, in principle, allows similar behavior (Type 2 dynamics) to be achieved in numerous ways, essentially by increasing any K_{v} conductance. This might be important for the robustness and homeostasis of this behavior (Marder and Goaillard 2006). Moreover, the dynamics of the individual neurons matter because the population-level coding is closely connected to the prevailing operating mode of the constituent neurons of the network (Cardin et al. 2009; Hong et al. 2012; Sohal et al. 2009).

#### Magnitude of injected conductances.

An important methodological question is whether the injected currents used in this study are within physiological limits. Estimating conductance levels from Hodgkin-Huxley modeling in the somata of pyramidal cells has yielded values between 0.032 μS and 1.2 μS for Na_{v} channels and between 0.24 μS and 0.76 μS for K_{v} channels (Table 2). In experiment-based Hodgkin-Huxley models of cortical pyramidal cells, different values for the somatic density of the conductances have been used (Baranauskas and Martina 2006; Colbert and Pan 2002; Mainen and Sejnowski 1996; McCormick et al. 2007). The behavior of these models is within the range of real cortical pyramidal cells (McCormick et al. 2007). Table 3 summarizes a few estimates. The conductance is many times higher if the contribution of the proximal dendrites is considered. In our study, the average capacitance was found to be 111 pF (methods), corresponding to a membrane area of 11,100 μm^{2}, assuming the often-used value for membrane capacitance of 1 μF/cm^{2} (Hodgkin and Huxley 1952). This area is approximately seven times larger than the soma and a fifth of the whole dendritic area (Mainen et al. 1995). Assuming a K_{v} conductance density of 0.15 nS/μm^{2} (McCormick et al. 2007), the measured membrane area would have a total K conductance of ∼1.5 μS. However, because the density is considerably higher in the soma than in the dendrites (Korngreen and Sakmann 2000), the total *ḡ*_{K} is most likely ∼0.5–1 μS. Therefore, the shift from Type 1 to Type 2 dynamics observed in our study (∼1 μS) lies in the upper physiological range.

#### Conclusions.

In conclusion, we have confirmed the channel density hypothesis experimentally. The spiking pattern of neurons can be switched between Type 1 and Type 2 by up- or downregulating channel densities. Furthermore, the studies suggest that the dynamics of neocortical pyramidal cells differ in some respects from previously studied neurons. Our modeling suggests that an intermediate-level increase of slow K_{v} channel density transforms the dynamics into a specific form of Type 2 dynamics, characterized mathematically by Andronov-Hopf bifurcations in a three-equilibrium point phase region not known for other neurons. The functional relevance of this is unknown.

The possibility of modifying the threshold dynamics of neurons by altering channel densities suggests novel scenarios for the action of channel-active drugs such as general anesthetics, implying mechanisms where selective blocking ion channels in critical neurons induces a switch from one brain state (e.g., associated with wakefulness) characterized by certain frequency patterns to another state (e.g., associated with general anesthesia) characterized by other frequency patterns. Network modeling has shown that such ideas are feasible (Halnes et al. 2007).

The possibility of channel density-induced modifications of neuronal dynamics also suggests novel pathophysiological scenarios. K_{v}2.1 channel density is subject to modulation by changes in neural activity (Misonou et al. 2004), by ischemia (Misonou et al. 2005a), and by muscarinic receptor activation (Mohapatra and Trimmer 2006). The mechanisms of synchronization at the network level are still not well understood, but the mechanisms at a cellular level have been extensively studied and a tight connection between the bifurcation structure and the phase response curve has been established (Gouwens et al. 2010; Zeberg et al. 2012). Type 2 neurons synchronize more easily than Type 1 neurons (Izhikevich 2010). Type 2 interneurons have recently been shown to account for the cortical gamma oscillations (20–80 Hz) (Cardin et al. 2009), which are considered to provide a temporal structure for information processing in the brain. The possibility of increased synchrony in brain networks due to a switch of integrators into resonators because of upregulated K_{v} channels may be relevant in the explanation of some forms of epilepsy.

## GRANTS

This work was supported by grants from Swedish Research Council (VR) Grant No. 15083, the Karolinska Institutet Faculty funds, the Biotechnology and Biological Sciences Research Council (BBSRC), the European Commission, and the Swedish Society of Medicine. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: H.Z., H.P.C.R., and P.Å. conception and design of research; H.Z. and H.P.C.R. performed experiments; H.Z., H.P.C.R., and P.Å. analyzed data; H.Z., H.P.C.R., and P.Å. interpreted results of experiments; H.Z. prepared figures; H.Z. and P.Å. drafted manuscript; H.Z., H.P.C.R., and P.Å. edited and revised manuscript; H.Z., H.P.C.R., and P.Å. approved final version of manuscript.

## Appendix

The existence of a limit cycle in a constrained region of the plane with an unstable fixed point is ensured by the Poincaré-Bendixon theorem (Bendixon 1901); since trajectories from the fixed points must remain on the plane, they will eventually approach a limit cycle. With more fixed points the analysis becomes more complex. There are several possibilities with three fixed points (one by necessity being a saddle). A common scenario is that one acts as an attractor and the other as a repeller, and when the attractor merges with the saddle the system is left oscillating around the repeller. A less well-studied scenario is oscillations around all three fixed points. This is allowed by index theory (Guckenheimer and Holmes 1983)—a limit cycle must always encapsulate a value of +1; an attractor/repeller has the value of +1 and a saddle the value of −1, bringing the sum to +1. Oscillations around two fixed points are therefore forbidden. An example of a limit cycle around three fixed points is depicted in Fig. A1.

The limit cycle is separated by an inner separatrix in the form of a homo-clinic orbit including the saddle, hence the name “big saddle homo-clinic bifurcation” (a minor variant is also possible). There are three fixed points in the plane represented by dots: a stable focus (λ = *a* ± *bi*), a saddle (λ_{1} = −*a*_{1},λ_{2} = *a*_{2}), and an unstable node (λ_{1} = *a*_{1},λ_{2} = *a*_{2}). There is a “ridge” between the saddle and the unstable node, separating trajectories going to the left and to the right. However, both sides of the ridge are included in the attractor domain of the stable focus defined by the homo-clinic orbit. Trajectories going to the stable focus originate from the eigenvectors of the two other fixed points, and the lines crossing the phase plane are the nullclines. Before the bifurcation takes place, the stable focus is a global attractor.

- Copyright © 2015 the American Physiological Society