1. Introduction

Scientists have been always fascinated by the functioning of the human brain and always attempted to understand its complexity.

Only in 1899 Santiago Ramon y Cajal, exploiting the experimental techniques developed by Camillo Golgi, was able to discover that the nervous system is made by individual cells, later called neurons [31]. His studies laid the foundations for the so-called ''Neuron Doctrine'' and gave him the Nobel Prize in Physiology and Medicine in 1906 [11]. Although some scientists suggest to rethink the Neuron Doctrine [6], it remains the pillar of modern neuroscience.

It is worth observing that in each individual there is not an only type of neuron, but several ones. Nevertheless, they share many common properties. From the cell body (called soma) starts a number of ramifying branches called dendrites. These structures constitute the input pole of a neuron. From the soma originates also a long fiber called the axon. It is considered as the output line of a neuron since through the axons terminals, called synapses, the exchange of information with other neurons takes place [33,35].

Indeed, the basic elements of the communication among the neurons are pulsed electric signals called action potentials or spikes. In fact, the neuronal cell is surrounded by a membrane, across with there is a difference in electrical charge (called also potential), that depends on the different concentrations of ions, especially Sodium (Na$^+$), Potassium (K$^+$) and Calcium (Ca$^{++}$), inside and outside the neuron. If the membrane potential exceeds a certain threshold, then the neuron generates a brief electrical pulse, that propagates along the axon. Finally the synapses transfer this electrical signal to the other surrounding neurons [19,21,33].

Moreover, a neuron can exhibit rich dynamical behaviors, such as resting, excitable, periodic spiking, and bursting activities. In particular, the ability of periodic firing has been recorded in isolated neurons since the 1930s [37] and Hodgkin [17,19] was the first to propose a classification of neural excitability, depending on the frequency of the action potentials generated by applying external currents.

In literature several nonlinear dynamical systems have been proposed to suitably model the dynamics of the electrical activity observed in a single neuron.

However, the paper by Hodgkin and Huxley [18] on the physiology of the squid giant axon remains a milestone in the science of nervous system, and the model proposed therein has been extensively studied, due to its richness.

It is known [13,15,32] that, depending on the value of the external current stimuli $I$, the Hodgkin-Huxley (HH) model exhibits periodic behaviors. Moreover, for a certain range of $I$, it behaves as an hard oscillator [26,29], that is there is a coexistence of a stable equilibrium and a stable limit cycle. Thus, this implies the existence of unstable limit cycles in order to separate the two basins of attraction. Few authors have been able to characterize these periodic solutions, how they emerge and disappear depending on the intensity of the external current [13,15,32].

In general, it is not an easy task to detect a periodic solution of a nonlinear dynamical system. Several methods are exploited to predict the existence of limit cycles and to study their stability, both in time and frequency domain [2,4,5,23,25,27,28,30]. In particular, for the HH model, this is even more difficult due to the high nonlinear structure of the system.

Our aim is to propose a technique based on the harmonic balance method [2,27] in order to characterize the periodic solutions exhibited by a nonlinear dynamical system. In particular, for the HH model we will show how this technique is more efficient than the previous approaches based either on finite differences, collocation or shooting methods [9,13,14,15,32]. Indeed, by applying the harmonic balance method we obtain an approximated but analytical expression of the periodic solutions, therefore we get a good approximation of the monodromy matrix and the Floquet analysis is straightforward. Moreover, it is worth noting that the harmonic balance method has an exponential convergence, while the collocation methods exploited in the previous papers [9,13,14,15,32], even with a mesh adaptation, have only a polynomial one. Finally, only from 2 to 50 harmonics are needed to correctly approximate the periodic solutions, so the linear system to solve is low-dimensional.

In the end, the application of the harmonic balance method to the Hodgkin-Huxley model displays several advantages: a straightforward implementation without mesh adaptation, exponential convergence rate, an analytical representation of the stable and unstable periodic solutions. Moreover, this method can be coupled with a simple continuation method based on a selective parameterization instead of the exploitation of the Keller arc-length parameterization [7].

The paper is structured as follows: in Section 2 the basics of collocation and harmonic balance methods are briefly recalled and our technique, that we call multi-base harmonic balance method is introduced. In Section 3 firstly we present the structure of the HH model and its main characteristics. We show how it is possible to obtain the bifurcation diagram of the HH model by mainly exploiting the multi-base harmonic balance method. Furthermore, all the bifurcations are analyzed via the harmonic balance method and Floquet analysis. Finally, concluding remarks are offered in Section 4.

---

2. Computation of periodic solutions

