Firing activities of a fractional-order FitzHugh-Rinzel bursting neuron model and its coupled dynamics

Fractional-order dynamics of excitable systems can be physically described as a memory dependent phenomenon. It can produce diverse and fascinating oscillatory patterns for certain types of neuron models. To address these characteristics, we consider a nonlinear fast-slow FitzHugh-Rinzel (FH-R) model that exhibits elliptic bursting at a fixed set of parameters with a constant input current. The generalization of this classical order model provides a wide range of neuronal responses (regular spiking, fast-spiking, bursting, mixed-mode oscillations, etc.) in understanding the single neuron dynamics. So far, it is not completely understood to what extent the fractional-order dynamics may redesign the firing properties of excitable systems. We investigate how the classical order system changes its complex dynamics and how the bursting changes to different oscillations with stability and bifurcation analysis depending on the fractional exponent (0 < α ≤ 1). This occurs due to the memory trace of the fractional-order dynamics. The firing frequency of the fractional-order FH-R model is less than the classical order model, although the first spike latency exists there. Further, we investigate the responses of coupled FH-R neurons with small coupling strengths that synchronize at specific fractional-orders. The interesting dynamical characteristics suggest various neurocomputational features that can be induced in this fractional-order system which enriches the functional neuronal mechanisms.

Collective oscillatory dynamics and synchronous activity are the fundamental phenomena in dynamical systems 1,2 . It has both theoretical importance and biophysical significance in computational neuroscience. Mathematical biophysical models [3][4][5][6] are the primary tools to characterize the nervous system. The foremost exciting step in neural dynamics is to understand the system architecture of individual neurons in terms of mathematical models of membrane potential. Various types of spiking and bursting are the dynamical responses of excitable cells 7,8 . Such type of research analyzes the chaotic behavior of excitable systems. When mathematical models are described as a single neuron or network, then the nonlinear dynamical techniques are applied to study the emerging oscillatory patterns and the synchronization phenomena. The classical order dynamical models depend on the immediate previous response, however the fractional-order derivative depends on all the previous responses, so it has a memory effect [9][10][11] . It can produce a different kinds of multiple time scale neuronal dynamics 12 . It [13][14][15][16] provides a wide range in understanding the rich dynamical and neuronal responses. Fractional-order calculus originated from a letter written to Leibnitz by L'Hospital [13][14][15] . Now, it has become a promising and reliable mathematical tool that includes hereditary properties or memory dependence phenomena [13][14][15][16][17] . The discussion of multiple time scale dynamics has been studied in some previous articles 18,19 which have the potential importance in signal processing. It has been studied that the firing rate of neocortical pyramidal neurons with the injected applied sinusoidal current can be well approximated with fractional-order derivative 20 .
Many researchers have worked on the fractional-order dynamical systems 21,22 . Results show that it follows power-law dynamics [23][24][25] . In human memory, the power-law dynamics was investigated earlier 26,27 , the accuracy of memory dependence decays at a rate nearly equal to t −α where α < < 0 1 . The power-law adaptation helps in describing some dynamical behavior of biological systems 12,28 . In recent years, fractional-order derivative has become very useful in modeling biological phenomena [13][14][15][16]29 , viscoelastic properties of tissue 30 , tissue electrode interface 31 , the kinetic property of drug delivery 32,33 , diffusion process [34][35][36] , biophysical neuron models and neural networks [37][38][39][40][41][42] . It has been found that cognitive behaviour can be modelled using fractional dynamics 43 . It was observed that fractional-order dynamics is used in vestibular oculomotor system 44 and fly motion of sensitive neurons H1 28 . It can include the mechanism of synapses 28 and the geometrical properties of excitable cells 44,45 . Single neuron models are analyzed by using fractional-order dynamics such as Hodgkin-Huxley (H-H), Morris-Lecar (M-L), FitzHugh-Nagumo (FHN), Hindmarsh-Rose (H-R) models, etc 10,39,[46][47][48][49] . In this study, it has been demonstrated that the fractional-order dynamics of a solitary nerve membrane can be analyzed by using a suitable biophysical model that exhibits elliptic bursting. It reflects the rate of change of information through membrane voltage that leads to previous history-dependent activities. Different parameter regimes corresponding to qualitatively different dynamical properties are analyzed. This article investigates the FH-R neuron [50][51][52] in a fractional-order domain to get a general idea about the spike patterns and spike frequency with stability and bifurcation scenarios. The relation between spiking and bursting is a significant question as well as a fascinating phenomenon in mathematical neuroscience, especially in neural coding. Bursting presents a recurrent transition between repetitive spiking and quiescent state. The switching phases depend on the strength of the slowly changing current stimulus to the dendrite. An exciting feature of the elliptic bursting is that the frequency of emerging spiking activity and ceasing the spiking is nonzero; at that time, the amplitude of the oscillations may be small 53 . It was experimentally studied that this type of bursting can be found in trifacial nerves controlling the jaw movement of rodents 54 .
It has been previously found that the slow variable in the 2D FHN model creates mathematical complexity that allows various dynamics for the membrane voltage of the neuron model including chaos. Therefore, fractional-order FH-R single neuron model has an excellent qualitative feature that exhibits many diverse oscillations of action potentials. Necessary and sufficient conditions are investigated for asymptotic stability analysis of the fractional-order commensurate FH-R model. Bifurcation shows the qualitative changes between the quiescent state and the oscillatory state 5,6 . The fractional-order FH-R model is investigated with a certain fixed-parameter sets, and extreme numerical computations are derived for examining the dynamical characteristics with analytical analysis using the fractional exponent as a predominant parameter. As a consequence, this generalization of the classical order model can produce biophysical variability. We present the effect of the fractional-order dynamics on the synchronization criterion in different ensembles for coupled oscillators. We observed different dynamical behavior with various fractional-orders that were not present in the classical order model.

