Hyperpolarization-Activated Current Induces Period-Doubling Cascades and Chaos in a Cold Thermoreceptor Model

In this article, we describe and analyze the chaotic behavior of a conductance-based neuronal bursting model. This is a model with a reduced number of variables, yet it retains biophysical plausibility. Inspired by the activity of cold thermoreceptors, the model contains a persistent Sodium current, a Calcium-activated Potassium current and a hyperpolarization-activated current (Ih) that drive a slow subthreshold oscillation. Driven by this oscillation, a fast subsystem (fast Sodium and Potassium currents) fires action potentials in a periodic fashion. Depending on the parameters, this model can generate a variety of firing patterns that includes bursting, regular tonic and polymodal firing. Here we show that the transitions between different firing patterns are often accompanied by a range of chaotic firing, as suggested by an irregular, non-periodic firing pattern. To confirm this, we measure the maximum Lyapunov exponent of the voltage trajectories, and the Lyapunov exponent and Lempel-Ziv's complexity of the ISI time series. The four-variable slow system (without spiking) also generates chaotic behavior, and bifurcation analysis shows that this is often originated by period doubling cascades. Either with or without spikes, chaos is no longer generated when the Ih is removed from the system. As the model is biologically plausible with biophysically meaningful parameters, we propose it as a useful tool to understand chaotic dynamics in neurons.

In this article, we describe and analyze the chaotic behavior of a conductance-based neuronal bursting model. This is a model with a reduced number of variables, yet it retains biophysical plausibility. Inspired by the activity of cold thermoreceptors, the model contains a persistent Sodium current, a Calcium-activated Potassium current and a hyperpolarization-activated current (I h ) that drive a slow subthreshold oscillation. Driven by this oscillation, a fast subsystem (fast Sodium and Potassium currents) fires action potentials in a periodic fashion. Depending on the parameters, this model can generate a variety of firing patterns that includes bursting, regular tonic and polymodal firing. Here we show that the transitions between different firing patterns are often accompanied by a range of chaotic firing, as suggested by an irregular, non-periodic firing pattern. To confirm this, we measure the maximum Lyapunov exponent of the voltage trajectories, and the Lyapunov exponent and Lempel-Ziv's complexity of the ISI time series. The four-variable slow system (without spiking) also generates chaotic behavior, and bifurcation analysis shows that this is often originated by period doubling cascades. Either with or without spikes, chaos is no longer generated when the I h is removed from the system. As the model is biologically plausible with biophysically meaningful parameters, we propose it as a useful tool to understand chaotic dynamics in neurons.