There are several methods that are exploited to predict the existence of limit cycles in nonlinear systems and to study their stability. Classical methods, based on integration schemes, such as the simple shooting method, are generally sensitive to the stability properties of solutions, so they cannot be exploited to detect unstable periodic solutions. In this section, first of all, we review two methods that overcome this problem and belong to a class of spectral methods. Indeed, we present collocation and harmonic balance methods, whose main idea is the approximation of the exact solution by the projection on a finite-dimensional subspace. It has been shown that the harmonic balance method is more efficient than the collocation method [20], but is computationally onerous when the periodic solution contains high harmonics. Therefore, in this paper we propose a new technique, a multi-base harmonic balance method, that permits to obtain a good approximation of a periodic function even if the nonlinear dynamical system under study is highly nonlinear.

---

2.1. Collocation methods

Let us consider an autonomous dynamical system

$\dot x = f(x)$

(1)

where $f$ is a vector field defined on $\mathbb{R}^n$, $n\geq 1$, and $x\in \mathbb{R}^n$. A solution $x = X$ of a continuous-time system is periodic if $X(t+T) = X(t)$ $\forall \;t$, for some $T>0$. The minimal $T$ for which this equality holds is called the period of the solution. This periodic solution $X$ of least finite period $T>0$ of the system corresponds to a closed orbit $\Gamma$ in $\mathbb{R}^n$. On this orbit each initial time corresponds to a location $x = x_0$.

It is interesting to notice that searching a periodic solution of an ODE is equivalent to the resolution of a boundary value problem (BVP). In fact, if $x = X$ is a $T$-periodic solution of (1), then it is solution of the following BVP:

$\begin{cases} \dfrac{dx}{dt} = f(x)\\ x(0) = X(0) = X(T). \end{cases}$

(2)

System (2) belongs to the general class of nonlinear boundary value problems and several methods have been proposed in literature to solve it [2]. The more intuitive one is probably the simple shooting method [24]. The idea is to find an initial condition $X_0$ and a period $T$ such that $X(0) = X_0 = X(T),$ where the unknown is the couple $(X_0, T)$, with $X_0 = (x_1, ..., x_n) \in \mathbb{R}^n$ and $T \in \mathbb{R}^{+}$. If we define the functional $G$ as

$G(X_0, T ) =\varphi(X_0, T ) - X_0,$

(3)

where $\varphi(X_0, T )$ is the solution of the Cauchy problem $\dot x =f (x)$ with the initial condition $x(0) =X_0$, then it is straightforward to see that a periodic solution of (2) is a zero of $G$. The nonlinear system

$G(X_0, T ) =\varphi(X_0, T ) - X_0 = 0,$

(4)

has $n$ equations for $n+1$ unknowns $X_0 = (x_1, ..., x_n)$ and $T$. Therefore, one fixes the first component $x_1$ of vector $X_0$ to be on the trajectory (see [32]) and consequently, solves (4) by Newton's method for the unknowns $x_2, ..., x_n, T$. Unfortunately, this method is sensitive to the initial conditions. Thus, it fails to detect unstable limit cycles. In order to avoid this flaw, the multiple shooting method has been proposed [2].

On the contrary, the collocation method is independent on the stability of the periodic solutions under consideration. We briefly recall the main properties of this method [34].

First of all, since $T$ is usually unknown, with a simple normalization

of the time scale, it is possible to write (2) on the interval $[0, 1]$:

$\begin{cases} \dfrac{d u}{d \tau} = Tf(u)\\ u(0) = u(1), \end{cases}$

(5)

where $\tau = \frac{t}{T}$ is the new time variable. Clearly, a solution $u(\tau)$ of (5) corresponds to a $T$-periodic solution of (2). However, the boundary condition in (5) does not define a unique periodic solution. Indeed, any time shift of a solution to the periodic BVP (5) is still a solution. Thus, an additional condition has to be appended to problem (5) in order to select a solution among all those corresponding to the cycle.

The idea of the collocation methods for BVPs is to approximate the analytical solution by a piecewise polynomial vectorial function $P(t)$ belonging to $\mathbb{R}^n$ that satisfies the boundary conditions and the original problem on selected points, called collocation points.

Let us consider the partition $0 = t_{0} < t_{1} < \dots < t_{N} = 1$, then the approximated solution $P(t)$ has the following form:

$P(t)/_{[{t_i}, {t_{i+1}}]} = P_{i}(t), \, \, \, \, \, \, \, \, i = 0, \dots, N-1,$

where $P_{i}$ is a polynomial of degree $m$ for all $i = 0, \dots, N-1$, and $P_{i}(t_{i+1}) = P_{i+1}(t_{i+1})$, in order to have a continuous polynomial on the whole interval $[0, 1]$. On each sub-interval $[{t_i}, {t_{i+1}}]$ we introduce the collocation points

$t_{ij} =t_i + \rho_{j}(t_{i+1}-t_{i}), ~i = 0, ~\dots, ~N-1, ~ j =1, ~\dots, ~m,$

