Document [original]

Chaos ARTICLE pubs.aip.org/aip/cha

Role of coupling delay in oscillatory activity in

autonomous networks of excitable neurons with

dissipation

Cite as: Chaos 33, 073114 (2023); doi: 10.1063/5.0147883

Submitted: 27 February 2023 ·Accepted: 19 June 2023 ·

Published Online: 5 July 2023

View Online

Export Citation

CrossMark

A. V. Bukh,1I. A. Shepelev,1,2 E. M. Elizarov,1S. S. Muni,3E. Schöll,4and G. I. Strelkova1,a)

AFFILIATIONS

1Institute of Physics, Saratov State University, 83 Astrakhanskaya Street, Saratov 410012, Russia

2Almetyevsk State Petroleum Institute, 2 Lenina Street, Almetyevsk 423458, Tatarstan, Russia

3Department of Physical Sciences, Indian Institute of Science Education and Research Kolkata, Campus Road, Mohanpur,

West Bengal 741246, India

4Institut für Theoretische Physik, Technische Universität Berlin, 10623 Berlin, Germany; Bernstein Center for Computational

Neuroscience Berlin, 10115 Berlin, Germany; and Potsdam Institute for Climate Impact Research, 14473 Potsdam, Germany

Note: Control of self-organizing nonlinear systems.

a)Author to whom correspondence should be addressed: strelko[email protected]

ABSTRACT

We study numerically effects of time delay in networks of delay-coupled excitable FitzHugh–Nagumo systems with dissipation. Generation

of periodic self-sustained oscillations and its threshold are analyzed depending on the dissipation of a single neuron, the delay time, and

random initial conditions. The peculiarities of spatiotemporal dynamics of time-delayed bidirectional ring-structured FitzHugh–Nagumo

neuronal systems are investigated in cases of local and nonlocal coupling topology between the nodes, and a first-order nonequilibrium phase

transition to synchrony is established. It is shown that the emergence of an oscillatory activity in delay-coupled FitzHugh–Nagumo neurons is

observed for smaller values of the coupling strength as the dissipation parameter decreases. This can provide the possibility of controlling the

spatiotemporal behavior of the considered neuronal networks. The observed effects are quantified by plotting distributions of the maximal

Lyapunov exponent and the global order parameter in terms of delay and coupling strength.

Published under an exclusive license by AIP Publishing. https://doi.org/10.1063/5.0147883

Excitability is a common property of many physical and biolog-

ical systems. Systems consisting of coupled excitable elements

have been widely studied as models for many natural phenomena,

such as the communication between neurons in the brain. Since

the work of Hodgkin and Huxley1and the development of the

basic mathematical model by FitzHugh2and Nagumo et al.,3the

reported research on the subject has grown enormously. Delays

are inherent in neuronal networks due to finite conduction veloc-

ities and synaptic transmission. Recently, many researchers have

investigated the effects of time delays in neuronal networks and

found many delay-induced phenomena. In particular, it has been

established that time-delayed coupling can be used as a pow-

erful tool for stabilizing various complex spatiotemporal pat-

terns and for controlling different types of synchronization in

one- and multilayer networks. In the present work, we consider

the FitzHugh–Nagumo oscillator that represents a paradigmatic

model for neuronal excitability. The considered two-variable sys-

tem of equations includes an additional parameter, which takes

into account the dissipation of a neuron and, thus, differs from

the simplified FitzHugh–Nagumo model. Starting from a pair of

delay-coupled FitzHugh–Nagumo oscillators and proceeding via

simple ring neuronal networks with delayed coupling, we analyze

how the spatiotemporal dynamics of the networks depends on

the parameter of dissipation, the delay time, the coupling param-

eters, and randomly distributed initial conditions. We provide

linear stability analyses and plot distributions of the maximal

Lyapunov exponent and the global order parameter depending

on the delay time and the coupling strength. Furthermore, we

establish nonequilibrium phase transitions from multiple clus-

ters to full synchronization. The presented results are compared

with the previously obtained findings for networks of simplified

FitzHugh–Nagumo models.

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-1

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

I. INTRODUCTION

Excitable systems play important roles in understanding and

modeling various natural phenomena, such as, for example, the

transmission of impulses between neurons in the brain, the car-

diac arrhythmia, the appearance of organized structures in the

cortex of egg cells, etc.4–6Generation of a single spike in the elec-

trical potential across the neuron membrane is a typical example

of excitable behavior. Such excitable units usually appear as con-

stitutive elements of complex systems and can transmit excitations

between them. The phenomenological FitzHugh–Nagumo set of

equations, which is a two-dimensional simplification of the four-

variable Hodgkin–Huxley model, has been proven to be a successful

model for the description of the spike generation of the neuronal

axon.2,3A biological neuron is a dissipative object (unit) in which

a single spike activity decays rather quickly due to high electrical

resistance. Accordingly, the higher the dissipation in the neuron,

the greater the energy required to excite the spike. The dissipa-

tive nature of spike generation in neurons was described in Refs. 7

and 8for some neuronal models. The presence of dissipation in the

FitzHugh–Nagumo neuron is an essential condition for the occur-

rence of self-sustained oscillations in the neuron, i.e., one of the

main modes of operation of the neuron. Effects related to dissipa-

tion in the FitzHugh–Nagumo neuron were considered in Refs. 9

and 10.

Time delays are a fundamental part of almost all biological phe-

nomena. The finite propagation speed of the action potential along

the axons of neurons and time lapses in information transmission

(synaptic process) and reception (dendritic process) between neu-

rons produces time delays in real neurons and their networks.11

The time required for neuronal communication can be signifi-

cantly prolonged due to the physical distance between sending and

receiving units,12,13 a finite velocity of signal transmission,14 mor-

phology of dendrites and axons,15,16 and information processing

time of the cell.17 Thus, time delays can play a crucial role in the

dynamics of neuronal networks and should be taken into account

in mathematical modeling and analysis.18,19 The role of delay in

signal transmission in brain circuits is also worth to be noted.20

Neglecting realistic time delays in mathematical models can lead

to discrepancies between theoretical and experimental findings and,

thus, prevents insight into relevant physiological mechanisms.

Many recent studies have been devoted to delay-induced

phenomena21–25 and effects of time delays on the synchroniza-

tion dynamics of neuronal networks.26–37 It is a well-known, and

often used, fact that time delay can destabilize a stationary point

and introduce oscillatory behavior. The case of two delay-coupled

FitzHugh–Nagumo systems has been examined in detail showing

that stable periodic oscillations may coexist with a stable steady

state.23,38–40 The bifurcation phenomenon is fully induced by the

delay τand represents a new form of oscillatory synchronization

exhibiting a period close to 2τ. It has been reported that delay-

enhanced synchronization may be essential for information trans-

mission in neuronal networks.17,41 It has also been revealed that

coupling delays present in the electrical or chemical synaptic con-

nections can influence the synchronization of neuronal firing.42,43 In

Ref. 44, the synchronization phenomenon in a time-delayed chaotic

system with unknown and uncertain parameters was studied, and

an intermittent adaptive control scheme was developed to guarantee

synchronization between neurons. It has been investigated how the

phase lag synchronization between the neurons of different brain

regions is governed by the spatiotemporal organization of the brain

by using self-sustained time-delayed chaotic oscillators.19

With the discovery of special patterns of partial synchroniza-

tion, called chimera states,45,46 in networks of nonlocally coupled

systems, a remarkable part of research was addressed to the effects

of time delays on birth, stability, and control of chimera struc-

tures in various complex networks.47–57 It has been shown that in

networks with fractal connectivity, desired spatiotemporal patterns

can be stabilized by varying the time-delayed coupling between the

nodes.53,57 It has been demonstrated that delay allows one to con-

trol amplitude chimeras51 and coherence-resonance chimeras by

adjusting delay time and feedback strength.52 The interlayer delay

in multilayer networks might act as a general organizing force to

synchronize spatiotemporal patterns between layers. Recently, it

has been shown that relay synchronization of chimeras in multi-

plex (triplex) networks can be controlled by time-delayed intra- and

interlayer coupling.54,55

In our work, we use the paradigmatic example of the

FitzHugh–Nagumo system in the form and for the parameter range

when the system displays excitable behavior. The form of the

used equations differs from the simplified and usually considered

