Nonlinear Dynamics, Chaos and Control of the Hindmarsh-Rose Neuron Model

Mathematics has changed over time to comprise interdisciplinary fields of research, and considering this, biomathematics has arisen as an interface study. In this work, we analyze the dynamical behavior of the Hindmarsh-Rose neuron model, which describes the neuronal bursting in a single neuron. A stability study through the Lyapunov exponents method is proposed and evidence of a chaotic dynamics is presented. This chaotic behavior is biologically comparable to an individual undergoing an epileptic seizure, in which the application of an efficient controller represents a proposal for preventing epilepsy from happening. Therefore, a control design based on the State-Dependent Riccati Equation is proposed aiming to reduce the oscillation of the system to a desired orbit. The results show that the controller is efficient and robust as a method for preventing epileptic seizures.


Introduction
There is a great motivation, from the bioengineering point of view, in studying the Hindmarsh-Rose neuron model, because it represents well a biological neuron and is, thus, capable of simulating several behaviors of a real neuron; for example: periodic, aperiodic and chaotic behaviors. In the literature, the chaotic behavior present in the Hindmarsh-Rose neuron model is biologically comparable to an individual, or patient, with epilepsy.
All conscious and unconscious activities are controlled by the nervous system, which is formed by billions of nerve cells receiving information that comes from inside and outside of the human body. Excitability and conductibility are physiological properties characteristic of nerve cells and also muscle cells; these cells are thus able to react to a given stimulus and transmit it as electrochemical impulses through their membranes [2]. The membrane of nerve cells has electrical properties that are defined by the molecular organization of their components, and these properties may affect the capacity of electrically excitable cells in terms of information delivery. When a neuron receives a stimulus, whether it is chemical, mechanical or electrical, there may be a change in the membrane permeability, thus allowing for a high input of sodium and low output of potassium in the cell. This occurs along with an inversion of the charges around the membrane, making it depolarized, which in turn leads to an action potential. This difference in potential varies according to the number of ions in the intracellular medium and has a steady value of -70 mV.
In the literature, the neuron membrane is often compared to an electrical circuit. Mathematical models were proposed to study the behavior of this membrane in its quiescent state and the flow of nerve impulses. Among them, we highlight the Hodgkin-Huxley model [9], which contains a mathematical 2 R. S. Lima and F. R. Chavarette description for the behavior of the nerve membrane, fiber conduction and excitation, and the Hindmarsh-Rose model [7,8], which presents a modification of the classical model of FitzHugh and Nagumo [6,14], whose aim is to describe neuronal bursting generation.
In mechanical engineering, there is a great interest in controlling chaotic behaviors. In practice, this implies in developing state control laws, such as feedback control [11], in order to stabilize the system around the unstable equilibrium points. Considering this, many control techniques have been developed with the purpose of controlling chaotic systems, and each of these techniques has its peculiarities and different results. To control such chaotic dynamics, a control design using the State-Dependent Riccati Equation (SDRE) is applied. The objective is to simulate the chaotic behavior of the Hindmarsh-Rose model [7] using a 4 a order Runge-Kutta numerical integrator [10], and ensure, through the Lyapunov exponents method, that a chaotic behavior exists.

Mathematical Model
The nerve impulse is characterized by a propagation of such depolarization through the neuron and thus, after the impulse passes, the membrane undergoes a repolarization, recovering its normal quiescent state and ending impulse transmission [2,19].  The Myelin Sheath acts as a protection for the Axon and helps accelerate the conduction of nervous impulses. Figure 1 (b) presents the electrical Hodgkin-Huxley model, where E k , E N a and E Cl , represent the Nernst Potential 1 for ions of potassium, sodium and chloride, respectively. Based on the squid giant axon, Hodgkin and Huxley [9] used values of -12, 115 and 10.6 mV for E k , E N a and E Cl , respectively (it must be stressed that they had assumed the membrane difference in potential to be measured relatively on the inside of the cell, and thus the sign convention was the opposite to that adopted since the beginning of their work). C = 1µF cm 2 represents the membrane capacitance. The parameters g k , g N a and g Cl represent the conductance for potassium, sodium and chloride, respectively. These parameters correspond to the resistive elements, and because g k and g N a are tension-dependent, they were defined as variable resistors.
Bursting generation has been extensively studied in the context of the Hindmarsh-Rose model, which establishes a dimensionless state variable for the membrane potential x(t), and other two (also dimensionless) variables, namely y(t), associated with the fast ion flows generated by the transfer of Na + e K + and z(t), which describes the dynamics of other low channels. These variables are nonlinearly related in the composition of the membrane potential, which can be translated by the following dynamics [7]: where a, b, c, d, s, r, x r and Iare parameters of the system, which, depending on their values, can simulate a wide range of dynamical behaviors that are topologically equivalent to the ones observed experimentally. This makes the Hindmarsh-Rose model one of the most emblematic models concerning the qualitative study of neuronal bursting [1]. The system presents many behaviors, and one of them is a typical chaotic behavior of the membrane potential, that is, a behavior that is aperiodic and sensitive to the initial conditions adopting the external stimulation (variable I) as a control parameter.