where $0\leq \rho_{1}< \dots < \rho_{m} \leq 1$.

Then, the request that $P(t)$ satisfies the BVP (5) on the collocations points leads to the following non-linear algebraic system of equations:

$\label{eq:colloc} \frac{1}{T} \dot P(t_{ij}) = f(P (t_{ij} )), \, \, i =0, \dots, N-1, j =1, \dots, m,$

(6)

with the boundary conditions

$\label{eq:cond} P (0) =P (1).$

(7)

Thus, the state vector is given by

$U = (P(0), P(t_{ij})_{0\leq i \leq N-1, 1 \leq j \leq m}, P(1), T)\in \mathbb{R}^{q},$

where $n$ is the dimension of $X$ in (5), $m\cdot N$ is the number of unknowns in (6), and $q = m \cdot N \cdot n +2n+1$.

Moreover, a further condition is necessary to determine the unknown parameter $T$. This is the so-called phase condition. For example, a condition based on the Poincaré section has been used in [32], while an integral condition can be found in [8]. In this work, since over one period the derivative of a periodic function is null in at least one point, we use the condition

$\dot u_1(0) = 0,$

where $u_1$ is the first component of $u\in \mathbb{R}^n$.

Several solvers using collocation methods have been proposed in the literature. For example, COLSYS/COLNEW [1,3] and AUTO [8] which use collocation methods with Gaussian points, or bvp4c, bvp5c and bvp6c [22,34] that are routines with Lobatto points.

In particular, the bvp4c method exploits 3 Lobatto points on each subinterval. In this case $P(t)$ is a piecewise vectorial cubic polynomial, $P(t)\in C^{1}([0, ~1]^n)$ [22,34].

The bvp4c methods controls the error $||u(t)-P(t)||$ indirectly, by minimizing the norm of the residue $r(t) = ||\dot{P}(t)-F(P(t))||$. Thus, $P(t)$ is considered as the exact solution of the following perturbed problem

$\dot{P}(t) = F(P(t))+r(t), ~~H(P(0), P(1)) = \delta,$

where $H$ is a function of $P(0)$ and $P(1)$ that describes the boundary conditions, and $\delta$ is the associated residue (in our case, $\delta = (P(1)-P(0))$).

The basic idea of this technique is to minimize the residue over each sub-interval $[x_i, ~x_{i+1}]$, to adjust the mesh as one goes along, such that the norm of the residue $r$ and $\delta$ tend to zero. This assures that the approximated solution $P(t)$ converges to the exact solution of the problem. It is worth noting that this mesh adaptation method is able to obtain the convergence even in case of bad initial conditions and with an approximation of order 4, i.e. $||u(t)-P(t)||< C h^4,$ where $C$ is a constant and $h$ is the maximum mesh step. For further details, see [34].

---

2.2. Harmonic balance (HB) method

It is worth noting that any periodic smooth function of period $T$ can be represented as an infinite Fourier series

$X(t) =A_0 + \sum\limits_{k = 1}^{\infty} \left(A_k \cos \left(k \dfrac{2\pi}{T} t\right) + B_k \sin \left(k \dfrac{2\pi}{T} t\right)\right),$

(8)

where

$\begin{align}\label{eq:Four_coeff} A_0& = \dfrac{1}{T}\displaystyle\int_{0}^{T} X(t) dt, \nonumber\\ A_k & = \dfrac{2}{T}\displaystyle\int_{0}^{T} X(t) \cos \left(k \dfrac{2\pi}{T} t\right)dt, ~~~~k = 1, ..., \infty.\\ B_k& = \dfrac{2}{T}\displaystyle\int_{0}^{T} X(t) \sin \left(k \dfrac{2\pi}{T} t\right)dt, \nonumber \end{align}$

(9)

The idea of the harmonic balance method [27] is to search for an approximation of the solution of (1) as a truncated series

${X}_{K}(t) =A_0 +\sum\limits_{k = 1}^{K} \left(A_k \cos \left(k \dfrac{2\pi}{T} t\right) + B_k \sin \left(k \dfrac{2\pi}{T} t\right)\right),$

(10)

where $K$ is the number of harmonics taken into account. If the periodic solution $X(t)$ is smooth, then the truncated series ${X}_{K}(t)$ converges to $X(t)$ rapidly, without exhibiting the Gibbs phenomenon [36]. In fact, if the unknown periodic solution is smooth, then it is possible to show that its Fourier series converges rapidly and uniformly. If $X$ has continuous derivatives up to order $n \geq 1$, and $X^{(n)}(t)$ is integrable, we can integrate $n$ times by parts the coefficients $A_k$ and $B_k$ in (9) and obtain