INTRODUCTION
Chaotic behavior in neural systems has been observed for many years. Experimental observations of non-periodic responses range from molluscan neurons (Aihara et al., 1984) to rat sciatic nerves (Gu, 2013), including lobster CPGs (Abarbanel et al., 1996) and fish's Mauthner cells (Faure et al., 2000) (for a review, see Korn and Faure, 2003). In addition, chaotic behavior has been analyzed in detail in several models of neural excitability. Notable examples include the Plant model for the R15 bursting cell in Aplysia (Plant and Kim, 1976) that shows a chaotic regime between bursting and beating firing modes (Canavier et al., 1990); the Chay model of pancreatic β cells (Chay and Rinzel, 1985); and the Huber & Braun (H&B) model of cold thermoreceptors (Braun et al., 1998;Feudel et al., 2000).
Most of the models reported to show chaotic activity present different types of bursting oscillations, that arise from the interaction between fast membrane voltage dynamics and a slower current or intracellular mechanism, such as Calcium dynamics (Chay and Rinzel, 1985;Canavier et al., 1990;Falcke et al., 2000). Because of this evidence, the interaction of different time scales has been proposed as the origin of the chaotic behavior. A simpler model, of only 3 variables, that produces burst firing and is known as the Hindmarsh-Rose model (Hindmarsh and Rose, 1984), presents chaotic dynamics for some parameter combinations. Despite its simplicity, it shows a variety of aperiodic behaviors that are being actively studied by several groups (Holden and Fan, 1992;Abarbanel et al., 1996;Barrio and Shilnikov, 2011;Barrio et al., 2014), giving useful insight into the intrinsic mathematical mechanisms that drive bursting dynamics. However, a drawback of the Hindmarsh-Rose model is that some of its equations and parameters lack actual biophysical meaning, and thus its usefulness to interpret biological data is somewhat limited.
Instead, here we report and analyze the chaotic behavior of a bursting model inspired on the temperature-dependent firing patterns observed in cold thermoreceptors (Braun et al., 1980). Our model is derived from the H&B model (Braun et al., 1998), in which a slow membrane oscillation is driven by a mixed Sodium/Calcium current and a Calcium-activated Potassium current. Fast Sodium and Potassium currents produce action potentials, and the usual effect of temperature on channel kinetics makes the model to display different spiking patterns such as bursting, tonic, skipping and chaotic. Recently, a hyperpolarization-activated current (I h ) was added to the model in order to reproduce and explain experimental results with genetic and pharmacological suppression of I h (Orio et al., 2012). This extended model will be referred here as HB+Ih and it is the main object of study in this paper. We present here a parameter sweeping approach (Barrio and Shilnikov, 2011;Barrio et al., 2014) to explore the regions of chaotic behavior and its dependence on certain model parameters.
The original H&B model shows chaotic behavior at very low temperatures (<10 • C), thus limiting the possible biological interpretations of this behavior. In contrast, we found that the new HB+Ih model displays chaotic behavior at physiological temperatures, namely in the 32-38 • C range. Moreover, the chaotic behavior is highly dependent on the presence of I h and its associated activation parameters; this is one of the most relevant features of the extended model. Importantly, we also show that chaos relies only in the slow oscillation subsystem, as chaos persists in the absence of the fast conductances that cause spiking. In addition, a bifurcation analysis shows that the transition from periodic to aperiodic behaviour-that is, from simple to chaotic dynamics-is organized by period doubling cascades (Guckenheimer and Holmes, 1983).
The manuscript is organized as follows: in Section 2 we describe the HB+Ih model and the numerical methods employed to analyze chaotic behavior. In Section 3, we present the results that include numerical simulations of the full system where we calculate different measures of chaos as we vary the system's parameters. Then we switch to the slow subsystem, doing systematic parameter explorations and bifurcation analysis. In Section 4 we summarize and discuss our findings.

Mathematical Model
The basis of our model is the Orio et al. (2012) model that reproduces the static firing patterns of cold thermoreceptors. The equation for the membrane voltage V is: where C m is the membrane capacitance; I d , I r , I sd , I sr are depolarizing (Na V ), repolarizing (K dr ), slow depolarizing (Na P / Ca T ) and slow repolarizing (K Ca ) currents, respectively. I h stands for hyperpolarization-activated current, and lastly I l represents the leak current. Currents are defined as: r, sd, h, l; (2) where a i is an activation term that represents the open probability of the channels (a l ≡ 1), with the exception of a sr that represents intracellular Calcium concentration. Parameter g i is the maximal conductance density, E i is the reversal potential and the function ρ(T) is a temperature-dependent scale factor for the current. The activation terms a r , a sd , and a h follow the differential equations: where On the other hand, a sr follows Finally, The function φ(T) in Equations (4) and (6) is a temperature factor for channel kinetics. The temperature-dependent functions for conductance ρ(T) in Equations (2-3), and for kinetics φ(T) in Equations (4) and (6) are given, respectively, by: Unless stated otherwise, the parameters used are given in Table 1.
Note that our set of equations presents some modifications from the original model in Orio et al. (2012). The first difference is that we do not consider the cold-inhibited trek current. This potassium current also contributes to the cold response as its inhibition produces depolarization of the cell (Viana et al., 2002;Noël et al., 2009). As the temperature dependence of the model is secondary to our objective, for simplicity reasons we omitted it. Secondly, and partially compensating the absence of the trek current, the leak current I l is now temperature-dependent (as the rest of ionic currents) due to the ρ(T) temperature factor. We also want to stress a departure from the original H&B model, namely the introduction of a saturable a sr -dependent expression in equation 3. This modification, already introduced by Orio et al. (2012), is necessary in order to perform a meaningful parameter exploration. In the original H&B formulation, a sr can grow far above 1 when I sd is high, making the g sr parameter no longer to be the maximal sr conductance.

