Impact of external applied currents in BVP model

The understanding of the activity of neurons in the brain has been modeled as nonlinear systems using mathematical modeling for decades. Nonlinearity in brain dynamics is complex structure to do mathematically but computational techniques make this area of research quite interesting and easy to study the dynamics. With advancement of new technology, mathematical and computational studies are more preferable to understand the behavior of neurons in a single cell to global cognitive process. In the present study, the impacts of different externally applied currents on the behavior of neurons in a simple BVP model (Bonhoeffer-Vander Pol Model) are analyzed thoroughly. The results of BVP model are similar to the characteristics of neurons shown by the Hodgkin-Huxley Model. In the BPV model, when system is stable, neurons are in resting-state. Unlike Hodgkin-Huxley model which follows all-or-none law, the BVP model does not follow this all-or-none rule. In the BVP model, there is an intermediate phase where no spike forms, but when sufficiently large input applied then spikes emerge. On applying constant current in BVP model, system is stable while it exhibits oscillatory behavior when current is applied externally above threshold value of it. If sinusoidal, continuous wavelet, and har wavelet form of external applied currents are injected then continuous firing emerges which have several interesting dynamics. Numerical simulations have been performed to understand the bifurcation analysis of the BVP model. Oneparameter and two-parameter bifurcation diagrams have been drawn in which threshold current values are discussed.


Introduction
The brain is a major organ that is responsible for thought, feelings, learning, and memory, etc. Brain controls many other functions of body parts, breathing, heartbeat, etc. The neuron is called the fundamental, structural, and functional component of the brain [1]. The figure 1 shows how neurons communicate with each other. Neurons receive input as a signal and send as output to other neurons which causes transmission of informations [1]. Neuron activity is due to a change in the membrane potential of the cell [2]. The cell has an excess of negative charges inside as compared to outside of the membrane, and it causes differences in membrane potential [3,4]. This change of membrane potential generates electrical signals. Neuron communicates with other neurons by transmitting and exchanging electrical (Fig. 2) via synapse which is also called an action potential [1,4].
An action potential can be fired by neurons only when neurons receive sufficiently large input otherwise not [1]. The Resting-state of a neuron is the condition of the membrane when the resulting ionic current is zero [5], and the potential is called resting potential. Hodgkin-Huxley [1,4] proposed the famous physiological model named Hodgkin-Huxley model which describes the generation of action potential experimentally as well as computationally in the axon of giant squid [6][7][8]. The experimental data and dynamics of neurons described in this study became the basis for several studies related to neurons behavior and expansion of different concepts and ideas [9] through mathematical modeling. The main aim of these studies is to resolve the complexity of neuron's activity and develop a simple model that mimics the properties of previous mathematical models. Hodgkin-Huxley model is a four-dimensional mathematical model which is difficult to analyze. FitzHugh extracted minimum informations from Hodgkin-Huxley equations and had formed a model named Bonhoeffer-Vander Pol Model (BVP Model). This is a simple model in two dimensions with three parameters which gives maximum information as similar to the Hodgkin-Huxley model [6]. Tasaki and Hgiwara [10] had shown that after injecting tetracthylanmonium chloride, there is a formation of the action potential [6]. FitzHugh combined the equations of the Vander Pol Model [4] and Bonhoeffer Model [5,7] to a single minimum-parameter model named as BVP (see 2) since this model exhibits equivalent geometrical properties [6].

