Controlling Equilibrium and Synchrony in Arrays of FitzHugh– Nagumo Type Oscillators

Abstract

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.

1. Introduction

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 [6]. 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 [7]. 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 [24]. 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 [25].

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 [30]. 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 R_n can be implemented by means of a negative impedance converter [34]. Typical train of spikes from its output is presented in Figure 2.

Figure 1.

Circuit diagram of the electronic analog of spiking neuron. R_n is a negative resistance.

Figure 2.

Voltage spikes from the circuit in Figure 1, generated by means of Electronics Workbench Professional software. R_n = − 680 Ω, D is a semiconductor diode (BAV99 type) with saturation current I_S = 10 nA (δ = 10⁻⁵), L = 100 mH, C = 100 nF, (ρ = 1 kΩ), r = 50 Ω (β = 0.05), I = 1 mA (γ = 1).

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:

$$x=\frac{V_C}{V^{\ast }},y=\frac{\rho I_L}{V^{\ast }},t\to \frac{t}{\sqrt{\mathit{LC}}},\alpha =\frac{\rho }{\left|R_n\right|},\beta =\frac{r}{\rho },\gamma =\frac{\mathit{I\rho }}{V^{\ast }},\delta =\frac{I_S\rho }{V^{\ast }},\mu =\frac{{\mathit{qV}}^{\ast }}{{\mathit{nk}}_{B}T},\rho =\sqrt{\frac{L}{C}},$$E1

where V^* = 1 V, k_B is the Boltzmann constant, T is the absolute temperature (in K), q is the elementary charge, k_BT/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:

$$\begin{array}{l}\dot{x}=F\left(x\right)-y-\gamma ,\\ \dot{y}=x-\mathit{\beta y}.\end{array}$$E2

Activation function F(x) in Eq. (2) is a strongly asymmetric one (Figure 3):

$$F\left(x\right)=\mathit{\alpha x}+\delta \left[\text{exp}\left(-\mathit{\mu x}\right)-1\right].$$E3

Figure 3.

Activation function F(x) from formula (3). α = 1.5, δ = 10⁻⁵, and μ = 20. Black dot on the curve marks the equilibrium coordinate x₀ = −0.12 from formula (5) at β = 0.1 and γ = 1.

F(x) essentially differs from the odd function F_FH(x) = x−x³/3, introduced by FitzHugh [35] and used in many later papers, e.g., in [28]. It also differs from the asymmetric three-segment [x < −1, −1≤ x ≤ 1, x > 1] piecewise linear function F_PL(x) = αx + d(x + 1)H(−x–1) + g(x–1)H(x–1) suggested in [36], 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 F_PL(x), the F(x) is a smooth function, and therefore it seems a more realistic option.

For αβ < 1 and

$$\gamma <<\frac{1-\mathit{\alpha \beta }}{\mathit{\mu \beta }}\text{ln}{\delta }^{-1}$$E4

the equilibrium solution of Eq. (2) is given by the fixed point coordinates

$$x_0=-\frac{\mathit{\beta \gamma }}{1-\mathit{\alpha \beta }},y_0=-\frac{\gamma }{1-\mathit{\alpha \beta }}.$$E5

Due to the exponent in the activation function F(x), strong inequality (4) practically can be replaced with a simple inequality:

$$\gamma \le \frac{1-\mathit{\alpha \beta }}{2\mathit{\mu \beta }}\text{ln}{\delta }^{-1}.$$E6

Note empiric factor 2 is in the denominator. Eqs. (2), linearized around the fixed point (5), read

$$\begin{array}{l}\dot{x}=\mathit{\alpha x}-y,\\ \dot{y}=x-\mathit{\beta y}.\end{array}$$E7

The corresponding characteristic equation is

$${\lambda }^2-\left(\alpha -\beta \right)\lambda +1-\mathit{\alpha \beta }=0.$$E8

It has two eigenvalues

$${\lambda }_{1,2}=\frac{\left(\alpha -\beta \right)}{2}\pm \sqrt{\frac{{\left(\alpha -\beta \right)}^2}{4}-\left(1-\mathit{\alpha \beta }\right)}=\frac{\left(\alpha -\beta \right)}{2}\pm \sqrt{\frac{{\left(\alpha +\beta \right)}^2}{4}-1}.$$E9

