A STUDY ON A NEW BIFURCATION PARAMETER IN A MODIFIED HUBER-BRAUN MODEL

Original scientific paper Thermally sensitive neurons represent a bursting-spiking activity that is indicated by rapid repetitive spiking trains of action potentials pursued by dormant periods. Synchronization of such behavior in a network of coupled spiking neurons such as the epileptogenic zone in the brain may cause some neurological disorders such as epileptic seizures. This paper introduces a new approach for predicting the seizure onset in a model of an epileptic neuron. The parameters which are used for simulations have been selected in such a way that they would be potentially applicable in the non-invasive brain stimulation approaches such as repetitive transcranial magnetic stimulation (rTMS). In this regard, a modified Huber-Braun model of a thermally sensitive neuron exposed to external rTMS-induced voltages is presented. Applying the chaos theory, the bifurcation diagram of a modified Huber-Braun model with a new bifurcation parameter is used to estimate the time at which the bifurcation takes place whereby allowing a more accurate prediction of the seizure onset based on the modified model.


Introduction
Epilepsy affects more than 60 million individuals worldwide [1]. This disease is caused by unexpected ephemeral discharges in neurons of the brain. An epileptic seizure is a synchronous neuronal action. The duration of such behavior may vary from a few seconds to several minutes [1]. For approximately two thirds of the patients, seizures can be suppressed by anticonvulsive drugs [2]. Epilepsy surgery may treat 7 % to 8 % of patients [1]. But for approximately 25 % of epilepsy patients, the seizures are resistant to these accessible control approaches and treatments. The epilepsy of this group of patients is called refractory epilepsy [2]. Recently there have been considerable efforts to predict and control the seizures in the patients with refractory epilepsy.
Scientists have studied different seizure onset patterns and have tried to detect and even control the seizures by extracting the features related to the seizure onset based on responses of the neural models of the epileptogenic zone or the EEG records of the epilepsy patients [2,3].
One of the strong mathematical methods that has been used for pattern recognition and feature extraction in the studies on the neurological disorders is the bifurcation theory [4,5]. In this theory, the changes that happen in the qualitative behavior or topological structure of a given family, such as solutions of a family of differential equations, are mathematically studied. These changes are produced by varying some parameters of the systems. Those parameters are called bifurcation parameters. The seizure detection and control methods based on the features related to bifurcation theory have been presented in many studies. For example, the correlation dimension, the largest Lyapunov exponent [4], the detection of bifurcation in in-vitro neuronal models and in the EEG records of epilepsy patients [5].
The seizure control strategies are based on stimulating the epileptogenic zone in focal epilepsy. These stimulation-based therapies are divided into several categories such as invasive and non-invasive as well as responsive and non-responsive. The non-responsive or open loop stimulation means applying the stimulus in a certain predefined time interval without any seizure detection. But during a non-responsive or closed loop stimulation-based therapy, the seizures are detected as they start and then bursts of stimuli are triggered by the stimulator [6].
Several invasive responsive and non-responsive stimulation techniques as well as neuro-stimulators such as responsive neuro-stimulation (RNS) [1], nonresponsive Vagus nerve stimulation (VNS) [7] and deep brain stimulation (DBS) [8] have been proposed. These invasive stimulations are applied to the epileptogenic zone in the form of a focal electrical stimulation [9], focal drug delivery [4], or focal cooling [10].
The responsive non-invasive stimulation-based therapies are considerably painless in comparison with the invasive stimulation approaches [11]. Therefore, they have attached attention recently. Researchers are trying to introduce more sensitive seizure detection methods and to design novel responsive neuro-stimulators. Among the non-invasive stimulation methods, transcranial magnetic stimulation (TMS) and repetitive transcranial magnetic stimulation (rTMS) have been used experimentally to suppress seizures in epilepsy patients [11,12]. It has been shown that applying the low-frequency high intensity rTMS on the focus of an epileptic seizure has a meaningful antiepileptic result on patients [13]. Besides, the rTMS therapy can also decrease the frequency of the epileptiform discharges, whereby improving the psychological status of these patients [13].
The sensitivity of the aforementioned non-invasive stimulation-based therapies has not been acceptable for clinical applications yet [11]. That is because of the lack of in-vivo models to examine the efficacy of these stimulations and neuro-modulatory therapies. Hence, more investigations are still needed for modification of the existing in-vitro models and to present more applicable and reliable control systems to control the seizures in these models. Furthermore, better features for prediction algorithms should be extracted.
One of the neuronal models which is used to simulate the bursting-spiking behavior of an epileptic neuron is the Huber-Braun model of thermally sensitive hypothalamic neurons [3]. This model exhibits a wide range of different impulse sequences as a function of the temperature [14]. Several studies have been performed in order to find the bifurcation parameters in this model [15], or to control a network of Huber-Braun neurons using an invasive stimulation protocol [3]. In order to use this model to study the non-invasive stimulation techniques, some modifications are needed. For instance, the induced electromagnetic forces and induced potentials due to the magnetic stimulation should be applied on the membrane of the neuron and their effects on the responses of the model should be considered.
This paper presents an extension of the Huber-Braun model, simulates a Huber-Braun neuron exposed to magnetic stimulation induced voltages, introduces a new bifurcation parameter in the modified model and rebuilds the resulted bifurcation diagram.
The rest of this paper is organized as follows.
The proposed modified Huber-Braun model is presented in Section 2. Section 3 deals with the results and discussions and finally the conclusion is presented in Section 4.

