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

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, *k*_{B} is the Boltzmann constant, *T* is the absolute temperature (in K), *q* is the elementary charge, *k*_{B}*T*/*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:

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

*F*(*x*) essentially differs from the odd function *F*_{FH}(*x*) = *x*−*x*^{3}/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

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 (*x*_{0}, *y*_{0}) 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.

## 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 *C*_{0}, *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*_{0}.

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

Here,

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

Here,

There is a problem in Eq. (17) with the term <*βy*_{0}>*βy*_{0}> ≠ <*β*> < *y*_{0}>. 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*_{0}> ≈ <*β*> <*y*_{0}>. 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

where

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.001*i*, and *ωf* = 0.1, the threshold is *k*_{1} = 1.34.

The second minor Δ_{2} is more cumbersome and yields quadratic equation:

where

Eq. (26) has an analytical solution

which provides two different values. For *α* = 1.5, *βi* = 0.05 + 0.001*i*, and *ωf* = 0.1, the values are *k*_{2} = 1.44 and *k*_{3} = −12.3. Eventually, we evaluate the threshold *k*_{th} = max(*k*_{1}, *k*_{2}, *k*_{3}) = 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 Δ_{2} > 0, the inequality for the third minor Δ_{3} > 0 can be replaced simply with *h*_{0} > 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.

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

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

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

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

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* = −2*kN*, 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.

## 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}/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.