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

1. 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.

2. Methods

2.1. 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:

$\begin{array}{ll}{C}_m\frac{dV}{dt}=-I_{sd}-I_{sr}-I_h-I_d-I_r-I_l,& (1)\end{array}$

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:

$\begin{array}{ll}{I}_i=\rho (T)g_i{a}_i(V-E_i)i=d,r,sd,h,l;& (2)\end{array}$

$\begin{array}{ll}{I}_{sr}=\rho (T)g_{sr}\frac{a_{sr}^2}{a_{sr}^2+0.4^2}(V-E_{sr}),& (3)\end{array}$

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:

$\begin{array}{ll}\frac{d{a}_i}{dt}=\varphi (T)\frac{a_i^{\infty }(V)-a_i}{{\tau }_i}i=r,sd,h,& (4)\end{array}$

where

$\begin{array}{ll}{a}_i^{\infty }\left(V\right)=\frac{1}{1+exp\left(-s_i\left(V-V_i^0\right)\right)}.& (5)\end{array}$

On the other hand, a_sr follows

$\begin{array}{ll}\frac{d{a}_{sr}}{dt}=\varphi (T)\frac{-\eta I_{sd}-\kappa a_{sr}}{{\tau }_{sr}}.& (6)\end{array}$

Finally,

$\begin{array}{ll}{a}_d=a_d^{\infty }=\frac{1}{1+exp\left(-s_d\left(V-V_d^0\right)\right)}.& \left(7\right)\end{array}$

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:

$\begin{array}{ll}\rho (T)=1.3^{\frac{T-25}{10}}\varphi (T)=3^{\frac{T-25}{10}}.& (8)\end{array}$

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.

TABLE 1

Table 1. Parameters of the HB+Ih model.

2.2. Numerical Estimation of Chaotic Behavior

2.2.1. 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).

2.2.2. Calculation of Lyapunov Exponent from Interval Time Series

In order to determine the Lyapunov exponent (LE) of the inter-spike 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:

$\begin{array}{l}{P}_k=\left[IS{I}_k,IS{I}_{k+1},IS{I}_{k+2},\dots ,IS{I}_{k+m-1}\right],k=1,\dots ,n-m+1.\end{array}$

For a given state point P_i ∈ R_m, we select all its neighbors $\left\{P_i^{*}\right\}$ within a certain vicinity of radius ϵ and we measure the mean Euclidean distance d₀ from P_i to the elements in $\left\{P_i^{*}\right\}$. 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₁, d₂, …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 $\left\{P_i^{*}\right\}$, namely the sets $\left\{P_{i+1}^{*}\right\},\left\{P_{i+2}^{*}\right\},\dots ,\left\{P_{i+r}^{*}\right\}$. The procedure is repeated for every single point in the series, thus obtaining a large number of distances d₀, d₁, …, d_r.

Finally, we take the averages of the distances over all points 〈d₀〉, 〈d₁〉, …, 〈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.

2.2.3. 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₁s₂ ⋯ 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₂n and c(n) counts the number of steps necessary to reconstruct the sequence s₁s₂ ⋯ 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₁ 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+1s_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.

2.3. 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.

2.4. 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.

3. Results

3.1. 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 1

Figure 1. Firing patterns observed in the HB+Ih model. (A–E) Typical voltage time series (left), Inter-spike intervals (ISIs, middle) and ISI histograms (right) for the model at the shown temperatures. In (A–D) the histograms consider around 1100 spikes from 150 s simulations. In (E), the histogram was built from 2977 spikes obtained in a 1000 s simulation.

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).

FIGURE 2

Figure 2. (A) ISI bifurcation plot of the HB+Ih model as temperature is changed. Each temperature was independently simulated (from 0 to 40 °C in 0.01 °C intervals), and the spikes were considered to occur each time the voltage crossed a −15 mV threshold. The color of the dots represents the Lyapunov exponent calculated from the ISI series. Two sections are shown below in expanded horizontal scale. (B) ISI bifurcation plot of the model in the absence of I_h (g_h = 0).

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.

3.2. 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.

FIGURE 3

Figure 3. Chaotic behavior of the model at different combination of the slow conductance densities g_sd, g_sr and g_h. (A) Firing rate in spikes/s. (B) MLE of the voltage trajectories. (C) Lempel-Ziv complexity. (D) Firing pattern. 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.

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.