$A_k = \frac{(-1)^{(n-1)}}{r^n}\frac{2}{T} \int_{0}^{T} X^{(n)}(t)\sin(rt)dt$

and

$B_k = \frac{(-1)^{(n)}}{ r^n} \frac{2}{T}\int_{0}^{T} X^{(n)}(t)\cos(rt)dt,$

with $r = k \dfrac{2\pi}{T}$. Therefore, it is clear, by applying the Riemann-Lebesgue lemma, that $A_k$ and $B_k$ are $O(\frac{1}{k^n})$, and when $n$ increases to infinity, the convergence is uniform and exponential. For more details, see [12].

Analogously, the function $f({X}_{K}(t))$ will be still periodic of period $T$, therefore also this term can be expanded in a truncated Fourier series as

$\label{n-Fourier} f({X}_{K}(t)) = N_0+\sum\limits_{k = 1}^{K}\left(M_k\cos\left(k \dfrac{2\pi}{T} t\right)+ N_k\sin\left(k \dfrac{2\pi}{T} t\right)\right),$

(11)

wherein each coefficient $M_k$, $N_k$ will be function of all the coefficients $A_k$, $B_k$ of ${X}_{K}(t)$. This dependence can be expressed by the relations:

$\begin{align} \label{coeff-Fourier} N_0 & =\frac{1}{T}\int_0^T f({X}_{K}(t))\, dt, \notag \\ M_k & =\frac{2}{T}\int_0^T f({X}_{K}(t))\cos\left(k\dfrac{2\pi}{T} t\right)\, dt, \\ N_k & =\frac{2}{T}\int_0^T f({X}_{K}(t))\sin\left(k\dfrac{2\pi}{T} t\right)\, dt. \notag \end{align}$

(12)

We note $\overline{X}_K = (A_0, A_1, B_1, ..., A_K, B_K)$ the vector of Fourier coefficients, and

$e_{j}(t) =\begin{cases} 1~~~\text{if}~~j = 0, \\ \cos \left(j \dfrac{2\pi}{T} t\right)~~~\text{if}~~j = 2k, ~~~~~~k = 1, ...K.\\ \sin \left(j \dfrac{2\pi}{T} t\right)~~~\text{if}~~j = 2k+1, \end{cases}$

the Fourier base. Introducing (10) and (12) in (1), and taking into account the orthogonality of the elements of the Fourier base, we obtain a nonlinear algebraic system in the unknowns $\bar X_K$ and $T$.

Unfortunately, if the function $f$ is highly nonlinear and when the harmonic balance is applied with an high number of harmonics, it is not easy to find an analytical expression of coefficients (12) and the numerical resolution of the nonlinear algebraic system could be computationally too expensive [23].

Remark 1. It is worth recalling the drawbacks of the methods presented before:

$\bullet$ The simple shooting method fails when one wants to detect an unstable limit cycle.

$\bullet$ The harmonic balance method is more efficient and has a higher rate of convergence than the collocation methods [20], but it is not exploitable when a high number of harmonics is needed or the system is highly nonlinear.

In the following we present a technique mixing Fourier series and trigonometric Lagrange interpolation in the harmonic balance method that permits to avoid such inconveniences.

---

2.3. A multi-base harmonic balance method

If $X(t)$ is a $T$-periodic continue function, then the $n^{th}$ trigonometric Lagrange interpolation polynomial of $X(t)$ with equally spaced nodes is the following [38]

$L_n(X(t), t) = \sum\limits_{j = 0}^{n} X(t_{j})l_{j}(t),$

where

$l_{j}(t) = \dfrac{\sin\left(\left(n+\dfrac{1}{2}\right)\left(t\dfrac{T}{2\pi}-t_j\right)\right)}{\sin\left(\dfrac{1}{2}\left(t\dfrac{T}{2\pi}-t_j\right)\right)}, \, \, \, \, t_j = \dfrac{jT}{2n+1}, \, \, \, \, j = 0, 1, ..., 2n.$

Then, the functions $(l_{j})_{j = 0, \dots, 2n}$ constitute an Hermitian base of the periodic functions space.

Let us suppose $f$ in (2) to be a polynomial of degree $d$, and let us consider $n = d\times K$. Let $Y_{F}(\overline{X}_K)$ be the coefficients of $f({X}_{K}(t))$ in the Fourier base $(e_{j})_{j = 0, \dots, 2K}$ and $Y_{L}(\overline{X}_K)$ the components of $f(X_{K}(t))$ in the Lagrange base $(l_{j})_{j = 0, \dots, 2n}$. It is easy to see that

$(Y_{L}(\overline{X}_K))_{j} = f(X_{K}(t_{j})).$

Moreover, it is possible to deduce $Y_{F}(\overline{X}_K)$:

$Y_{F}(\overline{X}_K) =P\Gamma^{-1} Y_{L}(\overline{X}_K),$