FitzHugh–Nagumo system23,39 since it takes into account a param-

eter that is responsible for the dissipation of a neuron. The larger

this parameter, the less the dissipation. In our numerical simulation

of delay-coupled FitzHugh–Nagumo oscillators, we choose two dif-

ferent values of the dissipation parameter, one of which is close to

the boundary of bistability and the other one is close to the bor-

derline of the self-sustained oscillatory region. We first consider

two delay-coupled neurons and then extend our numerical simu-

lation to a ring network with local and nonlocal coupling topology

between the neurons. We study the interplay between dissipation in

the single unit, the delay time, the coupling parameters, and initial

conditions in forming the spatiotemporal dynamics of the networks.

The numerical simulations performed and the results obtained can

be straightforwardly generalized to the dynamics of delay-coupled

simplified FitzHugh–Nagumo models.23,39,58

II. SYSTEM UNDER STUDY

As an object of our study, we consider the FitzHugh–Nagumo

oscillator, which represents one of the simplest neuron models and

is widely used in numerical simulations. This two-variable system

is a paradigmatic model for neural excitability and is described as

follows:

εdx

dt =x−x3/3−y,

dt =γx−y+β,

(1)

where xis the fast variable (activator) and represents the volt-

age across the cell membrane and yis the slow recovery variable

(inhibitor), which corresponds to the recovery state of the resting

membrane of a neuron. All the control parameters are dimension-

less. The small parameter ε > 0 is the ratio of the activator to

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-2

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

FIG. 1. Diagram of dynamical regimes of the single neuron (1) in the (γ,β)

parameter plane for ε=0.01. Regions I–III correspond to the bistable (gray

color), self-sustained oscillatory (pink color), and excitable regions (white color),

respectively. In our paper, we consider two sets of parameters of the single neu-

ron marked by red dots, which correspond to the excitable regime (γ=0.5 and

γ=0.7 with β= −0.5).

inhibitor time scales. The parameter βdetermines the asymmetry,

and the parameter γis responsible for dissipation in the neuron. In

an analog circuit of the FitzHugh–Nagumo oscillator (1), the param-

eter γis directly proportional to a ratio of two resistances R0/R6,

where R6is the resistance of the oscillatory circuit.59 The value of

R6determines the amount of dissipation in the circuit. If the R6

value increases, the dissipation increases, and hence, the parame-

ter γdecreases and vice versa. Both parameters βand γalso define

the dynamical behavior of the neuron. Note that the form of the

FitzHugh–Nagumo equation (1) differs from the original form,2,60,61

but it is also investigated in a number of works.62–64 An analog elec-

trical FitzHugh–Nagumo neuron, which is described by the system

of equations, such as (1), was used in Ref. 59 to analyze its spiking

responses on pulse stimulation.

The FitzHugh–Nagumo neuron model (1) enables us to con-

sider a variety of different dynamical behaviors, namely, excitable,

self-sustained oscillatory, and bistable regimes. The ranges of exis-

tence of each dynamical regime are highlighted in the bifurcation

diagram in the (γ,β) parameter plane shown in Fig. 1. The bista-

bility region is labeled as I, the self-sustained oscillatory region is II,

and the excitability region is denoted as III. The solid lines in the dia-

gram in Fig. 1 correspond to the fold bifurcation of the equilibrium

points for γ.0.72 and to the Hopf bifurcation for γ&0.72. The

relevant values for fold bifurcations can be estimated analytically

approximately as follows:65

β= ± 2

3p1−γ(1−γ).

The location of the bifurcation lines in Fig. 1 is independent of the

parameter ε. Inside region III (Fig. 1), the neuron settles to the

excitable regime when there are no self-sustained oscillations with-

out an external force or perturbation. A single stimulus of a certain

intensity can excite the neuron activity, but it quickly decays and the

system returns to the equilibrium.65

In the following, in order to study how the parameter of dissi-

pation can affect the dynamics of delay-coupled FitzHugh–Nagumo

FIG. 2. Phase portraits (a) and (b) and time series (c) and (d) for the system (1)

at γ=0.5 (a) and (c) and γ=0.7 (b) and (d) with ε=0.01 and β= −0.5.

The xand ynull-isoclines are depicted as red dashed and green dotted curves,

respectively. Several phase trajectories approaching the stable node Nfor differ-

ent initial conditions are plotted in black. The arrows show the direction of motion.

The vertical red lines in the time series plots (c) and (d) show the time when the

maximum spike amplitude is achieved.

neurons, we choose two different values of γmarked by red dots in

Fig. 1. One of them, γ=0.5, is close to the boundary of bistability

(Fig. 1) and the other one, γ=0.7, is closer to the boundary of the

self-sustained oscillatory regime. We recall that the second case cor-

responds to smaller dissipation in a single neuron. Throughout the

paper, we fix the values of ε=0.01 and β= −0.5. This corresponds

to the excitable regime observed in a single FitzHugh–Nagumo

oscillator within a rather wide range of γ∈[0.2, 1.1] (see Ref. 65).

Figures 2(a) and 2(b) show phase trajectories approaching the sta-

ble node Nfrom a set of different initial conditions for the system

(1) for fixed βand εand for the two chosen values of γ. As can

be seen, in the case of smaller dissipation [Fig. 2(b)], the phase tra-

jectories starting with different initial states lie more closely and

densely to each other as compared with the case of larger dissipa-

tion [Fig. 2(a)]. Some of the trajectories shown in Fig. 2 move along

the big excitability loop in the directions indicated by the arrows.

As can be seen from the time series plotted in Figs. 2(c) and 2(d),

the time for achieving a maximum spike amplitude is smaller when

dissipation is low (tmax =1.23 for γ=0.7). As a consequence, the

trajectories converge toward the stable node Nmore slowly. In the

case of a higher level of dissipation [Fig. 2(c)], a larger time is taken

to show the maximum spike amplitude (tmax =1.55 for γ=0.5),

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-3

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

and the trajectories approach the stable node faster than compared

with a smaller dissipation.

III. TWO DELAY-COUPLED FITZHUGH–NAGUMO

NEURONS

Since time delay is inevitable due to the finite propagating

speed in the signal transmission between the neurons, we start

with considering two identical FitzHugh–Nagumo neurons with

time-delayed coupling and explore how the self-sustained oscilla-

tory activity of this system depends on the delay time, the coupling

strength, and the parameter of dissipation of the single neuron.

For bidirectional delayed coupling between the two FitzHugh–

Nagumo oscillators, we have the following system of equations:

εdx1

dt =x1−x3

1/3−y1+σ(x2(t−τ) −x1),

dy1

dt =γx1−y1+β,

εdx2

dt =x2−x3

2/3−y2+σ(x1(t−τ) −x2),

dy2

dt =γx2−y2+β,

(2)

where xiand yi,i=1, 2 are the dynamical variables defined in

Eq. (1) for a single neuron. The coupling strength between the

oscillators is determined by the parameter σ, and τ > 0 represents

the time delay in signal transmission. This system with the simpli-

fied FitzHugh–Nagumo oscillator form ˙

yi=xi+βhas already been

studied in Refs. 23,39, and 56.

A. Delay-induced oscillations

The dynamics of the system (2) is explored for the parame-

ter values indicated in Sec. I. For numerical simulation within this

work, we use the Runge–Kutta method with time step h=0.005.

The maximum Lyapunov exponent is defined here as

3(t0,T):=1

Tln kξ(t0+T)k

kξ(t0)k(3)

and represents the finite-time growth rate of a small pertur-

bation ξ(t0)introduced in the variational differential equations.

The initial conditions are given in the form of two random

values in both neurons. The initial states (history function) of

the coordinates are chosen as follows: x1,2(t∈[−τ; 0])=x0

1,2 and

y1,2(t∈[−τ; 0])=y0

1,2, where x0

1,2 ∈[1, 2] and y0

1,2 ∈[0, 1] to induce

the in-phase oscillatory regime and from the ranges x0

1∈[1, 2],

1∈[0, 1], x0

2∈[−2, −1], y0

2∈[−1, 0] to obtain the anti-phase

oscillatory regime. Without delay and for any values of the coupling

strength σ, the two neurons demonstrate oscillations during only a

very short transient time, after which they switch to the same equi-

librium point without any oscillations. When the delayed coupling