Numerical Calculation of Maximal Lyapunov Exponent for Ordinary Differential Equations
The Lyapunov exponents give a measure of the exponential separation of nearby trajectories in a given direction (Guckenheimer and Holmes, 1983;Liu, 2010;Strogatz, 2014). In particular, a maximal Lyapunov exponent (MLE) greater than zero indicates sensitive dependence to initial conditions and, hence, is widely used as an indicator of chaos. We calculated MLEs from trajectories in the full variable space, following a standard numerical method based on that of Sprott (2003) (see also Jones et al., 2009).

Calculation of Lyapunov Exponent from Interval Time Series
In order to determine the Lyapunov exponent (LE) of the interspike interval (ISI) series, we proceeded as described in Kantz and Schreiber (2004). The method is based on Takens reconstruction theorem (Broer and Takens, 2011). Briefly, an ISI time series of length n is transformed into an m-dimensional reconstructed R m phase space, in which every k-th state point is specified by a vector with m elements, each one of them taken from the original ISI time series: For a given state point P i ∈ R m , we select all its neighbors {P * i } within a certain vicinity of radius ǫ and we measure the mean Euclidean distance d 0 from P i to the elements in {P * i }. The value of ǫ is chosen and constantly adjusted so that a maximum of 0.05% of points in R m fall within the vicinity. Next, the distances d 1 , d 2 , . . . d r are calculated from the following points in the series P i + 1 , P i + 2 , . . . , P i + r to the corresponding points that follow the elements in {P * i }, namely the sets The procedure is repeated for every single point in the series, thus obtaining a large number of distances d 0 , d 1 , . . . , d r .
Finally, we take the averages of the distances over all points d 0 , d 1 , ..., d r and the LE is taken to be the slope of the plot log( d i ) vs. i. We employed r = 6 and calculated LE for m (reconstructed dimension) = 7, 9 and 11. If for any value of m the regression yielded a p-value lower than 0.05 for the slope being different to 0, then we averaged the corresponding LE values. The number r, steps taken into account to advance the ISI series, and m, the embedding dimension, were empirically chosen by looking at the consistency of the results.

Lempel-Ziv Complexity Estimation
Lempel-Ziv complexity estimation method is an approximation to the Kolmogorov and Martin-Löf definition (Lempel and Ziv, 1976). This uses the idea that a computer program-as it scans an n-word string S = s 1 s 2 · · · s n from left to right-adds a new word to its memory (or "vocabulary") every time it discovers a sub-string of consecutive digits not previously encountered. The size of the vocabulary encountered and the rate at which new words are found along S are used in the Lempel-Ziv complexity measure. In this paper we are interested in the analysis of spike trains, thus to generate a binary sequence for a given spike-train it is necessary to divide the complete interval of measurement analysis in small sub intervals of size less than the minimum ISI and put one if there is a spike in the interval and zero if not.
Roughly speaking, the calculation of complexity is given by c(n)/b(n) where b(n) = n/ log 2 n and c(n) counts the number of steps necessary to reconstruct the sequence s 1 s 2 · · · s n of size n. The procedure to find c(n) can be explained with the diagram given in Kaspar and Schuster (1987) and summarized as follows: the first digit s 1 is always inserted to the vocabulary. Then, let s r be the last digit of the sequence S that has been reconstructed. We consider Q = s r + 1 and ask if Q is contained in the vocabulary of S. If s r + 1 can be obtained by repeating elements from the vocabulary, we define a new Q = s r + 1 s r + 2 and ask if it is in the vocabulary of S and so on until Q becomes so large that it cannot be obtained by copying a word from the vocabulary of SQπ (the operator π discards the last string added to SQ). Then, a new word is inserted into the vocabulary. c(n) is the number c of production steps to create a string, being the steps the vocabulary elements plus any repetition operation.