fractional-order fH-R Model
FitzHugh and Rinzel introduced FH-R model (1976, in an unpublished article) 50,52,53,55 , which is the modification of the classical FHN neuron model. The 2D FHN model 3,4 illustrates a geometrical explanation of interesting biophysical phenomena that are relevant to neuronal excitabilities and spike generation. It exhibits continuous spiking with a specific external stimulus. However, it is not capable of generating various fascinating firing patterns produced in cortical neurons. FH-R neuron model which is the improved version of the FHN model, can produce abundant firing activities for some parameters when it is varied in a specific fixed range. The fast-slow subsystems describe the model; the fast subsystem consists of classical FHN equation 3 . The slow subsystem is one dimensional. It is biologically plausible and computationally efficient single neuron model. The commensurate fractional-order FH-R model is described as where v, w and y represent the membrane voltage, recovery variable and slow modulation of the current respectively. I measures the constant magnitude of external stimulus current, and α is the fractional exponent which ranges in the interval α < ≤ (0 1). a, b, c, d, δ and μ are the system parameters. The system reduces to the original classical order system when α = 1. μ indicates a small parameter that determines the pace of the slow system variable, y. The fast subsystem (v-w) presents a relaxation oscillator in the phase plane where δ is a small parameter. v is expressed in mV (millivolt) scale. Time t is in ms (millisecond) scale 4,6,11 . It exhibits tonic spiking or quiescent state depending on the parameter sets for a fixed value of I. The parameter a in the 2D FHN model corresponds to the parameter c of the FH-R neuron model 53,55 . If we decrease the value of a, it causes longer intervals between two burstings, however there exists a relatively fixed time of bursting duration. With the increasing of a, the interburst intervals become shorter and periodic bursting changes to tonic spiking.
The relation between injected current stimulus and membrane potential to generate an action potential, i.e., spike 10 , was previously introduced. The ideal resistor-capacitor theory describes the passive cell membrane dynamical analysis and the non-ideal resistor-capacitor circuit diagrams can characterize the oscillatory behavior 10,16,39,56,57 . The theory preserves the membrane voltage behavior. It plays a significant role in analyze the dielectric behavior of the cell membranes 16,39,57 . It was observed in experimental results that fractional-order dynamics follow a general power-law relation 10,56 . In the electrical activities of neurons, it was shown that the power-law dynamics follows α = .
0 76 and 0.86 for warm and cold frog sciatic neurons respectively 39 . The non-ideal capacitor theory for current-voltage relation was described by a fractional-order derivative as follows, It follows a power-law dynamics and preserves the memory effects in the variations of the membrane voltage [9][10][11]16,39 . We consider this contrast in the fractional-order condition. The membrane voltages and specific membrane potential changes may instigate seizure-like activity in epilepsy 55 . It may cause reactions in the muscles for the specific strength of the stimulus. This type of bursting phenomenon can be explored in a more general way so that it may span in different research areas 55 . Let us study the fractional-order fast-slow system that contributes Method numerical solution scheme. To examine the fractional dynamics of FH-R model, we consider the most familiar definition of the fractional derivative in Caputo sense [13][14][15] . Consider the fractional-order derivative of a variable x(t) for the fractional exponent α ∈ (0, 1) as folllows using the definition, we have where gamma function is defined as 1 . An additional advantage of Caputo order derivative is that the derivative of a constant is zero. It is efficient to integrate all the previous activities of the function weighted by a function that follows power-law dynamics. Now, applying the L1 scheme 9,41,42,58 on Eq. (3), approximating the fractional-order derivative as and combining Eqs (2) and (4), the numerical solution of Eq. (2) can be formulated as where, t k represents the k th time step and = ∆ t k t k .
The variable x is considered as ≡ v w y x ( , , ) in our numerical results. Approximation of the fractional-order derivative for the membrane voltage (v(t)) is given by Similarly, we can derive numerically the expressions for other two variables (w and y) of Eq. (1). Hence the numerical solution of Eq. (2) can be summarized as the difference between the markov term weighted by the gamma function and the memory trace. Memory trace has the main functional role in the fractional-order system as it integrates all the past activities. The markov term weighted by the gamma function is given by ) . The memory trace has no effect for α = 1 and the fractional-order system behaves like classical order model. The nonlinearity in the memory trace increases as we decrease the fractional-order α from 1 and the system dynamics depends on time. The fractional-order FH-R system is numerically integrated by using this scheme. We have considered different sets of parameters as follows 53 I 0 3125 and remaining parameters are similar as above. We perform the analysis of FH-R model with these parameter sets. The system shows different firing patterns like elliptic bursting, tonic spiking/regular spiking, fast-spiking and high amplitude single spike with small amplitude oscillations. The different firing activities together with the mode transitions are investigated for different parameter regimes corresponding to qualitatively various dynamical behavior of a nerve cell. the characteristics of the fractional-order biophysical model. Stability analysis. The fixed points 3 3 respectively. Depending on the nature of the discriminant of the cubic polynomial