is introduced between the neurons, periodic self-sustained oscil-

lations are excited in the system starting with a certain threshold

of the coupling strength (σ≈0.2 for γ=0.5, τ&1 and σ≈0.1

for γ=0.7, τ&1). Depending on initial conditions, the neurons

demonstrate both complete in-phase synchronization and complete

anti-phase synchronization with each other.

The neuronal dynamics is illustrated in Fig. 3 for two different

values of the dissipation parameter (γ=0.5 and γ=0.7), the time

delay (τ=1 and τ=5), and for the anti-phase states. In this case,

the initial conditions are chosen in such a way that only one neu-

ron is excited (x0

1∈[1, 2], y0

1∈[0, 1], x0

2∈[−2, −1], y0

2∈[−1, 0]).

Panels (a1)–(b2) in Fig. 3 show results for γ=0.5, and (c1)–(d2)

illustrate regimes for γ=0.7. Time series for variables xi(t),i=1, 2

(black and orange curves, respectively) are plotted in panels with

index 1, and the corresponding projections of phase portraits for

both neurons are pictured in panels with index 2. It is seen that

despite the excitable regimes of the isolated neurons and the absence

of any external forces, the delayed coupling leads to the emergence

of periodic self-sustained oscillations (the maximal Lyapunov expo-

nent 3is equal to 0). This is a well-known result and has been

reviewed in Ref. 23. The temporal behavior for both neurons reflects

anti-phase spike oscillations [Figs. 3(a1)–3(d1)]. They are charac-

terized by a slow motion around the equilibrium (in the isolated

neuron) and a fast motion along a limit cycle that is illustrated

in Figs. 3(a2)–3(d2) for different values of the parameters. Simi-

lar results for two delay-coupled FitzHugh–Nagumo neurons in the

simplified form have been obtained numerically23 and analytically.56

Besides, the oscillations of the first and second neurons correspond

to the same limit cycle. For a larger value of γ=0.7 [Figs. 3(c1),

3(c2),3(d1), and 3(d2)], the size of the limit cycle extends in the

yaxis, while it remains almost unchanged along the xaxis. It is

related to the slope of the y-nullcline as γgrows. Hence, the x- and

y-nullclines intersect for larger values of the yvariable, and the oscil-

lation amplitude in this variable becomes larger (compare Figs. 3

and 2).

We now analyze how the neuronal dynamics changes if we vary

the time delay τin the system (2). Comparison of the time series

plotted in Figs. 3(a1) and 3(c1) for τ=1 and in Figs. 3(b1) and 3(d1)

for τ=5 gives evidence that the oscillation period increases when

the time delay becomes larger. As can be seen from the plots, this

result is independent of the parameter of dissipation. Our calcula-

tions show that the oscillation period Tstrongly depends on the

delay τas shown analytically in the simplified FitzHugh–Nagumo

model.23,56 Moreover, we reveal that the period is always equal to

2τfor the chosen initial conditions. This effect is rather easy to

explain.29,39 Since the initial conditions are chosen randomly, only

a single neuron is excited, while the other one is quiescent after a

short time. After the time delay, the second (quiescent) neuron gets

a signal from the first (excited) neuron and is also excited. In turn,

the first neuron switches to the quiescent regime until it receives a

signal from the second neuron via the coupling. Hence, the oscil-

lations in the system are anti-phase to each other, and the period

is strictly equal to the doubled time delay. It should be noted that

the emergence of delay-induced oscillations has a threshold charac-

ter. This means that the oscillations occur only when the coupling

strength σexceeds a certain threshold. It is reasoned by the follow-

ing fact. The intensity of a signal that a neuron receives from its

neighbor is equal to σA, where Ais the oscillation amplitude. This

signal can be considered an external periodic force. In this case, the

neuron demonstrates its activity only if the external force amplitude

is larger than a certain threshold level. It is worth noting that the

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-4

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

FIG. 3. Anti-phase oscillations for the two delay-coupled FitzHugh–Nagumo neu-

rons (2) for γ=0.5 [panels (a1), (a2), (b1), and (b2)] and γ=0.7 [panels

(c1), (c2), (d1), and (d2)] and for the time delay τ=1 [panels (a1), (a2), (c1),

and (c2)] and τ=5 [panels (b1), (b2), (d1), and (d2)]. Time series for the first

(i=1, black curves) and second (i=2, orange curves) neurons are shown in

panels (a1), (b1), (c1), and (d1). The corresponding phase portrait projections are

illustrated in panels (a2), (b2), (c2), and (d2). Initial conditions correspond to the

case when only one neuron is excited. Other parameters: σ=0.3, ε=0.01,

and β= −0.5.

dynamics shown in Fig. 3 is reminiscent of the reverberation of the

activity between the two nodes. An interesting application could be

the relation with the propagation of the signals in layered neuronal

networks explored in Ref. 66.

The next question is what changes will happen if both neurons

are initially excited? In this case, the neurons begin to oscillate in-

phase starting with a certain value of σ. Examples of the time series

xi(t),i=1, 2, are depicted in Fig. 4 for two values of dissipation

γ=0.5 [Figs. 4(a) and 4(b)] and γ=0.7 [Figs. 4(c) and 4(d)] and

for the time delay τ=1 [Figs. 4(a) and 4(c)] and τ=5 [Figs. 4(b)

and 4(d)]. It is seen that now both neurons oscillate synchronously

in-phase. In the case of smaller dissipation [Figs. 4(c) and 4(d)],

the form of spikes in the time series becomes more pronounced

and sharper for both time delay values as compared with the corre-

sponding spike sequences for larger dissipation [Figs. 4(a) and 4(b)].

Comparing the time series obtained for the same time delay in

Figs. 3(a1) and 3(c1) (anti-phase oscillations) and Figs. 4(a) and 4(b)

(in-phase oscillations) indicates that the frequency is doubled for

the anti-phase oscillations. At the same time, the projections of

phase portraits are completely equal to the previous cases depicted

in Figs. 3(a2) and 3(b2). The observed differences occur due to a dif-

ferent mechanism of delay-induced oscillations. For the first choice

FIG. 4. In-phase oscillations for the two delay-coupled FitzHugh–Nagumo neu-

rons (2) for γ=0.5 (a) and (b) and γ=0.7 (c) and (d). Time series for the first

(i=1, black curves) and second (i=2, orange curves) neurons are plotted for

τ=1 (a) and (c) and τ=5 (b) and (d). Initial conditions correspond to the case

when both neurons are initially excited. Other parameters: ε=0.01, σ=0.5,

and β= −0.5.

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-5

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

FIG. 5. Dependences of the period Ti(a) and frequency fi(b) for two neurons

(i=1, 2) on the time delay τfor anti-phase oscillations (red curves) and in-phase

oscillations (black curves) in the system (2) for γ=0.5. Functions T(τ) =τand

f(τ) =1/τ are shown by squares and T(τ) =2τand f(τ) =1/2τby circles.

Other parameters: ε=0.01, σ=0.5, and β= −0.5.

of initial conditions, the neurons are excited consistently, one after

another, with the oscillation period T=2τ. For the second case of

initial conditions, both neurons are activated simultaneously with

the period being equal to the delay. For the general discussion of

such exchange of pulses, including delayed feedbacks with different

delay times, see Ref. 29.

We now evaluate how the period and the frequency of oscilla-

tions depend on the delay for both types of initial conditions men-

tioned above. The calculated dependences are plotted in Fig. 5 for

γ=0.5 and for anti-phase (red curves) and in-phase (black curves)

oscillations. As follows from the graphs in Fig. 5(a), the period Ti

of oscillations in both neurons labeled with i=1, 2 is strictly equal

to the doubled time delay 2τfor the anti-phase oscillations and to

τfor the in-phase oscillatory regime. In both cases, the period Ti

increases monotonically and linearly as the delay time τgrows. In

turn, the frequency fipresents the function inverse of the time delay

1/(2τ) for the anti-phase oscillations and 1/τ for the in-phase oscil-

lations [Fig. 5(b)]. This is in line with the results in Refs. 23,29,

39, and 56. Moreover, our calculations show that both the period

and the frequency of oscillations are defined only by the delay and

are independent of the coupling strength σand the parameter γ.

However, there is a certain threshold with respect to σwhen the