BVP Model
The BVP model of a cell membrane produce maximum informations of neuron's activity and henceforth, it is important to study. Parameter estimation has been performed by many researchers for FHN model (FitzHugh-Nagumo model) and its modified models [11]. Using the numerical simulations technique, the sensitivity of parameters and dynamical behavior of the BVP model are analyzed and discussed here. BVP model can be formed with the help of following system of differential equations: where x is membrane potential, y is the recovery variable, I(t) is the magnitude of externally applied stimulus, and a=0.7, b=0.8, c=3.0 are constants Table 1. The dynamical qualitative 2.5) to system 1 to observe the behavior of the system and pattern of the relation between membrane potential and recovery variable. After that, we have performed a numerical simulation to understand bifurcation analysis (see 2. 3) using XPPAut. One-parameter bifurcation diagrams (see 2.3.1 and two-parameter bifurcation diagrams (see 2.3.2) have been performed to understand that how and when the qualitative structure of the system 1 changes as well as the relation between the parameters and its influence on the system respectively. Finally the results have been be concluded in discussion section (see 3).

Phase Space Analysis
In the phase-space analysis, with initial condition, the solutions of the differential equations are given by trajectories or paths [2] which may be useful to understand the dynamics of neurons. Firstly, we varied the initial values of membrane potential to positive and negative with recovery variable to zero value. We observed that for positive membrane potential neurons fire action potential but not for negative values. However, in both the situation neurons go to resting-state as time passes. This shows that there is some threshold value of membrane potential above which action potential is fired. But there is some intermediate region between an excited neuron and those which are not excited which is called quasi-threshold. In the case of the Hodgkin-Huxley model, there is all-or-none law [1]. That implies that either action potential will fire or not and there is no intermediate case. By observing phase-space some interesting outputs have been observed which is quite different from the Hodgkin-Huxley model.  Firstly the BVP model contains a quasi-threshold phenomenon and secondly stable limit cycle even though system 1 is potentially unstable [2].
We can observe from the figure (Fig. 3) that the trajectories approach the stable closed path. This graphical approach shows that there is a formation of spikes in the neural activity described by the Hodgkin-Huxley. x and y are considered as membrane potential and recovery variables respectively. The equilibrium point is shown by the intersection of x and y nullcline. We can observe that for parameters value I(t)=0, c=3, b=0.8, and a=0.7, system 1 is unstable and approaching towards a curve (Fig. 3). But as we increase the external applied current I(t)=0.5 then all trajectories approach to a stable limit cycle (Fig. 4). It is seen from the phase plane that the middle position of x-nullcline is very unstable. Therefore a sufficient amount of current is required to neurons for the firing of the action potential. For I(t)=0.5, system 1 generates a limit cycle but for I(t)=0 there is no limit (Fig. 4). The formation of the limit cycle in the system indicates the spike formation in the neural activity.

Time Series Analysis
Time series analysis is one of the important tools to understand the dynamics of any mathematical model. It is a powerful technique to observe the behavior of the system 1. Here time series analysis has been performed for different external applied currents and the relationship between variables of the system has been shown to understand the dynamics of neurons.  No activity can be observed in the behavior of the membrane potential when there is no external applied current in system 1 ( Fig. 5 & Fig. 6 ). Here membrane potential shows refractory stage where impulse is not transformed to action potential.

When constant external current is applied to neuron
When constant external current is applied to the system 1, the different behavior of membrane potential and relation between membrane potential and recovery variable has been observed for the different values of the current.
Neurons go to resting-state after some perturbation in the case of I(t)=1 µA/cm 2 (Fig. 7). This perturbation is also not sufficient to generate action potential. If we increase the value of externally applied current, then the periodic behavior of neurons has been obtained. Here the impulse is sufficient to generate action potential. If we slightly change the current to I(t)=1.1 µA/cm 2 we have observed that there is periodic firing in the behavior of membrane potential (Fig. 9). Also from the phase portrait, it is clear that the system exhibits a limit cycle (Fig. 10). When we increase the current to I(t)=1.5 µA/cm 2 we can see clearly from figure 9 that periodic behavior continues and the system again produces a limit cycle (Fig. 10). Again for I(t)=5 µA/cm 2 neuron shows a resting state after a single spike (Fig. 21). It has been observed that for I(t)> 1 µA/cm 2 and I(t)≤ 4.5 µA/cm 2 system shows periodic behavior and forI(t)> 4.5 µA/cm 2 system shows stable behavior (Fig. 7-22). Now we have applied varying current to analyze the dynamics of neurons behavior. When we have applied sinusoidal current, we have observed only the periodic behavior of the neurons. But when we have applied wavelet current then we have obtained stable as well as periodic behavior of the neurons. Parameters are very sensitive in the case of wavelet currents. Even small perturbation has changed the dynamics of the neurons.

