Open access peer-reviewed chapter

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

By Elena Adomaitienė, Skaidra Bumelienė and Arūnas Tamaševičius

Submitted: December 5th 2017Reviewed: January 24th 2018Published: February 16th 2018

DOI: 10.5772/intechopen.74337

Downloaded: 245

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.

Keywords

  • nonlinear dynamics
  • spiking neuron model
  • FitzHugh−Nagumo oscillator
  • arrays of coupled oscillators
  • equilibrium states
  • synchronization
  • control methods

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 Rn 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. Rn is a negative resistance.

Figure 2.

Voltage spikes from the circuit in Figure 1, generated by means of Electronics Workbench Professional software. Rn = − 680 Ω, D is a semiconductor diode (BAV99 type) with saturation current IS = 10 nA (δ = 10−5), 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=VCV,y=ρILV,ttLC,α=ρRn,β=rρ,γ=V,δ=ISρV,μ=qVnkBT,ρ=LC,E1

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:

ẋ=Fxyγ,ẏ=xβy.E2

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

Fx=αx+δexpμx1.E3

Figure 3.

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

F(x) essentially differs from the odd function FFH(x) = xx3/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 FPL(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 FPL(x), the F(x) is a smooth function, and therefore it seems a more realistic option.

For αβ < 1 and

γ<<1αβμβlnδ1E4

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

x0=βγ1αβ,y0=γ1αβ.E5

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

γ1αβ2μβlnδ1.E6

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

ẋ=αxy,ẏ=xβy.E7

The corresponding characteristic equation is

λ2αβλ+1αβ=0.E8

It has two eigenvalues

λ1,2=αβ2±αβ241αβ=αβ2±α+β241.E9

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.

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−5, 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

ẋi=Fxiyiγ,ẏi=xiβiyi,E10
Fxi=αxi+δexpμxi1.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):

ẋi=Fxiyiγ+kxxi,ẏi=xiβiyi,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=1Ni=1Nxi.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−5, 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−5, 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:

ẋi=Fxiyiγ+kzxi,ẏi=xiβiyi,ż=ωfxz.E14

Figure 8.

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

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:

ẋ=Fyγ+kzx,ẏ=xβy,ż=ωfxz.E15

Here,

x=1Ni=1Nxi,y=1Ni=1Nyi,F=1Ni=1NFxi,βy=1Ni=1Nβiyi,ωf=NLCRC0.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

0=x0y0γ+kz0x0,0=x0βy0,0=x0z0.E17

Here,

x0=1Ni=1Nx0i,y0=1Ni=1Ny0i,βy0=1Ni=1Nβiy0i,β=1Ni=1Nβi.E18

There is a problem in Eq. (17) with the term <βy0>β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

x0=βγ1αβ,y0=γ1αβ,z0=x0.E19

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

ẋ=αxy+kzx,ẏ=xβy,ż=ωfxz.E20

The corresponding characteristic equation is

λ3+h2λ2+h1λ+h0=0,E21

where

h2=α+β+k+ωf,h1=1αβ+βkαβωf,h0=1αβω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 kth = 1.44, where the largest eigenvalues become negative.

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

H=h2h001h100h2h0.E23

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

Δ1=h2>0,Δ2=h2h1h0>0,Δ3=h0Δ2>0.E24

The first minor is Δ1 > 0, if

k>k1=αβωf.E25

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:

βk2+dk+g=0.E26

where

d=12αβ+β2α2βωf,g=αβ1αβαβωf+ωf2.E27

Eq. (26) has an analytical solution

k2,3=d2β±d24β2gβ,E28

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

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−5, μ = 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):

ẋi=Fxiyiγ+kxxi,ẏi=xiβiyi

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

x=1Ni=1Nxi0.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:

ẋi=Fxiyiγ+k0xi,ẏi=xiβiyi.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−5, μ = 20, and k = 1. Control is switched on at t = 300 (xCN = <x> in the coupling term is replaced with xCN = 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 Rn.

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

i=1NkxixCNGxCN=0,E31
xCN=VCNV,G=ρRn.E32

Eq. (31) yields:

xCN=kNx/kN+G.E33

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:

ẋi=Fxiyiγ+kxxi,ẏi=xiβiyi.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−5, μ = 20, and k = 1. Control is switched on at t = 300 (xCN = <x> in the coupling term is replaced with xCN = − < 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 xx3/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.

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Elena Adomaitienė, Skaidra Bumelienė and Arūnas Tamaševičius (February 16th 2018). Controlling Equilibrium and Synchrony in Arrays of FitzHugh– Nagumo Type Oscillators, Nonlinear Systems - Modeling, Estimation, and Stability, Mahmut Reyhanoglu, IntechOpen, DOI: 10.5772/intechopen.74337. Available from:

chapter statistics

245total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Appell-Gibbs Approach in Dynamics of Non-Holonomic Systems

By Jiří Náprstek and Cyril Fischer

Related Book

First chapter

On Nonoscillatory Solutions of Two-Dimensional Nonlinear Dynamical Systems

By Elvan Akın and Özkan Öztürk

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us