FIGURE 4

Figure 4. (A) ISI bifurcation of the model as the g_sd parameter is changed and g_h is fixed to 0.2 mS/cm². Numbers denote different firing modes between which chaotic regions are found. Inset shows a detail of the ISI-doubling events to the right of region (1), from g_sd = 0.217–0.219. (small black bar in the large plot) (B) Sample voltage trajectories showing the firing modes. (1), (2), and (3) correspond to different skipping patterns; (4) is bursting with skipping; (5) is regular tonic; and (6) is regular bursting. In all traces the scale bar is 250 ms.

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.

FIGURE 5

Figure 5. Absence of chaos in the model when g_h = 0. (A–D) are as described in Figure 3.

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.

FIGURE 6

Figure 6. Chaotic behavior dependency on the time constants of slow conductances τ_sd, τ_sr and τ_h. (A) ISI bifurcation diagram as τ_h is changed. Color of ISIs represents the LE of the ISI sequence. At the top, the line depicts the MLE of the voltage trajectories for the same values of τ_h. (B) MLE of voltage trajectories for different combinations of τ_h and τ_sd. (C) MLE for different combinations of τ_h and τ_sr. In (B,C), the white line indicates the corresponding region of the parameter space that is explored in (A).

3.3. 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.

FIGURE 7

Figure 7. Chaotic behavior in the slow subsystem of the model. (A) Voltage trajectory of the model with g_d = g_r = 0, and $g_{sd}=0.222(mS/c{m}^2)$. The rest of parameters are as given in 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.

FIGURE 8

Figure 8. Bifurcation analysis of the slow subsystem. (A) Intervals between V_eq crossings at different g_sd values. The colors represent the LE of the interval sequence. (B) Bifurcation of the system, as calculated by continuation methods in XPPAUT. Red lines represent stable fixed points; black line is an unstable fixed point. Green lines are maxima and minima of stable periodic atractors, blue are maxima and minima of unstable periodic attractors. The inset shows only the stable/unstable fixed point in a wider g_sd range. HB = Hopf bifurcation, LP=Limit point. (C), zoom of the C1 and C2 regions shown in (B). In C1, the arrows show period doubling events. In this Figure, g_h = 0.4

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.

FIGURE 9

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.

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.

FIGURE 10

Figure 10. Bifurcation analysis of the slow subsystem in the absence of I_h (g_h = 0). (A) Intervals between V_eq crossings at different g_sd values. The colors represent the LE of the interval sequence. (B) Bifurcation of the system, as calculated by continuation methods in XPPAUT. Red lines represent stable fixed points; black line is a unstable fixed point. Green is maxima and minima of stable periodic atractors. The inset shows only the stable/unstable fixed point in a wider g_sd scale. HB = Hopf bifurcation, LP = Limit point.

