We present a case study of the FitzHugh–Nagumo (FHN) type model with a strongly asymmetric activation function. The proposed model is an electronically rather than a biologically inspired approach. The asymmetric exponential model imitates the shape of spikes in real neurons better than the classical FHN model with a cubic van der Pol activation function. An array of mean-field coupled non-identical FHN type oscillators is considered. The effect of mutual synchronization (phase locking) of units, originally oscillating at their individual frequencies, is demonstrated. Several feedback control methods, including stable tracking filter technique, mean field nullifying, and repulsive coupling are shown either to stabilize unstable equilibrium states or to suppress synchrony of the coupled FHN oscillators. The stability of the equilibrium states is analyzed by employing the eigenvalues, obtained from the characteristic equation, and by using the diagonal minors of the Routh–Hurwitz matrix. Nonlinear differential equations are solved numerically.
- nonlinear dynamics
- spiking neuron model
- FitzHugh−Nagumo oscillator
- arrays of coupled oscillators
- equilibrium states
- control methods
The stability of any either natural or artificial system is a valuable and desired property. Therefore, the control of dynamical systems, in particular stabilization of their unstable equilibrium (UEQ) states, is an important problem in basic science and engineering applications, if periodic or chaotic oscillations are unacceptable behaviors. Usual control methods, based on proportional feedback control [1, 2] require knowledge of a mathematical model of a dynamical system or at least the exact coordinates of the UEQ in the phase space for the reference point. However, in many real complex systems, especially in biology, physiology, economics, sociology, and chemistry neither the full reliable models nor the exact coordinates of the UEQ are a priori known. Moreover, the position of the UEQ may change with time because of external unknown and unpredictable forces. In these cases, adaptive, that is model-independent and reference-free methods, automatically tracing and stabilizing unknown UEQ, can be helpful [3, 4, 5].
Synchronization is a universal and very common phenomenon, widely observed in nature, science, engineering, and social life . Coupled oscillators and their arrays, exhibiting synchrony, range from pendulum clocks to various biological populations. In many cases, synchronization plays a positive role. However, sometimes, it has an unfavorable impact. Strong synchronization of neurons in human brain is an example. It is assumed that synchrony of spiking neurons in a neuronal population causes the symptoms of the Parkinson’s disease and essential tremor . Therefore, development of the methods and practical techniques for controlling, more specifically, for suppressing synchrony of coupled oscillators, in general, and particularly with possible application to neuronal arrays, is of great importance [8, 9, 10].
A variety of adaptive feedback methods for stabilizing UEQ of nonlinear dynamical systems have been described in literature. Here, we mention only some of them, e.g., derivative control technique [11, 12, 13], stable filter technique [3, 4, 14, 15, 16, 17], unstable filter technique [18, 19, 20], and combined filters techniques [21, 22, 23]. A comprehensive list and an overview of control methods developed to stabilize UEQ states can be found in . We note that the above mentioned techniques deal with single unstable dynamical systems. Stabilization of a network of coupled oscillators has been considered in a recent paper .
Suppression of synchrony in arrays of oscillators by means of feedback methods has been described in many papers [7, 8, 9, 10, 26, 27, 28, 29]. More publications and discussion on the feedback techniques for control of synchrony are presented in [24, 25].
Another way to avoid synchrony in arrays of oscillators is a non-feedback method using external periodic drive at relatively high frequency (much higher than the natural frequency of the oscillators). In neurology, it is known as deep brain stimulation (DBS), applying about 150 Hz periodic pulses to certain brain areas . It is a clinically approved therapy for patients with the Parkinson’s disease symptoms. However, mechanism of the DBS is not fully understood. There are several papers considering the Hodgkin–Huxley and the FitzHugh–Nagumo models and demonstrating that high frequency forcing can stabilize the UEQ of the neuronal oscillators and thus inhibit spiking cells [31, 32, 33].
In this chapter, we present a case study of the FitzHugh–Nagumo (FHN) type model with a strongly asymmetric activation function (Section 2). An array of mean-field coupled non-identical FHN type oscillators is considered in Section 3. The effect of mutual synchronization (phase locking) of units originally oscillating at their individual frequencies is demonstrated. Several feedback control methods, including stable tracking filter technique (Section 4), mean field nullifying (Section 5), and repulsive coupling (Section 6) are shown either to stabilize UEQ states or to suppress synchrony of the coupled FHN type oscillators.
2. Single FHN type oscillator
An extremely simple electrical circuit, imitating a single spiking neuron, is sketched in Figure 1. The negative resistance Rn can be implemented by means of a negative impedance converter . Typical train of spikes from its output is presented in Figure 2.
We apply the Kirchhoff’s laws to electrical circuit in Figure 1, use the Shockley current–voltage characteristic for the diode, and introduce the following dimensionless quantities:
where V* = 1 V, kB is the Boltzmann constant, T is the absolute temperature (in K), q is the elementary charge, kBT/q is the thermal potential (≈ 25 mV at room temperature, T = 293 K), n is a diode ideality factor, sometimes called emission coefficient (assumed value n = 2). Then, differential equations, convenient for analysis and numerical integration, are derived:
F(x) essentially differs from the odd function FFH(x) = x−x3/3, introduced by FitzHugh  and used in many later papers, e.g., in . It also differs from the asymmetric three-segment [x < −1, −1≤ x ≤ 1, x > 1] piecewise linear function FPL(x) = αx + d(x + 1)H(−x–1) + g(x–1)H(x–1) suggested in , where d >> g and H(u) is the Heaviside unit step function, i.e., H(u > 0) = 1, H(u ≤ 0) = 0. In contrast to the FPL(x), the F(x) is a smooth function, and therefore it seems a more realistic option.
For αβ < 1 and
the equilibrium solution of Eq. (2) is given by the fixed point coordinates
Due to the exponent in the activation function F(x), strong inequality (4) practically can be replaced with a simple inequality:
Note empiric factor 2 is in the denominator. Eqs. (2), linearized around the fixed point (5), read
The corresponding characteristic equation is
It has two eigenvalues
If α > β, then both real parts of the eigenvalues, Reλ1,2 are positive, proving that the equilibrium (x0, y0) is an unstable fixed point. If α + β > 2, it is a node, and if α + β < 2, it is a spiral.
3. Array of FHN type oscillators
An array of isolated (non-coupled) oscillators is given by
Here and elsewhere i = 1, 2, …, N. Note that the structure of function F(x) and parameters α, δ, and μ are the same for all oscillators, whereas the damping parameters βi in Eq. (10) are intentionally set different for each oscillator to make them slightly non-identical units.
Now we introduce interaction between oscillators. To be specific, we consider mean-field coupling, which is also called “star” coupling in electronics (Figure 5):
Here, k = ρ/R* is the strength of coupling and
Typical phase portraits for isolated and coupled (synchronized) oscillators are shown in Figure 6.
Intricate phase trajectories in Figure 6 (left) indicate that the oscillators are not synchronized, but oscillate at their individual frequencies, whereas simple closed loop in Figure 6 (right) proves that oscillators are in synchrony (phase-locked), i.e., oscillate at the same frequency. For synchronized oscillators, the phase difference is not necessarily zero (the phase portrait is not fine diagonal), but it does not change with time. The mean variable <x> for the two cases is shown in Figure 7. The amplitude of mean-field variable <x> is relatively low for isolated oscillators (k = 0), but becomes large for synchronized state (k = 1).
4. Stabilizing equilibrium states in array of oscillators
When an external capacitor is applied to the coupling node CN (Figure 8), the overall system becomes (2 N + 1)-dimensional system:
Here, z is a dimensionless dynamical variable related to voltage across the external capacitor C0, z = VC0/V*, the mean <x> is given by formula (13), and ωf is the dimensionless cut-off frequency of the filter composed by R* and C0.
Analysis of the high-dimensional system is very complicated. Therefore, we consider a mean-field approach. We average all terms in Eq. (14) over all oscillators i = 1, 2, …, N:
Eq. (15) is not suitable to describe full dynamics of the system. However, we can exploit it to find equilibrium coordinates. If inequality (6) is valid for all oscillators with different βi, the steady-state equations read
There is a problem in Eq. (17) with the term <βy0>. In general, <βy0> ≠ <β> < y0>. However, if the ranges of the multiplicands βi and y0i in (18) are considerably different (in our case, the individual equilibrium coordinates y0i, in comparison with βi, much weaker depend on i), then <βy0> ≈ <β> <y0>. Similarly to a single oscillator, considered in Section 2, for αβi < 1, the equilibrium coordinates are
Linearization of Eqs. (15) around the equilibrium coordinates yields:
The corresponding characteristic equation is
Numerical solution of Eq. (21) is presented in Figure 9 for different values of the coupling parameter k. The equilibrium is stable, if the real parts of all three eigenvalues are negative, Reλ1,2,3 < 0.
Necessary and sufficient conditions of stability can be found analytically from the Hurwitz matrix
The Routh–Hurwitz stability criterion claims that the system is stable, if all diagonal minors of the matrix H are positive:
The first minor is Δ1 > 0, if
For α = 1.5, βi = 0.05 + 0.001i, and ωf = 0.1, the threshold is k1 = 1.34.
The second minor Δ2 is more cumbersome and yields quadratic equation:
Eq. (26) has an analytical solution
which provides two different values. For α = 1.5, βi = 0.05 + 0.001i, and ωf = 0.1, the values are k2 = 1.44 and k3 = −12.3. Eventually, we evaluate the threshold kth = max(k1, k2, k3) = 1.44. It is in a very good agreement with numerical value of kth obtained from Reλ1,2,3(k) in Figure 9.
Once Δ2 > 0, the inequality for the third minor Δ3 > 0 can be replaced simply with h0 > 0. This can be further simplified to (1–α<β>) > 0, since ωf > 0 by definition. Finally, we come to inequality α<β> <1, which satisfied by itself, because it was already used as an assumption to derive the equilibrium coordinates; see formulas (19).
5. Mean-field “nullifying” technique
A straightforward way to desynchronize the mean-field coupled oscillators is to “nullify” the mean field at the coupling node CN, i.e., to remove the reason of synchronization. The corresponding diagram is shown in Figure 11.
and emphasize that the mean-field value <x> by itself is not zero:
The control technique implicates that the mean-field variable <x> is not fully nullified, but its value at the coupling node CN, <x> CN is set zero:
The coupling node is simply grounded, as sketched in Figure 11. Numerical results are shown in Figure 12. Note that when the control is switched on, the value of actual mean-field variable <x> becomes relatively small, but is not zero.
6. Repulsive coupling technique
Voltage at the coupling node xCN is found from the Kirchhoff’s law for current:
Eq. (31) yields:
Evidently, for Rn = 0, the G → ∞ and xCN = 0, as expected. It is the case considered in previous Section 5. If Rn is negative, say it provides value of G = −2kN, then xCN = − <x>. It is the case of so-called repulsive coupling:
Numerical results are presented in Figure 14. Similarly to the mean-field “nullifying” technique, the mean <x> becomes small, which is a typical feature of either non-synchronized or antiphase synchronized oscillators.
A modification of the FitzHugh–Nagumo (FHN) model of a spiking neuron has been proposed. In the original model, developed by FitzHugh , the cubic activation function x–x3/3 has been replaced with a strongly asymmetric exponential one. This function provides more realistic shape of the membrane voltage spikes. Synchronization effect in an array of mean-field coupled non-identical FHN type oscillators has been demonstrated.
Three methods for controlling arrays of coupled FHN type oscillators have been described:
Stable filter technique aimed to damp spikes in coupled oscillators. It is based on replacing the mean variable <x> at the coupling node with its filtered version z.
Mean field nullifying technique, <x> = 0 (grounding the coupling node).
Repulsive coupling technique, following the idea described in  and shown for an array of Kuramoto 1D phase oscillators. It is based on replacing the mean-field variable <x> at the coupling node with the inverse version “– <x>.”
The above control techniques have different physical mechanisms behind, ranging from stabilization of the equilibrium states to desynchronization and antiphase synchronization. However, all of them ensure low value of mean-field variable in the array.