When external current is applied in the form of sinusoidal current to neuron
The behavior of neurons has been observed for no and constant external applied current and stable and periodic behavior has been observed. Now, an external current is being varied by a sine function (Eqn 2) and its behavior has been observed for the system 1. Continuous periodic behavior has been obtained when the sinusoidal external current has been applied.

I(t)=a sin(t)
( We have observed periodic behavior for all lower as well as the higher amplitude of the sinusoidal current (Fig. 23). The Limit cycle has been obtained in phase plane analysis shown in figure 24. Now some wavelets have been considered to explore the dynamics of neurons.
Here two types of wavelet functions have been considered. First one is continuous wavelet in the form of sine and second one is Har wavelet. Time series analysis has been done for both the functions that are treated as external applied curerrent to neurons. Wavelet theory is a mathematical technique that has been emerging as an important tool in the field of not only computer science but also signal processing, computational neuroscience, and other subjects. Neural network and time series estimation needs a modeling task. Neural networks have limitations for non-linear systems. Wavelet analysis becomes very helpful while dealing with the non-linear system. Here continuous wavelet has been considered as an external applied current. Time series for different values of parameters have been examined (Fig. 25-32). Relation between membrane potential and recovery variable have been shown with help of phase diagram (Fig. 25-32). Periodic firing can be observed clearly from the given figure (Fig. 25). An interesting pattern has been observed in phase diagram (Fig. 26). Voltage time series for all given values are almost similar but different pattern in phase diagram. In figure 26, trajectories are evolving from down to up for a=0.7, b=0.8, c=3 (Fig.14). But as the values are perturbed to a=0.7, b=1.5, c=3, and then a=0.7, b=2, c=3 trajectories start evolving from up to down (Fig. 28 & 30). For a=0.7, b=2, c=3 same pattern like figure 26 has been found (Fig. 31).

When external current is applied to neuron in the form of Har wavelet
The Haar wavelet is a simple wavelet and its mother wavelet (Eqn 4) is given by ϕ(t) on interval [0,1], as follows: Let us consider the externally applied current as a single pulse of the function t given by I(t) (Eqn 5) as follows (Fig. 33): This pulse is jumped function considered in the interval . Now we have analyzed the time series by using numerical computation to understand the behavior of system 1. When we have applied wavelet current as an externally applied current periodic behavior has been observed. Neurons show periodic firing for a=0.1, b=0.8, c=3 (Fig. 34) and also for a=0.1, b=1, c=2 there is late firing of action potential (Fig. 35). The Limit cycle has been observed in phase plane analysis which clearly shows the periodic behavior of neurons. We have found that when stimulus is applied in the form of continuous wavelet function, neurons behaviour is observed periodic for all values of given parameters (Fig. 36 & 37).

One-parameter Bifurcation Analysis
BVP model shows multistability since stable equilibrium and stable periodic orbit both exist. System (1) contains only odd power so it has symmetric behaviour. We have analyzed that membrane potential x is increasing slowly up to HB1 (Fig. 38) and then rapidly up to HB2 concerning external applied current I. System (1) shows subcritical Hopf bifurcations at both point HB1 (Fig. 39) (I(t)=0.3465) and HB2 (I(t)=1.404). Stable behavior is obtained for I(t) ≤ 0.3465 and I(t) > 1.404, also unstable for 0.3465 < I(t)< 1.404 (Fig. 38(left)). As we increase the external applied current I(t), system 1 shows unstable behavior switched from stable behavior. If we further increase the external applied current then again system shows stable behavior. Now we have analyzed the influence of externally applied current I(t) on recovery variable y. It has been observed that the recovery variable increases slowly below HB1 and increases rapidly from HB1 to HB2 (Fig. 40). The system is stable for I(t)<HB1 and I(t) >HB2, but it is unstable for HB1<I(t) <HB2. As we increase the external applied current the system shows unstable behavior and further increment again leads to stable behavior. When we have analyzed the relation between membrane potential with parameter a, we found that membrane potential decreases slowly up to HB1 and then decreases   rapidly up to HB2 with respect to parameter a. The system goes under sub-critical Hopf bifurcation HB1 (a=-0.4228) and HB2 (a=0.4228) which can be shown in figure 40(right). It has been observed that the system shows stable behavior for a < HB1 and a > HB2, but unstable behavior for HB1 < a < HB2. Figure 42 gives the relation between membrane potential x and parameter b. It has been observed that membrane potential decreases slowly with respect to parameter b. the system shows unstable behavior before HB1 (b=0.4153) and stable for b > HB1. One-parameter bifurcation analysis has been done and results have been shown in the figures (Fig. 38-42). Now two-parameter analysis will be discussed for all the three parameters of the system 1.

Two-parameter Bifurcation Analysis
We have observed the influence of other parameters on the system by two-parameter bifurcation analysis. We have observed the influence of other parameters on the system by two-parameter bifurcation analysis. Each parameter has been observed concerning external applied current to understand the dynamics of the system 1. Here, a stable region is shown by a grey shaded region, and the periodic region is shown by the black shaded region. For a Figure 43. Two-parameter bifurcation of parameter a concerning the external applied current with parameters a=0.8, b=0.7 & c=3. Stable state is shown by grey shaded region) and unstable by black shaded region  vs.I(t), a stable region is obtained for higher values of parameter a and lower values of external applied current as well as for lower values of parameter a and higher values of external applied current (Fig. 43). The periodic region is obtained for a short-range of parameters a and an external applied current. For b vs. I(t), a periodic region is obtained for lower values of parameter b, and a stable region is obtained for higher values of parameter b, preferably b> 1 (Fig. 44).
For c vs. I(t), the stable region is obtained for the lower values of external applied current and higher values of parameter c as well as higher values of external applied current whereas periodic region has been observed for a specific range of external applied current and parameter c (Fig. 45).