where $P$ is the projection matrix

$P = \begin{pmatrix} I_{2K+1}&0&\dots&0 \end{pmatrix}\in \mathbb{R}^{2K+1, 2n+1},$

$I_{2K+1}$ is the identity matrix of rank $2K+1$, and $\Gamma^{-1}$ is the transition matrix between the two bases $(l_{j})_{j = 0, \dots, 2n}$ and $(e_{j})_{j = 0, \dots, 2n}$. It is easy to notice that the elements of $\Gamma^{-1}$ can be determined as the values of the functions $(e_{j}(t))_{j = 0, \dots, 2n}$ in the nodes $t_j$:

$\Gamma^{-1} = \left(\begin{matrix} 1 &\cos(1\times t_0) &\sin(1\times t_0) &\cdots\cdots\cdots &\cos(n\times t_0) &\sin(n\times t_0)\\ \vdots &\vdots &\vdots &\cdots\cdots\cdots &\vdots &\vdots\\ \vdots &\vdots &\vdots &\cdots\cdots\cdots &\vdots &\vdots\\ 1 &\cos(1 \times t_{2n}) &\sin(1 \times t_{2n}) &\cdots\cdots\cdots &\cos(n\times t_{2n}) &\sin(n \times t_{2n})\\ \end{matrix} \right)$

It is worth observing that the choice of $n = dK$ and the utilization of the projection matrix $P$ permit to avoid a sort of aliasing phenomenon [16], that could take place if $n = K$. Moreover, this technique can be generalized to the case of a non-polynomial nonlinearity and, in this case, $n$ can be determined by looking at the convergence rate of the Fourier coefficients with respect to the number of considered harmonics.

As far as we know, this joint exploitation of Fourier and Lagrange basis in the harmonic balance method is an original approach. It is worth noting that the change of basis between the Fourier and the Lagrange ones permits to avoid to find directly the Fourier coefficients of $f({X}_{K}(t))$ by using (12). In particular, this is extremely useful in the case of highly nonlinear systems, such as the HH model, since formulas (12) would require huge calculations, while our technique permits to obtain the Fourier coefficients in a more efficient way.

Finally, from (10) we obtain

$\tilde D\overline{X}_K = Y_F(\overline{X}_K),$

(13)

where

$\tilde D = D\otimes I_{n}$, $D$ is the matrix of differential time operator in the base $(e_j)$:

$D = \left(\begin{matrix} 0 &\cdots &\cdots &\cdots &\cdots &0 \\ \vdots &D_1 &\ddots &\ddots &\ddots &\vdots \\ \vdots &\ddots &\ddots &\ddots &\ddots &\vdots \\ \vdots &\ddots &\ddots &D_k &\ddots &\vdots \\ \vdots &\ddots &\ddots &\ddots &\ddots &\vdots \\ 0 &\cdots &\cdots &\cdots &\cdots &D_K \\ \end{matrix} \right)$

and $D_k$ is the $2\times2$ matrix of differential time operator in the sub-base $(e_{2k}, ~ e_{2k+1}):$

$D_k =\left( \begin{matrix} 0 &k\\ -k &0 \end{matrix} \right).$

Remark 2. Generally, the choice of $K$ depends on the degree of nonlinearity. In order to detect when an accurate solution has been found, the following criterion is used:

$|| X_K-X_{K+1}||\leq \varepsilon,$

(14)

for a given $\varepsilon>0$. It is important to remark that, if we introduce the slightly different notation with respect to (10),

$\begin{align*} X_K& ={X}_{K}(t) =A_0^K +\sum\limits_{k = 1}^{K} \left(A_k^K \cos \left(k \dfrac{2\pi}{T} t\right) + B_k^K \sin \left(k \dfrac{2\pi}{T} t\right)\right)\\ {X}_{K+1}(t)& =A_0^{K+1} +\sum\limits_{k = 1}^{K+1} \left(A_k^{K+1} \cos \left(k \dfrac{2\pi}{T} t\right) + B_k^{K+1} \sin \left(k \dfrac{2\pi}{T} t\right)\right), \end{align*}$

then, in general, we have $A_k^K\neq A_k^{K+1}$ and $B_k^K\neq B_k^{K+1}$ for all $k$, and therefore (14) is not trivial.

---

2.4. Floquet multipliers

In this section we will show how the computation of the Floquet multipliers, and thus the stability analysis of the limit cycles under consideration, is straightforward when the limit cycles have been detected via the harmonic balance method [10].

Let $x(t)$ be a periodic solution of equation (1) with period $T$. The stability analysis of a periodic solution $x(t)$ can be carried out by computing the eigenvalues of the monodromy matrix $M = Z(T),$ where $Z(t)$ is the solution of the matrix differential equation