Materials and methods
Huber-Braun model has been developed originally from the Hodgkin-Huxley (HH) model to imitate temperature dependent alterations of static spike train designs recognized in electroreceptors of hypothalamic neurons of rat [5]. Huber-Braun neurons are described by noisy impulsive vibrations which are mirrored in the spike sequences. The deactivation terms in the classical HH model are ignored in the Huber-Braun model but instead, the primary model is lengthened by two extra slow currents. These slow currents are corresponding to the preliminary discoveries of bursting behavior in the biological neurons [14]. The equation describing this model is: where m and are the capacitance and potential of the membrane, respectively, and is the leakage current that can be calculated as and are fast Hodgkin-Huxley sodium and potassium currents, respectively, which account for spike generation.
represents the slow current that is mostly generated by sodium ions and to a less extent, calcium ions.
is the slow calcium-dependent potassium current. The last two currents explain the oscillatory behavior of the model. The amount of these ionic currents may be obtained using the following equation: where ( ) = 1 0 is a temperature dependent scaling factor, is the maximum conductance of membrane and stands for the Nernst potential for each ionic current. The variable represents the activation of each ionic channel of the membrane and for the currents and , is the solution of the following differential equation: where ( ) = 2 0 is the scaling factor for each activating variable, is the time constant and ∞ is the steady state activation variable which is governed by the following equation with being the sharpness and 0 being the halfactivating voltages. Since the sodium current changes are faster than other ionic currents, = ∞ . Therefore, The variable which is the activating variable of the current is described as where represents the +2 ions increased followed by the passage of current and shows the active elimination of intracellular calcium ions. The values used for the parameters of this model are the same as the ones in [15] and are listed in Tab. 1 for easier reference.
Several studies focused on investigating the behavior of the Huber-Braun model in response to changes in its parameters. For example, different simulations have been performed in both deterministic and noisy conditions [3] [14÷16]. The temperature attributes of the model also have been studied and a wide range of different temperature-depending spike patterns in the model responses has been observed [14]. Altering the temperature, the Homoclinic bifurcation was found in the classic model and the resulting bifurcation diagram was rebuilt [15]. Using the external current injected to the model via an invasive electrical stimulation, a period doubling bifurcation has been observed in the model and its bifurcation diagram for the external current has been produced [16]. In a recent study, a network of Huber-Braun neurons has been simulated and some control protocols using a feedback current have been investigated [3].

Modified Huber-Braun model
In order to modify this classic model to be applicable in the non-invasive stimulation techniques, an external voltage, , was added to the membrane voltage which represents the rTMS-induced voltage. Therefore, the nerve voltage is modeled as: Replacing by in the classic equations, yields the new set of modified Huber-Braun model equations.
The membrane potential is the most important variable through which the different behaviors of the classic model can be seen. The membrane voltage spikes are responsible for transmitting the neural massages. Thus, the inter-spike interval (ISI) which is the time between two successive spikes, is one of the important measures that should be calculated. Additionally, depending on the changes in temperature, the bifurcations occur in the values of inter-spike interval [15,16]. According to the studies on the effects of transcranial magnetic stimulation on spike trains of epilepsy [11] [17], we estimated the peak to peak amplitude of rTMSinduced voltage to be around 4 mV. This was based on the fact that the peak to peak amplitude of the motorevoked potentials during magnetic stimulation has been calculated to be equal to 2 mV [11]. As shown in the next section, this selected range is suitable for investigating the effects of different induced voltage magnitudes on the neural responses.

