Subcritical Hopf Bifurcation and Stochastic Resonance of Electrical Activities in Neuron under Electromagnetic Induction

Introduction

Memristor or memory resistor was supposed by Chua (1971) 46 years ago as the forth fundamental circuit device along with resistor, inductor and capacitor, and it was successfully realized by Stan William's group at HP Labs in 2008 for the first time (Strukov et al., 2008). Because of the huge storage potential and the complex nonlinearity, the memristor has recently attracted considerable attention in theoretical and applied neuroscience (Itoh and Chua, 2008; Yogesh and Stephen, 2009; Wen et al., 2013; Bao et al., 2015; Chen et al., 2017; Zha et al., 2017).

In modern society, human being or animals are inevitably more or less exposed in the electrical hazards of ubiquitous electromagnetic radiation, and the fact that this radiation can have severe consequence on biological rhythm and recognition has attracted much attention (World Health Organization, 2011), but how the electromagnetic radiation changes the behavior of neural systems is still unclear. Fortunately, the device of memristor emerges and it can act as a feasible tool for exploring the influence of electromagnetic radiation on neural system activities, since one can keep the consistency of physical dimension (or unit) when modeling the membrane potential and magnetic flux into coupling systems (Wu et al., 2016). Several investigations have been done in this regard. For example, Lv et al. (Lv and Ma, 2016) proposed a comprehensive modified Hindmarsh-Rose neuron model by introducing the magnetic flux as a fourth variable, and Lu et al. (2017) imposed different types of electrical stimulus impended with a high-low frequency current on this improved HR model to investigate mode selection in neural activity. Guo et al. (Ren et al., 2017) used memristor to discuss the polarization and magnetization in excitable neural model. Ma et al. (Wu et al., 2016; Ma et al., 2017) adopted a magnetic flux across the membrane to describe the electromagnetic induction and the spiral waves have been induced. Although the abundant firing patterns have been revealed, the underlying dynamical mechanism responsible for these patterns in these newly-built models has not been disclosed.

In nervous systems, noise not only has various origins but seldom acts as a trivial disturbance (Tanabe and Pakdaman, 2001; Hasegawa, 2004; Faisal et al., 2008; Shao and Kang, 2014; Sun and Shi, 2014). One of the anti-intuitive phenomena of noise is often termed as stochastic resonance (SR), where a suitable noise can amplify the external weak coherent signal under certain nonlinearity. In the absence of electromagnetic interference, many theoretical or experimental literatures have shown that living organisms can utilize noise as a benefit in detecting or transferring weak signal on both cellular and system levels (Mark et al., 2009). Kang et al. (2005) showed the existence of signal-to-noise ratio gain of SR based on the leaky integrate-and-fire neuron model, Jiao and Wang (2010) observed SR under the effect of synaptic transmission noise; Sun and Li (2016) demonstrated that the partial time delay can induce a stochastic multi-resonance in a Watts-Strogatz neuronal network. Nevertheless, to our knowledge the phenomenon of SR has not been explored in neural system under electromagnetic disturbance. Therefore, we naturally wonder whether a weak coherent oscillation in the electromagnetic disturbance can be amplified through SR.

After a modified FitzHugh–Nagumo (FHN) model with flux-controlled memristor is introduced in Section A modified FitzHugh–Nagumo neuron model, some analytical and numerical results on the bifurcation behavior of the system are derived in Section Analysis of bifurcation. And then, SR in the modified model with a weak periodic modulation is exhibited and explained in Section Stochastic resonance. Finally, conclusion and discussion are given in Section Conclusion and Discussion.

A Modified FitzHugh–Nagumo Neuron Model

Let us start with the conventional FHN neuron model (Fitzhugh, 1961; Nagumo et al., 1962).

