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.

Let us start with the conventional FHN neuron model (Fitzhugh, 1961; Nagumo et al., 1962). (1){v˙=v(v−a)(1−v)−w+Iextw˙=ε(v−dw) 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)=dφ(q)dq. In reverse, a flux-controlled memristor should be described by i(t) = W(φ(t))v(t) with W(φ)=dq(φ)dφ. 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(φ)=dq(φ)dφ=k(α+3βφ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 v˙=v(v−a)(1−v)−w+k(α+3βφ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, φ˙=k1v−k2φ+φ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. (2){v˙=v(v−a)(1−v)−w+k(α+3βφ2)vw˙=ε(v−dw)φ˙=k1v−k2φ+φext. 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.

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 {v0(v0−a)(1−v0)−w0+k(α+3βφ02)v0=0v0−dw0=0k1v0−k2φ0+φext=0 and thus v0[(v0−a)(1−v0)−1d+k(α+3β(k1v0+φextk2)2)]=0 with three solutions given by v01=0,v02=−B+B2−4AC2A,v03=−B−B2−4AC2A where A=3kk12βk22-1     B=6kk1βφextk22+1+a     C= 3kβφext2k22-a-1d+kα. Let E₀₁, E₀₂ and E₀₃ to be the resultant equilibrium points corresponding to v₀₁, v₀₂ and v₀₃, then E01=(0,0,φext/k2)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.

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,φ)=[−3v2+2(1+a)v−a+k(α+3βφ2)−16kβφvε−εd0k10−k2], One can easily obtain the characteristic determinant of system (2) at E₀₁ as |λI3−J(0,0,φextk2)| =|λ−A110−ελ+εd0−k10λ+k2|                                         =(λ+k2)[(λ−A1)(λ+εd)+ε]=0 with A1=-a+k(α+3β(φextk2)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 −ω2±(εd−A1)ω i−A1εd+ε=0 i.e., {−ω2−A1εd+ε=0εd−A1=0, which leads to A₁ = εd and ω=ε-ε2d2. Thus, the branch of the equilibrium point E₀₁ undergoes two Hopf bifurcations at parameter φext=±εd+a−kα3kβk2 which can be reduced to φext(1)=-2.381 and φ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 φext(1), and another hysteresis loop can be captured near φext(2) by some partial amplification., therefore the Hopf bifurcations at A and B are both subcritical.

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.

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 φext=r sin(2πft)+φext0+ξ(t) where rsin(2πft)+φext0 is a subthreshold signal of bias φext0, 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 (3){v˙=v(v−a)(1−v)−w+k(α+3βφ2)vw˙=ε(v−dw)φ˙=k1v−k2φ+r sin(2πft)+φext0+ξ(t) 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=10log10S(ω)N(ω) 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 φext0, 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.

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.

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 θ(D)=2π×BI(D)IBI(D)=2π×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.

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.

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.

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.