If α > β, then both real parts of the eigenvalues, Reλ_1,2 are positive, proving that the equilibrium (x₀, y₀) is an unstable fixed point. If α + β > 2, it is a node, and if α + β < 2, it is a spiral.

Numerical solution of nonlinear equation Eq. (2) is presented in Figure 4.

Figure 4.

Waveforms x(t) from Eq. (2) with α = 1.5, γ = 1, δ = 10⁻⁵, and μ = 20 for different damping β. (Top) β = 0.05 and (bottom) β = 0.1. Note, different inter-spike periods in the two plots.

3. Array of FHN type oscillators

An array of isolated (non-coupled) oscillators is given by

$$\begin{array}{l}{\dot{x}}_i=F\left(x_i\right)-y_i-\gamma ,\\ {\dot{y}}_i=x_i-{\beta }_i{y}_i,\end{array}$$E10

$$F\left(x_i\right)=\alpha x_i+\delta \left[\text{exp}\left(-\mu x_i\right)-1\right].$$E11

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):

$$\begin{array}{l}{\dot{x}}_i=F\left(x_i\right)-y_i-\gamma +k\left(⟨x⟩-x_i\right),\\ {\dot{y}}_i=x_i-{\beta }_i{y}_i,\end{array}$$E12

Figure 5.

Diagram of mean-field coupled oscillators. R* are coupling resistors and CN is a coupling node.

Here, k = ρ/R* is the strength of coupling and

$$⟨x⟩=\frac{1}{N}\sum _{i=1}^N{x}_i.$$E13

Typical phase portraits for isolated and coupled (synchronized) oscillators are shown in Figure 6.

Figure 6.

Phase portraits. N = 24, α = 1.5, β _i = 0.05 + 0.001i, γ = 1, δ = 10⁻⁵, and μ = 20. (Left) Isolated oscillators either from Eq. (10) or Eq. (12) with k = 0, and (right) coupled oscillators from Eq. (12) with k = 1.

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).

Figure 7.

Waveform of the mean-field variable <x> from Eq. (12). N = 24, α = 1.5, β _i = 0.05 + 0.001i, γ = 1, δ = 10⁻⁵, and μ = 20. Coupling (k = 1) is switched on at t = 100.

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:

$$\begin{array}{l}{\dot{x}}_i=F\left(x_i\right)-y_i-\gamma +k\left(z-x_i\right),\\ {\dot{y}}_i=x_i-{\beta }_i{y}_i,\\ \dot{z}={\omega }_f\left(⟨x⟩-z\right).\end{array}$$E14

Figure 8.

Diagram of mean-field coupled oscillators with a stabilizing capacitor C₀. Stable RC filter is composed of coupling resistors R* and capacitor C₀ (see formulas (16)).

Here, z is a dimensionless dynamical variable related to voltage across the external capacitor C₀, z = V_C0/V*, the mean <x> is given by formula (13), and ω_f is the dimensionless cut-off frequency of the filter composed by R* and C₀.

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:

$$\begin{array}{c}⟨\dot{x}⟩=⟨F⟩-⟨y⟩-\gamma +k\left(z-⟨x⟩\right),\\ ⟨\dot{y}⟩=⟨x⟩-⟨\mathit{\beta y}⟩,\\ \dot{z}={\omega }_f\left(⟨x⟩-z\right).\end{array}$$E15

Here,

$$⟨x⟩=\frac{1}{N}\sum _{i=1}^N{x}_i,⟨y⟩=\frac{1}{N}\sum _{i=1}^N{y}_i,⟨F⟩=\frac{1}{N}\sum _{i=1}^{N}F\left(x_i\right),⟨\mathit{\beta y}⟩=\frac{1}{N}\sum _{i=1}^N{\beta }_i{y}_i,{\omega }_f=\frac{N\sqrt{\mathit{LC}}}{R^{\ast }{C}_0}.$$E16

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

$$\begin{array}{l}0=⟨x_0⟩-⟨y_0⟩-\gamma +k\left(z_0-⟨x_0⟩\right),\\ 0=⟨x_0⟩-⟨\beta y_0⟩,\\ 0=⟨x_0⟩-z_0.\end{array}$$E17