$\begin{array}{cc}\{\begin{array}{l}\dot{v}=v(v-a)(1-v)-w+I_{ext}\hfill \\ \dot{w}=\epsilon (v-dw)\hfill \end{array}& \left(1\right)\end{array}$

where the fast-varying trans-membrane potential v and the slow current variable w are treated as dimensionless. In the model (1), the nonlinear term v(v − a)(1 − v) stands for total trans-membrane ionic currents per unit area, I_ext is the external forcing current. Since our purpose is to investigate the influence of electromagnetic induction, the external forcing current is set to zero. The parameters ε = 0.02, d = 1 and a = 0.5 are fixed such that the dynamical evolution of the model (1) starting from any initial state can asymptotically approach a resting equilibrium state when I_ext = 0. We keep this zero external input throughout the context.

In order to consider the effect of electromagnetic induction on membrane potentials of neuron, with the help of Strukov et al. (2008) we employ the memristor to realize the coupling and modulation on membrane potential from magnetic flux for maintaining the consistency of physical meaning (Wu et al., 2016). Note that the memristor characterizes the relation between charge and magnetic flux, so if let q be charge and φ the magnetic flux, then a voltage across a charge-controlled memristor can be modeled as v(t) = M(q(t))i(t) with $M(q)=\frac{d\phi (q)}{dq}$. In reverse, a flux-controlled memristor should be described by i(t) = W(φ(t))v(t) with $W(\phi )=\frac{dq(\phi )}{d\phi }$. Here the explicit forms of M(q) and W(φ) should depend on the design of the device of memristor.

We take a flux-controlled memristor of $W(\phi )=\frac{dq(\phi )}{d\phi }=k(\alpha +3\beta {\phi }^2)$ (Itoh and Chua, 2008; Bao et al., 2010a; Wu et al., 2016) to modify the conventional FHN model. We choose α = 0.1 and β = 0.02 to generate complex dynamical behavior (Bao et al., 2010a,b; Wu et al., 2016). Thus, the time evolution equation of the membrane potential v becomes

$\dot{v}=v(v-a)(1-v)-w+k(\alpha +3\beta {\phi }^2)v.$

According to the Faraday's law, the change of the magnetic flux φ is dominated by varying voltage and magnetic flux, we can suppose that the time derivative of φ is a linear function of v and φ, that is to say,

$\dot{\phi }=k_{1}v-k_2\phi +{\phi }_{ext}$

where φ_ext, as a bifurcation parameter, is used to describe the bias in external forcing magnetic field.

Therefore, the improved FHN model, which takes the effect of electromagnetic induction into consideration, has the following form of a set of three-variable nonlinear ordinary differential equations.

$\begin{array}{cc}\{\begin{array}{l}\dot{v}=v(v-a)(1-v)-w+k(\alpha +3\beta {\phi }^2)v\hfill \\ \dot{w}=\epsilon (v-dw)\hfill \\ \dot{\phi }=k_{1}v-k_2\phi +{\phi }_{ext}.\hfill \end{array}& \left(2\right)\end{array}$

Although the original FNH model has an exclusively stable asymptotic state for the given parameters of ε, d, and a, the introduction of electromagnetic induction can induce complex periodic or bursting firing patterns in the modified model. We will explore the involving dynamical mechanism in the absence of noise and in the presence of noise in Sections Analysis of Bifurcation and Stochastic Resonance, respectively.

Analysis of Bifurcation

The concept of bifurcation in nonlinear dynamical theory can be categorized into static bifurcation and dynamic bifurcation. Usually, the former means the change in number or the stability of equilibrium points, while the latter refers to the similar changes relating to limit cycles (Zhang, 2005; Xie et al., 2008a,b). If an equilibrium state is stable, the affiliating system will evolve closer and closer to it if the initial state falls within its basin of attraction; otherwise, the system will leave it forever if the initial state is not exact on it. In neuron model, the stable equilibrium usually stands for a resting state, and the limit cycle corresponds to the repetitive firing state. Bifurcation has been scrutinized analytically and numerically by many researchers in computational neuronal science. Quantities of literatures have dedicated to the bifurcation behavior of neuron models during the past decades (Hassard, 1978; Guckenheimer and Labourian, 1993; Eugene, 2000; Xie et al., 2008b; Jia and Gu, 2017), but the involving investigations have been extended to the improved models with electromagnetic induction, even if the abundant firing patterns have been revealed (Lv and Ma, 2016; Wu et al., 2016; Lu et al., 2017). In this section, we aim to disclose the underlying dynamical mechanism responsible for the emergence of the firing pattern in the improved system (2).

At first, let us pick out the constant solution, namely the equilibrium points of the system (2), where the variables v, w and φ stay there forever if there is no external perturbation. Suppose E₀(v₀, w₀, φ₀) is one of the equilibrium points, then there holds

$\{\begin{array}{l}{v}_0(v_0-a)(1-v_0)-w_0+k(\alpha +3\beta {\phi }_0{}^2)v_0=0\hfill \\ v_0-d{w}_0=0\hfill \\ k_1{v}_0-k_2{\phi }_0+{\phi }_{ext}=0\hfill \end{array}$

and thus

$v_0\left[(v_0-a)(1-v_0)-\frac{1}{d}+k\left(\alpha +3\beta {(\frac{k_1{v}_0+{\phi }_{ext}}{k_2})}^2\right)\right]=0$

with three solutions given by

$v_{01}=0,v_{02}=\frac{-B+\sqrt{B^2-4AC}}{2A},v_{03}=\frac{-B-\sqrt{B^2-4AC}}{2A}$

where $A=\frac{3k{k}_1^2\beta }{k_2^2}-1B=\frac{6k{k}_1\beta {\phi }_{ext}}{k_2^2}+1+aC=\frac{3k\beta {\phi }_{ext}^2}{k_2^2}-a-\frac{1}{d}+k\alpha$. Let E₀₁, E₀₂ and E₀₃ to be the resultant equilibrium points corresponding to v₀₁, v₀₂ and v₀₃, then $E_{01}={(0,0,{\phi }_{ext}/k_2)}^T$.

Next, let us explore the stability of the three equilibrium points with φ_ext as bifurcation parameter. Note that the stability of an equilibrium point is determined by the eigenvalues of its Jacob matrix. That is to say, if all the eigenvalues are of negative real part, then the equilibrium point is stable, otherwise it might be marginally stable or unstable. According to the distinction criteria (Zhang, 2005), branches of the equilibrium points and their stability on the φ_ext − v plane are shown in Figure 1, where the solid curves represent the stable branches, the dash lines indicate the unstable ones, and evidently the points A~H indicate the occurrence of bifurcation.

FIGURE 1

Figure 1. Schema of the branches of equilibrium points for trans-membrane potential and their stability vs. the bifurcation parameter φ_ext. The parameters are set as k = 1, k₁ = 0.5, and k₂ = 0.9. The branches locate at lower left and right stand for E₀₂ and the branches situate in upper left and right are E₀₃. The flat line is E₀₁. In figures, the solid curves stand for stable branches, and the dash curves denote unstable branches. Remark: The figure is plotted by the software of xppaut.

In fact, we can further distinguish the bifurcation types of the bifurcation points A~H. In general, the appearance of a pair of pure imaginary eigenvalues signifies Hopf bifurcation, and the emergence of a zero eigenvalue predicates fork-type or saddle-node bifurcation (Zhang, 2005). In this paper, our emphasis is put on the identification of Hopf bifurcation. For the system (2), if the resultant 3 × 3 Jacobi matrix has one eigenvalue of negative real part and two zero real parts at some critical value of the bifurcation parameter, then as usual we say that Hopf bifurcation occurs, through which a constant membrane potential solution becomes unstable, but a stable periodically oscillatory action potential solution appears.

Let us take the equilibrium point E₀₁ as example to demonstrate how Hopf bifurcation occurs. From the linearization Jacobian matrix of the system (2)

$J(v,w,\phi )=\left[\begin{array}{ccc}-3v^2+2(1+a)v-a+k(\alpha +3\beta {\phi }^2)& -1& 6k\beta \phi v\\ \epsilon & -\epsilon d& 0\\ k_1& 0& -k_2\end{array}\right],$

One can easily obtain the characteristic determinant of system (2) at E₀₁ as

$\begin{array}{l}|\lambda I_3-J(0,0,\frac{{\phi }_{ext}}{k_2})|=\left|\begin{array}{ccc}\lambda -A_1& 1& 0\\ -\epsilon & \lambda +\epsilon d& 0\\ -k_1& 0& \lambda +k_2\end{array}\right|\\ =(\lambda +k_2)[(\lambda -A_1)(\lambda +\epsilon d)+\epsilon ]=0\end{array}$

with $A_1=-a+k(\alpha +3\beta {(\frac{{\phi }_{ext}}{k_2})}^2)$. Obviously, there is an eigenvalue λ = −k₂ satisfying Re(k₂) ≠ 0, and thus if a pair of pure imaginary eigenvalues λ = ±iω exists, there must hold true

$-{\omega }^2\pm (\epsilon d-A_1)\omega i-A_1\epsilon d+\epsilon =0$

i.e.,

$\{\begin{array}{l}-{\omega }^2-A_1\epsilon d+\epsilon =0\hfill \\ \epsilon d-A_1=0,\hfill \end{array}$

which leads to A₁ = εd and $\omega =\sqrt{\epsilon -{\epsilon }^2{d}^2}$. Thus, the branch of the equilibrium point E₀₁ undergoes two Hopf bifurcations at parameter

${\phi }_{ext}=\pm \sqrt{\frac{\epsilon d+a-k\alpha }{3k\beta }}{k}_2$

which can be reduced to ${\phi }_{ext}^{(1)}=-2.381$ and ${\phi }_{ext}^{(2)}=2.381$ for given parameters as denoted in the caption of Figure 1. Clearly, they correspond to the bifurcation points A and B on Figure 1.

Similar to the above analysis, we also can know that the branch of the equilibrium point E₀₂ undergoes Hopf bifurcations at C and D, and the branch of the equilibrium point E₀₃ undergoes Hopf bifurcations at points E and F, with the bifurcation parameter φ_ext being −5.386(C), 5.512(D), −4.113(E) and 3.236(F), respectively. Further, we can confirm but skip the details that G and H, as the intersections of the branches of E₀₁ and E₀₂, are saddle-knot bifurcation points with bifurcation parameters φ_ext = ±4.347, which essentially. Here we stress that all of the numerical errors do not exceed 0.001.

The occurrence of Hopf bifurcation has been verified by the existence of a pair of pure imaginary eigenvalues, but Hopf bifurcation can further distinguished into two types: subcritical Hopf bifurcation and supercritical Hopf bifurcation by checking dependence of the membrane potential difference over a sufficiently large time span via bifurcation parameter. In general, in subcritical Hopf bifurcation the oscillatory solution occurs before the bifurcation point, but in supercritical Hopf bifurcation the oscillatory solution emerges after the bifurcation point, and whether a hysteresis loop exists is a typical distinction criteria. In fact, if for some given parameter (v_max − v_min) ≈ 0, we can regard that the system approaches an equilibrium state; otherwise, we can pick out an oscillatory solution, which can correspond to the limit cycle in the phase space if the occurrence of the oscillation is due to a bifurcation induced by an equilibrium point. For the system (2), the difference (v_max − v_min) vs. the varying φ_ext is depicted in Figure 2. From the picture, a hysteresis loop can be clearly observed near ${\phi }_{ext}^{(1)}$, and another hysteresis loop can be captured near ${\phi }_{ext}^{(2)}$ by some partial amplification., therefore the Hopf bifurcations at A and B are both subcritical.

FIGURE 2

Figure 2. Dependence of the trans-membrane potential difference on the bifurcation parameter φ_ext (A). The green solid branch shows how the trans-membrane potential difference changes as the parameter φ_ext increases, while the blue dash branch gives the change trend as the parameter φ_ext decreases. We can see that within certain parameter range there are two hysteresis loops, the left-hand one is evident and the magnification of the right-hand one (B) has been inserted in the figure. The existence of the two hysteresis loops show the bifurcation at ±2.381 both are subcritical Hopf bifurcation.

For an intuitive understanding, we depict the time series of the trans-membrane potential and the phase diagram when φ_ext > 0 in Figure 3. As Figure 3 shows, the membrane potential will stay at the resting level when the bifurcation parameter is less than the critical bifurcation (Figures 3A,B), but it will evolve according to a periodic motion as the bifurcation parameter increases (Figures 3C–F). More precisely, the Figures 3C–F respectively correspond to subthreshold oscillation and superthreshold oscillation (impulsive discharge) of neurons. Moreover, Figures 3G,H reveal that the model (2) will attain another equilibrium state of a high asymptotical membrane potential, which should be morbid for neuronal activity.

FIGURE 3

Figure 3. Time evolution of the trans-membrane potential and the phase diagram of the model (2). From up to down, the bifurcation parameter φ_ext is equal to 2.25 (A,B), 2.3805 (C,D), 3.0 (E,F), and 3.4 (G,H), respectively.

Stochastic Resonance

It has been extensively proven that weak periodic signal can be amplified by random fluctuation in many nonlinear systems through the principle of SR (Heneghan et al., 1996). This phenomenon has also been well-documented in neural systems in the absence of electromagnetic induction (André, 1993; Kang et al., 2005; Gosak et al., 2007). Here what we interest in is to check whether a weak subthreshold oscillation in the electromagnetic field can be amplified by a suitable amount of fluctuation.

Assume that the external forcing magnetic field consists of a subthreshold signal and noise, i.e., the form of φ_ext is described as follows

${\phi }_{ext}=r\sin (2\pi ft)+{\phi }_{ext}^0+\xi (t)$

where $rsin(2\pi ft)+{\phi }_{ext}^0$ is a subthreshold signal of bias ${\phi }_{ext}^0$, that is to say, the particle will not jump from one state to the other if it is only driven by this signal, and ξ(t) is Gaussian white noise satisfying < ξ(t)) > = 0 and < ξ(t)ξ(s) > = 2Dδ(t − s). With these factors taken into account, the system (2) can be rewritten into

$\begin{array}{cc}\{\begin{array}{l}\dot{v}=v(v-a)(1-v)-w+k(\alpha +3\beta {\phi }^2)v\hfill \\ \dot{w}=\epsilon (v-dw)\hfill \\ \dot{\phi }=k_{1}v-k_2\phi +r\sin (2\pi ft)+{\phi }_{ext}^0+\xi (t)\hfill \end{array}& \left(3\right)\end{array}$

In this section, the same parameters as in Figure 1 are used, unless otherwise stated.

Many physical measure indexes are suitable for quantifying the phenomenon of SR, such as spectral amplification factor, mutual information, resident time distribution and signal-to-noise ratio (SNR), and all these indexes can reflect the beneficial role of noise from different viewpoints. Here, we choose the SNR defined by

$SNR=10{\log }_{10}\frac{S(\omega )}{N(\omega )}$

to characterize SR. Here S(ω) represents the height of the peak at the signal frequency and N(ω) is the power spectrum density of background noise. Figure 4 shows the dependence of SNR via noise intensity D under different signal amplitude r and biased intensity of electromagnetic field ${\phi }_{ext}^0$, and these non-monotonic dependence curves exhibits the occurrence of SR in the model (3). This observation confirms that SR can occur in the presence of electromagnetic induction.

FIGURE 4

Figure 4. Dependence of signal to noise ratio on noisy intensity under different amplitude (A) or biased intensity of electromagnetic field (B) at ε = 0.005 and 2πf = 0.001. We fix b = 2.2328 in (A) and r = 0.28 in (B). It is shown from (A) that the width of the resonant peak evidently enhances as the r increases, and from (B) that the resonant peak is dramatically shifted to larger noise intensity as b away from the critical bifurcation point (2.3383).

Here let us explain why the phenomenon of SR occurs in the system (3) from the perspective of energy. As shown in Figure 5, the circle represents the power of the system offered by the subthreshold signal. If the power at some time is greater (less) than the fixed power as shown in the figure, it means the periodic signal is supplying (consuming) energy to the system. Obviously, when the noise is absent, the power from the subthreshold signal cannot drive the membrane potential to cross over the firing threshold, but with the help of the energy of a suitable noise, the membrane potential can cross the threshold over a time interval corresponding to an arc with a central angle θ. This to a certain extent is the dynamical mechanism underlying the occurrence of SR.

FIGURE 5

Figure 5. Qualitative analysis for stochastic resonance from the viewpoint of power. Every point on the circle represents the power offered only by the subthreshold signal ($rsin(2\pi ft)+{\phi }_{ext}^0$) to the system at different time. The height of center of the circle represents the fixed power offered by ${\phi }_{ext}^0$. The power provided by the noise raises the height of the circle. From the schemata, there should be suitable noise level D = D₀ such that the membrane potential transits between the state 1(resting state) to state 2(burst state) almost periodically, with a period corresponding to an arc with central angle θ(D).

Now we turn to explain why the width of resonant peak becomes large as the signal amplitude enhances and why the resonant peak shifts to smaller noise intensity as the bifurcation parameter closes to the critical value, as observed in Figure 4. From Figure 5, it is clear that as the signal amplitude increases the radius of the circle becomes larger. This means that the power supplied by the periodic signal has more dominated capability in the interplay of noise and signal. As a result, the width of resonant peak for stronger signal will wider than that for the weaker signal. Similarity, when the bias intensity increases, the “Fixed power” in Figure 5 will be boosted and the required energy for the system to reach the threshold will become smaller. Hence, the resonant peak will appear at smaller noise intensity.

In order to explain the meaning of θ(D) in Figure 5, let us resort to the statistics of the time history of membrane potential (shown in Figure 6). It is easy to see that the firing pattern in Figure 6 belongs to bursting. Discarding the transient evolution, we can calculate the interburst intervals (IBIs), burst interval (BI) and resting interval (RI). Here the IBI is referred to as the time interval between adjacent bursts (Gritsun et al., 2011), and we define the BI as burst interval, namely the duration interval of a single bursting and the RI as the time interval from the end of the first burst to the beginning of the next one. For example, as shown in Figure 6B, the IBI, BI and RI correspond to t_C − t_A, t_B − t_A and t_C − t_B, respectively. Denoting IBI(D) = 〈IBI〉, BI(D) = 〈BI〉 and RI(D) = 〈RI〉, then θ(D) in Figure 5, proportional to time, can be calculated by

$\theta (D)=2\pi \times \frac{BI(D)}{IBI(D)}=2\pi \times \frac{BI(D)}{BI(D)+RI(D)}.$

By means of the evolution curves of these average intervals via the noise intensity D in Figure 7, θ(D) is less than π when the noise intensity is under the critical value at the intersection, and it is larger than π when the noise intensity is between this critical value and 0.7.

FIGURE 6

Figure 6. Time history of trans-membrane potential under different noise intensity D. The parameters ε = 0.005, r = 0.28, f = 0.001/2π, and ${\phi }_{ext}^0=0.328$, and from (A–D), the noise intensity D is 0, 0.25, 0.5, 0.75, respectively.

FIGURE 7

Figure 7. The evolution curves of IBI(D) (red solid), BI(D) (blue circle), and RI(D) (green dot) via noise intensity at ε = 0.005, r = 0.28, f = 0.001/2π, and ${\phi }_{ext}^0=0.328$. Note that the IBI(D) cannot be defined when there is no noise, since the long time of the deterministic system (3) is static, but once the noise is induced and if the noise intensity falls within the range (0.07, 0.72) in the figure, the IBI(D) will fluctuate around the period of the sinusoidal drive (black dash), and this suggests that noise should be utilized by the model (3) in weak signal detecting. Here we emphasize that none of IBI(D)), RI(D), and BI(D) can be measured in absence of noise or in the case of too much noise, since the absence of noise leads the membrane potential to a static level, while too much noise causes the periodic bursting like patterns mix into one single noisy bursting.

Noting that in the phenomenon of SR noise helps detect weak signal, so let us demonstrate how the weak subthreshold signal is detected in the system (3). It is clear that when the noise is absent, the membrane potential eventually stays at the resting level (Figure 6A), but as the noise intensity increases until attains a suitable range, the membrane potential will evolve according to an approximate periodic motion with fluctuations (Figures 6B–D), where the relevant approximate period is basically the signal period. Therefore, if we analyze the time history within this suitable noise range, the period of the weak signal will be identified. In fact, as shown in Figure 7, the mean IBI is nearly equal to the signal drive, although this quantity cannot be defined in the deterministic case.

Viewing from the history of SR, the phenomenon was named after the term of “resonance” to some extent. In the symmetrical overdamped bistable system, the occurrence of SR is relating to a match relation between a half of the signal period and the mean first passage time, while in the underdamped oscillator system, there exists a coincidence between the driving frequency and the noise-tuned inherent frequency peak (Kang et al., 2003). These match relations suggest that there exists certain frequency interval such that the phenomenon of SR in the system (3) cannot occur if the signal frequency falls outside of this range. Figure 8 just verifies this point. For instance, when 2πf = 0.005, as the figure shows that the SNR as function of noise intensity monotonically decreases, which indicates no SR for the given signal amplitude.

FIGURE 8

Figure 8. Dependence of signal to noise ratio on noisy intensity under different signal frequency at ε = 0.005, r = 0.28, and ${\phi }_{ext}^0=0.328$. From the figure, it is clear that the SNR as a function of noise is monotonic, so SR can occur if 2πf ≤ 0.001, but if 2πf = 0.008, the function becomes monotonically decreasing, and as a result, there is no SR in this case. This figure tells that the phenomenon of SR has obvious frequency dependence, that is to say, the phenomenon can only occur within certain frequency range.

Conclusion and Discussion

By taking the electromagnetic induction into account, an improved FHN model is proposed by means of the flux-controlled memristor. With the technique of linear stability analysis, the bifurcation diagram of the deterministic autonomous model is obtained and especially the subcritical Hopf bifurcation points are identified. For the noisy weak signal modulated system, we have observed the phenomenon of SR near the subcritical Hopf bifurcation point. From the viewpoint of energy we give an explanation of the occurring mechanism. By defining the several mean intervals relating to burst, we also discuss how to detect a weak signal based on the principle of SR in this model. Our investigations once again demonstrate that electromagnetic radiation can induce electrical activity in neurons. Moreover, our investigations suggest that the phenomenon of SR could be utilized by neurons in detecting weak signal in the presence of electromagnetic induction.

We would like to have some discussion on the phenomenon of SR disclosed in this investigation from weak signal detection. It is well known that neurons communicate with each other through action potential, and the timing of action potential trains often contains more significant information than their shape, thus it should be more inspiring to disclose the phenomenon of SR both from the noise dependence of the SNR and the statistical quantities of trains of interspike intervals (ISIs). Indeed, as pointed out in the introduction section, the experimental and theoretical investigations have suggested that noise helps in detecting or transmitting weak signal (Tanabe and Pakdaman, 2001; Hasegawa, 2004; Kang et al., 2005; Faisal et al., 2008; Mark et al., 2009; Jiao and Wang, 2010; Shao and Kang, 2014; Sun and Shi, 2014; Sun and Li, 2016), however, most of the existing investigations only disclosed SR by showing an optimal noise level which could make the system have a better SNR than other noise level, but failed to try to find the relation between the maximal SNR and the mean of ISIs. Different from these existing investigations, in this paper we find the relation of the maximal SNR and the evolution of the mean of IBIs (〈IBI〉), namely, the optimal noise level for the SNR is an interval and over the same interval the 〈IBI〉 is almost equal to the signal period. Although we are not sure whether this relation is universal in general neural systems, we infer it might be a predictable conclusion for SR occurring near subcritical Hopf bifurcation point. If this inference can also be confirmed in other neuron models, undoubtedly it will provide a biologically feasible scheme for weak signal detection.

Author Contributions

Y-XF designed the model, did most of the numerical calculations and finished the initial draft. Y-MK guided in elementary methods about linear stability analysis and stochastic resonance, and she rewrote the whole presentation. YX drew Figure 1 with the software of xppaut.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The work was financially supported by National Natural Science Foundation with Grant Nos. 11772241 and 11672219.

References

André, L. (1993). Stochastic resonance in neuron models. J. Stat. Phys. 70, 309–327. doi: 10.1007/BF01053970

CrossRef Full Text | Google Scholar

Bao, B., Liu, Z., and Xu, J. (2010a). Steady periodic memristor oscillator with transient chaotic behaviours. IEEE Electron. Lett. 46, 228–230. doi: 10.1049/el.2010.3114

CrossRef Full Text | Google Scholar

Bao, B., Xu, J., and Liu, Z. (2010b). Initial state dependent dynamical behaviors in a memristor based chaotic circuit. Chin. Phys. Lett. 27:070504. doi: 10.1088/0256-307X/27/7/070504

CrossRef Full Text | Google Scholar

Bao, H., Park, J. H., and Cao, J. (2015). Adaptive synchronization of fractional-order memristor-based neural networks with time delay. Nonlinear Dyn. 82, 1343–1354. doi: 10.1007/s11071-015-2242-7

CrossRef Full Text | Google Scholar

Chen, C., Li, L., Peng, H., and Yang, Y. (2017). Fixed time projective synchronization of memristor-based neural networks with discrete delay. Sci China Inform. Sci. 60:032201.

Chua, L. O. (1971). Memristor—the missing circuit element. IEEE Trans. Circuit Theory. 18, 507–519. doi: 10.1109/TCT.1971.1083337

CrossRef Full Text | Google Scholar

Eugene, M. I. (2000). Neural excitability, spiking and bursting. Int. J. Bifurcat. Chaos 10, 1171–1266. doi: 10.1142/S0218127400000840

CrossRef Full Text | Google Scholar

Faisal, A. A., Selen, L. P., and Wolpert, D. M. (2008). Noise in the nervous system. Nat. Rev. Neurosci. 9, 292–303. doi: 10.1038/nrn2258

PubMed Abstract | CrossRef Full Text | Google Scholar

Fitzhugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466. doi: 10.1016/S0006-3495(61)86902-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Gosak, M., Marhl, M., and Perc, M. (2007). Spatial coherence resonance in excitable biochemical media induced by internal noise. Biophys. Chem. 128, 210–214. doi: 10.1016/j.bpc.2007.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Gritsun, T., Le feber, J., Stegenga, J., and Rutten, W. L. (2011). Experimental analysis and computational modeling of interburst intervals in spontaneous activity of cortical neuronal culture. BiolCybern 105, 197–210. doi: 10.1007/s00422-011-0457-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Guckenheimer, J., and Labourian, I. S. (1993). Bifurcation of the Hodgkin and Huxley equations: a new twist. Bull. Math. Biol. 55, 937–952. doi: 10.1007/BF02460693

CrossRef Full Text | Google Scholar

Hasegawa, H. (2004). Dynamic mean-field approximation to small-world networks of spiking neurons: form local to global and/or from regular to random couplings. Phys. Rev. E 70, 066107–066111. doi: 10.1103/PhysRevE.70.066107

PubMed Abstract | CrossRef Full Text | Google Scholar

Hassard, B. (1978). Bifurcation of periodic solutions of the Hodgkin-Huxley model for the squid giant axon. J. Theoret. Biol. 71, 401–420. doi: 10.1016/0022-5193(78)90168-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Heneghan, C., Chow, C. C., Collins, J. J., and Thomas, T. I. (1996). S B Lowen, Malvin C Teich. Information measures quantifying aperiodic stochastic resonance. Phys. Rev. E 53, 2228–2231. doi: 10.1103/PhysRevE.54.R2228

CrossRef Full Text

Itoh, M., and Chua, L. O. (2008). Memristor oscillators. Int. J. Bifurcat. Chaos. 18, 3183–3206. doi: 10.1142/S0218127408022354

CrossRef Full Text | Google Scholar

Jia, B., and Gu, H. (2017). Dynamics and physiological roles of stochastic firing patterns near bifurcation points. Int. J. Bifurcat. Chaos 27:1750113. doi: 10.1142/s0218127417501139

CrossRef Full Text | Google Scholar

Jiao, X., and Wang, Y. (2010). Stochastic resonance of an integral-and-fire neuron model with threshold driven by synaptic noise. J. Dyn. Control 8, 273–276. doi: 10.3969/j.issn.1672-6553.2010.03.017. (In Chinese).

CrossRef Full Text

Kang, Y. M., Xu, J. X., and Xie, Y. (2003). Observing stochastic resonance in an underdamped bistable Duffing oscillator by the method of moments. Phys. Rev. E 68:036123. doi: 10.1103/PhysRevE.68.036123

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, Y. M., Xu, J. X., and Xie, Y. (2005). Signal-to-noise ratio gain of a noisy neuron that transmits subthreshold periodic spike trains. Phys. Rev. E 72:021902. doi: 10.1103/PhysRevE.72.021902

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, L., Jia, Y., Liu, W., and Yang, L. (2017). Mixed stimulus-induced mode selection in neural activity driven by high and low frequency current under electromagnetic radiation. Complexity 2017:7628537. doi: 10.1155/2017/7628537

CrossRef Full Text | Google Scholar

Lv, M., and Ma, J. (2016). Multiple modes of electrical activities in a new neuron model under electromagnetic radiation. Neurocomputing 205, 375–381. doi: 10.1016/j.neucom.2016.05.004

CrossRef Full Text | Google Scholar

Ma, J., Wang, Y., Wang, C., Xu, Y., and Ren, G. (2017). Mode selection in electrical activities of myocardial cell exposed to electromagnetic radiation. Chaos Solitons Fractals 99, 219–225. doi: 10.1016/j.chaos.2017.04.016

CrossRef Full Text | Google Scholar

Mark, D., McDonnell, and Derek, A. (2009). What is stochastic resonance? Definitions, misconceptions, debates, and its relevance to biology. PLoS Comput. Biol. 5:e1000348. doi: 10.1371/journal.pcbi.1000348

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagumo, J. S., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission like simulating nerve axon. Proc. IRE 50, 2061–2070. doi: 10.1109/JRPROC.1962.288235

CrossRef Full Text | Google Scholar

Ren, G., Xu, Y., and Wang, C. (2017). Synchronization behavior of coupled neuron circuits composed of memristors. Nonlin. Dyn. 88, 893–901. doi: 10.1007/s11071-016-3283-2

CrossRef Full Text | Google Scholar

Shao, Y., and Kang, Y. (2014). Effect of spatially correlated noise on stochastic synchronization in globally coupled FizHugh-Nagumo neuron systems. Theor. Appl. Mech. Lett. 4:013006. doi: 10.1063/2.1401306

CrossRef Full Text | Google Scholar

Strukov, D. B., Snider, G. S., Stewart, D. R., and Williams, R. S. (2008). The missing memristor found. Nature 453, 80–83. doi: 10.1038/nature06932

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., and Li, G. (2016). Stochastic multi-resonance induced by partial time delay in a Watts-Strogatz small-world neuronal network. Acta Phys. 65:120502. doi: 10.7498/aps.65.120502. (In Chinese).

CrossRef Full Text

Sun, X., and Shi, X. (2014). Effects of channel blocks on the spiking regularity in clustered neuronal networks. Technol. Sci. 57, 879–884. doi: 10.1007/s11431-014-5529-x

CrossRef Full Text | Google Scholar

Tanabe, and Pakdaman, K. (2001). Dynamics of moments of FitzHugh-Nagumo neuronal models and stochastic bifurcations. Phys. Rev. E 63:031911. doi: 10.1103/PhysRevE.63.031911

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, S., Zeng, Z., Huang, T., and Chen, Y. (2013). Fuzzy modeling and synchronization of different memristor-based chaotic circuits. Phys. Lett. A. 377, 2016–2021. doi: 10.1016/j.physleta.2013.05.046

CrossRef Full Text | Google Scholar

World Health Organization (2011). Electromagnetic Fields and Public Health Mobile Phones. Issues.

Wu, F., Wang, C., Xu, Y., and Ma, J. (2016). Model of electrical activity in cardiac tissue under electromagnetic induction. Sci. Rep. 6:28. doi: 10.1038/s41598-016-0031-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y., Chen, L., Kang, Y., and Aihara, K. (2008a). Controlling the onset of Hopf bifurcation in the Hodgkin-Huxley model. Phys. Rev. E 77:061921. doi: 10.1103/PhysRevE.77.061921

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y., Aihara, K., and Kang, Y. (2008b). Change in types of neuronal excitability via bifurcation control. Phys. Rev. E 77:021917. doi: 10.1103/PhysRevE.77.021917

PubMed Abstract | CrossRef Full Text | Google Scholar

Yogesh, N. J., and Stephen, J. W. (2009). The elusive memristor, properties of basic electrical circuits. Eur. J. Phys. 30, 661–675. doi: 10.1088/0143-0807/30/4/001

CrossRef Full Text | Google Scholar

Zha, J., Huang, H., Huang, T., and Cao, J. (2017). A general memristor model and its application in programmable analog circuits. Neurocomputing 267, 134–140. doi: 10.1016/j.neucom.2017.04.057

CrossRef Full Text | Google Scholar

Zhang, Q. (2005). Bifurcation and Chaos Theory And Application. Tianjin: Tianjin University Press. (In Chinese).