and the memory trace is given by
, the system (1) can have maximum three equilibrium states. Throughout this study, the assumption (A) < + bd d b holds (based on the numerical data). Proposition I. The cubic function F(v*) is strictly increasing and there exists only one branch of equilibrium state = The discriminant of F′ is given by . Using the assumption (A), we obtain ′ < D F ( ) 0 that implies ′ > ⁎ F v ( ) 0 and the function F is strictly increasing (and invertible) on . Thus, it has only one real root 1 . The Jacobian of the system (1)  www.nature.com/scientificreports www.nature.com/scientificreports/ From assumption (A) and , we obtain that at least one of the roots of the characteristic polynomial Q(λ) is negative. Considering the value of the parameter = d 1 (which is constant and fixed for all the parameter sets), we have , 0). We will discuss the case when μ δ < b for analytical treatment and proceeding in the similar way, we can also derive the analytical results for the case when μ δ > b . The system changes its stability through Hopf bifurcation and it occurs when the trace of the Jacobian matrix vanishes i.e., . v H denotes the system variable where Hopf bifurcation occurs.
Proposition II. The equilibrium state E(q) of system (1) is asymptotically stable (independent of the fractional exponent, α) for any

. It can be obtained that in both the cases
and other two roots satisfy respectively. From the above discussion, we can conclude that the roots lie on the negative real axis, so the equilibrium state E(q) is asymptotically stable and independent of the fractional exponent.
Proposition III.
Numerical results. Assumption (A) holds for all the above mentioned sets of parameters, hence there exists only one branch of equilibrium state E(q) for the system (1). We can also find γ = − − . * . − . = − . = . F(0 9674) 4 5331 for parameter sets I and II respectively. From the proposition II, E(q) is asymptotically stable for any γ ≤ q F( ) 1 or γ ≥ q F( ) 2 , or equivalently, for any < .
I 0 1390 or > . I 3 1610, independent of α. When ∈ . . I (0 1390, 3 1610), the equilibrium state E(q) becomes unstable for some values of I and α, and the stability criteria is given by the proposition III for this range. Now, suppose the value of applied stimulus = .
I 0 4, the equilibrium point is stable for α < . 0 6951 (set II). The condition μ δ > b holds for the parameter set III and the stability of the equilibrium solution is given by the same condition (7). Thus, the equilibrium point is asymptotically stable for all α < . 0 95665. The equilibrium point for the set IV is .
. . ( 0 00028055, 0 0613089, 0 576231) 3 respectively. Therefore, the equilibrium point E is a saddle point. Further, the equilibrium point is stable for α < . 0 956455 at parameter set V. The chaotic behavior of the system (1) can be obtained for the different parameter sets using the necessary condition α π λ λ > | | − (2/ )tan ( Im( ) /Re( )) 1 60 . We consider the parameter set I. The system has one real equilibrium point − .
− . . E( 0 885098, 0 231373, 0 110098) and the eigenvalues at this equilibrium point is given by λ 1 , λ 2 and λ = − . . ± . i ( 0 000196427, 0 076349 0 245811 ) 3 . The equilibrium point E is a saddle point with index 2 60,61 . Now, using the above condition we obtain that the system (1) exhibits chaos for α > . www.nature.com/scientificreports www.nature.com/scientificreports/ Bifurcation analysis. The bifurcation analysis for the classical order FH-R model is performed using the MATCONT software. The applied stimulus (I) is treated as the predominant parameter and other parameters are fixed at their base values with μ = .
I 0 138716 (HB1) and = . I 3 161277 (HB2) respectively (see Fig. 1). The thick blue line in the figure indicates the stable equilibrium state while the dotted blue line indicates the unstable equilibrium state. The system has stable focus node for < .
I 0 138716 and > . I 3 161277 respectively. The system has saddle focus for . < ≤ . I 0 13872 0 6 and . < ≤ . I 2 6 316128 respectively. The FH-R model shows elliptic bursting at = . I 0 3125 (set I) 53 , however as we increase the value of I ( . < ≤ . I 0 4 0 6), it shows co-existence of regular spiking with bursting. Further, with the increase of I, the system exhibits regular spiking. The system has saddle point for . < ≤ . I 0 6 2 6 and in this region the system first shows regular spiking and then it shows the first spike latency with the increase of I. For . < ≤ . . I (0 1390, 3 1610). Figure 2(a,b) show the stable/unstable region in the α I ( , ) -plane for equilibrium state E(q) with parameter sets I and II  I 0 3125, then the critical value is α .
. I [0 7, 2 6] for the numerical analysis as > D Q I ( ( )) 0. The equilibrium state is asymptotically stable for most of the values of the fractional exponent α and unstable for very few values of α. For = .
I 0 4, the critical value for the Hopf bifucation is α . = .
Excitatory responses of the single fractional-order FH-R model. Now, we study the fractional-order FH-R neuron model to investigate various firing activities. We evaluate how the classical order dynamics changes its neuronal behavior and how bursting changes to different firing patterns with respect to stability and bifurcation analysis for different fractional exponents. We used the time step ∆ = .
t 0 1 for numerical results. We considered the initial conditions as small random perturbations around the fixed points for all the numerical simulations. Here the random perturbation is taken from a uniformly distributed random number in the interval (0, 1). An interesting feature is that the integer-order FH-R model produces elliptic bursting at the parameter set I. It generates decay and growth of small amplitude oscillations during the silent phase of bursting, and it is not damped rapidly. The firing pattern has multiple numbers of spikes in each burst having some subthreshold oscillations between two bursts (see Fig. 3(a)). The active and silent phases of bursting change as we slowly decrease the fractional exponent from classical/integer-order one α < ≤ (0 1). At α = . 0 98, it shows fast-spiking and after some time, it generates mixed-mode oscillations with high amplitude single spiking and low amplitude oscillations (see Fig. 3(b)). When it is further decreased to α = .
0 92, it displays a transition from tonic spiking into different spiking pattern with high amplitude and low amplitude oscillations (see Fig. 3(f)). The system shows mixed-mode oscillations at α = .
0 85 (see Fig. 3(g)). There is first spike latency when the system is in this transition mode, and the firing frequency decreases. Then, the system exhibits a Hopf bifurcation at α = .
0 99, it shows regular low amplitude spikes relative to classical order model, and it has first spike latency (see Fig. 3(j)). Further, it converges to quiescent state at α = .
0 95 (see Fig. 3(k)). Here, the system also converges to a stable fixed point ( = . ⁎ v 0 891229). Next, we consider the parameter set IV. The classical order model shows fast-spiking (see Fig. 3(l)). Then, it produces first mixed-mode oscillations then regular spiking at α = .
0 85 (see Fig. 3(m)). The system has first spike latency with the decrease of fractional-orders and firing frequency decreases. The first spike latency in the system increases with the decrease of fractional exponent α = .
0 80 (see Fig. 3(n)). Finally, the parameter set V has been considered. The classical order model exhibits single high amplitude spikes with small amplitude oscillations not decaying to completely silent phase or oscillation death (see Fig. 3(o)). When the fractional-order is decreased, the spike frequency is decreased and the period of small amplitude oscillations increases, i.e., it is growing with larger time duration with α = .
0 95 i.e.; the system converges to the stable fixed point ( = − . ⁎ v 0 948702) (see Fig. 3(q)). We characterize these diverse neuronal responses with stability and bifurcation analysis. The results show that the v-memory trace is zero at α = 1. The memory trace does not affect the dynamics of the system at α = 1 (see Fig. 4(a,d)). When fractional-order is decreased from classical order, new dynamical responses emerge. The voltage memory trace displays oscillations, i.e., it has major effects on membrane voltage dynamics and also membrane voltage affects the memory trace (see Fig. 4(b,e)). The fractional-order system is in steady-state at lower fractional-orders, α for all the parameter sets, i.e., the memory trace becomes too small, and it cannot significantly affect the membrane voltage dynamics to evoke a spike (see Fig. 4(c,f)). Therefore, the fractional-order excitable system changes to different oscillations as α increases above a threshold value for a fixed set of parameters.

