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:

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.

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.

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.

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.

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⁻⁵mS/cm²) (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.

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.

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., 2014, 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 1-parameter 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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.