oscillations are excited, and this value depends on the delay and the

control parameters; see analytical results in Refs. 29 and 56. Our cal-

culations show that the dissipation parameter γcan influence the

threshold for the neural activity in the system (2). Dependences of

the σthreshold value on γare plotted in Fig. 6 for fixed delay time

τ=5 and for the in-phase and anti-phase oscillations. The parame-

ter γis varied within the region of excitable dynamics (region III in

Fig. 1). It is seen that both dependences coincide, and the threshold

value gradually decreases as the parameter γincreases, i.e., as the

dissipation level in the neuron decreases.

To evaluate the region of existence of delay-induced oscilla-

tions in the system (2), we calculate the maximal Lyapunov exponent

3and plot its values in the (τ,σ) parameter plane for two different

FIG. 6. Dependence of the threshold for the oscillation excitation σth vs the dis-

sipation parameter γin the system of two delay-coupled neurons for in-phase

(red curve with dots) and anti-phase (green dashed curve with dots) oscillations

at τ=5. Other parameters: ε=0.01 and β= −0.5.

values of γ=0.5 and γ=0.7. The corresponding diagrams are

presented in Figs. 7(a) and 7(b) for anti-phase oscillations and in

Figs. 7(c) and 7(d) for in-phase oscillations.

Obviously, the values of 3must be negative when there are

no oscillations and zero or positive (only for chaos) for the oscil-

latory regime. In our case, the delay-induced oscillations are always

non-chaotic. Hence, the maximal Lyapunov exponent is always zero.

Thus, the distribution of 3in the (τ,σ ) parameter plane enables

one to distinguish the quiescent (red color, negative values of 3)

and the oscillatory (yellow color, 3=0) regimes and, thus, to define

the boundary between them. The diagrams show that there is a

threshold for the oscillation excitation with respect to both the cou-

pling strength σand the delay τ. When τ=0, neither anti-phase

nor in-phase stable oscillations can be excited even for very strong

coupling. Apparently, the shortest time delay is determined by a

minimal duration of a spike impulse in the neuron under excita-

tion by the delayed coupling. In return, the threshold with respect

to σis caused by a minimal amplitude of the influence of the neigh-

boring neuron through the delayed coupling. The smaller the delay,

the stronger the coupling strength must be, up to a certain value

of τ; after that, the threshold with respect to σdoes not change. It

should be noted that, as is clearly seen from Fig. 7 and is verified by

our calculations, the σthreshold is substantially lower for larger val-

ues of γ. Increasing γdecreases the neuron dissipation, and thus, a

lower amplitude of excitation is needed to induce oscillations in the

FitzHugh–Nagumo neuron. The diagrams of the maximal Lyapunov

exponent in Fig. 7 give evidence that the threshold for the oscillatory

regime with respect to the delay is larger for in-phase oscillations

than that for anti-phase ones, while the σthreshold is practically the

same for both cases. Note that negative values of the maximal Lya-

punov exponent in the quiescent region decrease with an increase in

the duration of the delay τ.

B. Linear stability analysis of the equilibrium

Let the unique fixed point be x∗=(x∗

1,y∗

1,x∗

2,y∗

2). The fixed

point is given by solving simultaneously two nonlinear equations for

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-6

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

FIG. 7. Distributions of the maximal Lyapunov exponent 3in the (τ,σ) parameter plane for anti-phase oscillations (upper row) and in-phase oscillations (lower row) at

γ=0.5 (a) and (c) and γ=0.7 (b) and (d). Other parameters: ε=0.01, β= −0.5.

each pair (x∗

1,y∗

1)=(x∗

2,y∗

2),

x1−x3

3−y1=0,

γx1−y1+β=0,

x2−x3

3−y2=0,

γx2−y2+β=0.

(4)

The equilibrium point in the excitable regime is given by

x∗

1=x∗

2and y∗

1=y∗

2. Linearizing system (2) around the fixed point

x∗by setting x(t)=x∗+δx(t), one obtains

δ˙x=1

εAδx(t)+σ

εBδx(t−τ), (5)

where

A=





ξ−1 0 0

γ ε −ε0 0

0 0 ξ−1

0 0 γ ε −ε





, (6)

B=





0 0 1 0

0 0 0 0

1 0 0 0

0 0 0 0





, (7)

where we have introduced ξ:=1−x∗2

1−σ=1−x∗2

2−σ < 0

since x∗

2>1. The ansatz

δx(t)=exp(λt)u, (8)

where uis an eigenvector of the Jacobian matrix A, leads to the fol-

lowing characteristic equation for the eigenvalues λ, given by setting

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-7

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

the determinant

det A

ε+σ

εBexp(−λτ) −λIN×N=0.

Expanding the matrices inside the determinant, we get



ε−λ−1

εσe−λτ 0

γ−1−λ0 0

σe−λτ 0ξ

ε−λ−1

0 0 γ−1−λ



=0. (9)

Because of the symmetry between 1 and 2, the determinant can be

factorized in analogy with the simplified model23 into

ξ

ε−λ(−1−λ) +γ

ε2

−σ

εe−λτ (−1−λ)2

=0, (10)

which gives the characteristic equation

λ2−λξ

ε−1−ξ

ε+γ

ε±(λ +1)σ

εe−λτ =0. (11)

This transcendental equation has infinitely many complex solutions

λ.Figure 8 shows the largest real part of λvs τ. Obtaining ξfrom

the above equation does not admit a simple form as with Ref. 23.

For all values of τ, the real part is negative and the modulus of the

real part of the eigenvalue decreases monotonically as τincreases.

Asymptotically for large τ, the eigenvalues are given by the quadratic

characteristic equation

λ2+λ|ξ|

ε+1+|ξ| + γ

ε=0,

where we have used that ξ < 0. Its solutions

λ1,2 = − 1

2|ξ|

ε+1 1∓s1−4ε|ξ| + γ

(|ξ| + ε)2!(12)

have always negative real parts. Hence, the equilibrium is always sta-

ble. For both γ=0.5 [in Fig. 8(a)] and γ=0.7 [in Fig. 8(b)], we

FIG. 8. Real part of eigenvalues λvs the delay time τfor γ=0.5 (a) and

γ=0.7 (b). The equilibrium point of the system (4) for two coupled

FitzHugh–Nagumo units is stable for all parameter values shown. The curve shifts

slightly upward with increasing σ(see blow-up in the insets). The parameters are

set as ε=0.01, β= −0.5.

FIG. 9. Phase projections in the (x1,x2) coordinate plane for the system (2) at

γ=0.5 (a) and γ=0.7 (b) with ε=0.01 and β= −0.5. In the planes, the

stable equilibrium E, the in-phase (dotted red line), and the anti-phase (green

curve) periodic orbits are shown.

plot Re(λ) vs τ. Increasing γdoes not seem to affect the stability of

the equilibrium point with respect to τ. Thus, the equilibrium point

coexists with limit cycle oscillations, and there is multistability in the

two coupled FitzHugh–Nagumo neurons. This property is verified

numerically and illustrated in Fig. 9 for two values of the parameter

γ. It is seen that there are three stable regimes in the phase plane,

i.e., the equilibrium point Nand two limit cycles, one correspond-

ing to in-phase oscillations (dotted red line) and the other one to

anti-phase oscillations (green curve).

IV. DYNAMICS OF A RING OF FITZHUGH–NAGUMO

NEURONS WITH TIME-DELAYED COUPLING

Knowing the features of the emergence of oscillations in the

two FitzHugh–Nagumo neurons with delayed coupling, we extend

our studies to a ring of excitable FitzHugh–Nagumo neurons with

delayed coupling. This network is described by the following system

of equations:

εdxi(t)

dt =xi(t)−x3

i(t)/3−yi(t)+σ

i+P

j=i−P

[xj(t−τ) −xi(t)],

(13)

dyi(t)

dt =γxi(t)−yi(t)+β,i=1, ...,N,

xi+N(t)=xi(t),yi+N(t)=yi(t),

where the index iof the dynamical variables xiand yi,i=1, ...,N

determines the position in the network and N=50 is the total num-

ber of elements. The parameter Pis the coupling range and denotes

the number of neighboring nodes where the ith oscillator is coupled

with from each side. Here, we consider the case when all the oscil-

lators are identical, and each of them interacts (through the variable

x) with one node from the left and one from the right. Thus, P=1