Synchronization in coupled fractional FH-R models.
To study the dynamics of the coupled FH-R neurons, we consider two synaptically coupled FH-R neurons in the fractional domain. Generally, synaptically coupled excitable cells produce in-phase or anti-phase synchronous activity depending on coupling structure and strengths. We show the transitions to complete synchronization 37,61,63 (CS) regime of the coupled fractional-order excitable neurons. We use bidirectional coupling, i.e., electrical coupling between two FH-R neurons, which is biophysically efficient synaptic coupling mechanism.  www.nature.com/scientificreports www.nature.com/scientificreports/ where g e is the coupling strength, = i 1, 2 and = j 2, 1. The above system is bidirectionally coupled via membrane voltage. The synchronization regimes and its stability are examined by similarity functions 64 . Complete synchronization of these two coupled systems indicates the stability of zero solutions of the error system. The desired synchronization state is achieved by using suitable coupling strengths and appropriate fractional exponents at different parameter sets (see Fig. 5). Now, we introduce a statistical measure known as similarity function to estimate the synchronization error between the coupled neuronal oscillators to produce CS. The function is defined www.nature.com/scientificreports www.nature.com/scientificreports/ 1 /2 and S(γ) measures the phase lag between the coupled excitable systems. The smaller value of S(0) shows a high correlation between driver and response oscillators. The functional value S(0) confirms CS (see Fig. 6) regime as it converges to zero with different initial conditions. This verified the coupling scheme and effectiveness of the method for synchronization.