Here,

$$⟨x_0⟩=\frac{1}{N}\sum _{i=1}^N{x}_{0i},⟨y_0⟩=\frac{1}{N}\sum _{i=1}^N{y}_{0i},⟨\beta y_0⟩=\frac{1}{N}\sum _{i=1}^N{\beta }_i{y}_{0i},⟨\beta ⟩=\frac{1}{N}\sum _{i=1}^N{\beta }_i.$$E18

There is a problem in Eq. (17) with the term <βy₀>$⟨\beta y_0⟩$. In general, <βy₀> ≠ <β> < y₀>. However, if the ranges of the multiplicands β_i and y_0i in (18) are considerably different (in our case, the individual equilibrium coordinates y_0i, in comparison with β_i, much weaker depend on i), then <βy₀> ≈ <β> <y₀>. Similarly to a single oscillator, considered in Section 2, for αβ_i < 1, the equilibrium coordinates are

$$⟨x_0⟩=-\frac{⟨\beta ⟩\gamma }{1-\alpha ⟨\beta ⟩},⟨y_0⟩=-\frac{\gamma }{1-\alpha ⟨\beta ⟩},z_0=⟨x_0⟩.$$E19

Linearization of Eqs. (15) around the equilibrium coordinates yields:

$$\begin{array}{c}⟨\dot{x}⟩=\alpha ⟨x⟩-⟨y⟩+k\left(z-⟨x⟩\right),\\ ⟨\dot{y}⟩=⟨x⟩-⟨\beta ⟩⟨y⟩,\\ \dot{z}={\omega }_f\left(⟨x⟩-z\right).\end{array}$$E20

The corresponding characteristic equation is

$${\lambda }^3+h_2{\lambda }^2+h_1\lambda +h_0=0,$$E21

where

$$h_2=-\alpha +⟨\beta ⟩+k+{\omega }_f,h_1=1-\alpha ⟨\beta ⟩+⟨\beta ⟩k-\left(\alpha -⟨\beta ⟩\right){\omega }_f,h_0=\left(1-\alpha ⟨\beta ⟩\right){\omega }_f.$$E22

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.

Figure 9.

Real parts of the eigenvalues from Eq. (21). N = 24, α = 1.5, β_i = 0.05 + 0.001i, <β> = 0.0625, and ω_f = 0.1. Arrow in the plot indicates the threshold coupling parameter k_th = 1.44, where the largest eigenvalues become negative.

Necessary and sufficient conditions of stability can be found analytically from the Hurwitz matrix

$$H=\left(\begin{array}{lll}{h}_2& h_0& 0\\ 1& h_1& 0\\ 0& h_2& h_0\end{array}\right).$$E23

The Routh–Hurwitz stability criterion claims that the system is stable, if all diagonal minors of the matrix H are positive:

$${\mathrm{\Delta }}_1=h_2>0,{\mathrm{\Delta }}_2=h_2{h}_1-h_0>0,{\mathrm{\Delta }}_3=h_0{\mathrm{\Delta }}_2>0.$$E24

The first minor is Δ₁ > 0, if

$$k>k_1=\alpha -⟨\beta ⟩-{\omega }_f.$$E25

For α = 1.5, β_i = 0.05 + 0.001i, and ω_f = 0.1, the threshold is k₁ = 1.34.

The second minor Δ₂ is more cumbersome and yields quadratic equation:

$$⟨\beta ⟩k^2+\mathit{dk}+g=0.$$E26

where

$$\begin{array}{l}d=1-2\alpha ⟨\beta ⟩+{⟨\beta ⟩}^2-\left(\alpha -2⟨\beta ⟩\right){\omega }_f,\\ g=-\left(\alpha -⟨\beta ⟩\right)\left[1-\alpha ⟨\beta ⟩-\left(\alpha -⟨\beta ⟩\right){\omega }_f+{\omega }_f^2\right].\end{array}$$E27

Eq. (26) has an analytical solution

$$k_{2,3}=-\frac{d}{2⟨\beta ⟩}\pm \sqrt{\frac{d^2}{4{⟨\beta ⟩}^2}-\frac{g}{⟨\beta ⟩}},$$E28