4. 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 (INa, P) and Calcium-activated Potassium channels (IK, 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 INa, 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 period-doubling (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 period-doubling 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 period-doubling 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 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 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, 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.

Author Contributions

KX, JM, MC, and PO performed numerical simulations and analysis. DQ and PO performed bifurcation continuation analysis. KX, JM, PA, and PO wrote the manuscript. KX, JM, MC, DQ, PA, and PO revised and approved the manuscript.

Funding

This work was funded by Grants Fondecyt 1130862, ACT-1113, and the Advanced Center for Electrical and Electronic Engineering (FB0008 Conicyt, Chile). The Centro Interdisciplinario de Neurociencia de Valparaíso (CINV) is a Millennium Institute supported by the Millennium Scientific Initiative of the Ministerio de Economía (Chile). PA was partially funded by Proyecto Fondecyt Iniciación 11150306.

Conflict of Interest Statement

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.

References

Abarbanel, H. D., Huerta, R., Rabinovich, M. I., Rulkov, N. F., Rowat, P. F., and Selverston, A. I. (1996). Synchronized action of synaptically coupled chaotic model neurons. Neural comput. 8, 1567–1602.

PubMed Abstract | Google Scholar

Aguirre, P. (2015). Bifurcations of two-dimensional global invariant manifolds near a non-central saddle-node homoclinic orbit. SIAM J. Appl. Dyn. Sys. 14, 1600–1644. doi: 10.1137/151004367

CrossRef Full Text | Google Scholar

Aguirre, P., Krauskopf, B., and Osinga, H. (2013). Global invariant manifolds near homoclinic orbits to a real saddle: (non)orientability and flip bifurcation. SIAM J. Appl. Dyn. Sys. 12, 1803–1846. doi: 10.1137/130912542

CrossRef Full Text | Google Scholar

Aguirre, P., Krauskopf, B., and Osinga, H. (2014). Global invariant manifolds near a Shilnikov homoclinic bifurcation. J. Comput. Dyn. 1, 1–38. doi: 10.3934/jcd.2014.1.1

CrossRef Full Text | Google Scholar

Aihara, K., Matsumoto, G., and Ikegaya, Y. (1984). Periodic and non-periodic responses of a periodically forced Hodgkin-Huxley oscillator. J. Theor. Biol. 109, 249–269.

PubMed Abstract | Google Scholar

Arnol'd, V. I., Afrajmovich, V., Il'yashenko, Y. S., and Shil'nikov, L. (1994). Dynamical Systems V: Bifurcation Theory and Catastrophe Theory, Vol. 5. Berlin; Heilderberg: Springer Science & Business Media. Available online at: http://www.springer.com/us/book/9783540181736

Barrio, R., Angeles Martínez, M., Serrano, S., and Shilnikov, A. (2014). Macro- and micro-chaotic structures in the Hindmarsh-Rose model of bursting neurons. Chaos 24:023128. doi: 10.1063/1.4882171

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrio, R., Ibáñez, S., and Pérez, L. (2017). Hindmarsh–rose model: Close and far to the singular limit. Phys Lett. A 381, 597–603. doi: 10.1016/j.physleta.2016.12.027

CrossRef Full Text | Google Scholar

Barrio, R., and Shilnikov, A. (2011). Parameter-sweeping techniques for temporal dynamics of neuronal systems : case study of Hindmarsh-Rose model. J. Math. Neurosci. 1:6. doi: 10.1186/2190-8567-1-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Biel, M., Wahl-Schott, C., Michalakis, S., and Zong, X. (2009). Hyperpolarization-Activated cation channels: from genes to function. Physiol. Rev. 89, 847–885. doi: 10.1152/physrev.00029.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Braun, H. A., Bade, H., and Hensel, H. (1980). Static and dynamic discharge patterns of bursting cold fibers related to hypothetical receptor mechanisms. Pflüg. Archiv. Eur. J. Physiol. 386, 1–9.

PubMed Abstract | Google Scholar

Braun, H. A., Dewald, M., Schäfer, K., Voigt, K., Pei, X., Dolan, K., et al. (1999a). Low-dimensional dynamics in sensory biology 2: facial cold receptors of the rat. J. Comput. Neurosci. 7, 17–32.

PubMed Abstract | Google Scholar

Braun, H. A., Dewald, M., Voigt, K., Huber, M., Pei, X., and Moss, F. (1999b). Finding unstable periodic orbits in electroreceptors, cold receptors and hypothalamic neurons. Neurocomputing 26–27, 79–86.

Google Scholar

Braun, H. A., Huber, M. T., Dewald, M., Schäfer, K., and Voigt, K. (1998). Computer Simulations of Neuronal Signal Transduction: The Role of Nonlinear Dynamics and Noise. Int. J. Bifurc. Chaos 08, 881–889.

Google Scholar

Broer, H., and Takens, F. (2011). Dynamical Systems and Chaos. New York, NY: Springer. Available online at: http://www.springer.com/us/book/9781441968692

Google Scholar

Canavier, C. C., Clark, J. W., and Byrne, J. H. (1990). Routes to chaos in a model of a bursting neuron. Biophys. J. 57, 1245–1251.

PubMed Abstract | Google Scholar

Chay, T. R. (1985). Chaos in a three-variable model of an excitable cell. Phys. D Nonlinear Phenomena 16, 233–242.

Google Scholar

Chay, T. R., and Rinzel, J. (1985). Bursting, beating, and chaos in an excitable membrane model. Biophys. J. 47, 357–366.

PubMed Abstract | Google Scholar

Doi, S., and Kumagai, S. (2005). Generation of very slow neuronal rhythms and chaos near the hopf bifurcation in single neuron models. J. Comput. Neurosci. 19, 325–356. doi: 10.1007/s10827-005-2895-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Doi, S., Nabetani, S., and Kumagai, S. (2001). Complex nonlinear dynamics of the Hodgkin–Huxley equations induced by time scale changes. Biol. Cybern. 85, 51–64. doi: 10.1007/PL00007996

PubMed Abstract | CrossRef Full Text | Google Scholar

Falcke, M., Huerta, R., Rabinovich, M. I., Abarbanel, H. D. I., Elson, R. C., and Selverston, A. I. (2000). Modeling observed chaotic oscillations in bursting neurons: the role of calcium dynamics and IP3. Biol. Cybern. 82, 517–527. doi: 10.1007/s004220050604

PubMed Abstract | CrossRef Full Text | Google Scholar

Faure, P., Kaplan, D., and Korn, H. (2000). Synaptic efficacy and the transmission of complex firing patterns between neurons. J. Neurophysiol. 84, 3010–3025.

PubMed Abstract | Google Scholar

Feudel, U., Neiman, A., Pei, X., Wojtenek, W., Braun, H., Huber, M., et al. (2000). Homoclinic bifurcation in a Hodgkin-Huxley model of thermally sensitive neurons. Chaos 10, 231–239. doi: 10.1063/1.166488

PubMed Abstract | CrossRef Full Text | Google Scholar

Frank, M., and Stengos, T. (1988). Chaotic dynamics in economic time-series. J. Economic Surveys 2, 103–133.

Google Scholar

Gu, H. (2013). Biological experimental observations of an unnoticed chaos as simulated by the Hindmarsh-Rose model. PLoS ONE 8:e81759. doi: 10.1371/journal.pone.0081759

PubMed Abstract | CrossRef Full Text | Google Scholar

Guckenheimer, J., and Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. New York, NY: Springer. Available online at: http://www.springer.com/br/book/9780387908199

Google Scholar

He, C., Chen, F., Li, B., and Hu, Z. (2014). Neurophysiology of HCN channels: From cellular functions to multiple regulations. Prog. Neurobiol. 112, 1–23. doi: 10.1016/j.pneurobio.2013.10.001. Available online at: http://jn.physiology.org/content/86/3/1351.long

PubMed Abstract | CrossRef Full Text | Google Scholar

Herzog, R. I., Cummins, T. R., and Waxman, S. G. (2001). Persistent TTX-resistant Na+ current affects resting potential and response to depolarization in simulated spinal sensory neurons. J. Neurophysiol. 86, 1351–1364.

PubMed Abstract | Google Scholar

Hindmarsh, J. L., and Rose, R. M. (1984). A model of neuronal bursting using three coupled first order differential equations. Proc. R. Soc. B Biol. Sci. 221, 87–102.

PubMed Abstract | Google Scholar

Hines, M. L., and Carnevale, N. T. (2001). NEURON: a tool for neuroscientists. Neuroscientist 7, 123–135. doi: 10.1177/107385840100700207

PubMed Abstract | CrossRef Full Text | Google Scholar

Hines, M. L., Davison, A. P., and Muller, E. (2009). NEURON and Python. Front. Neuroinformat. 3:1. doi: 10.3389/neuro.11.001.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Holden, A. V., and Fan, Y. (1992). Crisis-induced chaos in the Rose-Hindmarsh model for neuronal activity. Chaos Solitons Fract. 2, 583–595.

Google Scholar

Jones, D. S., Plank, M., and Sleeman, B. D. (2009). Differential Equations and Mathematical Biology. Boca Raton, FL: CRC Press.

Google Scholar

Kantz, H., and Schreiber, T. (2004). Nonlinear Time Series Analysis. New York, NY: Cambridge University Press.

Google Scholar

Kaspar, F., and Schuster, H. G. (1987). Easily calculable measure for the complexity of spatiotemporal patterns. Phys. Rev. A 36, 842–848.

PubMed Abstract | Google Scholar

Korn, H., and Faure, P. (2003). Is there chaos in the brain? II. Experimental evidence and related models. Comptes Rendus Biol. 326, 787–840. doi: 10.1016/j.crvi.2003.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lempel, A., and Ziv, J. (1976). On the complexity of finite sequences. IEEE Trans. Inform. Theory 22, 75–81.

Google Scholar

Linaro, D., Champneys, A., Desroches, M., and Storace, M. (2012). Codimension-two homoclinic bifurcations underlying spike adding in the hindmarsh–rose burster. SIAM J. Appl. Dyn. Sys. 11, 939–962. doi: 10.1137/110848931

CrossRef Full Text | Google Scholar

Liu, Z. (2010). Chaotic time series analysis. Math. Prob. Eng. 2010:720190. doi: 10.1155/2010/720190

CrossRef Full Text | Google Scholar

Llinás, R. R. (1988). The intrinsic electrophysiological properties of mammalian neurons: insights into central nervous system function. Science 242, 1654–1664.

PubMed Abstract | Google Scholar

Lu, Q., Gu, H., Yang, Z., Shi, X., Duan, L., and Zheng, Y. (2008). Dynamics of firing patterns, synchronization and resonances in neuronal electrical activities: experiments and analysis. Acta Mechan. Sin. 24, 593–628. doi: 10.1007/s10409-008-0204-8

CrossRef Full Text | Google Scholar

Medvedev, G. S. (2006). Transition to bursting via deterministic chaos. Phys. Rev. Lett. 97:048102. doi: 10.1103/PhysRevLett.97.048102

PubMed Abstract | CrossRef Full Text | Google Scholar

Noël, J., Zimmermann, K., Busserolles, J., Deval, E., Alloui, A., Diochot, S., et al. (2009). The mechano-activated K+ channels TRAAK and TREK-1 control both warm and cold perception. EMBO J. 28, 1308–1318. doi: 10.1038/emboj.2009.57

PubMed Abstract | CrossRef Full Text | Google Scholar

Orio, P., Parra, A., Madrid, R., González, O., Belmonte, C., Viana, F., et al. (2012). Role of Ih in the firing pattern of mammalian cold thermoreceptor endings. J. Neurophysiol. 108, 3009–3023. doi: 10.1152/jn.01033.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Plant, R. E., and Kim, M. (1976). Mathematical description of a bursting pacemaker neuron by a modification of the Hodgkin-Huxley equations. Biophys. J. 16, 227–244.

PubMed Abstract | Google Scholar

Rinzel, J., and Ermentrout, G. B. (1998). “Analysis of neural excitability and oscillations,” in Methods in Neuronal Modeling: From Synapses to Networks, 2 Edn., eds C. Koch and I. Segev (Cambridge MA: MIT Press), 251–291.

Google Scholar

Sanhueza, M., and Bacigalupo, J. (2005). Intrinsic subthreshold oscillations of the membrane potential in pyramidal neurons of the olfactory amygdala. Euro. J. Neurosci. 22, 1618–1626. doi: 10.1111/j.1460-9568.2005.04341.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shilnikov, A. (2012). Complete dynamical analysis of a neuron model. Nonlinear Dyn. 68, 305–328. doi: 10.1007/s11071-011-0046-y

CrossRef Full Text | Google Scholar

Shilnikov, A., and Cymbalyuk, G. (2005). Transition between tonic spiking and bursting in a neuron model via the blue-sky catastrophe. Phys. Rev. Lett. 94, 2–5. doi: 10.1103/PhysRevLett.94.048101

PubMed Abstract | CrossRef Full Text | Google Scholar

Sprott, J. C. (2003). Chaos and Time-series Analysis, Vol. 69. Oxford: Oxford University Press.

Google Scholar

Strogatz, S. H. (2014). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Boulder: Westview Press.

Google Scholar

Terman, D. (1992). The transition from bursting to continuous spiking in excitable membrane models. J. Nonlinear Sci. 2, 135–182.

Google Scholar

Viana, F., de la Peña, E., and Belmonte, C. (2002). Specificity of cold thermotransduction is determined by differential ionic channel expression. Nat. Neurosci. 5, 254–260. doi: 10.1038/nn809

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X.-J. (1993). Genesis of bursting oscillations in the hindmarsh-rose model and homoclinicity to a chaotic saddle. Phys. D Nonlinear Phenom. 62, 263–274.

Google Scholar

Xu, J., Liu, Z.-R., Liu, R., and Yang, Q.-F. (1997). Information transmission in human cerebral cortex. Phys. D Nonlinear Phenom. 106, 363–374.

Google Scholar

Yang, M., Liu, Z., Li, L., Xu, Y., Liu, H., Gu, H., et al. (2009). Identifying distinct stochastic dynamics from chaos: a study on multimodal neural firing patterns. Int. J. Bifurc. Chaos 19, 453–485. doi: 10.1142/S0218127409023135

CrossRef Full Text | Google Scholar