Discussion and Results
Hodgkin-Huxley's proposed model was historic research in 1952 in the field of neuroscience to understand the transmission of the nerve impulse. Although the model is complex to handle mathematically it provides many characteristics about neurons computationally. After the Hodgkin-Huxley classic model, various models have been established later based on it. BVP model was a simple model as it contains only two variables and three parameters, but mimics important properties given by the Hodgkin-Huxley model. In this manuscript, we have analyzed the impact of different externally applied currents on system 1. Results observed after analyzing all the perspective of time series, bifurcation, and numerical simulation are given as follows: i. Neurons are resting when the system is stable or it is not stimulated.
ii. There is an intermediate stage where no spike forms. But sufficiently large input generates spikes.
iii. When we applied constant current then the system shows stable and periodic behavior for specific current values.
iv. When the sinusoidal form of external applied current is injected then we have observed the continuous firing of the action potential. The threshold value of the current is required to pass the information from one neuron to another neuron.
v. When continuous and Har wavelet are introduced in the form of external applied current, then we have found periodic signals with a different pattern. Also, an interesting pattern in the phase diagram has been observed in the form of emerging trajectories in upward and downward directions for different values of parameters.
vi. With the help of bifurcation analysis threshold value has been investigated and bifurcation points have been examined. One-parameter and two-parameter bifurcation analysis helped to understand the influence of parameters on the given system. This study will make an understanding of how the external current plays role in the transmission of neuron signals. The brain is a very complex part of the body and there are still so many things to explore regarding brain dynamics.

Acknowledgement
First author would like to acknowledge DST/INSPIRE/03/2016/000597 for financial support.