which provides two different values. For α = 1.5, β_i = 0.05 + 0.001i, and ω_f = 0.1, the values are k₂ = 1.44 and k₃ = −12.3. Eventually, we evaluate the threshold k_th = max(k₁, k₂, k₃) = 1.44. It is in a very good agreement with numerical value of k_th obtained from Reλ_1,2,3(k) in Figure 9.

Once Δ₂ > 0, the inequality for the third minor Δ₃ > 0 can be replaced simply with h₀ > 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).

Numerical results from Eqs. (14), demonstrating dynamics of equilibrium stabilization, are presented in Figure 10.

Figure 10.

Waveforms from Eq. (14). N = 24, α = 1.5, β_i = 0.05 + 0.001i, γ = 1, δ = 10⁻⁵, μ = 20, ω_f = 0.1, and k = 1.6. (Top) Mean-field variable <x>, (bottom) control term z – <x>. Control is switched on at t = 300 (<x> in the coupling term is replaced with z).

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.

Figure 11.

Diagram of mean-field coupled oscillators with the coupling node CN grounded.

We repeat here Eq. (12) from Section 3 for clarity and for comparison with Eq. (30):

$$\begin{array}{l}{\dot{x}}_i=F\left(x_i\right)-y_i-\gamma +k\left(⟨x⟩-x_i\right),\\ {\dot{y}}_i=x_i-{\beta }_i{y}_i\end{array}$$

and emphasize that the mean-field value <x> by itself is not zero:

$$⟨x⟩=\frac{1}{N}\sum _{i=1}^N{x}_i\ne 0.$$E29

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:

$$\begin{array}{l}{\dot{x}}_i=F\left(x_i\right)-y_i-\gamma +k\left(0-x_i\right),\\ {\dot{y}}_i=x_i-{\beta }_i{y}_i.\end{array}$$E30

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.

Figure 12.

Waveform of the mean-field variable <x> from Eq. (30). N = 24, α = 1.5, β_i = 0.05 + 0.001i, γ = 1, δ = 10⁻⁵, μ = 20, and k = 1. Control is switched on at t = 300 (x_CN = <x> in the coupling term is replaced with x_CN = 0).

6. Repulsive coupling technique

An alternative method of desynchronization of coupled oscillators is the repulsive coupling, also called “repulsive synchronization” technique [26]. Diagram is sketched in Figure 13.

Figure 13.

Diagram of mean-field coupled oscillators with coupling node CN, grounded via resistor R_n.

Voltage at the coupling node x_CN is found from the Kirchhoff’s law for current:

$$\sum _{i=1}^{N}k\left(x_i-x_{\mathit{CN}}\right)-{\mathit{Gx}}_{\mathit{CN}}=0,$$E31

$$x_{\mathit{CN}}=\frac{V_{\mathit{CN}}}{V^{\ast }},G=\frac{\rho }{R_n}.$$E32

Eq. (31) yields:

$$x_{\mathit{CN}}=\mathit{kN}⟨x⟩/\left(\mathit{kN}+G\right).$$E33

Evidently, for Rn = 0, the G → ∞ and x_CN = 0, as expected. It is the case considered in previous Section 5. If R_n is negative, say it provides value of G = −2kN, then x_CN = − <x>. It is the case of so-called repulsive coupling:

$$\begin{array}{l}{\dot{x}}_i=F\left(x_i\right)-y_i-\gamma +k\left(-⟨x⟩-x_i\right),\\ {\dot{y}}_i=x_i-{\beta }_i{y}_i.\end{array}$$E34

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.

Figure 14.

Waveform of the mean-field variable <x> from Eq. (34). N = 24, α = 1.5, β_i = 0.05 + 0.001i, γ = 1, δ = 10⁻⁵, μ = 20, and k = 1. Control is switched on at t = 300 (x_CN = <x> in the coupling term is replaced with x_CN = − < x>).

7. Conclusions

A modification of the FitzHugh–Nagumo (FHN) model of a spiking neuron has been proposed. In the original model, developed by FitzHugh [35], the cubic activation function x–x³/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 [26] 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.