Discussion
Fractional-order dynamics has more advantages to real-world applications. It may produce complex dynamics, such as switching to different stabilities, periodic nature, and chaotic behavior. Our nonlinear fractional-order biophysical model shows such types of complex dynamics. The theoretical analysis and numerical results reveal some interesting neuronal responses that can be useful for further investigation of the fractional-order excitable . 1, 0 90, 0 69 for parameter sets I and II respectively. For α = 1, the v-memory trace has no effect on the fractional dynamics. The nonlinearity in the system increases as we decrease the fractional-order (α) and the system shows no oscillation at lower values of α. systems. The characteristics of a fast-slow FH-R model has been introduced in this article by using a commensurate fractional-order derivative. It has been examined that how fractional exponents influence the dynamics of the system and make it different from classical-order FH-R model that exhibits elliptic bursting. It shows different types of oscillations; spike frequencies based on different sets of parameters. The fractional exponent plays a significant role in generating and destroying bursting. It changes the nature of the system dynamics. It also makes us to understand the information processing in coupled systems 11,20,37,39 . Synchronization of fractional-order excitable systems, especially chaotic systems has potential applications to control secure communication 63 . We observed that spikes also produce for small fractional-orders in the fast-slow neuron model. The transition states for various firing modes, including the quiescent states are discussed for different fractional-orders with various sets of parameters. The significance of our work is that we consider biologically relevant electrical coupling, i.e., bidirectional coupling for two neurons showing different types of oscillations and establish the CS criterion in fractional-order coupled systems. It can be extended to a network of neurons with such type of fractional-order neurons for a bidirectional or gap junction coupling scheme 11,42 . It has become a challenging task to select suitable neuron model with appropriate parameter sets that exhibits different dynamical behavior when we introduce the fractional-order component in the system. This type of study of excitable biophysical systems is limited as different complex mathematical solutions arise; however, some techniques have already been developed to investigate the fractional-order dynamics. It may play a significant role in understanding the signal processing dynamics, noise-induced electrical activity, the synaptic mechanism in neuronal populations, different neurocomputational features, properties of different types of neural networks for complex brain functioning in healthy and diseased state conditions such as neurological disorders. Further study with fractional-order dynamics is needed to investigate the excitable single neuron model with its coupled nature and the dynamics of the different structured excitatory-inhibitory neuronal population. 0 99 and (c,f) α = . 0 98 at parameter sets I and III respectively. The synchronization error between the coupled neuronal oscillators converges to zero as we increase the coupling strength (g e ) in a fixed range which confirms CS between the oscillators. However, with further decrease of α, it does not show CS.