corresponds to local coupling between the nodes. The boundary

conditions are periodic, and the initial states are fixed for each neu-

ron as follows: xi(t∈[−τ; 0])=x0

iand yi(t∈[−τ; 0])=y0

i, where

iand y0

iare randomly uniformly distributed within the intervals:

xi(0)∈[−2, 2], yi(0)∈[−1, 1]. The other parameters have the same

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-8

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

meaning as for the system (2). Delay-coupled networks of simpli-

fied FitzHugh–Nagumo systems have been studied numerically and

analytically for a ring and for more general topologies in Refs. 30,

56, and 67.

A. Delay-induced oscillations

We vary the coupling strength σand the delay time τin the

coupling term as in the case of the two coupled FitzHugh–Nagumo

neurons described above. Interestingly, the delay-induced oscilla-

tions in the ring equation (13) can only be in-phase synchronized

for any initial conditions under consideration, and stable anti-phase

oscillations are not observed.

We start with analyzing how the threshold of the oscillation

excitation depends on the coupling strength and the delay. We also

calculate distributions of the maximal Lyapunov exponent 3in

the (τ,σ) parameter plane, which are presented in Fig. 10 for two

different values of the dissipation parameter γ.

Our simulations show that all delay-induced oscillations in the

neural network (2) are regular [3=0 within the whole region of

oscillation existence (Fig. 10)]. Comparing the diagrams constructed

for the two delay-coupled FitzHugh–Nagumo neurons (Fig. 7) and

the neural network (Fig. 10) with delayed coupling gives evidence

of a similarity in the emergence of oscillations with respect to

the coupling strength and the delay in both cases. The maximal

Lyapunov exponent 3decreases with decreasing delay in the non-

oscillatory regime as in the case of two coupled neurons. However,

the threshold with respect to σis slightly larger for the network.

Thus, compared to the two coupled neurons, the region of self-

sustained oscillations decreases both with respect to the dissipation

parameter γand delay τ. Since the scenario of oscillation excitation

is very similar for γ=0.5 and γ=0.7, we continue our numerical

studies of the neural network dynamics for the case of γ=0.5.

We now fix the delay time to τ=5 and explore the evolu-

tion of the ring dynamics as the coupling strength σincreases.

When σis sufficiently small, the initial fluctuations cannot induce

oscillations in the network (13) (red and orange colors in the dis-

tribution in Fig. 10). Delay-induced oscillations occur when 3≈0

(light-yellow colors in the diagrams in Fig. 10). However, unlike the

case of two coupled neurons, only a part of the network begins to

demonstrate oscillations, while the rest remains quiescent. This is

similar to a bump state in neuroscience.68,69 The transient process

for σ=0.3 is illustrated by a space-time plot in Fig. 11(a1). It is

seen that the transient process lasts a sufficiently long time within

several dozens of periods. There are one wide and three narrow clus-

ters, which demonstrate the oscillatory dynamics at the initial stage.

However, in the course of time, the size of the clusters changes and

the wide cluster significantly narrows, and one of the small clusters

completely disappears. As a result, only three clusters are observed

asymptotically, which is shown in the space-time plot in Fig. 11(a2),

where the vertical time scale is expanded for better visualization.

Note that the neurons in the two clusters are excited in-phase with

each other, while their instantaneous phases do not coincide with

those for the neurons in the third cluster. Apparently, the oscilla-

tion period is not strictly equal to the time delay, T=τ, but is very

close to it, T≈τ. Furthermore, the period can slightly change in

different clusters independently from each other. This leads to an

advance or delay of the instantaneous phases over a long observation

time. The rest of the ring oscillators remain quiescent. Moreover,

the oscillators of the clusters with oscillating dynamics demonstrate

their activity during a sufficiently short time and are at rest most of

the time between pulses caused by the delay, which is characteristic

for the FitzHugh–Nagumo oscillator dynamics. The instantaneous

spatial profile at the moment of oscillation excitation is shown in

Fig. 11(a3). The neurons are not excited simultaneously but with

a certain time lag. The neurons belonging to the quiescent cluster

and interacting with the edge elements of the oscillatory clusters are

also perturbed during the pulse excitation, but this perturbation is

still insufficient to excite the neurons. Apparently, increasing the

coupling strength σcan lead to the excitation of other neurons.

FIG. 10. Distributions of the maximal Lyapunov exponent 3in the (τ,σ) parameter plane for two different values of the dissipation parameter γ=0.5 (a) and γ=0.7

(b) for random initial conditions in the network (13). Other parameters: ε=0.01, β= −0.5, P=1, N=50.

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-9

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

FIG. 11. Delay-induced oscillations observed within the region with 3≈0 in the diagram in Fig. 10(a) for three increasing values of the coupling strength: σ=0.3

(a1)–(a3), σ=0.4 (b1)–(b3), and σ=0.5 (c1)–(c3). The evolution of the network (13) dynamics is illustrated by space-time plots for the transient process within t∈[0, 250]

[left column, (a1)–(c1)] and the asymptotic process [middle column, (a2)–(c2)] within t∈[2475, 2500] and by snapshots for the xivariables [right column, (a3)–(c3), red curves].

The dotted line in (c3) corresponds to the quiescent regime. Note that the time scales in the left and middle columns are different. Other parameters: τ=5, γ=0.5,

ε=0.01, β= −0.5, P=1, N=50.

We increase the coupling strength to σ=0.4. The correspond-

ing space-time plots and the snapshot of the ring dynamics are

presented in Figs. 11(b1)–11(b3). As can be seen from the plots,

increasing σcauses the change in the transient time of the process.

The initial spatial distribution of the clusters is very similar to the

case of σ=0.3 [Fig. 11(a1)]. However, the sizes of the oscillatory

clusters do not practically change in time, and the transient process

becomes sufficiently short. At the same time, the oscillation period

in different clusters continues to change slightly around the delay

time τ, and the neurons are not activated simultaneously. The stable

regime is exemplified by the space-time plot pictured in Fig. 11(b2).

The width of the oscillatory cluster remains the same as at the ini-

tial stage of oscillations. The instantaneous phases of oscillations in

different clusters are shifted with respect to each other, and their

phase differences are not constant in time, as in the case of σ=0.3.

The snapshot of the system state taken for the stable regime at the

moment of oscillation excitation [Fig. 11(b3)] gives evidence that

the neurons are also excited with a certain lag as in the previous

case.

We continue to increase the coupling strength between the

FitzHugh–Nagumo neurons in (13) and now consider the case of

σ=0.5. As follows from the space-time plot for the transient pro-

cess presented in Fig. 11(c1), the oscillations originate in the same

way as in the two previously considered cases due to the same set

of initial conditions. However, all the oscillatory clusters expand

over time. This effect is related to the fact that the edge neurons

of the oscillatory clusters interact with the quiescent neurons. Since

the coupling strength is sufficiently large, this interaction excites the

quiescent neurons, and thus, they join with the oscillatory clusters.

In turn, these neurons subsequently excite the adjacent elements.

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-10

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

Thus, all the neurons of the network (13) begin to demonstrate the

oscillatory dynamics, which is observed over a sufficiently long time.

This asymptotically stable regime is illustrated by the space-time plot

shown in Fig. 11(c2). The neurons are not excited simultaneously,

and there is a short time lag at the initial state of excitation. The oscil-

lation period is strictly equal to T=τ. The instantaneous phases of

all the nodes are close to each other but not equal, which shows a

snapshot of the system state in Fig. 11(c3). The instantaneous pro-

file shows a certain spatial distribution of the xivalues at the moment

of neuron excitation.

Our numerical simulations show that the duration of the tran-

sient process can depend on the number of initially excited neurons.

This fact is illustrated in Fig. 12 where the space-time plots are

depicted for the cases when one, two, random (about 30 nodes),

and all the network neurons (N=50) are initially excited. The

numerical results testify that there is a critical number of initially

excited neighboring neurons, which allows the ring to demonstrate

the oscillatory behavior. All the network neurons remain in the qui-

escent state if only a single neuron is initially excited [Fig. 12(a)]

and show a gradual propagation of oscillation activity when at least

two neighboring neurons are initially excited [Fig. 12(b)]. How-

ever, in this case, the transient process is rather long, as is seen in