Results and discussion
By numerically simulating the suggested modified Huber-Braun model with MATLAB, the voltage intervals of this output were calculated. Then, by altering the value of the electromagnetically induced voltage as a bifurcation parameter, different behaviors of the membrane potential were explored. These behaviors are tonic firing, bursting-spiking discharges, chaotic, subthreshold oscillations and steady state which are related to different states of the neuron. Some examples of these behaviors are illustrated in Fig. 1. Moreover, the bifurcation diagram was successfully produced for the induced potential as shown in Fig. 2. Depending on the magnitude of the induced voltage, three different areas can be observed in the bifurcation diagram. As can be seen in Fig. 2 and Fig. 3(a), when is smaller than -0,72 mV, the neuron experiences a tonic firing behavior which means that it produces a single spike per each period. When increases up to 0,7 mV, the number of spikes per each period increases and for the induced voltages from 0,7 to 1,385 mV the number of spikes of membrane potential decreases. In this region, the interval between every two successive spikes remains constant. In other words, the trajectory of V moves twice along its path at its previous speed and the neuron represents a bursting-spiking behavior. This kind of bifurcation is called period doubling (see Fig. 2 and Fig.  3) and it lasts until = 1,385 mV.
At the other end of the zone, for the magnitudes above 1,385 mV, the neuron stops firing and no more spikes will be generated.
The modified membrane voltage starts to represent sub-threshold oscillations and soon afterwards it reaches a steady state.

Figure 4 The largest Lyapunov exponents
During the time that the induced voltage is diminished, at first we can see a distinct spike during each period then the bursting-spiking discharges come into view and the number of spikes grows. Between these two areas, a chaotic behavior is observed.  Fig. 4 provides information about the largest Lyapunov exponents (LLE) of the system responses with respect to the electromagnetically induced voltages on the membrane surface. As indicated in Fig. 4, after = 1,385 mV, the LLEs of the system become positive. This means that the neuron experiences a chaotic behavior at this zone. As can be seen in Fig. 5, when = 1,385 mV, even very small variations of will lead to perceptible qualitative variations in the dynamic behavior of the system. In other words, the membrane response patterns will vary tremendously from representing a single spike per period to producing sub-threshold oscillations. Since the period doubling bifurcation occurs before the chaotic behavior and synchronization of neuronal discharges in epilepsy, it seems that detecting the starting time of bifurcation could be a cue that a seizure will happen. Furthermore, by introducing as an external controllable bifurcation parameter, it could be possible to desynchronize and control the chaotic behavior of the neuron using this parameter and cause the bifurcation point to move to an unreachable point at the same time.

Conclusion
In this study, the Huber-Braun model of thermally sensitive hypothalamic neurons of rat was modified by adding the rTMS-induced voltage to it. Then, different behaviors of the model response were studied based on alteration of this parameter and a bifurcation diagram of this parameter was successfully generated. Since the modified model's behaviors are similar to pathological cell's action potentials which are caused by diseases like epilepsy, the results may be used in prediction and control of epileptic seizures by a responsive non-invasive neurostimulator. In other words, the results of this study can be used in designing an epileptic seizure onset detector as well as a magnetic neuro-stimulator to control the epileptic seizures. The onset of the bifurcation of the membrane voltage could be considered as the proper time for a closed loop neuro-stimulator to trigger a feedback control signal to suppress the synchronization of the neural activity in the epileptogenic zone. For future works, the parameters of cable equation of a Huber-Braun nerve fiber exposed to a magnetic stimulator could be studied. Furthermore, in order to obtain more precise and reliable results, the distribution of the electromagnetically induced voltages with respect to time, position and current in the stimulating coil of a magnetic stimulator should be considered.