Numerical Simulations
As the differential equations that describe the model (2.1) are simple, involving only polynomials with the variable x, they can be simulated using analogical electronic circuits that are composed of tension multipliers and dividers, and integrators. Such a model may be called an analogical computer or an electronic neuron (EN). ENs based on model (2.1) were used with the purpose of restoring the functioning of a damaged biological neural circuit. These damages may be abnormal electric brain discharges, which characterizes a different behavior compared to its ideal condition, such as epileptic seizures, which present the same nature as nerve impulses [5,3,15].
The fixed parameters used were b = 3, c = 1, d = 5, s = 4, x r = −1.56, I = 3; the varying ones, in turn, were a = 0-2 and r = 0-0.2; also, the initial conditions were x = 0.1, y = 0.2 and z = 0.1 [7], leading to the diagram in Figure 2, which represents the stability between the variation of membrane potential and the fast ionic flows generated by the transfer of Na + and K + , and Figure 3, which represents the stability associated with the fast ionic flows generated by the transfer of Na + and K + and the dynamical description of slow channels.

Control Design via State-Dependent Riccati Equation
The SDRE nonlinear regulator has the same structure as the linear quadratic regulator (LQR), except that all the matrices are state-dependent. The SDRE approach for obtaining a suboptimal solution of the control problem has the following procedure [4] and [13]: 1. Represent the model in state-space form. Use direct parametrization to bring the nonlinear dynamicsẋ = f (x) + g(x) to the state-dependent coefficient (SDC) form, as follows: x ∈ R n is the state vector, u ∈ R m is the control law, y ∈ R s is the output vector.
2. Define the initial conditions x(0) = x 0 , and choose the coefficients of the positive definite weighting matrices Q(x) and R(x), which determine the relative importance of state x(t) and control effort u(t), respectively.

Solve the state-dependent Riccati equation given by:
4. Construct the nonlinear feedback control via: R. S. Lima and F. R. Chavarette

The control law (4.3) is calculated so that the performance index given by
In the multivariable case, there always exists an infinite number of SDC parameterizations. Therefore, the choice of the matrix A(x) is not unique [4].
The pair {A(x), B(x)} is a controllable parametrization of the nonlinear system in a region Ω if {A(x), B(x)} is pointwise controllable in the linear sense for all x ∈ Ω. Therefore, the choice of A(x) must be such that the state-dependent controllability matrix [B(x) A(x)B(x) ... A n−1 (x)B(x)] has full rank [13].
Applying the above procedure in the nonlinear system [7], we obtain: Assuming C(x), Q(x) and R(x) are constant matrices, we have: The responses of the controlled and non-controlled system are shown in Figures 6 and 7.

Results and Discussion
Firstly, the structural stability of the system was analyzed based on the eigenvalues of the Jacobian matrix [12,10] by varying the parameters a, r and s, and considering the remaining parameters fixed, which led to the stability diagrams. The points (•) in Figures 2 and 3 represent the linear stability of the system, and the points ( * ) represent the instability of the system. A random pair (a = 1 and r = 0.15) was selected among the stability points to carry out a thorough analysis about the behavior of an epileptic seizure, because it is essentially abnormal brain discharges that demonstrate how the nervous system works under different unstable conditions. Afterwards, by analyzing the dynamics of the Lyapunov exponents by the Wolf Method [18], we verified whether the system is indeed unstable. Figures 4 and 5 depict how the exponents behave over the sampling dimensionless time.   The Lyapunov exponents obtained were approximately, according to Table 1, λ 1 = 0.010639, λ 2 = 0.0086279. A system is said to be stable if none of the Lyapunov exponents is positive. The existence of a positive Lyapunov exponent characterizes the chaotic behavior of the system. Furthermore, considering that two positive Lyapunov exponents were found, the Hindmarsh-Rose model is considered hyperchaotic [16]. Based on the results presented, the presence of chaos for all the variation of the membrane potential was verified, since the depolarization process until repolarization, representing an individual undergoing an epileptic seizure; this fact highlights the need of applying a control law to stabilize the whole process. Figure 7: Phase Portrait -(a) x vs y (non-controlled system) ; (b) x vs y (controlled system); (c) x vs z (non-controlled system); (d) x vs z (controlled system); (e) y vs z (non-controlled system); and (f) y vs z (controlled system).

Conclusions
In this work, the dynamics of a model called the Hindmarsh-Rose System with chaotic behavior is analyzed through the Lyapunov exponents. A positive exponent was found, which demonstrates the existence of a chaotic behavior and characterizes an epileptic seizure. To stabilize such chaotic behavior, a control method based on the State-Dependent Riccati equation was implemented as a strategy for reducing the chaotic motion of the system to a desired small periodic orbit. The choice for a SDRE controller is due to the fact that it leads to less computational cost, as well as a faster convergence, if compared to feedback control. Figures 6 and 7 depict the effectiveness of the control strategy applied to solve these problems. The potential results of this work may be used, in future research, on mechanisms of action for new antiepileptic drugs.
It must be stressed that the results presented in our work can be used in similar or correlated areas, such as Artificial Intelligence, Pharmacology and Electronic Circuits Technology, for example, with the purpose of developing microchips and new medicines for preventing epilepsy.