Fig. 12(b), while all the neurons are excited more rapidly when ran-

domly chosen neighboring neurons or all the neurons are initially

excited [Figs. 12(c) and 12(d)].

In order to demonstrate how the coupling strength σaffects

the number of excited (firing) neurons, we calculate distributions of

the ratio of the number of firing nodes Nfto the whole number of

neurons in the network N[Figs. 13(a) and 13(d)], the global order

parameter

r=*

ei2i+t

where 2i=arctan(yi/xi)denotes the geometric phase of the ith

unit,49 in the (τ,σ ) parameter plane [Figs. 13(b) and 13(e)], and

spatial distributions of minimal values of xivariables in the (i,σ)

plane [Figs. 13(c) and 13(f)]. The notation h...itmeans averaging

over time. The calculated dependences are plotted in Fig. 13 for

two values of the dissipation parameter γ=0.5 [Figs. 13(a)–13(c)]

and γ=0.7 [Figs. 13(d)–13(f)]. The distributions Nf/N(τ,σ ) are

obtained using ten different sets of random initial conditions for

averaging the results.

As follows from the diagrams shown in Figs. 13(a) and 13(d)

and Figs. 13(b) and 13(e), two different threshold values with respect

to the coupling strength σcan be distinguished. Only a small part of

the neurons are excited when the first threshold, σth, min, is exceeded.

Below this value, all the neurons are quiescent and the order param-

eter ris very close to 1 [Figs. 13(b) and 13(e)]. As is seen from

the diagrams, the first threshold value is rather large for the case

of zero and very small delay times (τ < 2), rapidly decreases up to

a certain level as the delay time increases (τ > 4) and then remains

unchanged and independent of the delay time. Note that this thresh-

old value is higher for larger dissipation [it is about 0.21 for γ=0.5,

Figs. 13(a) and 13(b)] and lower for smaller dissipation [it is about

0.1 for γ=0.7, Figs. 13(d) and 13(e)].

When σexceeds the first threshold value σth, min, the number of

firing neurons Nfbegins to increase. This process can be divided into

FIG. 12. Delay-induced oscillations observed within the region with 3≈0 in the diagram in Fig. 10(a) for one (a), two (b), random (about 30 nodes) (c), and 50 (d) initially

excited neurons with σ=0.5. The evolution of the network (13) dynamics is illustrated by space-time plots. Other parameters: τ=5, γ=0.5, ε=0.01, β= −0.5,

P=1, and N=50.

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-11

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

FIG. 13. Diagrams for the ratio of the number of firing nodes Nfto the whole number of neurons (a) and (d) and for the global order parameter r(b) and (e) in the (τ,σ)

parameter plane and spatial distributions of minimal values of variable x(quiescent: red, oscillatory: yellow) (c) and (f) for the delay time τ=5 in the network (13) for two

different values of the dissipation parameter γ=0.5 (a)–(c) and γ=0.7 (d)–(f). Other parameters: ε=0.01, β= −0.5, P=1, and N=50.

two stages: first, the number Nfincreases slightly (first stage) and

then very abruptly (second stage). The first stage takes part within

a certain finite range of the coupling strength σand is related to

partial synchronization [light green color in Figs. 13(b) and 13(e)].

This range is short in the case of τ=1, while it expands for a larger

delay [Figs. 13(a),13(b),13(d), and 13(e)]. Moreover, for τ > 2, this

σ-range is twice wider for larger dissipation [Figs. 13(a) and 13(b)]

than in the case of smaller dissipation [Figs. 13(d) and 13(e)]. A sin-

gle or several clusters of neurons with the oscillating synchronous

dynamics are formed within this σ-range. The process of cluster

formation is well illustrated by the spatial distributions of mini-

mal values of the xcoordinate shown in Figs. 13(c) and 13(f) for

fixed τ=5. It is seen that the synchronous clusters are distributed

randomly along the network space and have varying widths. Our cal-

culations have indicated (not shown in this paper) that this property

depends on random initial condition distributions.

The σ-range corresponding to the first stage is bounded by the

second threshold value, σth, all, where all the neurons are excited and

Nf=N[full synchronization in Figs. 13(b) and 13(e)]. Note that

this transition occurs abruptly, and the second threshold is inde-

pendent of τstarting with τ > 1 [Figs. 13(a) and 13(d)]. Similar as

for the first threshold value described above, the value of σth, all is

larger for larger dissipation [it is about 0.48 for τ > 1 and γ=0.5,

Figs. 13(a) and 13(b)] and smaller for smaller dissipation [it is about

0.19 for τ > 1 and γ=0.7, Figs. 13(d) and 13(e)]. Panels (c) and

(f) in Fig. 13 show that this two-stage process describes a first-order

nonequilibrium phase transition from the quiescent state to full syn-

chronization via partial cluster synchronization with increasing σ

for fixed τ=5. This is similar to nucleation phenomena found

in equilibrium and nonequilibrium systems, which have recently

become a focus of research on synchronization of networks.70 The

final step is an abrupt single-step transition to full synchrony.

In order to get more insight into the effect of dissipation on

the dynamics of the FitzHugh–Nagumo neuron network, we cal-

culate dependences of the two threshold values σth,min and σth,all vs

the dissipation parameter γwhen γis varied within the region of

excitable dynamics (region III in Fig. 1). The numerical results are

plotted in Fig. 14. It is seen that increasing γ(decreasing the dissipa-

tion) leads to a gradual and nonlinear decrease in the thresholds of

oscillation excitation. Besides, these values are visibly different for

0.4 < γ < 0.9, but when the dissipation parameter becomes very

large, γ > 0.9 (very low dissipation), and approaches the boundary

of the self-sustained oscillatory regime, the threshold values almost

coincide and vanish.

B. Linear stability analysis of the equilibrium

The variational equation for a ring network of delay-coupled

FitzHugh–Nagumo systems is given by

δ˙x=1

εAδx(t)+σ

εBδx(t−τ), (14)

where x=(x1,y1,x2,y2,...,xN,yN)and Nis the size of the ring

network considered. Here, Nis set to be 50.

The analysis is similar to the one performed for the simpli-

fied FitzHugh–Nagumo model.30 For the network (14), we get the

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-12

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

FIG. 14. Dependences of the threshold values σth,min (red curve with dots) and

σth,all (green dashed curve with dots) vs the dissipation parameter γat τ=5

for randomly distributed initial conditions in the network (13). Other parameters:

ε=0.01, β= −0.5, P=1, N=50.

FIG. 15. Real part of the eigenvalues 3vs delay time τfor γ=0.5 (a) and

γ=0.7 (b). The equilibrium point of system (14) for the ring of coupled

FitzHugh–Nagumo units is stable for all parameter values shown. The curve shifts

slightly upward with increasing σ(see blow-up in the insets). The parameters are

set as ε=0.01, β= −0.5.

(a) (b)

FIG. 16. Distributions of the maximal Lyapunov exponent 3in the (τ,σ) parameter plane for the network (13) for four different values of the coupling range P: (a) P=2,

(b) P=3, (c) P=4, and (d) P=5. Other parameters: ε=0.01, β= −0.5, γ=0.5, and N=50.

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-13

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

following expressions for the matrices Aand B:







ξ−1 0 0 ... 0 0

γ ε −ε0 0 ... 0 0

0 0 ξ−1... 0 0

0 0 γ ε −ε . . . 0 0

....... .

0... 0 0 . . . ξ −1

0... 0 0 . . . γ ε −ε







, (15)







0 0 1 0 0 ... 1 0

0 0 0 0 0 ... 0 0

1 0 0 0 1 ... 0 0

0 0 0 0 0 ... 0 0

....... .

1 0 ... 1 0 0 0

0 0 0 0 0 ... 0 0







, (16)

where ξ=1−x∗

2−σ. We apply the ansatz δx(t)=exp(λt)uto

(14), which gives the following eigenvalue problem:

A

ε+σB

2εexp(−λτ)u=λu. (17)

To obtain the values of λ, we solve the determinant



ε+σB

2εe−λτ −λIN×N

=0.

We also observe a similar trend in the Re(λ) vs τplot for the ring of

FitzHigh–Nagumo units as in Sec. III B. For all values of τ, the real

part is negative and the modulus of the real part of the eigenvalue

decreases as τincreases and then tends to a fixed negative value