Bifurcation Analysis
Equilibrium states and periodic solutions of a dynamical system may undergo critical transitions under parameter variation. These re-arrangements may result in drastic changes -known as bifurcations-of the global dynamics including the onset of chaos (Guckenheimer and Holmes, 1983;Broer and Takens, 2011;Strogatz, 2014). Among the most simple bifurcation phenomena that one can find are saddle-node or limit point (LP) bifurcations -characterized by the sudden birth or disappearance of two equilibrium points-; and a Hopf bifurcation (HB) -where a periodic orbit is born from an equilibrium.
For the purposes of this work, the so-called period doubling or flip bifurcation plays a crucial role to understand the transition to chaos. This bifurcation is characterized by the loss of stability of a periodic orbit of period, say T, and the simultaneous birth of a secondary periodic orbit with period ≈ 2T. This process may repeat itself many times under parameter variation -within a relatively small range of parameter values-in a phenomenon knows as a period doubling cascade. At each occurrence of a period doubling bifurcation within the cascade a new orbit emerges with approximately twice the period of the one that had been born at the previous bifurcation. The consequence of this mechanism is the existence of aperiodic (chaotic) dynamics for a range of parameter values at one end of the cascade bifurcation values; see (Guckenheimer and Holmes, 1983;Broer and Takens, 2011) and the references therein for more details. Today, different software packages allow one to detect and continue a given bifurcation in one or two control parameters. In this paper, we make use of XPPAUT (and the numerical routines within it) to carry on a careful, detailed computational bifurcation analysis of equilibria and periodic orbits.

Numerical Simulation and Analysis
For long simulations, to obtain large ISI sequences, the model was implemented and run in the Neuron simulation environment (Hines and Carnevale, 2001) and run from Python scripts (Hines et al., 2009). Typically, ISI sequences were obtained from 1,000 s of simulation after 30 s of equilibration that were discarded to remove transient behaviors. For MLE calculation, simulations were solved with a fourth-order Runge-Kutta scheme written in Python. No detectable differences were found between Python and Neuron simulations. Data analysis and plotting was performed with Python and the libraries Numpy, Scipy, and Matplotlib.

Chaotic spiking in the HB+Ih Model
The HB+Ih model (Equations 1-8) studied here is an extension of the Huber & Braun (H&B) model of cold thermoreceptor (Braun et al., 1998). To this model, a hyperpolarization-activated current (I h ) was added in order to agree with experimental data obtained with I h blockers and HCN1 knock-out mice (Orio et al., 2012). Like the original model, and reproducing the behavior of cold thermoreceptors under static temperature conditions, this model shows a variety of firing patterns as the temperature is changed. Figure 1 shows typical time series (voltage trace) of deterministic simulations of the model at five different temperatures. At 20, 24.76, and 26 • C the model displays a periodic bursting pattern, decreasing the number of spikes per burst as temperature increases. At 33 • C, periodic tonic firing is observed, and at 36.3 • C, the pattern becomes irregular with "skipping, " i.e., some oscillations do not generate a spike and thus the intervals are distributed in a polymodal fashion. The irregular firing, evidenced in the ISI plot and the ISI histogram at 36.3 • C, suggest a typical chaotic dynamic. Figure 2A shows an ISI bifurcation plot against temperature. Visual inspection shows a high multiplicity of ISI values at almost every transition between firing modes: around 35 • C when skipping patterns appear (right zoom), around 29 • C when bursting occurs (left zoom) and then each time a new spike is added to the bursting pattern. This multiplicity of intervals suggests an irregular firing which is characteristic of chaotic behavior. However, calculation of the Lyapunov Exponent (LE) from the ISI time series (color code in Figure 2) reveals that not all the spike patterns that have a large number of ISI values are chaotic. Some of them, like the pattern at 24.76 • C in Figure 1B, have a large number of ISI values but still are highly repetitive and thus display a LE value near 0. We designate these firing patterns as "complex" but not chaotic. There are also some chaotic firing patterns around 10 • C, near the transition to the tonic firing behavior, like the original H&B model, which display chaos only between 7 and 12 • C (Feudel et al., 2000). However, the H&B model does not display chaotic dynamics at higher temperatures, where the model becomes physiologically more relevant. We suspect that chaos at high temperatures is introduced by the presence of the I h , and this is confirmed in Figure 2B where our model is simulated in the absence of I h (g h = 0). This diagram looks different to what has been described for the original H&B model (Feudel et al., 2000), because of some differences in parameters and the use of a saturable function of a sr in the I sr expression (equation 3). However, important qualitative features are conserved: no chaos (nor complex firing patterns) is observed above 10 • C, and chaotic firing patterns are only seen near the transition to tonic firing at low temperatures. This means that chaos at high temperatures is mainly introduced by I h and not by the other minor modifications that were made to the model (see Section 2.1). Although the HB+Ih model was inspired in the firing patterns of cold thermoreceptors at different temperatures, the combination of ion currents in this model is not exclusive to sensory neurons. The dynamics of this model can resemble other neurons in the CNS, where temperature changes are of lesser importance. Thus, we decided to study how the chaotic dynamics of the model depends on other parameters and the rest of the work presented here was performed with the temperature fixed at 36 • C, very close to the physiological value.

