Saddle-node loop bifurcations ubiquitously occur in conductance-based model neurons and promote spike-based coding and synchronization

Information processing in the brain crucially depends on coding properties of single neurons shaped by their intrinsic dynamics. Our theoretical study shows that neuronal excitability, spike-based coding, and synchronization change drastically around a special dynamical regime: the codimension-two saddle-node loop bifurcation. This bifurcation maximizes the entrainment range and turns out to be more relevant for spike-based processing than the typically considered Bogdanov-Takens point. We show that the saddle-node loop bifurcation occurs ubiquitously in a generic class of neuron models (planar type-I). Our theory, applied to optical stimulation techniques, reveals that they could manipulate a variety of information processing characteristics in nerve cells beyond pure stimulation.

The dynamical properties of a neuron shape the way its input is encoded into a series of spike times [1][2][3]. While the precise dynamics depend on a multitude of cellintrinsic parameters, the qualitative behavior of a neuron is determined by the resulting bifurcation structure [4]. Here, we concentrate on the bifurcation at spike onset and discuss the saddle-node loop (SNL) bifurcation with its far-reaching implications for spike-based coding and synchronization.
As codimension-two bifurcation, the SNL bifurcation is reached by tuning two system parameters. We focus on input and capacitance, which are both fundamental parameters in single-neuron models [5], and are, for example, relevant for the experimental technique of infrared neural stimulation [6]. This technique excites neurons by a membrane current arising from a transient increase in capacitance [7]. We show, however, that changes in capacitance could also push the neurons closer to an SNL bifurcation. Our results demonstrate that, in this case, a single neuron's coding efficiency and its ability to synchronize in a network are significantly altered.
A capacitance increase drives neurons towards an SNL bifurcation. In conductance-based neuron models, the dynamics of the membrane voltage v are modeled as a balance of input, ionic, and capacitive currents, where C m is the membrane capacitance. The termĊ m v allows for the transient depolarizing current during infrared neural stimulation [7]. Beyond these transient dynamics, we show that in a relaxed state, i.e.,Ċ m = 0, C m has a rich impact on neuronal dynamics. Common neuron models show two basic behaviors: They are either at rest or spike repetitively. For low I in , such a model relaxes to a fixed point (the resting state), while it shows limit cycle dynamics (spiking) for high I in . In our analysis, we consider models in which the limit cycle is created through a saddle-node on an invariant cycle (SNIC) bifurcation occurring at some I in = I SN1 . For example, in the Wang-Buzsaki model [8] with I in = I SN1 , we observe that an increase in C m leads to an SNL bifurcation (later denoted small SNL) as illustrated in Fig. 1(a). Generally, this bifurcation replaces the invariant cycle of the saddle-node bifurcation by the loop of a saddle homoclinic orbit (HOM) bifurcation [9]. For both SNIC and SNL bifurcation, the limit cycle arises from a loop connected to a saddle-node fixed point. For the SNIC bifurcation, the orbit passes (infinitely slowly) through the saddle-node fixed point along its slow, semi-stable manifold. At the SNL bifurcation, the orbit is deformed: While the exit of the saddle-node still follows the semi-stable manifold, the approach of the saddle-node switches to the faster, strongly-stable manifold, see Fig. 1(a). This switch accelerates the dynamics along the approach, while the time required for the exit remains unchanged, resulting in a lower period once the limit cycle is detached, exemplified for the Wang-Buzsaki model in Fig. 2(a), see also [10].
The SNL bifurcation is accompanied by a change in PRC symmetry. For the SNIC bifurcation, the velocities of entry and exit dynamics at the saddle-node fixed point are similar. This symmetry is broken at the SNL bifurcation. As a consequence, the phase response curve (PRC) and spike-based coding properties (because of their sensitivity to PRC symmetry) are changed, as we outline in the following.
In order to relate neuronal dynamics to spike-based coding, we introduce a phase variable [11], φ ∈ R, where the spike times (maximum voltage) are mapped to φ(t spike ) = k for k ∈ Z. In absence of perturbations, the phase evolves along the limit cycle with period T aṡ φ = 1/T . The PRC, Z, is a periodic function that maps the phase of occurrence of a small perturbation to the phase shift of the next spike time, Z : φ → ∆φ. For convenience, we refer to the PRC of a limit cycle that detaches from a bifurcation (SNIC or SNL) as the PRC at this bifurcation.
The PRC of the Wang-Buzsaki model at a nondegenerated SNIC bifurcation has a (reflection) symmetric shape shown in Fig. 1(b), middle panel, comparable to the canonical PRC derived by normal form theory, Z SNIC (φ) = 1 − cos(2πφ) [12]. Towards the SNL bifurcation, the PRC develops an asymmetric component as shown in Fig. 1(b), top panel, accompanied by an increase in Fourier components.
This PRC asymmetry is a direct consequence of the broken symmetry in the dynamics at the SNL bifurcation. Because limit cycle dynamics are maximally slowed when passing through the semi-stable manifold near the saddle-node (more precisely, the ghost of the former saddle-node), a perturbation that propels the dynamics over the saddle-node will maximally advance the next spike, i.e., the maximum of the PRC aligns with the phase of the saddle-node passage. For a SNIC with its symmetric dynamics, the saddle-node is reached at halfperiod, i.e., the PRC maximum lies at φ = 0.5. In comparison, the PRC maximum is shifted to the left at the SNL bifurcation, because the faster approach along the strongly-stable manifold advances the saddle-node to earlier phases. The shift of the maximum away from the center destroys the symmetric shape of the PRC. More precisely, while the SNIC flow approaches the saddle-node along the semi-stable manifold, the SNL flow is almost instantaneously injected at the saddle-node: Subsequent to the saddle-node, the SNL flow aligns with the second half of the SNIC flow. Thus, also the PRC at the SNL bifurcation should align with the second part of the PRC at the SNIC bifurcation, Z SNL (φ) ≈ c Z SNIC (φ/2 + 0.5), with some constant c. This intuition is supported by numerical continuation ( Fig. 1(b)).
Mathematically, the SNIC bifurcation normal form with its parabolic dynamics,ẋ = I in +x 2 results in a point symmetric flow, x(t) = −x(−t), and a point symmetric Jacobian, J(x) = −J(−x). Inserting these two quantities into the adjoint equation,Ż = J T Z, directly leads to a reflection symmetric PRC, Z SNIC (t) = Z SNIC (−t). In contrast, the asymmetric dynamics at the SNL bifurcation lead to an asymmetric PRC.
The symmetry breaking generalizes beyond the SNL bifurcation: For a normal saddle with different stabilities along stable and unstable manifold, a homoclinic orbit shows an asymmetric PRC, compare [12]. Because the symmetry breaking in the dynamics immediately follows from the orbit flip bifurcation (as defined in [13]) at a saddle-node fixed point [9,13], the consequential symmetry breaking in the PRC is a general property of the SNL bifurcation.
Consequences of passing the SNL bifurcation for spike-based coding and synchronization. The symmetry breaking in the PRC has important consequences for spike-based coding and synchronization summarized in Fig. 2. Neurons close to the SNL point behave radically different from what is expected for SNIC type-I neurons: (i) The model's linear response function passes considerably higher frequencies at the SNL point, see Fig. 1(c). The increased cutoff frequency of the filter results from the larger amount of high Fourier components in the PRC, compared to the single Fourier component of Z SNIC (and also compared to the canonical type-II PRC). (ii) As a consequence, the lower bound of the information rate transmitted per spike about a time-varying stimulus is increased at the SNL point, as shown in Fig. 2(c) [3]. (iii) Models around the SNL point show enhanced reliability, i.e., the locking to a time-varying white noise stimulus is reached earlier. This manifests in a smaller Lyapunov exponent of the locking state ( Fig. 2(b)), given by λ = − 1 0 (Z ′ (φ)) 2 dφ, where ′ denotes the derivative with respect to φ [14][15][16][17]. (iv) In an all-to-all coupled Kuramoto network with delta synapses, the frequency range for which the population still synchronizes peaks around the SNL point, compare Fig. 2(b) [11].
Together, this suggests that passing the SNL point could affect a variety of neural coding schemes: Regarding spike-timing codes, the shape of the PRC at the SNL bifurcation is similar to the optimal PRC for synchronization by common Poisson input [18]. Furthermore, the lower bound of the information rate is increased, because, at the SNL point, higher frequencies emerge in the linear response function. Regarding phase-codes [19], a consequence of the PRC at the SNL bifurcation is a temporally more precise locking of spikes to a reference signal. This prevents spurious spikes, i.e., fewer phase slips occur, reducing the decoding error. Regarding rate codes, at the SNL bifurcation, the lower bound of the decoding error (Cramer-Rao bound) is decreased, which is connected to the slope of the firing rate as a function of the input I in [20]. This slope is known to be larger for the spike onset via a HOM bifurcation compared to the spike onset via a SNIC bifurcation [4], and hence rises around the SNL point.
Generic occurrence of the SNL bifurcation. When systematically changing the capacitance in planar conductance-based models, two SNL bifurcations enclose the SNIC bifurcation, compare Fig. 3(a). Decreasing the capacitance, an SNL point is passed before the Bogdanov-Takens (BT) bifurcation is reached, at which the classical switch from type-I to type-II occurs [21].
To derive the bifurcation structure, we augment Eq. 1 with a second dimension, n, for the gating kinetics .
We prove under the following assumptions that this bifurcation structure is ubiquitous for planar conductancebased neuron models: (A1) We consider a generic class of type-I non-bursting conductance-based neuron models, for which at some capacitance value C SNIC , the spike onset occurs at a non-degenerated SNIC bifurcation at some I in = I SN1 . (A2) At C m = C SNIC , stable dynamics for I in < I SN1 are only given by a single stable fixed point, the unique resting state, and limit cycle dynamics terminate in a bifurcation denoted excitation block. (A3) The nullcline of the voltage has an inverted N -shape, and the nullcline of the second variable is strictly monotonically increasing with voltage. (A4) To fulfill the assumptions in Ref. [5], we require that n ∞ (v) is a positive, bounded, twice differentiable function that becomes sufficiently flat in the limit v → ±∞, lim v→±∞ v ∂ v n ∞ (v) = 0, which is fulfilled in all typical neuron models.
The location of the nullclines is given by I cap (v, n) = 0 and n ∞ (v) − n = 0; nullclines and fixed points are hence independent of C m . The voltage nullcline is shifted upwards by an increase in input I in . Depending on I in , we observe the following fixed points: For low I in , the model has a single, stable fixed point, which we call P rest . A fold bifurcation occurs at some I in = I SN0 , at which a saddle, P saddle , and a node, P block , are born. P block is unstable so that P rest can be the only stable fixed point at C m = C SNIC as requested by (A2). In a second fold bifurcation at I in = I SN1 , P rest and P saddle collide.
Lemma 1: A variation in C m generically breaks homoclinic orbits to hyperbolic fixed points. A variation in C m breaks homoclinic orbits if the corresponding Melnikov integral, M , is non-zero [13]. Following the presentation in Equation (6.12) in [22], we evaluate the Melnikov integral, M , with respect to C m for a homoclinic orbit, given by its flow h(t), as where K(t) = exp − t 0 div F (h(s)) ds . With ∀t : 0 < K(t), and, to observe spiking, ∃t : 0 < I 2 cap (h(t)), the Melnikov integral is strictly positive, 0 < M . With that, following a HOM branch in one direction, the corresponding input values will either continuously increase or decrease, i.e., the branch in the bifurcation diagram cannot "fold back" in I in .
Lemma 2: A HOM branch is stable when it connects to a non-degenerated fold bifurcation of a stable node and a saddle. The collision of a stable node and a saddle results in a saddle-node fixed point with a zero and a negative eigenvalue. The sum of both eigenvalues results in a negative saddle-quantity for the saddle-node, and, thus, the homoclinic orbit arising from this saddle-node fixed point is stable.
Observation 1: Anchoring of the bifurcation diagram in the limit C m → 0 yields I big HOM < I Hopf . In this limit, the conductance-based model is transformed into a relaxation oscillator with voltage as fast variable, as sketched in Fig. 3(c) [23]. For sufficiently small C m , an increase in I in results in a sequence of bifurcations given in Theorem 2 of Ref. [24]. Relevant for our consideration is I big HOM < I Hopf < I SN1 , where a big HOM bifurcation occurs at I big HOM and a subcritical Hopf bifurcation destabilizes P rest at I Hopf .
Observation 2: The big HOM branch continues in the direction of increasing I in . As stated in Ref. [24], the big HOM branch, I big HOM (C m ), is an increasing function of C m for sufficiently small C m , and, with Lemma 1, the big HOM branch continues in the direction of increasing I in for larger C m .
Lemma 3: For I in = I SN1 , increasing C m from zero passes first a BT point, and then an SNL point, before a non-degenerated SNIC bifurcation occurs, C BT1 < C big SNL < C SNIC . With Observation 2 and the fact that in a planar system limit cycle branches cannot cross each other, the ordering in I in immediately implies an ordering in C m . In particular, Observation 2 implies that the big HOM branch connects as stable limit cycle (Lemma 2) at some C big SNL to the fold branch at I SN1 . We identify the point (I SN1 , C big SNL ) as big SNL bifurcation. Ref. [24] ensures that the Hopf branch does not return to the limit C m → 0, and as it cannot cross the big HOM branch it connects at some C BT1 to the fold branch at I SN1 in a BT point. Furthermore, because no stable orbit exists by (A2) for I in < I SN1 at C m = C SNIC , we conclude C big SNL < C SNIC .
These arguments have proven the bifurcation sequence in the lower part of the bifurcation diagram founded on the limit C m → 0. In the following, we show that the upper part of the bifurcation diagram results from the unfolding of a second BT point.
Observation 3: The bifurcation diagram contains exactly two BT points. Our system shows two fold bifurcations, which are independent of C m , and for each of those we find exactly one BT point, one at (I SN1 , C BT1 ), the other at (I SN0 , C BT0 ) [5]. In particular, this ensures that no second BT point occurs for I in = I SN1 .
Observation 4: The BT point at (I SN0 , C BT0 ) is related to the excitation block and occurs at C BT0 > C SNIC . The Hopf bifurcation that arises from the BT point at (I SN0 , C BT0 ) stabilizes P block and must hence be related to the excitation block that occurs by (A2) after the SNIC bifurcation with a sufficient increase in input. Observe that the branch of the Hopf bifurcation can neither pass the big HOM branch identified before, nor the branch of SNIC bifurcations that occurs between C big SNL and C SNIC . Furthermore, fixing C m to C SNIC , (A2) implies that the stabilization of P block by the Hopf bifurcation cannot pass the line spanned by C = C SNIC for I in < I SN1 . Thus, C BT0 > C SNIC .
Lemma 4: A second SNL bifurcation occurs at some C small SNL > C SNIC . From the BT point at (I SN0 , C BT0 ) arises by normal form theory a branch of a small HOM bifurcation in the direction of increasing input I in , and, with Lemma 1, it continues in this direction, in analogy with Observation 2. Hence, we find some C m = C small SNL for which the small HOM branch, stable due to Lemma 2, connects to the fold bifurcation at I SN1 . We identify the point (I SN1 , C small SNL ) as a small SNL bifurcation. From the BT point at (I SN0 , C BT0 ), a Hopf branch arises next to the small HOM, and a limit cycle exists between both. As long as the small HOM bifurcation exists, the Hopf bifurcation cannot terminate the SNIC bifurcation as demanded by (A2). Consequently, the SNL point occurs at some C small SNL > C SNIC .
In summary, we have shown that C BT1 < C big SNL < C SNIC < C small SNL . This generic bifurcation structure occurs with the membrane capacitance C m as bifurcation parameter at I in = I SN1 . For a model starting at a SNIC bifurcation, a variation in the capacitance will thus pass an SNL bifurcation before a BT point is reached.
Implications of the SNL bifurcation. We have shown that the SNL bifurcation entails interesting changes in spike-based coding properties and synchronization and that, for planar conductance-based models with a SNIC bifurcation, an SNL bifurcation is the first bifurcation reached for lowered or increased capacitance, respectively. With the three bifurcation parameters capacitance, input, and leak conductance, the identified sequence of bifurcations collapses into a codimensionthree cusp BT point [5]. This potentially generalizes the described bifurcation structure beyond the planar case. While our mathematical arguments are based on the membrane capacitance, the SNL bifurcation is also reached by tuning other parameters, including maximal gating conductances, tonic inhibition [5], neuromodulators [25], and gating time constants.
Drastic changes in spike-based coding can only be expected at a bifurcation that affects not only the fixed points, but also the stable limit cycle, such as the SNL bifurcation. This distinguishes the SNL bifurcation from the BT point where the transition of type-I to type-II dynamics occurs. The BT-associated Hopf bifurcation in relevant neuron models is typically subcritical, such that the arising limit cycle is unstable. This Hopf bifurcation hence affects the resting fixed point and changes subthreshold dynamics and filtering [5], but has only mild implications for spiking dynamics. In particular, the negative component of the PRC at the subcritical Hopf bifurcation is not derived from the subcritical Hopf itself, but likely reflects a combination of the PRCs at the limit cycle bifurcations responsible for limit cycle creation and excitation block.
On the practical side, our results demand for a careful interpretation of infrared neural stimulation results. For a neuron close to an SNL bifurcation, potentially, the effect of the capacitance increase caused by infrared neural stimulation goes beyond the reported transient capacitive current. We showed that the capacitance increase by itself can alter neuronal excitability, spike-based coding, and synchronization.
Passing the SNL point causes a region of bistability of stable limit cycle and resting state P rest , see Fig. 3(ab). The bistability leads to hysteresis in the response of stimulated neurons. With regard to optical stimulation, such hysteresis can result in spiking that persists beyond the duration of the transient current stimulation. The respective statistics of switching between the two stable states (i.e., between silent periods and periodic firing) and their dependence on the timing of excitation and inhibition may hence be characteristic of the different SNL types and allow for a distinction between these dynamical regimes, see also [4]. Altogether, we predict that infrared neuronal stimulation may result in unexpected "side effects", switching neurons into a significantly different dynamical regime.
Finally, we remark that the phase description, on which our results are partly based, is only valid for stimulus perturbations whose effect on the amplitude of the limit cycle decays sufficiently within one period. In our case, the relative decay rate, as given by the ratio of the magnitude of the Floquet exponent and the period ( Fig. 2(a)), is small throughout the capacitance interval enclosed by the SNL points. This means that the phase description in our region of interest is also valid for relatively large inputs [11].
In summary, the general occurrence of the SNL bifurcation next to a SNIC bifurcation, together with the derived symmetry changes, highlights the SNL bifurcation as a hitherto underestimated bifurcation with prominent importance for neuronal dynamics. The homeostatic regulation of intrinsic dynamics towards the SNL point may even constitute a useful strategy for those neurons that process particularly high frequencies.
Funded by the German Federal Ministry of Education and Research, Germany (01GQ0901, 01GQ1403).   (c) Schematic illustration for the limit Cm → 0, in which the system corresponds to a relaxation oscillator: The solid line with inverted N-shape represents the voltage nullcline, the dashed line represents the gating nullcline. At some Iin, a big HOM orbit around all fixed points is created and the resting state looses stability.