$\dot{Z} = D_x f(x(t)) Z, ~~Z(0) = I,$

(15)

$D_x f(x(t))$ is the Jacobian matrix evaluated in $x(t)$ and $I$ is the identity matrix [10,24]. These eigenvalues are also called Floquet multipliers. It is easy to see that the monodromy matrix $M$ always has an eigenvalue equal to $1$. The other $n-1$ eigenvalues of $M$ determine the stability and bifurcation behavior of the periodic orbits. In fact, $x(t)$ is stable if all the eigenvalues are inside the unit circle in the complex plane, unstable if there is at least an eigenvalue outside, and a bifurcation occurs when one eigenvalue crosses the unit circle.

Since, in general, it is not easy to solve (15), the monodromy matrix $M$ can be approximated as proposed in [10]. Let us define $A(t) =D_x f(x(t))$ and for sake of simplicity let us assume that $A(t)$ is real. Let us divide the interval $[0, T]$ into $n$ equal parts by $0 =t_0 < t_1 < \dots < t_n =T$, i.e. $h =t_{k+1} - t_k =T/n$. Notice that for periodicity the nodes $t_0 =0$ and $t_n =T$ are equivalent.Let us replace the matrix $A(t)$ with the piece-wise constant matrix $A_h(t)$ defined

$A_h(t) = A_{h, k}, ~~t \in [t_k, ~t_{k+1}), (k = 0, \dots, n-1),$

where $A_{h, k}$ is a constant matrix, satisfying $\min A(t)\leq A_{hk}\leq \max A(t),$ $t\in [t_k, ~t_{k+1})$. Usually, the constant value assumed in the left extreme of the subinterval, namely $A_{hk} =A(t_k)$, is chosen.

Thus, it is easy to see that the solution $Z_h(t)$ of the approximated system

$\dot{Z_h} = A_h(t) Z_h, ~~Z_h(0) = I$

is given by

$Z_h(t) = \exp ((t-t_k)A_{hk})\exp(hA_{h, k-1})\cdots \exp(hA_{h0}), \, \, \, \, \, t\in [t_k, ~t_{k+1}],$

so, the following approximation of the monodromy matrix $M$ is obtained:

${M_h} ={Z_h}(T) =\prod\limits_{k \in \{ 0, ..N - 1\} } {\exp } (h{A_{h, k}}).$

For more details, see [10].

We remark that when the periodic solution is found via the harmonic balance method, we have an approximated but analytical solution in the form (10). Therefore, the matrix $A_h(t)$ can be computed for any interval subdivision $0 = t_0<...<t_k<...t_{n-1}<t_n = T$ and the Floquet multipliers can be easily calculated. On the contrary, if other numerical methods, such as the collocation method, are exploited, the calculation of the Floquet multipliers is less simple, since first of all a numerical interpolation of the periodic solution is needed.

---

3. Numerical results

In this section we show how our technique based on the HB method permits to efficiently characterize the periodic solutions of the Hodgkin-Huxley (HH) model. First of all, the structure of the HH model and its main characteristics are briefly presented.

---

3.1. The Hodgkin-Huxley model

The Hodgkin-Huxley model for a neuron consists in a set of four nonlinear ordinary differential equations in the four variables $X = (V, m, h, n)$, where $V$ is the membrane potential, $m$ and $h$ are the activation and inactivation variables of the sodium channel and $n$ is the activation variable of the potassium current. The corresponding equations are the following [15,18]:

$\left\{ \begin{array}{l} C\dfrac{dV}{dt} = -I-[(V-E_{Na})\overline{g}_{Na} m^3 h+(V-E_{K})\overline{g}_{K} n^4+(V-E_{L})\overline{g}_{L}], \\ \\ \dfrac{dn}{dt} =\alpha_{n}(V)(1-n)-\beta_{n}(V)n, \\ \\ \dfrac{dh}{dt} =\alpha_{h}(V)(1-h)-\beta_{h}(V)h, \\ \\ \dfrac{dm}{dt} =\alpha_{m}(V)(1-m)-\beta_{m}(V)m, \end{array} \right.$

(16)

where $I$ is the external current stimulus, $C$ is the membrane conductance, $\overline{g}_{i}$ are the shifted Nernst equilibrium potentials, $E_i$ are the maximal conductances, $\alpha(V)$ and $\beta(V)$ are functions of $V$, as follows:

$\begin{array}{lcl} \alpha_{n}(V) = 0.1 \text{expc}(0.1 (10+V) ),&& \beta_{n}(V) = \exp(V/80)/8, \\ \alpha_{h}(V) =0.07\exp(V/20), && \beta_{h}(V) =1/( 1+\exp( 0.1(30+V)) ), \\ \alpha_{m}(V) = \text{expc}(0.1(25+V)),&& \beta_{m}(V) = 4 \exp(V/18), \end{array}$

and $\text{expc}(x)$ is given by [15]

$\text{expc}(x) = \left\{ \begin{array}{ll} \dfrac{x}{\exp(x)-1}&\text{if}\, \, \, \, x\neq 0\\\\ 1&\text{if}\, \, \, \, x = 0. \end{array} \right.$

Finally, the typical values for the other parameters are [13,15]:

$E_{K} = 12 \;mV, ~~E_{Na} = -115\; mV, ~~E_{L} = -10.599 \;mV$

$\overline{g}_{K} = 36 \;mS/cm^2, \overline{g}_{Na} = 120 \;mS/cm^2, \overline{g}_{L} = 0.3 \;mS/cm^2.$

For small values of the current stimulus $I$ the system exhibits a stable equilibrium point. If $I$ is increased, then a stable periodic solution with large amplitude appears, while the equilibrium point remains stable. This means that necessarily unstable solutions are present, in order to separate the two basins of attraction. In [15], the author shows that, depending on the value of $I$ the HH model presents from one to three unstable limit cycles. Moreover, in a certain range of values of $I$ the equilibrium point becomes unstable, but finally the stable periodic solution disappears through a Hopf bifurcation and the equilibrium point regains its stability.

At present, few works about the detection of the periodic solutions of the HH model and their related bifurcations exist in literature, for example [9,13,14,15,32], because of the high dimension of the system and its high nonlinearity. Furthermore, the existing works approach the problem by exploiting different methods (finite differences, collocation or shooting methods), that are not so simple to handle with.

We show that through our technique based on the harmonic balance method [2,27] it is possible to efficiently characterize and numerically approximate the periodic solutions exhibited by the HH model (16) and the related bifurcations, depending on the intensity of the external current stimuli $I$.

In particular, we point out how the harmonic balance method can be efficiently exploited in the route-to-chaos region found by [13], that is in the region of the Hodgkin-Huxley bifurcation diagram where the dynamics is more awkward and thus more interesting. This method permits to obtain analytically the stable and unstable periodic solutions of the Hodgkin-Huxley model, therefore we get a better characterization of the unstable chaos and of the Hopf bifurcations. In fact, the Floquet multipliers can be easily computed and the flip bifurcation effortlessly detected. In particular, we show that in the regions close to the Hopf bifurcations we are able to obtain the unstable periodic solution with just two harmonics (thus 12 unknowns variables).

---

3.2. Continuation technique and Newton's algorithm

In this work, we are interested in detecting periodic solutions of HH system, depending on the values of the external current $I$. We will use a continuation method where at each iteration we fix the parameter $I$ and determine the unknown periodic solution $X(I)$ and its period $T(I)$, starting from the solution found at the previous step.

Therefore, the numerical computation of a periodic solution is based on the resolution of the nonlinear algebraic system

$\label{eq:nonlinear} F(X, T, I) = 0,$

(17)

issue of one of the methods presented above (mainly HB method), where $X$ is the state vector, $T$ is the period approximation, and $I$ is the parameter.

The continuation method implementation exhibits two main difficulties: on the one hand, the construction of a right initial solution, and on the other hand the progression of the algorithm for critical values of parameter $I$, that is the turning points. It is possible to ride out the first one by starting from a stable limit cycle branch, that can be suitably numerically approximated with a fourth-order Runge-Kutta method. For the second one, Chan and Keller [7] proposed the arc-length continuation method. In our case, since in a neighborhood of the turning points the function $T(I)$ is monotone [32], we can use $T$ as parameter in order to follow the branch continuation.

Let us remark that, in our case, the branches are locally linear, so, for suitable choices of $T$ and $I$ in our algorithm, the next point is obtained via a correction and by using Newton method [7].

Our codes are available at: http://lmah.univ-lehavre.fr/codes/codes.html

---

3.3. Bifurcation diagram and branches of periodic solutions

Both collocation and harmonic balance methods are appropriate for detecting all the periodic solutions, both the stable and unstable ones. The multi-base harmonic balance method uses Fourier series, therefore we have an analytical expression of the solution and the convergence rate is of exponential order. In our case, in general 50 harmonics are enough to get the desired accuracy (see Remark 2), so we have at most only 405 unknowns variables.

In Fig. 1 the bifurcation diagram for the HH model has been obtained by jointly and optimally exploiting the three methods presented above, that is shooting, collocation and multi-base HB methods. The stability analysis of the detected limit cycles has been carried out by the calculation of the Floquet multipliers, by applying the numerical algorithm proposed in [10] to the approximated solution.

It is possible to see that the dynamical behavior of HH system can be decomposed in two main regions, depending on the value of the external current $I$. For $I > I_{2} = 9.73749234$, there is one equilibrium point and one stable periodic solution, that disappears through a Hopf bifurcation for $I = I_1 = 154.500$. Moreover, there is a second region of the bifurcation diagram, for $I \in~[0, ~I_{2}]$, that is more interesting, since its dynamical behavior is more complex and rather less understood. A zoom of the diagram for $I \in~[0, ~I_{2}]$ can be found in Fig. 2. It is possible to see that in this second part of the diagram the system undergoes three saddle-node of cycles bifurcations at $I_3 = 7.92198549$, $I_4 = 7.84654752$ and $I_5 = 6.26490316$. They correspond to the knees of the bifurcation diagram and they consist in the collision and disappearance of two periodic solutions. Moreover, the system exhibits a period-doubling bifurcation [32] at $I_6 = 7.92197768$ that can be detected by the joint application of the harmonic balance method and the Floquet analysis.

For $I$ close to $I_5$, the periodic solutions detected by the multi-base HB method exhibit a sort of Gibbs phenomenon [36], as it can be seen in Fig. 3, and this does not permit to accurately detect the saddle-node of cycles bifurcation. Therefore only in this region shooting and collocation methods have been used in order to find the stable and unstable periodic solutions, respectively. The remaining of the diagram has been found via the multi-base HB method, by choosing and controlling the minimal number of required harmonics (see Remark 2). In particular, for $I \in~[0, ~I_{2}]$, 50 harmonics have been considered for the approximation of the stable limit cycle with higher amplitude, while for $I>7$ the unstable limit cycles required only 30 harmonics, and for $I > 8$ the number of harmonics can be gradually reduced since the unstable limit cycle becomes more and more regular. Finally, close to the Hopf bifurcation at $I_2$ only one harmonics is sufficient to get the best approximation of the unstable periodic solution.

Therefore, we can conclude that multi-base HB method works very well in the region between $I_4$ and $I_2$, that is in the richest and most interesting part of the bifurcation diagram. Indeed, it permits to carry out bifurcation analysis in a more straightforward way.

In the following, we analyze more accurately the various limit cycles bifurcations that take place.

---

3.4. Analysis of the limit cycles bifurcations

Hopf bifurcations. In this paragraph, we are interested in Hopf bifurcations, that take place at $I = I_{1} = 154.500$ and $I_{2} = 9.73749234$. A view into $(V, n, I)$ space projection of stable and unstable periodic solutions, in a neighborhood of the two Hopf bifurcations, are shown in Fig. 4. These results suitably match with the theoretical results proved in [13,15].

It is worth noting that in both cases over a large interval of $I$ close to the Hopf bifurcations, the periodic solutions are almost sinusoidal (see Fig. 5). Therefore, only one or two harmonics are needed to conveniently approximate this solution via the harmonic balance method.

Saddle node of cycles bifurcations. In our case, there are two types of saddle node of cycles bifurcation: for $I =I_5 = 6.26490316$ we have a simultaneous appearance of two limit cycles (one stable and the other unstable), while at $I =I_3 = 7.92198549$ and $I = I_4 = 7.84654752$ we have the collision of two unstable periodic solutions (see Figs. 6, 7 and 8). For detecting such bifurcations, we use the Floquet analysis, by searching when an additional Floquet multiplier crosses the unit circle in $+1$.

The Floquet multipliers for these three cases are represented in Fig. 9. It is possible to see that a multiplier leaves or enters in the unit circle through $+1$.

Period doubling bifurcation. Finally, in this section, we consider the period-doubling bifurcation. By exploiting the Floquet analysis, we can easily detect this bifurcation since in this case a Floquet multiplier crosses the unit circle through $-1$, as it is shown in Fig. 10. Table 1 shows the values of the Floquet multipliers for different values of $I$. As $I$ tends to $I_6 = 7.92197768$, the second most negative Floquet multiplier tends to $-1$.

---

4. Conclusion

In 1952 Hodgkin and Huxley developed the pioneer and still up-to-date mathematical model for describing the activity of the squid giant axon. Depending on the value of the external current stimuli, this fourth-order nonlinear dynamical system exhibits many complex behaviors, such as multiple periodic solutions (both stable and unstable) and chaos.

Previous works have treated this problem by using several numerical methods, such as shooting and finite difference methods, that are not so simple to handle with. In this paper, we propose a multi-base HB method, a technique based on the harmonic balance method, permitting to detect the stable and unstable periodic solutions and the associated bifurcations of a nonlinear dynamical system. In particular, we have shown how our multi-base harmonic balance method is extremely handy, permits to obtain the bifurcation diagram of the HH model, and works very well in the most complex part of such diagram. Furthermore, harmonic balance and Floquet analysis have permitted to suitably detect the period-doubling bifurcation that entails a route-to-chaos in the HH model.

---