Chaotic Behavior in the Full System and the Role of Slow Conductances
We simulated the model with different combinations of the slow currents maximal densities g sd , g sr , and g h , and calculated the Lyapunov exponent of the ISI time series and the maximum Lyapunov exponent (MLE) from the voltage trajectories. As an alternative measure of chaotic behavior, we calculated the Lempel-Ziv complexity of the ISI data which has been used previously to prove the existence of chaos in neural models (Xu et al., 1997;Lu et al., 2008;Yang et al., 2009) and other disciplines (Frank and Stengos, 1988;Lu et al., 2008). MLE and complexity measures are shown in Figures 3B,C, respectively, together with the mean firing rate (average spikes per second) in Figure 3A and the firing pattern (bursting, tonic, skipping or no firing) in Figure 3D. In particular, Figure 3A shows that most of the explored region in parameter space corresponds to firing rates below 10 spikes/s. Though MLE ( Figure 3B) and Complexity ( Figure 3C) results are not completely overlapping, a high degree of correspondence can be seen on the results. Moreover, LE from ISIs mostly agrees with the results of MLE (not shown in Figure 3 but see for instance Figure 6A ). Chaotic behavior is concentrated in regions with low (<5) firing rate mostly, where skipping or polymodal firing pattern occurs. There is also chaotic firing at higher firing rates, but in narrower regions of the parameter space.
By looking at the center column of Figure 3 (g sd vs. g h ), we note that as g sd increases the system alternates between several firing patterns and it exhibits chaotic firing in a vicinity of almost every transition. To illustrate this better, we selected g h = 0.2 (white dashed line in B and D) and performed an ISI bifurcation diagram on the parameter g sd . Figure 4A shows that the model displays several firing modes and most (if not all) of the transitions between them imply a chaotic behavior. We find noteworthy the transitions between different skipping (polymodal) firing patterns -three of which are shown in Figure 4B-, and a transition between "skip-bursting" and regular tonic modes (denoted as (4) and (5), respectively). Figure 4 also shows how these chaotic transitions that separate different firing modes are created. Take for instance firing mode (1) where the ISI value is kept almost constant for a relatively large range of g sd values. As parameter g sd becomes greater than a certain critical value (inset), there are two possible ISI values. This duplication of ISI values continues as g sd is further increased, giving rise to what can be effectively understood as a cascade of "ISI doublings" -very much like in a period-doubling scenario. We confirmed this by repeating the plot at a much higher resolution (i.e. more values of g sd , separated by 10 −5 mS/cm 2 ) (inset). It is important to recall that each value of g sd was independently simulated, so there is no possible effect of the direction of parameter change; we also inspected the critical ISI time series to check that there were no transients involved in the results, thus ruling out an artifact due to initial conditions. The limit behaviour where g h = 0 can readily be seen in the center and right columns of Figure 3. This visual inspection suggests that in the absence of I h there is no chaotic behavior. To test this idea, we explored again the (g sd , g sr ) parameter subspace but now considering g h = 0 in Figure 5. The calculations clearly show that, in the absence of I h , no chaotic behavior is detected, even though similar firing rates and firing patterns are produced.
In order to better characterize the involvement of I h in the chaotic behavior of the model, we explored how the system depends on the time constant for this slow current, namely the parameter τ h . Figure 6A shows an ISI bifurcation diagram against τ h , showing that indeed the chaotic behavior depends on this parameter. In particular, the chaotic features disappear when the time constant is above 210 ms. In this Figure, we also show the good correspondence between the MLE calculated from voltage traces (top) and the LE calculated from the ISI sequences (see color code). Figures 6B,C show that the chaotic behavior depends on all the slow time constants, τ h , τ sr and τ sd .