for larger τ. Hence, the equilibrium point of the system is stable.

With an increase in σas shown in different colors in Fig. 15, the

curve slightly shifts upward. With an increase in γto 0.7, we see that

the shape of the curve is similar. In conclusion, the quiescent equi-

librium is always stable. Therefore, the equilibrium point coexists

FIG. 17. Delay-induced oscillations in the case of nonlocal coupling shown in asymptotic spatiotemporal diagrams for a system (13) at (a) P=2, σ=0.42; (b) P=2,

σ=0.48; (c) P=5, σ=0.46; and (d) P=5, σ=0.48. Other parameters: ε=0.01, β= −0.5, and N=50.

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-14

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

with limit cycle oscillations in regions where the oscillatory regime

is observed, and there is multistability.

V. INFLUENCE OF NONLOCAL COUPLING ON

DELAY-INDUCED OSCILLATIONS

To determine the effects of nonlocal coupling, we study the

oscillations in the oscillatory ring network (13) when each neuron

interacts with Pneighbors from the left and from the right. We vary

the coupling range within P∈[2, 5] and calculate distributions of

the maximal Lyapunov exponent 3in the (τ,σ) parameter plane

within the same intervals of τand σvariations as in Fig. 10. Numer-

ical results are shown in Fig. 16 for the fixed value of γ=0.5.

If we compare these diagrams with those for the local coupling

[Fig. 10(a)], we can see that introducing nonlocal coupling leads to a

significant increase in the threshold values for inducing oscillations

with respect to both the delay time τand the coupling strength σ.

Interestingly, there is a large difference between the threshold level

for the cases P=1 [Fig. 10(a)] and P=2 [Fig. 16(a)]. However, a

further growth of the coupling range Pleads only to a weak increase

in the thresholds with respect to σand τ.

At the same time, the dynamical regimes do not undergo any

qualitative changes in comparison with the case of local coupling

P=1 illustrated in Fig. 11. The pulse duration and spatial features

remain very similar. The spatiotemporal dynamics of these regimes

are illustrated by space-time plots for P=2 and P=5 in Figs. 17(a)

and 17(b) and 17(c) and 17(d), respectively. The left column in

Fig. 17 demonstrates the case when only a part of the ring is fir-

ing (lower values of the coupling strength), while the right column

corresponds to the case when the whole system fires. Thus, the first

case is similar to that one shown in Fig. 11(b2) and the second col-

umn looks like the regime in Fig. 11(c2) for the system with local

coupling P=1. As was shown,58 it is typical for nonlocal coupling

that multiple clusters are found for a small coupling range and single

clusters for a larger coupling range.

VI. CONCLUSIONS

Our detailed numerical analysis of the dynamics of the net-

works of delay-coupled FitzHugh–Nagumo oscillators with dis-

sipation has shown qualitative and, in some cases, quantita-

tive similarities with the dynamics of delay-coupled simplified

FitzHugh–Nagumo models.23,29,39,56 However, a number of specific

differences have been established, which are due to the dissipation

parameter γin the considered neuron model.

Our numerical simulations of the dynamics of two delay-

coupled FitzHugh–Nagumo oscillators have demonstrated that the

threshold of inducing both in-phase and anti-phase periodic oscil-

lations in the system with respect to the coupling strength σsub-

stantially depends on the dissipation parameter γwhen the delay

time τincreases. The larger γand, thus, the smaller the dissipa-

tion in the neuron, the less the σthreshold. The smaller dissipation

means that a substantially lower amplitude of excitation is needed

to induce oscillations in the FitzHugh–Nagumo neurons. Starting

with a certain value of τ, the σthreshold remains unchanged as

the delay time increases. The threshold for the oscillatory regime

with respect to the delay time τdoes not depend on the dissipation

parameter, and as our calculations have shown, it is larger for induc-

ing in-phase oscillations than that for anti-phase oscillations. The

same peculiarity in the σthreshold is observed for the ring network

of locally delay-coupled FitzHugh–Nagumo oscillators. However, in

this case, the dissipation parameter γalso has a significant effect on

the τthreshold value for the occurrence of oscillations. When the

dissipation is low (γis large), the delay-induced oscillations emerge

in the network for smaller values of τas compared with the case

of large dissipation. A sufficiently high level of dissipation prevents

neurons to switch to the self-oscillatory mode; i.e., stronger coupling

is needed to excite neurons.

We have analyzed in detail the transition of locally delay-

coupled FitzHugh–Nagumo neurons from the quiescent to the oscil-

latory state for two different values of the dissipation parameter γ

when the coupling strength σand the delay time τare varied. For

this purpose, we have calculated and constructed the distributions

of the global order parameter and the ratio of the firing neurons to

the whole number of the network nodes. Both two-parameter dia-

grams have indicated the presence of two different threshold values

with respect to σ. The first value corresponds to the case when only

a small part of the neurons are excited, i.e., synchronous clusters

are formed, and the second one—when all the neurons demon-

strate self-sustained oscillations, i.e., full synchrony is reached. It

has been established that both threshold values of σare substantially

smaller (at least twice) for weaker dissipation than for stronger dis-

sipation. The intermediate region with respect to σ(starting with

a small delay time τ≈2) where the network neurons are conse-

quently excited is three times narrower for large γas compared with

the case of smaller γ. With this, we have established a first-order

nonequilibrium phase transition from the quiescent state to the fully

synchronous state via cluster formation with increasing coupling

strength. The second, abrupt transition from clusters to full syn-

chrony is a single-step transition where multiple clusters abruptly

merge into one at a certain threshold coupling strength; it is sim-

ilar to nucleation phenomena recently observed in transitions to

synchrony in power grids71 and adaptive neuronal networks.70

We believe that the results presented in this work for networks

of delay-coupled FitzHugh–Nagumo neurons with dissipation will

contribute to the insight into synchronization transitions by general-

izing previous results obtained for the simplified FitzHugh–Nagumo

networks.

ACKNOWLEDGMENTS

The reported study was funded by the Russian Science Foun-

dation (Project No. 20-12-00119). E.S. acknowledges financial sup-

port from the German Science Foundation (DFG-Projektnummer

163436311–SFB 910, 429685422, and 440145547). A.V.B. acknowl-

edges the financial support provided by The Council for Grants of

the President of the Russian Federation (Project No. SP-774.2022.5).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-15

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11

Chaos ARTICLE pubs.aip.org/aip/cha

Author Contributions

A. V. Bukh: Investigation (equal); Methodology (equal); Soft-

ware (equal); Supervision (equal); Validation (equal); Visualiza-

tion (equal); Writing – review & editing (equal). I. A. Shepelev:

Conceptualization (equal); Funding acquisition (equal); Investiga-

tion (equal); Methodology (equal); Project administration (equal);

Supervision (equal); Writing – original draft (equal); Writing –

review & editing (equal). E. M. Elizarov: Conceptualization (equal);

Investigation (equal); Methodology (equal); Validation (equal);

Writing – original draft (equal). S. S. Muni: Formal analysis (equal);

Methodology (equal); Writing – review & editing (equal). E. Schöll:

Conceptualization (equal); Formal analysis (equal); Methodol-

ogy (equal); Project administration (equal); Supervision (equal);

Writing – review & editing (equal). G. I. Strelkova: Conceptualiza-

tion (equal); Formal analysis (equal); Funding acquisition (equal);

Methodology (equal); Project administration (equal); Supervision

(equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available

from the corresponding author upon reasonable request.

REFERENCES

1A. L. Hodgkin and A. F. Huxley, J. Physiol. 117, 500 (1952).

2R. FitzHugh, Biophys. J. 1, 445 (1961).

3J. Nagumo, S. Arimoto, and S. Yoshizawa, Proc. IRE 50, 2061 (1962).

4J. Keener and J. Sneyd, Mathematical Physiology (Springer, New York, 1998).

5“Computational cell biology,” edited by C. P. Fall, E. S. Marland, J. M. Wagner,

and J. J. Tyson (Springer New York, 2002).

6G. B. Ermentrout and D. H. Terman, Mathematical Foundations of Neuroscience

(Springer, New York, 2010).

7R. Appali, U. van Rienen, and T. Heimburg, in Advances in Planar Lipid Bilay-

ers and Liposomes, Advances in Planar Lipid Bilayers and Liposomes Vol. 16

(Elsevier, 2012), pp. 275–299.

8B. Lindner, Phys. Rev. Lett. 129, 198101 (2022).

9J. G. Freire and J. A. C. Gallas, Phys. Lett. A 375, 1097 (2011).

10Y. Yao and J. Ma, Eur. Phys. J. Plus 137, 1214 (2022).

11M. M. Asl, A. Valizadeh, and P. A. Tass, Front. Physiol. 9, 1849 (2018).

12A. Knoblauch and F. T. Sommer, Neurocomputing 52–54, 301 (2003).

13A. Knoblauch and F. T. Sommer, Neurocomputing 58–60, 185 (2004).

14J. E. Desmedt and G. Cheron, Electroencephalogr. Clin. Neurophysiol. 50, 382

(1980).

15Y. Manor, C. Koch, and I. Segev, Biophys. J. 60, 1424 (1991).

16S. Boudkkazi, E. Carlier, N. Ankri, O. Caillard, P. Giraud, L. Fronzaroli-

Molinieres, and D. Debanne, Neuron 56, 1048 (2007).

17Q. Wang, M. Perc, Z. Duan, and G. Chen, Phys. Rev. E 80, 026206 (2009).

18G. Stepan, Philos. Trans. R. Soc. A 367, 1059 (2009).

19S. Petkoski and V. K. Jirsa, Philos. Trans. R. Soc. A 377, 20180132 (2019).

20A. Pariz, I. Fischer, A. Valizadeh, and C. Mirasso, PLoS Comput. Biol. 17,

e1008129 (2021).

21A. Balanov, N. Janson, and E. Schöll, Physica D 199, 1 (2004).

22A. G. Balanov, V. Beato, N. B. Janson, H. Engel, and E. Schöll, Phys. Rev. E 74,

1539–3755 (2006).

23E. Schöll, G. Hiller, P. Hövel, and M. A. Dahlem, Philos. Trans. R. Soc. A 367,

1079 (2009).

24O. V. Popovych, S. Yanchuk, and P. A. Tass, Phys. Rev. Lett. 107, 228102

(2011).

25M. Kantner, E. Schöll, and S. Yanchuk, Sci. Rep. 5, 8522 (2015).

26C.-U. Choe, T. Dahms, P. Hövel, and E. Schöll, Phys. Rev. E 81, 025205 (2010).

27Y. N. Kyrychko, K. B. Blyuss, and E. Schöll, Eur. Phys. J. B 84, 307 (2011).

28J. Lehnert, T. Dahms, P. Hövel, and E. Schöll, EPL 96, 60013 (2011).

29A. Panchuk, D. P. Rosin, P. Hövel, and E. Schöll, Int. J. Bifurcat. Chaos 23,

1330039 (2013).

30S. A. Plotnikov, J. Lehnert, A. L. Fradkov, and E. Schöll, Phys. Rev. E 94, 012203

(2016).

31C. Wille, J. Lehnert, and E. Schöll, Phys. Rev. E 90, 032908 (2014).

32Z. G. Esfahani and A. Valizadeh, PLoS One 9, e112688 (2014).

33Y. N. Kyrychko, K. B. Blyuss, and E. Schöll, Chaos 24, 043117 (2014).

34A. Gjurchinovski, A. Zakharova, and E. Schöll, Phys. Rev. E 89, 032915 (2014).

35Z. G. Esfahani, L. L. Gollo, and A. Valizadeh, Sci. Rep. 6, 23471 (2016).

36A. Pariz, Z. G. Esfahani, S. S. Parsi, A. Valizadeh, S. Canals, and C. R. Mirasso,

NeuroImage 166, 349 (2018).

37A. Ziaeemehr, M. Zarei, A. Valizadeh, and C. R. Mirasso, Neural Netw. 132, 155

(2020).

38N. Buri´

c and D. Todorovi´

c, Phys. Rev. E 67, 066222 (2003).

39M. A. Dahlem, G. Hiller, A. Panchuk, and E. Schöll, Int. J. Bifurc. Chaos 19, 745

(2009).

40O. Vallès-Codina, R. Möbius, S. Rüdiger, and L. Schimansky-Geier, Phys. Rev. E

83, 036209 (2011).

41J. Tang, J. Ma, M. Yi, H. Xia, and X. Yang, Phys. Rev. E 83, 046207 (2011).

42M. Dhamala, V. K. Jirsa, and M. Ding, Phys. Rev. Lett. 92, 074104 (2004).

43X. Yang, H. Li, and Z. Sun, PLoS One 12, e0177918 (2017).

44Y. Wang, X. Zhang, L. Yang, and H. Huang, Int. J. Control Autom. Syst. 18, 696

(2020).

45Y. Kuramoto and D. Battogtokh, Nonlinear Phenom. Complex Syst. 5, 380

(2002).

46D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. 93, 174102 (2004).

47L. Larger, B. Penkovsky, and Y. Maistrenko, Phys. Rev. Lett. 111, 054103 (2013).

48V. Semenov, A. Zakharova, Y. Maistrenko, and E. Schöll, EPL 115, 10005

(2016).

49E. Schöll, Eur. Phys. J. Spec. Top. 225, 891 (2016).

50S. Ghosh, A. Kumar, A. Zakharova, and S. Jalan, EPL 115, 60005 (2016).

51A. Gjurchinovski, E. Schöll, and A. Zakharova, Phys. Rev. E 95, 042218 (2017).

52A. Zakharova, N. Semenova, V. Anishchenko, and E. Schöll, Chaos 27, 114320

(2017).

53J. Sawicki, I. Omelchenko, A. Zakharova, and E. Schöll, Eur. Phys. J. Spec. Top.

226, 1883 (2017).

54J. Sawicki, I. Omelchenko, A. Zakharova, and E. Schöll, Phys. Rev. E 98, 062224

(2018).

55J. Sawicki, S. Ghosh, S. Jalan, and A. Zakharova, Front. Appl. Math. Stat. 5, 19

(2019).

56D. Nikitin, I. Omelchenko, A. Zakharova, M. Avetyan, A. L. Fradkov, and

E. Schöll, Philos. Trans. R. Soc. A 377, 20180128 (2019).

57J. Sawicki, I. Omelchenko, A. Zakharova, and E. Schöll, Eur. Phys. J. B 92, 54

(2019).

58A. Zakharova, Chimera Patterns in Networks (Springer International Publish-

ing, 2020).

59A. S. Tchakoutio Nguetcho, S. Binczak, V. B. Kazantsev, S. Jacquir, and J.-M.

Bilbault, Nonlinear Dyn. 82, 1595 (2015).

60E. M. Izhikevich and R. FitzHugh, Scholarpedia 1, 1349 (2006).

61M. A. Dahlem, F. M. Schneider, and E. Schöll, J. Theor. Biol. 251, 202 (2008).

62V. I. Nekorkin, D. S. Shapin, A. S. Dmitrichev, V. B. Kazantsev, S. Binczak, and

J. M. Bilbault, Physica D 237, 2463 (2008).

63V. B. Kazantsev, Phys. Rev. E 64, 056210 (2001).

64V. I. Nekorkin, A. S. Dmitrichev, D. S. Shchapin, and V. B. Kazantsev, Mat.

Model. 17, 75 (2005).

65I. A. Shepelev, D. V. Shamshin, G. I. Strelkova, and T. E. Vadivasova, Chaos,

Solitons Fractals 104, 153 (2017).

66H. Rezaei, A. Aertsen, A. Kumar, and A. Valizadeh, PLoS Comput. Biol. 16,

e1008033 (2020).

67S. Lenhert, Biophys. J. 100, 507a (2011).

68C. R. Laing and O. Omel’chenko, Chaos 30, 043117 (2020).

69H. Schmidt and D. Avitabile, Chaos 30, 033133 (2020).

70J. Fialkowski, S. Yanchuk, I. M. Sokolov, E. Schöll, G. A. Gottwald, and

R. Berner, Phys. Rev. Lett. 130, 067402 (2023).

71L. Tumash, S. Olmi, and E. Schöll, EPL 123, 20001 (2018).

Chaos 33, 073114 (2023); doi: 10.1063/5.0147883 33, 073114-16

Published under an exclusive license by AIP Publishing

10 July 2023 06:42:11