The Chaotic Behavior in the Slow Subsystem
We further reduced the model by eliminating the fast spike mechanism and leaving only the slow oscillation mechanism. In other words, we take the instantaneous variable a d (t) ≡ 0 and the fast recovery variable a r (t) ≡ 0, which is the same as setting the parameters g d = g r = 0. In this way, the system now reduces, effectively, from five to four dimensions. As shown in Figure 7, the model still retains a chaotic behavior, showing a complex oscillatory pattern. Computations also show that -for certain parameter values-the solutions of interest converge in the long term to a strange attractor which is shown in Figure 7C in a projection onto the subspace of the a h , a sr , and a sd variables. To characterize this system, Figure 7B shows a return interval map (measured in ms) considering the voltage at the equilibrium point as a threshold (red line in Figure 7A). Figure 8A shows a bifurcation diagram that depicts the dependence of these return intervals on parameter g sd ; in addition, the color scale shows the corresponding LE measures. Note that chaotic regions are preceded -as g sd is increased-by doublings of return intervals very much like in Figure 4A. The color bar indicates: 0, no oscillations; 1, sub-threshold oscillations (no spikes); 2, oscillations and spikes with skipping; 3, regular tonic spiking; 4, burst firing (the shade represents the number of spikes per bursts); 5, tonic with firing rate between 20 and 50 spikes/s; 6, firing rate higher than 50 spikes/s. With this 4-dimensional reduced system we were also able to perform a bifurcation study of equilibrium points and periodic orbits, shown in Figure 8B for the voltage vs. the g sd parameter (the inset shows only the fixed points in a wider g sd range). This bifurcation diagram shows the birth of a stable periodic orbit at a Hopf bifurcation (labelled as HB). As g sd increases, this primary periodic orbit becomes the germ of a period doubling cascade (indicated by arrows in the enlargement in Figure 8C). This sequence of period doubling transitions coincides with the generation of chaotic regions identified in the return intervals plot. Note that only a limited number of period doubling branches were calculated and are shown here, because of space constraints. The period doubling events kept appearing as more branches were followed in the bifurcation. The chaotic regime is pulled back by a "reversed" period doubling mechanism as g sd is further increased; two of these bifurcations are illustrated here in inlet C2. Finally, for g sd > 0.3, a single stable periodic orbit exists. The oscillations eventually disappear at the g sd value corresponding to label LP which coincides with a saddle-node or limit point bifurcation of equilibria; this phenomenon is known as an infinite-period bifurcation or saddle-node homoclinic point (Aguirre, 2015). This combination of homoclinic phenomena and period-doubling cascades has been also reported and described as one of the mechanisms that produce chaos in the Hindmarsh-Rose equations (Linaro et al., 2012;Barrio et al., 2014Barrio et al., , 2017. A 2-parameter bifurcation diagram for the slow system is shown in Figure 9A. In this plot, we see how the period-doubling points, the Hopf bifurcation, and the Limit-Point bifurcation that ends the oscillation extend as bifurcation curves as both parameters g sd and g h are allowed to vary. As the maximal conductance g h is increased, more period doubling curves are added, enlarging the region of parameter values that allow chaos. Figure 9B shows the same bifurcation curves superimposed to the MLE values calculated from voltage trajectories. The resulting picture is revealing in that it shows how the regions of higher MLE fit perfectly into the predicted limits of chaotic dynamics generated by the period doubling phenomena. Hence, these findings emerge as another strong evidence to point out I h as the main responsible for the chaotic behavior of the model. This result is further enforced by the realization that -as with the full (spiking) system-there is no chaos in the absence of I h . Indeed, the period doubling bifurcation curves actually do not touch the g h = 0 axis. This fact becomes evident when the 1parameter bifurcation is done with g h = 0. Figure 10 shows that this system with g h = 0 indeed has no period doubling events, retaining only the Hopf and Limit-Point bifurcations.

DISCUSSION
In this article we have shown that a conductance-based model (denoted as HB+Ih model) displays chaotic behavior in many biologically plausible regions of the explored parameter space. These results are based on different numerical techniques to uncover and try to explain the regions of chaotic oscillations. To show the presence of chaos in a quantitative way, we have measured the maximal Lyapunov exponent of the system for different parameter scenarios, the Lyapunov exponent of the inter spike interval series and estimated the Lempel Ziv complexity. Bifurcation analysis on the reduced model without spikes allowed us to identify several period doubling cascades and propose them as the mathematical mechanism that originates chaos.
We discovered that the appearance of chaotic dynamics is related to the presence of a hyperpolarization-activated current (I h ), commonly found in neurons throughout the Central Nervous System (Biel et al., 2009;He et al., 2014) and that plays important roles in rhythmic firing, also controlling the membrane potential and regulating synaptic plasticity. The other distinctive elements of the model resemble persistent Sodium currents (I Na,P ) and Calcium-activated Potassium channels (I K,Ca ), also found in neurons that generate oscillatory rhythms (Llinás, 1988;Sanhueza and Bacigalupo, 2005). Thus, this model can serve as a general framework for studying how the interactions between different ion currents can generate chaotic oscillations.
The HB+Ih model offers the advantage of having a few number of variables while its equations and parameters are biophysically meaningful. In particular, chaos is observed in a 4-dimensional reduced model that considers the slow variables only. Unstable orbits and chaotic dynamics have been experimentally described in thermoreceptors (Braun et al., 1999a,b). The original H&B model, however, only exhibits such complex dynamics at around 8 • C, far from the physiological range (Feudel et al., 2000). The HB+Ih model, on the other hand, readily displays irregular and chaotic dynamics at temperatures above 30 • C, and thus it may be a more suitable system to explain the dynamics of cold thermoreceptors. Most of the parameter explorations presented here are based on the maximum conductance densities, that reflect the constantly changing levels of ion channel expression. Thus, the variety of firing patterns that we found here and the chaotic transitions between them are expected to be found under physiological conditions. Some adjustments incorporated by Orio et al. (2012) (compared to the original H&B model Braun et al., 1998) made a further shift toward biological plausibility. These changes include a slightly higher voltage dependence of the I Na,P -in accordance to what has been measured in somatosensory neurons (Herzog et al., 2001). Also, in the original H&B model the I sr current depends linearly on a sr , a variable that is not naturally bound by the model. Therefore, in a parameter space exploration where I sd can increase to high levels (because of a high g sd value), a sr would follow it to values much higher than 1, then losing its meaning as a channel open probability. In contrast, in the HB+Ih model, this term was replaced by a saturable binding term (see Equation 3) that allows the model to remain meaningful at high values of g sd .
Under the variation of parameters, even minimal 3-variable models of bursting neurons exhibit a rich variety of periodic and aperiodic dynamical patterns corresponding to different spiking and bursting regimes. The transitions between these patterns may contain complex dynamical structures such as perioddoubling (PD) cascades and deterministic chaos. For instance, it has been reported that transition mechanism from tonic spiking to bursting in a class of bursting neurons (square-wave bursters) is based on periodic spiking with a series of perioddoubling bifurcations followed by a homoclinic bifurcation of a saddle equilibrium (Terman, 1992;Wang, 1993;Feudel et al., 2000). Moreover, chaos has been proposed as a key signature for the transition between bursting and tonic firing (Chay, 1985;Rinzel and Ermentrout, 1989;Canavier et al., 1990;Terman, 1992;Feudel et al., 2000) (Terman (1992) gives a rigorous proof of the existence of Smale horseshoes). Then it has been shown that chaotic spiking can be generated close to the transition from spiking to bursting through perioddoubling cascades (Medvedev, 2006). In our model, chaotic regimes have been detected in regions at the transitions between different types of oscillations. The first region of chaos   Table 1. The red line depicts V eq , the value of V at the unstable singular point associated to the attractor. (B), Time intervals between successive crossings (with positive slope) of the V eq value. (C), 3D plot of a strange attractor in the a sd , a sr , a h sub-space (black trace). The 2-D projections onto the corresponding planes are shown in color.
appears through the transition from subthreshold oscillations (with no spikes) to irregular spiking (also called polymodal firing or skipping). This chaos appears through a cascade of Period Doublings, unrelated to the Hopf bifurcation that causes the appearance of subthreshold oscillations. In fact, when exploring different combination of g sd and g h , the range of g sd exhibiting subthreshold oscillation with no spikes become widened as g h increases, thus separating the Hopf bifurcation from the chaotic region. Even though at smaller g h the onset of chaos gets closer to the bifurcation, they seem to be unrelated. A second region of chaos emerges via periodic spikes to aperiodic spiking transition. These two kinds of transitions from periodic to aperiodic oscillations have been described as a mechanism that triggers chaotic spiking action potentials via period-doubling bifurcations. In the higher firing rate regimes, the HB+Ih model generates the third type of chaos at the transition from skipping to burst firing. In this case, the HB+Ih model exhibits chaotic bursting.
The period-doubling scenario was generally detected at one boundary of the chaotic regions (for instance, increasing g sd but not decreasing), being the other boundary of a different type. While other possible transitions to chaos have been found to be related to a boundary crisis (Arnol'd et al., 1994) as in the Hindmarsh-Rose model (Holden and Fan, 1992), here this question will remain the subject for future work.
The importance of time-scale differences in Hodgkin-Huxley type equations has already been investigated by some authors (Doi et al., 2001;Doi and Kumagai, 2005). In our case, we explored the contribution of different time scales and FIGURE 9 | (A) Two-D parameter bifurcation of the slow subsystem. The HB, LP and PD (period doubling) points identified in Figure 8B are followed as g sd and g h change. Note that the PD curves do not touch the g h = 0 axis. (B), The bifurcation curves are superimposed to the MLE calculated from voltage trajectories of the same system, to show that PD cascades delimit the chaotic regions. Dotted lines in (A,B) (black and white, respectively) refer to the parameter region explored in Figure 8. found that chaos appears only in certain ranges of slow time constants τ sd , τ sr and τ h . However, in contrast to the mentioned previous works, our model exhibits chaotic behavior within the biologically plausible values of these parameters.
Our findings emerge as an important contribution to the existing literature on the role of homoclinic bifurcations in concrete neuronal models (Feudel et al., 2000;Shilnikov and Cymbalyuk, 2005;Shilnikov, 2012). Future work on the chaotic behavior of the HB+Ih model can include a more detailed geometrical analysis of the chaotic attractors and a deeper investigation of bifurcations occurring at the onset of chaos. Indeed, chaotic dynamics can also be triggered by a wide range of global bifurcations such as homoclinic and heteroclinic phenomena (Aguirre et al., 2013(Aguirre et al., , 2014 which have yet to be analyzed in the HB+Ih model. The simplicity of this model and the fact that its equations and parameters maintain biophysical meaning, can make it a useful tool to understand how chaotic brain dynamics can be shaped by changes in ion channel expression or to give crucial insight to characterize properties of healthy and ill brains.