The Impact of Higher-Order Interactions on the Synchronization of Hindmarsh-Rose Neuron Maps under Different Coupling Functions

In network analysis, links depict the connections between each pair of network nodes. However, such pairwise connections fail to consider the interactions among more agents, which may be indirectly connected. Such non-pairwise or higher-order connections can be signified by involving simplicial complexes. The higher-order connections become even more noteworthy when it comes to neuronal network synchronization, an emerging phenomenon responsible for the many biological processes in real-world phenomena. However, involving higher-order interactions may considerably increase the computational costs. To confound this issue, map-based models are more suitable since they are faster, simpler, more flexible, and computationally more optimal. Therefore, this paper addresses the impact of pairwise and non-pairwise neuronal interactions on the synchronization state of 10 coupled memristive Hindmarsh-Rose neuron maps. To this aim, electrical, inner linking, and chemical synaptic functions are considered as 2- and 3-body interactions in three homogenous and two non-homogenous cases. The results show that through chemical pairwise and non-pairwise synapses, the neurons achieve synchrony with the weakest coupling strengths.

Keywords:

Subject: Physical Sciences - Theoretical Physics

1. Introduction

The work network refers to a set of nodes or agents interacting through links, which in fact, specify the configuration of the nodes’ connection. The study of the behavior of such connected nodes becomes more exciting when they have nonlinear dynamics. In mathematical neuroscience, the dynamics of each network node are defined by a neuronal model with the purpose of studying the brain’s function. As a result, many studies have been devoted to investigating neuronal collective behaviors or events that have real-world instances [1, 2]. Among such neuronal collective behaviors, synchronization has owned a dominant importance since this emergent phenomenon [3] itself includes a variety of subcategories, each of which is responsible for a biological process, disease, or function [4-6]. Complete synchronization [7], generalized synchronization [8, 9], phase or anti-phase synchronization [10, 11], lag synchronization [12], cluster synchronization [13], and chimera [14, 15] are well-known subcategories that have been examined analytically and/or numerically in literature. For instance, the necessary conditions for synchronizing the Hindmarsh-Rose (HR) neuron model via the diffusive coupling functions are given in [7]. The synchronization of two pre- and post-synaptic HR neurons is investigated in [16]. The synchronization of memristive HR (mHR) neurons with electrical and Field couplings is explored in [17]. The necessary conditions for the synchronization of the photosensitive FitzHugh-Nagumo (FHN) neurons are analytically and numerically studied in [18]. In another study carried out in [19], the synchronization of heterogeneous FHN neurons is studied. The effect of memristors as the neuronal synaptic pathways are studied for two HR in [20] and FHN in [21] neurons as well. The synchronization of the Morris-Lecar (ML) neurons with memristive autapse as the neurons’ self-feedback was taken into account in [22]. Some recent relevant studies focused on map-based neurons since it is believed that discrete-time neurons not only are able to mimic natural neuron behavior, such as spiking and bursting, but also they are more straightforward, faster, more flexible, and of less computational cost [23]. For illustration, the synchronization of the Rulkov neuron map under electrical [24, 25], inner linking [24], chemical [26, 27], hybrid [28, 29], and memristor [30] synapses are thoroughly investigated in the literature. Another synchronization study, reported in [31], was conducted on the mHR neuron map in a two-node structure network under different coupling functions, including bidirectional electrical, chemical, inner linking, and hybrid synaptic functions. The intra- and inter-layer synchronization of mHR neurons is numerically analyzed in [32].

In the literature, it is noticeable that many studies have paid attention to the pairwise interactions among neurons, and non-pairwise interactions have been neglected. Nonetheless, such non-pairwise or higher-order interactions have been proven to be existed not only among the interconnecting neuron population [33] but also among other coupled systems, including physical ones [34-36]. To nail the limitation of graph-based networks and to involve the multi-body interactions, the simplicial complexes can be considered to define the nodal interactions [36]. In this way, especially in neuronal network analysis, the connections that imply actual neuronal connectivity can be described more insightfully [37]. Consequently, some studies depict the effect of higher-order interactions on network synchronization. For instance, the synchronization of a higher-order network with HR neurons with 2- and 3-body interactions are investigated [38]. In this study, electrical and chemical higher-order interactions, as well as pairwise electrical connections, were studied, and the necessary conditions for the neurons to achieve synchrony are given analytically and numerically. In a similar study [39], the synchronization of

β

cells subjected to the 2-node and 3-node interactions was investigated. This study considered the higher-order chemical and electrical synapses, while the 2-node connections were assumed as a hybrid synapse. The impact of considering the degree of the higher levels of multi-node interactions was the objective of the study declared in [40]. This study focused on the dynamics of the higher-order network of the Rulkov maps with pairwise electrical and non-pairwise chemical synapses. The synchronization of a higher-order network of ML neurons with geometrical couplings was investigated in [41]. Besides the neuronal network analysis, higher-order interactions were studied on phase oscillators [42] and mathematical models [37].

Inspired by the aforementioned literature, the presented paper is devoted to investigating the synchronization of a higher-order network of mHR neuron maps subjected to different pairwise and non-pairwise synoptic conditions, including electrical, inner linking, and chemical coupling functions. The rest of the paper is organized as follows: the higher-order network is described in Section 2. The necessary conditions for synchronizing the mHR neuron under the assumed coupling schemes are analytically and numerically given in Section 3. Finally, Section 4 concludes the paper and sums up the important findings of the paper.

2. Higher-Order Network Model

The addition of simplicial complexes to the network model allows for considering higher-order interactions, including multi-body interaction, among the neurons involved in the network. A simplicial complex is a set of connected nodes building a topological structure [37]. For instance, 0-simplexes, 1-simplexes, and 2-simplex are respectively known as nodes, links, and triangles. Hence,

d

-complex structures can model the

d + 1

-body interactions, which are called higher-order interactions. In general, a map-based network with all possible higher-order interactions, by considering simplicial complexes in

1, \dots, d

dimensions, can be described as:

X_{i}^{n + 1} = F (X_{i}^{n}) + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) + σ_{2} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) + \dots + σ_{d} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(d)} H^{(d)} (X_{i}^{n}, X_{j_{1}}^{n}, . ., X_{j_{d}}^{n}),

(1)

where

X

is the state vector and

F (X)

is the dynamic vector of the system network.

N

is the network size,

G^{(d)} = {[G_{i j_{1} \dots j_{d}}^{(d)}]}_{N \times N \dots \times N \times d + 1}

is the adjacency tensor whose non-zero elements show nodes

i j_{1} \dots j_{d}

together form a

d

-simplex,

H^{(d)}

is the coupling function determining the relationships among the involved nodes in a

d

-dimensional simplicial structure, and

σ_{d}

is the coupling strength of

(d + 1)

-body interactions. Note that the superscript

n

shows the number of iterations, and subscript

i

indicates the node’s index.

Taking up to 2-simplexes, Network (1) can be rewritten in a more straightforward form below

X_{i}^{n + 1} = F (X_{i}^{n}) + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) + σ_{2} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) .

(2)

Here

G_{i j_{1}}^{(1)} = 1

shows there exists a link between two nodes and

G_{i j_{1} j_{2}}^{(2)} = 1

presents nodes

i j_{1} j_{2}

together construct a triangle. Figure 1a is a schematic representation of Network (2) with global couplings for

N = 10

as well as its adjacency matrix

G^{(1)}

(Figure 1 b) and adjacency tensor

G^{(2)}

(Figure 1 c).

Letting

F (X)

describes the dynamics of the mHR neuron map, and

G

determines the all-to-all network configuration for

N = 10

, this paper studies the effect of different pairwise and non-pairwise interactions (different

H^{(1)}

and

H^{(2)}

conditions) on the network synchronization. The mHR map is a three-dimensional neuron model proposed in [43] obtained by discretizing the flow-based model presented in [44]. According to the dynamics of mHR neuron map

F (X)

can be defined as

F (X) = \{\begin{matrix} f (x, y, ϕ) = x + ϵ (y - a x^{3} + b x^{2} - m \tanh (ϕ) x) \\ g (x, y, ϕ) = y + ϵ (c - d x^{2} - y) \\ h (x, y, ϕ) = ϕ - ϵ x \end{matrix},

(3)

where

x

is the membrane potential,

y

is the resting state, and

ϕ

is the magnetic flux with the strength of

m

. Other parameters are the constants affecting the dynamics of neurons’ spiking activity. Therefore,

a = 1

b = 3

c = 1

d = 5

ϵ = 0.1,

and

m = 1.4

is selected as the fixed parameter settings.

3. Results

Using the Master stability function (MSF) analysis, this section provides the necessary conditions for synchronizing the globally coupled mHR neurons with higher-order interactions under different pairwise and non-pairwise coupling conditions. First, we consider the cases wherein all interactions are homogeneous. As a result, electrical synapses, inner linking functions, and chemical synapses are considered as the 2-body and 3-body interactions separately. Thereafter, two non-homogeneous cases are taken into account wherein electrical and inner linking functions are considered as the 2-body connections while chemical 3-body interactions are maintained the same.

Furthermore, to approve the analytical results obtained through the MSF analysis, time-averaged synchronization error, henceforth called synchronization error, is regarded as the numerical assessment. The synchronization error is calculatable according to the following formula

E = \frac{1}{n (N - 1)} \sum_{k = 1}^{n} \sum_{\begin{matrix} j = 1 \\ j \neq i \end{matrix}}^{N} ‖X_{j}^{k} - X_{i}^{k}‖,

(4)

in which

‖\dots‖

symbolizes the Euclidean norm and

X = [x, y, ϕ]

3.1. Electrical Pairwise and Electrical Non-Pairwise Interactions

In the first homogenous case, both

2

-body and

3

-body interactions are assumed to be electrical. Therefore,

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) = [x_{j_{1}}^{n} - x_{i}^{n}, 0,0]

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) = [x_{j_{1}}^{n} + x_{j_{2}}^{n} - 2 x_{i}^{n}, 0,0]

. Thus, Network (2) can be updated as

X_{i}^{n + 1} = \{\begin{matrix} x_{i}^{n + 1} = f (X_{i}^{n}) + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} [x_{j_{1}}^{n} - x_{i}^{n}] + σ_{2} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [x_{j_{1}}^{n} + x_{j_{2}}^{n} - 2 x_{i}^{n}] \\ y_{i}^{n + 1} = g (X_{i}^{n}) \\ ϕ_{i}^{n + 1} = h (X_{i}^{n}) \end{matrix} .

(5)

According to the MSF formalism, a network can achieve synchrony when the synchronization manifold is stable. In the synchronization state, all neurons follow the same temporal pattern, i.e.,

X_{1}^{n} = X_{2}^{n} = \dots = X_{s}^{n}

. This leads to

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) \equiv 0

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) \equiv 0

. As a result, the synchronization manifold obeys the following relation

X_{s}^{n + 1} = F (X_{s}^{n}) = \{\begin{matrix} x_{s}^{n + 1} = f (X_{s}^{n}) \\ y_{s}^{n + 1} = g (X_{s}^{n}) \\ ϕ_{s}^{n + 1} = h (X_{s}^{n}) \end{matrix},

(6)

which is similar to the dynamics of an uncoupled neuron (Eq. 3). To investigate the stability of the synchronization manifold, a negligible perturbation is added to the synchronous states. Thus,

δ X_{i}^{n} = X_{s}^{n} - X_{i}^{n}

and the dynamics of

δ X_{i}^{n}

can be obtained through

δ X_{i}^{n + 1} = D F (X_{s}^{n}) δ X_{i}^{n} + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} [D H^{(1)} (X_{s}^{n}, X_{s}^{n}) δ X_{i} + D H^{(1)} (X_{s}^{n}, X_{s}^{n}) δ X_{j_{1}}] + σ_{2} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [D H^{(2)} (X_{s}^{n}, X_{s}^{n}, X_{s}^{n}) δ X_{i} + D H^{(2)} (X_{s}^{n}, X_{s}^{n}, X_{s}^{n}) δ X_{j_{1}} + D H^{(2)} (X_{s}^{n}, X_{s}^{n}, X_{s}^{n}) δ X_{j_{2}}],

(7)

where

D F (X_{s}^{n})

is the Jacobian of

F (X_{i}^{n})

in the synchronization manifold

X_{s}^{n}

that can be defined as

D F (X) = [\begin{matrix} \frac{\partial f (X)}{\partial x} & \frac{\partial f (X)}{\partial y} & \frac{\partial f (X)}{\partial ϕ} \\ \frac{\partial g (X)}{\partial x} & \frac{\partial g (X)}{\partial y} & \frac{\partial g (X)}{\partial ϕ} \\ \frac{\partial h (X)}{\partial x} & \frac{\partial h (X)}{\partial y} & \frac{\partial h (X)}{\partial ϕ} \end{matrix}] = [\begin{matrix} 1 - ϵ (3 a x^{2} - 2 b x + m \tanh (ϕ)) & ϵ & ϵ m x \\ - 2 d ϵ & 1 - ϵ & 0 \\ - ϵ & 0 & 1 \end{matrix}] .

(8)

Applying the assumptions, Eq. (7) becomes

δ X_{i}^{n + 1} = \{\begin{matrix} δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} [δ x_{j_{1}}^{n} - δ x_{i}^{n}] \\ + σ_{2} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [δ x_{j_{1}}^{n} + δ x_{j_{2}}^{n} - 2 δ x_{i}^{n}] \\ δ y_{i}^{n + 1} = D g (X_{s}^{n}) δ X_{i}^{n} \\ δ z_{i}^{n + 1} = D h (X_{s}^{n}) δ X_{i}^{n} \end{matrix} .

(9)

Letting

L^{(d)}

be the Laplacian matrix of

G^{(d)}

, then

L^{(d)} = D^{(d)} - G^{(d)}

, where

D^{(d)}

is the degree tensor whose elements are non-zero only on the main diagonal.

L^{(d)}

can be generally defined as

L^{(d)} = \{\begin{matrix} 0 & f o r i \neq j a n d G_{i j_{1}}^{(1)} = 0 \\ - (d - 1)! k_{i j_{1}}^{(d)} & f o r i \neq j a n d G_{i j_{1}}^{(1)} = 1 \\ d! k_{i}^{(d)} & f o r i = j \end{matrix},

(10)

where

k_{i j_{1}}^{(d)} = \frac{1}{(d - 1)!} \sum_{j_{2 = 1}}^{N} \dots \sum_{j_{d = 1}}^{N} G_{i j_{1} j_{2} \dots j_{N}}^{(d)}

Since the coupling function is only applied to the membrane potential,

δ x_{i}^{n + 1}

can be extended as

δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} + σ_{1} (\sum_{j_{1} = 1}^{N} D_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} - \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} - δ x_{i}^{n} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)}) + σ_{2} (\sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} D_{i j_{1} j_{2}}^{(2)} [δ x_{j_{1}}^{n} + δ x_{j_{2}}^{n}] - \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} L_{i j_{1} j_{2}}^{(2)} [δ x_{j_{1}}^{n} + δ x_{j_{2}}^{n}] - 2 δ x_{i}^{n} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)}) = D f (X_{s}^{n}) δ X_{i}^{n} - σ_{1} \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} - σ_{2} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} L_{i j_{1} j_{2}}^{(2)} [δ x_{j_{1}}^{n} + δ x_{j_{2}}^{n}] .

(11)

Since

\sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(2)} δ x_{j_{1}}^{n} = \sum_{j_{2} = 1}^{N} L_{i j_{2}}^{(2)} δ x_{j_{2}}^{n}

, we have

δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} - σ_{1} \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} - 2 σ_{2} \sum_{j_{1} = 1}^{N} L_{i j_{1} j_{2}}^{(2)} δ x_{j_{1}}^{n} .

(12)

In an all-to-all network configuration,

L^{(2)} = (N - 2) L^{(1)}

. Therefore,

δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} - (σ_{1} + 2 σ_{2} (N - 2)) \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} .

(13)

Consequently, Eq. (9) can be updated as

δ X_{i}^{n + 1} = \{\begin{matrix} δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} - (σ_{1} + 2 σ_{2} (N - 2)) \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} \\ δ y_{i}^{n + 1} = D g (X_{s}^{n}) δ X_{i}^{n} \\ δ z_{i}^{n + 1} = D h (X_{s}^{n}) δ X_{i}^{n} \end{matrix} .

(14)

Note that

D f (X_{s}^{n}) δ X_{i}^{n}

is block diagonal and

L^{(1)}

is diagonalizable. Considering

λ_{i}

, where

λ_{1} = 0, λ_{2} = \dots = λ_{N} = N

is the eigenvalues of

L^{(1)}

, and new variables

ζ

, the perturbation equations (Eq. (14)) can be projected to the linearized system below

ζ^{n + 1} = \{\begin{matrix} ζ_{x}^{n + 1} = D f (X_{s}^{n}) ζ^{n} - N (σ_{1} + 2 σ_{2} (N - 2)) ζ_{x}^{n} \\ ζ_{y}^{n + 1} = D g (X_{s}^{n}) ζ^{n} \\ ζ_{ϕ}^{n + 1} = D h (X_{s}^{n}) ζ^{n} \end{matrix} .

(15)

For a synchronization manifold to be stable, System (15) must be stable around the origin. According to the Lyapunov analysis, the non-positive values of the maximum Lyapunov exponent (

Λ

) show the synchronization manifold’s stability.

Figure 2a shows the values of

Λ

obtained for System (15) as a function of

0 \leq σ_{1} \leq 0.01

and

0 \leq σ_{2} \leq 0.0007

. The regions coded with purple spectra are the stability region for which

Λ \leq 0

. In the numeric approach, which is demonstrated in Figure 2b in the parameter plane

σ_{1}

σ_{2}

, the stability region coded in dark blue color with

E = 0

is the same as in Figure 2a. Overall, the synchronous and asynchronous regions can be distinguished by a linear line such that the more

σ_{2}

increases, the less

σ_{1}

is needed to synchronize the neurons and vice versa. Moreover, the results of pure 1-simplex (

σ_{2} = 0

) and pure 2-simplex (

σ_{1} = 0

) cases are presented in Figure 1c,d. According to Figure 1c,d, the synchronization is acquired for

σ_{1} \geq 0.0072

and

σ_{2} \geq 0.000455

. It can be seen that the neurons achieve synchrony in weaker strength of

σ_{2}

(higher-order case), compared to the

σ_{1}

(pairwise case) value needed to synchronize the neurons.

Figure 3 shows the neuron dynamics in the synchronization state (System (6)) using the phase diagram and time series. It should be noted that the neurons’ initial conditions are selected randomly around the origin.

3.2. Inner Linking Pairwise and Inner Linking Non-Pairwise Interactions

According to [24], an inner linking function is a more general and nonlinear form of the electrical synapse. Therefore, here another homogenous case is taken into account in which we have

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) = [f (X_{j_{1}}^{n}) - f (X_{i}^{n}), 0,0]

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) = [f (X_{j_{1}}^{n}) + f (X_{j_{2}}^{n}) - 2 f (X_{i}^{n}), 0,0]

. Thus, the network can be described as

X_{i}^{n + 1} = \{\begin{matrix} x_{i}^{n + 1} = f (X_{i}^{n}) + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} [f (X_{j_{1}}^{n}) - f (X_{i}^{n})] \\ + σ_{2} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [f (X_{j_{1}}^{n}) + f (X_{j_{2}}^{n}) - 2 f (X_{i}^{n})] \\ y_{i}^{n + 1} = g (X_{i}^{n}) \\ ϕ_{i}^{n + 1} = h (X_{i}^{n}) \end{matrix} .

(16)

When all neurons evolve synchronously,

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) \equiv 0

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) \equiv 0

. Therefore, the synchronization manifold is the same as in System (6) and demonstrated in Figure 3. To obtain the perturbation equations, similar to the previous case, a small perturbation is added to the synchronous neurons’ state and

δ X_{i}^{n} = X_{s}^{n} - X_{i}^{n}

. Using Eq. (7) and considering

L^{(d)} = D^{(d)} - G^{(d)}

, the perturbation system can be obtained through

δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} + σ_{1} D f (X_{s}^{n}) \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} [δ X_{j_{1}}^{n} - δ X_{i}^{n}] + σ_{2} D f (X_{s}^{n}) \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [δ X_{j_{1}}^{n} + δ X_{j_{2}}^{n} - 2 δ X_{i}^{n}] = D f (X_{s}^{n}) δ X_{i}^{n} + σ_{1} D f (X_{s}^{n}) (\sum_{j_{1} = 1}^{N} D_{i j_{1}}^{(1)} δ X_{j_{1}}^{n} - \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ X_{j_{1}}^{n} - δ X_{i}^{n} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)}) + σ_{2} D f (X_{s}^{n}) (\sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} D_{i j_{1} j_{2}}^{(2)} [δ X_{j_{1}}^{n} + δ X_{j_{2}}^{n}] - \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} L_{i j_{1} j_{2}}^{(2)} [δ X_{j_{1}}^{n} + δ X_{j_{2}}^{n}] - 2 δ X_{i}^{n} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)}) = D f (X_{s}^{n}) δ X_{i}^{n} - σ_{1} D f (X_{s}^{n}) \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ X_{j_{1}}^{n} - σ_{2} D f (X_{s}^{n}) \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} L_{i j_{1} j_{2}}^{(2)} [δ X_{j_{1}}^{n} + δ X_{j_{2}}^{n}] .

(17)

Considering

\sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(2)} δ x_{j_{1}}^{n} = \sum_{j_{2} = 1}^{N} L_{i j_{2}}^{(2)} δ x_{j_{2}}^{n}

and

L^{(2)} = (N - 2) L^{(1)}

for the globally coupled neurons, the perturbation system becomes

δ X_{i}^{n + 1} = \{\begin{matrix} δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} - D f (X_{s}^{n}) (σ_{1} + 2 σ_{2} (N - 2)) \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ X_{j_{1}}^{n} \\ δ y_{i}^{n + 1} = D g (X_{s}^{n}) δ X_{i}^{n} \\ δ z_{i}^{n + 1} = D h (X_{s}^{n}) δ X_{i}^{n} \end{matrix} .

(18)

Thereafter, the above-mentioned perturbation equations (Eq. (18)) can be stated in the linearized form using the new variable

ζ

ζ^{n + 1} = \{\begin{matrix} ζ_{x}^{n + 1} = D f (X_{s}^{n}) (1 - N (σ_{1} + 2 σ_{2} (N - 2))) ζ^{n} \\ ζ_{y}^{n + 1} = D g (X_{s}^{n}) ζ^{n} \\ ζ_{ϕ}^{n + 1} = D h (X_{s}^{n}) ζ^{n} \end{matrix} .

(19)

Similarly, the maximum Lyapunov exponent of System (19), shown in Figure 4a for

0 \leq σ_{1} \leq 0.012

and

0 \leq σ_{2} \leq 0.0009

can provide the necessary conditions to complete synchronization. Also, Figure 4b confirms the results obtained through the MSF analysis. It can be observed that the linear separator is maintained when the inner linking coupling functions are applied instead of the electrical synapses. However, the pure 1-simplex (Figure 4c) and 2-simplex (Figure 4d) cases better show that slightly stronger 2- and 3-body coupling strengths (

σ_{1} \geq 0.0095

and

σ_{2} \geq 0.0006

) required to synchronize the neurons. Thus, compared to Figure 2a,b, the separating line between the synchronous (in purple spectra in Figure 4a and dark blue in Figure 4b) and asynchronous regions is shifted towards the higher values of

σ_{1}

and

σ_{2}

. In the synchronous region, all neurons behave as shown in Figure 3 since the synchronization manifold remains the same as an isolated neuron’s dynamics.

3.3. Chemical Pairwise and Chemical Non-Pairwise Interactions

Electrical and inner linking functions are, in fact, more suitable to model the physical or the short-range neuronal pathway of information. Nonetheless, chemical synapses are proper to model either short- or long-range neuronal interactions [31]. Hence, in the last homogenous case, we consider

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) = [(v_{s} - x_{i}^{n}) Γ (x_{j_{1}}^{n}), 0,0]

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) = [(v_{s} - x_{i}^{n}) (Γ (x_{j_{1}}^{n}) + Γ (x_{j_{2}}^{n})), 0,0]

, where

v_{s} = - 1.4

is the reversal potential and

Γ (x) = \frac{1}{1 + e^{- k (x - θ)}}

with the slope of

k = 50

and the threshold of

θ = - 1.4

. As a consequence, the network can be expressed as

X_{i}^{n + 1} = \{\begin{matrix} x_{i}^{n + 1} = f (X_{i}^{n}) + σ_{1} (v_{s} - x_{i}^{n}) \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} Γ (x_{j_{1}}^{n}) \\ + σ_{2} (v_{s} - x_{i}^{n}) \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [Γ (x_{j_{1}}^{n}) + Γ (x_{j_{2}}^{n})] \\ y_{i}^{n + 1} = g (X_{i}^{n}) \\ ϕ_{i}^{n + 1} = h (X_{i}^{n}) \end{matrix} .

(20)

For a network with global couplings, we have

\sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} = (N - 1)

and

\sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} = (N - 1) (N - 2)

. Accordingly, in the synchronization state wherein

X_{1}^{n} = X_{2}^{n} = \dots = X_{s}^{n}

, and

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) \equiv σ_{1} (N - 1) (v_{s} - x_{s}^{n}) Γ (x_{s}^{n})

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) \equiv 2 σ_{2} (N - 1) (N - 2) (v_{s} - x_{s}^{n}) Γ (x_{s}^{n})

, the dynamics of neurons obey the following equations

X_{s}^{n + 1} = F (X_{s}^{n}) = \{\begin{matrix} x_{s}^{n + 1} = f (X_{s}^{n}) + (N - 1) (σ_{1} + 2 σ_{2} (N - 2)) (v_{s} - x_{s}^{n}) Γ (x_{s}^{n}) \\ y_{s}^{n + 1} = g (X_{s}^{n}) \\ ϕ_{s}^{n + 1} = h (X_{s}^{n}) \end{matrix} .

(21)

To analyze the stability of the synchronization state expressed in Eq. (21), the general Eq. (7) is used. Thus, letting

\sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} = (N - 1)

\sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} = (N - 1) (N - 2)

, and

L^{(2)} = (N - 2) L^{(1)}

, the purterbation eqauation

δ x_{i}^{n + 1}

reads

δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} [(v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) δ x_{j_{1}}^{n} - Γ (x_{s}^{n}) δ x_{i}^{n}] + σ_{2} D f (X_{s}^{n}) \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [(v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) δ x_{j_{1}}^{n} + (v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) δ x_{j_{2}}^{n} - 2 Γ (x_{s}^{n}) δ x_{i}^{n}] = D f (X_{s}^{n}) δ X_{i}^{n} + σ_{1} ((v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) (\sum_{j_{1} = 1}^{N} D_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} - \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n}) - Γ (x_{s}^{n}) δ x_{i}^{n} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)}) + σ_{2} ((v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) (\sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} D_{i j_{1} j_{2}}^{(2)} [δ x_{j_{1}}^{n} + δ x_{j_{2}}^{n}] - \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} L_{i j_{1} j_{2}}^{(2)} [δ x_{j_{1}}^{n} + δ x_{j_{2}}^{n}]) - 2 Γ (x_{s}^{n}) δ x_{i}^{n} \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)}) .

(22)

Finally, considering

\sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(2)} δ x_{j_{1}}^{n} = \sum_{j_{2} = 1}^{N} L_{i j_{2}}^{(2)} δ x_{j_{2}}^{n}

the perturbation equations can be obtained as

δ X_{i}^{n + 1} = \{\begin{matrix} δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} + (N - 1) (σ_{1} + 2 σ_{2} (N - 2)) \times \\ ((v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) - Γ (x_{s}^{n})) δ x_{i}^{n} \\ - (σ_{1} + 2 σ_{2} (N - 2)) (v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} \\ δ y_{i}^{n + 1} = D g (X_{s}^{n}) δ X_{i}^{n} \\ δ z_{i}^{n + 1} = D h (X_{s}^{n}) δ X_{i}^{n} \end{matrix} .

(23)

Afterward, the expression of the linearized system is

ζ^{n + 1} = \{\begin{matrix} ζ_{x}^{n + 1} = D f (X_{s}^{n}) ζ^{n} - (σ_{1} + 2 σ_{2} (N - 2)) \times \\ ((v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) + (N - 1) Γ (x_{s}^{n})) ζ_{x}^{n} \\ ζ_{y}^{n + 1} = D g (X_{s}^{n}) ζ^{n} \\ ζ_{ϕ}^{n + 1} = D h (X_{s}^{n}) ζ^{n} \end{matrix} .

(24)

The maximum Lyapunov exponent of System (24) and the synchronization error of Network (20) are demonstrated in Figure 5a and Figure 5b. From Figure 5a,b, it can be noticed that through the first- and second-order chemical interactions, the neurons synchronize for considerably weaker strength of the

σ_{1}

and

σ_{2}

. This drop is more remarkable for

σ_{1}

. Interestingly, despite the previous cases, here, several lines are needed to separate the synchronous and asynchronous zones. Although the stability region is presented in purple spectra in Figure 5a, the light purple color of the significant areas reveals that

Λ ≅ 0

. However, Figure 5b, wherein the dark blue color shows the stability region obtained in the numerical approach, manifests that such areas are stable if the initial conditions are appropriately selected. The pure 1-simplex (

σ_{2} = 0

) and 2-simplex (

σ_{1} = 0

) cases are also indicated in Figure 5c,d. Accordingly, different minor and major synchronous regions can be observed within the asynchronous zones. Nevertheless, for

σ_{1} > 0.00062

in the pure 1-simplex case (Figure 5c) and

σ_{2} > 0.00004

in the pure 2-simplex case (Figure 5d), no asynchronous areas can be observed within the synchronous region. Besides, as previously shown (System (21)), the dynamics of the neurons in the synchronization state depend on the value of the first- (

σ_{1}

) and second-order (

σ_{2}

) coupling strengths. For example, Figure 6 points out that the neurons have periodic bursting behavior of period-1 for

σ_{1} = 10

and

σ_{2} = 0.0001

and of period-2 for

σ_{1} = 10

and

σ_{2} = 0.00055

, both differ from the original chaotic dynamics of an uncoupled mHR neuron. The initial conditions are randomly chosen around the origin.

3.4. Electrical Pairwise and Chemical Non-Pairwise Interactions

As the first homogenous case, we consider the electrical synapses to model the short-range 2-body connections and chemical synapses as the long-range 3-body interactions. Hence, applying

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) = [x_{j_{1}}^{n} - x_{i}^{n}, 0,0]

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) = [(v_{s} - x_{i}^{n}) (Γ (x_{j_{1}}^{n}) + Γ (x_{j_{2}}^{n})), 0,0]

, the network dynamics can be obtained from

X_{i}^{n + 1} = \{\begin{matrix} x_{i}^{n + 1} = f (X_{i}^{n}) + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} [x_{j_{1}}^{n} - x_{i}^{n}] \\ + σ_{2} (v_{s} - x_{i}^{n}) \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [Γ (x_{j_{1}}^{n}) + Γ (x_{j_{2}}^{n})] \\ y_{i}^{n + 1} = g (X_{i}^{n}) \\ ϕ_{i}^{n + 1} = h (X_{i}^{n}) \end{matrix} .

(25)

Substituting

X_{1}^{n} = X_{2}^{n} = \dots = X_{s}^{n}

, we have and

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) \equiv 0

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) \equiv 2 σ_{2} (N - 1) (N - 2) (v_{s} - x_{s}^{n}) Γ (x_{s}^{n})

, in the synchronization state, the synchronization manifold can be acquired according to

X_{s}^{n + 1} = F (X_{s}^{n}) = \{\begin{matrix} x_{s}^{n + 1} = f (X_{s}^{n}) + 2 σ_{2} (N - 1) (N - 2) (v_{s} - x_{s}^{n}) Γ (x_{s}^{n}) \\ y_{s}^{n + 1} = g (X_{s}^{n}) \\ ϕ_{s}^{n + 1} = h (X_{s}^{n}) \end{matrix} .

(26)

System (26) shows that the behavior of the synchronous neurons depends on the value of the higher-order coupling strength (

σ_{2}

), which is here of chemical synaptic type. Looking more closely at Eq. (14) and Eq. (23), the perturbation equations needed to examine the stability of the synchronous state can be written as

δ X_{i}^{n + 1} = \{\begin{matrix} δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} + 2 σ_{2} (N - 2) (N - 1) \times \\ ((v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) - Γ (x_{s}^{n})) δ x_{i}^{n} \\ - (σ_{1} + 2 σ_{2} (N - 2) (v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n})) \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} \\ δ y_{i}^{n + 1} = D g (X_{s}^{n}) δ X_{i}^{n} \\ δ z_{i}^{n + 1} = D h (X_{s}^{n}) δ X_{i}^{n} \end{matrix} .

(27)

Consequently, the linearized system becomes

ζ^{n + 1} = \{\begin{matrix} ζ_{x}^{n + 1} = D f (X_{s}^{n}) ζ^{n} - σ_{1} N ζ_{x}^{n} \\ - 2 σ_{2} (N - 2) ((v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) + (N - 1) Γ (x_{s}^{n})) ζ_{x}^{n} \\ ζ_{y}^{n + 1} = D g (X_{s}^{n}) ζ^{n} \\ ζ_{ϕ}^{n + 1} = D h (X_{s}^{n}) \end{matrix} .

(28)

The results of the Lyapunov analysis of System (28) and the synchronization error of Network (25) are given in Figure 7a and Figure 7b. The regions coded in purple spectra (

Λ \leq 0

) in Figure 7a, and the dark blue regions (

E = 0

) in Figure 7b are the regions wherein the neurons achieve complete synchrony. Compared to homogenous cases, the edges between the synchronous and asynchronous regions are not linear. Moreover, two examples of one-dimensional cases for

σ_{2} = 0.00004

and

0 \leq σ_{1} \leq 0.01

(Figure 7c) and

σ_{1} = 0.008

and

0 \leq σ_{2} \leq 0.00005

(Figure 7d). Furthermore, As shown in Figure 8, in this case, the synchronous neurons can behave chaotically (for

σ_{1} = 0.01

and

σ_{2} = 0.00001

) and periodically (

σ_{1} = 0.001

and

σ_{2} = 0.00004

) based on the value of

σ_{2}

3.5. Inner Linking Pairwise and Chemical Non-Pairwise Interactions

In the final case of the study, the electrical function used as the pairwise neuronal interactions in the previous case is replaced by the inner linking function. Hence we get

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) = [f (X_{j_{1}}^{n}) - f (X_{i}^{n}), 0,0]

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) = [(v_{s} - x_{i}^{n}) (Γ (x_{j_{1}}^{n}) + Γ (x_{j_{2}}^{n})), 0,0]

. As a consequence, the Network (25) changes into

X_{i}^{n + 1} = \{\begin{matrix} x_{i}^{n + 1} = f (X_{i}^{n}) + σ_{1} \sum_{j_{1} = 1}^{N} G_{i j_{1}}^{(1)} [f (X_{j_{1}}^{n}) - f (X_{i}^{n})] \\ + σ_{2} (v_{s} - x_{i}^{n}) \sum_{j_{1} = 1}^{N} \sum_{j_{2} = 1}^{N} G_{i j_{1} j_{2}}^{(2)} [Γ (x_{j_{1}}^{n}) + Γ (x_{j_{2}}^{n})] \\ y_{i}^{n + 1} = g (X_{i}^{n}) \\ ϕ_{i}^{n + 1} = h (X_{i}^{n}) \end{matrix} .

(29)

Accordingly, due to the diffusive nature of the pairwise inner linking interactions, the coupling functions become

H^{(1)} (X_{i}^{n}, X_{j_{1}}^{n}) \equiv 0

and

H^{(2)} (X_{i}^{n}, X_{j_{1}}^{n}, X_{j_{2}}^{n}) \equiv 2 σ_{2} (N - 1) (N - 2) (v_{s} - x_{s}^{n}) Γ (x_{s}^{n})

, and thus, the synchronization manifold remains the same as in System (26). Inspired by Eq. (18) and Eq. (23), the stability of the synchronization manifold can be examined by performing the Lyapunov analysis on the perturbation system bellow

δ X_{i}^{n + 1} = \{\begin{matrix} δ x_{i}^{n + 1} = D f (X_{s}^{n}) δ X_{i}^{n} + 2 σ_{2} (N - 2) (N - 1) \times \\ ((v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) - Γ (x_{s}^{n})) δ x_{i}^{n} \\ - 2 σ_{2} (N - 2) (v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ x_{j_{1}}^{n} \\ - σ_{1} D f (X_{s}^{n}) \sum_{j_{1} = 1}^{N} L_{i j_{1}}^{(1)} δ X_{j_{1}}^{n} \\ δ y_{i}^{n + 1} = D g (X_{s}^{n}) δ X_{i}^{n} \\ δ z_{i}^{n + 1} = D h (X_{s}^{n}) δ X_{i}^{n} \end{matrix} .

(30)

Network (30) can then be projected to the linearized system as follows

ζ^{n + 1} = \{\begin{matrix} ζ_{x}^{n + 1} = D f (X_{s}^{n}) ζ^{n} (1 - σ_{1} N) ζ^{n} - 2 σ_{2} (N - 2) \times \\ ((v_{s} - x_{s}^{n}) Γ_{x} (x_{s}^{n}) + (N - 1) Γ (x_{s}^{n})) ζ_{x}^{n} \\ ζ_{y}^{n + 1} = D g (X_{s}^{n}) ζ^{n} \\ ζ_{ϕ}^{n + 1} = D h (X_{s}^{n}) ζ^{n} \end{matrix} .

(31)

The maximum Lyapunov exponents of System (31) are reported in Figure 9a for

0 \leq σ_{1} \leq 0.012

and

0 \leq σ_{2} \leq 0.00005

. In the same parameter intervals, the synchronization error of Network (29) is presented in Figure 9b. The purple in Figure 9a or dark blue regions in Figure 9b specifies the coupling strengths for which the neurons achieve synchrony. However, compared to Figure 7a,b, the stability region occupies a more significant area of the parameter plane

σ_{1}

σ_{2}

. Also, Figure 9c,d illustrates two one-dimensional examples for

σ_{2} = 0.000038

and

0 \leq σ_{1} \leq 0.012

(Figure 9c) and

σ_{1} = 0.0091

and

0 \leq σ_{2} \leq 0.00005

(Figure 9d). Note that, in the synchronous regions, the synchronization manifolds shown in Figure 8 can be observed of almost the same value of the coupling parameters

σ_{1}

and

σ_{2}

since the dynamics of the neurons in the synchronous state remain the same as System (26).

4. Conclusions

This paper investigated the impact of different 2- and 3-body interactions on the synchronization of mHR neuron maps configured in a higher-order network with global couplings. Therefore, a complete network of

N = 10

mHR neuron maps was considered in which the neurons were communicated through the homogenous and non-homogenous pairwise and non-pairwise coupling functions. The analysis of the stability of the synchronization state in each studied case was performed using the MSF formalism, which led to finding the necessary conditions for synchronization. Moreover, to approve the analytic results, the synchronization error of the corresponding network was calculated numerically. In homogenous cases, two- and three-neuron interactions were considered electrical, inner linking, and chemical, respectively. The results showed weaker pairwise and non-pairwise coupling strengths were needed to synchronize the mHR maps through chemical synapses. On the other hand, when neurons purely interacted through the inner linking functions, the synchronization occurred for the higher values of two-node and three-node coupling strengths. Interestingly, the synchronous and asynchronous regions were linearly separable in all homogenous cases, yet when chemical synapses were involved, multiple lines could be found between the regions. Two non-homogeneous cases were also taken into account, in both of which the three-node interactions were kept chemical since they are more suitable for long-range neuronal interactions. In the first case, two-node interactions were assumed to be electrical since they are more reasonable for short-range interactions. In the second case, the pairwise electrical synapse was replaced with the inner linking functions. The result indicated that when the inner linking function was considered to link each pair of neurons, the synchronous region occupied a significant part of the parameter plane compared to the pairwise electrical connections.

Author Contributions

Conceptualization, M.M. and A.A.; methodology, F.B. and A.H.J.; software, M.M. and A.A..; validation, S.J. and D.G.; Formal analysis, F.B. and A.H.J.; investigation, F.B.; resources, A.H.J.; visualization, A.A.; writing—original draft preparation, M.M., A.A., F.B. and A.H.J.; writing—review and editing, S.J. and D.G.; supervision, D.G.; All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article

Conflicts of Interest

The authors declare no conflict of interest.

References

Uzuntarla, M.; Barreto, E.; Torres, J.J. Inverse stochastic resonance in networks of spiking neurons. PLoS Comput. Biol. 2017, 13, e1005646. [Google Scholar] [CrossRef] [PubMed]
Uzuntarla, M.; Yilmaz, E.; Ozer, M. Vibrational resonance in a heterogeneous scale free network of neurons. Commun. Nonlinear Sci. Numer. Simul. 2015, 22, 367–374. [Google Scholar] [CrossRef]
Hong, H.; Choi, M.-Y.; Kim, B.J. Synchronization on small-world networks. Phys. Rev. E 2002, 65, 026139. [Google Scholar] [CrossRef] [PubMed]
Pecora, L.M.; Carroll, T.L. Synchronization in chaotic systems. Phys. Rev. Lett. 1990, 64, 821. [Google Scholar] [CrossRef] [PubMed]
Fell, J.; Axmacher, N. The role of phase synchronization in memory processes. Nat. Rev. Neurosci. 2011, 12, 105–118. [Google Scholar] [CrossRef] [PubMed]
Abreu, R.; da Silva, F.L.; Figueiredo, P. EEG synchronization measures predict epilepsy-related BOLD-fMRI fluctuations better than commonly used univariate metrics. Clin. Neurophysiol. 2018, 129, 618–635. [Google Scholar] [CrossRef] [PubMed]
Huang, L.; Chen, Q.; Pecora, L.M. Generic behavior of master-stability functions in coupled nonlinear dynamical systems. Phys. Rev. E 2009, 80, 036204. [Google Scholar] [CrossRef]
Acharyya, S.; Amritkar, R. Synchronization of nearly identical dynamical systems: Size instability. Phys. Rev. E 2015, 92, 052902. [Google Scholar] [CrossRef]
Koronovskii, A.; Moskalenko, O.; Hramov, A. Generalized synchronization in the action of a chaotic signal on a periodic system. Tech. Phys., 2014, 59, 629–636. [Google Scholar] [CrossRef]
Rosenblum, M.G.; Pikovsky, A.S.; Kurths, J. Phase synchronization of chaotic oscillators. Phys. Rev. Lett. 1996, 76, 1804. [Google Scholar] [CrossRef]
Yi, I.G.; Lee, H.K.; Kim, B.J. Antiphase synchronization of two nonidentical pendulums. Int. J. Bifurc. Chaos 2010, 20, 2179–2184. [Google Scholar] [CrossRef]
Shahverdiev, E.; Sivaprakasam, S.; Shore, K. Lag synchronization in time-delayed systems. Phys. Lett. A 2002, 292, 320–324. [Google Scholar] [CrossRef]
Sorrentino, F.; Pecora, L.M.; Roy, R. Complete characterization of the stability of cluster synchronization in complex dynamical networks. Sci. Adv. 2016, 2, e1501737. [Google Scholar] [CrossRef]
Parastesh, F.; et al. Chimeras. Phys. Rep. 2021, 898, 1–114. [Google Scholar] [CrossRef]
Andreev, A.; Frolov, N.; Hramov, A. Chimera state in complex networks of bistable Hodgkin-Huxley neurons. Phys. Rev. E 2019, 100, 022224. [Google Scholar] [CrossRef] [PubMed]
Malik, S.; Mir, A.H. Synchronization of hindmarsh rose neurons. Neural Netw. 2020, 123, 372–380. [Google Scholar] [CrossRef]
Usha, K.; Subha, P. Hindmarsh-Rose neuron model with memristors. Biosystems 2019, 178, 1–9. [Google Scholar] [CrossRef]
Hussain, I.; Jafari, S.; Perc, M. Synchronization and chimeras in a network of photosensitive FitzHugh–Nagumo neurons. Nonlinear Dyn. 2021, 104, 2711–2721. [Google Scholar] [CrossRef]
Plotnikov, S.A.; Fradkov, A.L. On synchronization in heterogeneous FitzHugh–Nagumo networks. Chaos Solitons Fractals 2019, 121, 85–91. [Google Scholar] [CrossRef]
Xu, Y.; Jia, Y.; Ahmad, B. Synchronization between neurons coupled by memristor. Chaos Solitons Fractals 2017, 104, 435–442. [Google Scholar] [CrossRef]
Li, R.; Wang, Z.; Dong, E. A new locally active memristive synapse-coupled neuron model. Nonlinear Dyn. 2021, 104, 4459–4475. [Google Scholar] [CrossRef]
Fan, W.; Chen, X.; Xu, Q. Firing patterns and synchronization of Morris-Lecar neuron model with memristive autapse. AEU - Int. J. Electron. Commun. 2023, 158, 154454. [Google Scholar] [CrossRef]
Ibarz, B.; Casado, J.M.; Sanjuán, M.A. Map-based models in neuronal dynamics. Phys. Rep. 2011, 501, 1–74. [Google Scholar] [CrossRef]
Sun, H.; Cao, H. Complete synchronization of coupled Rulkov neuron networks. Nonlinear Dyn. 2016, 84, 2423–2434. [Google Scholar] [CrossRef]
Wang, C.; Cao, H. Stability and chaos of Rulkov map-based neuron network with electrical synapse. Commun. Nonlinear Sci. Numer. Simul. 2015, 20, 536–545. [Google Scholar] [CrossRef]
Hu, D.; Cao, H. Stability and synchronization of coupled Rulkov map-based neurons with chemical synapses. Commun. Nonlinear Sci. Numer. Simul. 2016, 35, 105–122. [Google Scholar] [CrossRef]
Ge, P.; Cao, H. Synchronization of Rulkov neuron networks coupled by excitatory and inhibitory chemical synapses. Chaos 2019, 29, 023129. [Google Scholar] [CrossRef]
Rakshit, S.; Ray, A.; Ghosh, D. Synchronization and firing patterns of coupled Rulkov neuronal map. Nonlinear Dyn. 2018, 94, 785–805. [Google Scholar] [CrossRef]
Mehrabbeik, M.; et al. Synchronization and chimera states in the network of electrochemically coupled memristive Rulkov neuron maps. Math. Biosci. Eng. 2021, 18, 9394–9409. [Google Scholar] [CrossRef]
Li, K.; Bao, B.; Bao, H. Synchronization transitions in a discrete memristor-coupled bi-neuron model. Chaos, Solitons Fractals 2022, 165, 112861. [Google Scholar] [CrossRef]
Wang, S.; Wei, Z. Synchronization of coupled memristive Hindmarsh–Rose maps under different coupling conditions. AEU - Int. J. Electron. Commun. 2023, 161, 154561. [Google Scholar] [CrossRef]
Fan, W.; Wu, H.; Xu, Q. Synchronization and chimera in a multiplex network of Hindmarsh–Rose neuron map with flux-controlled memristor. Eur. Phys. J. Spec. Top. 2022, 231, 4131–4141. [Google Scholar] [CrossRef]
Ince, R.A.A.; Montani, F.; Arabzadeh, E.; Diamond, M.E.; Panzeri, S. On the presence of high-order interactions among somatosensory neurons and their effect on information transmission. J. Phys. Conf. Ser. 2009, 197, 012013. [Google Scholar] [CrossRef]
Böhle, T.; Kuehn, C.; Mulas, R.; Jost, J. Coupled hypergraph maps and chaotic cluster synchronization. Europhys. Lett. 2021, 136, 40005. [Google Scholar] [CrossRef]
Mulas, R.; Kuehn, C.; Jost, J. Coupled dynamics on hypergraphs: Master stability of steady states and synchronization. Phys. Rev. E 2020, 101, 062313. [Google Scholar] [CrossRef] [PubMed]
Carletti, T.; Fanelli, D.; Nicoletti, S. Dynamical systems on hypergraphs. J. Phys. Complexity 2020, 1, 035006. [Google Scholar] [CrossRef]
Gambuzza, L.V.; et al. Stability of synchronization in simplicial complexes. Nat. Commun. 2021, 12, 1255. [Google Scholar] [CrossRef]
Parastesh, F.; Mehrabbeik, M.; Perc, M. Synchronization in Hindmarsh–Rose neurons subject to higher-order interactions. Chaos 2022, 32, 013125. [Google Scholar] [CrossRef]
Anwar, M.S.; Ghosh, D. Stability of synchronization in simplicial complexes with multiple interaction layers. Phys. Rev. E 2022, 106, 034314. [Google Scholar] [CrossRef]
Mirzaei, S.; Mehrabbeik, M.; Chen, G. Synchronization of a higher-order network of Rulkov maps. Chaos 2022, 32, 123133. [Google Scholar] [CrossRef]
Tlaie, A.; Leyva, I.; Sendiña-Nadal, I. High-order couplings in geometric complex networks of neurons. Phys. Rev. E 2019, 100, 052305. [Google Scholar] [CrossRef] [PubMed]
Skardal, P.S.; Arenas, A. Higher order interactions in complex networks of phase oscillators promote abrupt synchronization switching. Commun. Phys. 2020, 3, 1–6. [Google Scholar] [CrossRef]
Bao, H.; Hua, Z.; Bao, B. Discrete memristive neuron model and its interspike interval-encoded application in image encryption. Sci. China Technol. Sci. 2021, 64, 2281–2291. [Google Scholar] [CrossRef]
Bao, H.; Hu, A.; Bao, B. Hidden bursting firings and bifurcation mechanisms in memristive neuron model with threshold electromagnetic induction. IEEE Trans. Neural Networks Learn. Syst. 2019, 31, 502–511. [Google Scholar] [CrossRef]

Figure 1. (a) A schematic representation of Network (2) with

N = 10

globally coupled nodes. Black nodes, green links, and the light yellow triangle represent 0-simplex, 1-simplex, and 2-simplex structures. Also, (b)

G^{(1)}

is the

N \times N

adjacency matrix, and (c)

G^{(2)}

is the

N \times N \times 2

adjacency tensors. Light yellow matrix elements indicate the existence of a relation among involved nodes, and green elements show otherwise. Thus,

G_{i j_{1}}^{(1)} = 1

shows nodes

i

and

j_{1}

are connected through a link, and nodes

G_{i j_{1}}^{(2)} = 1

shows nodes

i

j_{1}, and j_{2}

together construct a triangle.

Figure 1. (a) A schematic representation of Network (2) with

N = 10

globally coupled nodes. Black nodes, green links, and the light yellow triangle represent 0-simplex, 1-simplex, and 2-simplex structures. Also, (b)

G^{(1)}

is the

N \times N

adjacency matrix, and (c)

G^{(2)}

is the

N \times N \times 2

adjacency tensors. Light yellow matrix elements indicate the existence of a relation among involved nodes, and green elements show otherwise. Thus,

G_{i j_{1}}^{(1)} = 1

shows nodes

i

and

j_{1}

are connected through a link, and nodes

G_{i j_{1}}^{(2)} = 1

shows nodes

i

j_{1}, and j_{2}

together construct a triangle.

Figure 2. First row: (a) The maximum Lyapunov exponent of System (15) and (b) the synchronization error of Network (5) with

N = 10

for

0 \leq σ_{1} \leq 0.01

and

0 \leq σ_{2} \leq 0.0007

. The stability region for which

Λ \leq 0

and

E = 0

is coded in purples spectra and dark blue in the analytical and numerical approaches, respectively. Second row: The maximum Lyapunov exponent of System (15) (shown in orange) and the synchronization error of Network (5) (shown in navy blue) for the (c) pure 1-simplex (

σ_{2} = 0

) and (d) pure 2-simplex (

σ_{1} = 0

) cases.

Figure 2. First row: (a) The maximum Lyapunov exponent of System (15) and (b) the synchronization error of Network (5) with

N = 10

for

0 \leq σ_{1} \leq 0.01

and

0 \leq σ_{2} \leq 0.0007

. The stability region for which

Λ \leq 0

and

E = 0

σ_{2} = 0

) and (d) pure 2-simplex (

σ_{1} = 0

) cases.

Figure 3. (a) The phase diagram and (b) the time series of the mHR neuron maps in their synchronous state described in System (6) for

σ_{1} = 0.005

and

σ_{2} = 0.0005

. Other parameters are

a = 1

b = 3

c = 1

d = 5

ϵ = 0.1,

, and

m = 1.4

. The initial values are considered randomly around the origin.

Figure 3. (a) The phase diagram and (b) the time series of the mHR neuron maps in their synchronous state described in System (6) for

σ_{1} = 0.005

and

σ_{2} = 0.0005

. Other parameters are

a = 1

b = 3

c = 1

d = 5

ϵ = 0.1,

, and

m = 1.4

. The initial values are considered randomly around the origin.

Figure 4. First row: (a) The maximum Lyapunov exponent of System (19) and (b) the synchronization error of Network (16) with

N = 10

for

0 \leq σ_{1} \leq 0.012

and

0 \leq σ_{2} \leq 0.0009

. The stability region for which

Λ \leq 0

and

E = 0

is coded in purples spectra and dark blue in the analytical and numerical approaches, respectively. Second row: The maximum Lyapunov exponent of System (19) (shown in orange) and the synchronization error of Network (16) (shown in navy blue) for the (c) pure 1-simplex (

σ_{2} = 0

) and (d) pure 2-simplex (

σ_{1} = 0

) cases.

Figure 4. First row: (a) The maximum Lyapunov exponent of System (19) and (b) the synchronization error of Network (16) with

N = 10

for

0 \leq σ_{1} \leq 0.012

and

0 \leq σ_{2} \leq 0.0009

. The stability region for which

Λ \leq 0

and

E = 0

σ_{2} = 0

) and (d) pure 2-simplex (

σ_{1} = 0

) cases.

Figure 5. First row: (a) The maximum Lyapunov exponent of System (24) and (b) the synchronization error of Network (20) with

N = 10

for

0 \leq σ_{1} \leq 0.0007

and

0 \leq σ_{2} \leq 0.00005

. The stability region for which

Λ \leq 0

and

E = 0

is coded in purples spectra and dark blue in the analytical and numerical approaches, respectively. Second row: The maximum Lyapunov exponent of System (24) (shown in orange) and the synchronization error of Network (20) (shown in navy blue) for the (c) pure 1-simplex (

σ_{2} = 0

) and (d) pure 2-simplex (

σ_{1} = 0

) cases.

Figure 5. First row: (a) The maximum Lyapunov exponent of System (24) and (b) the synchronization error of Network (20) with

N = 10

for

0 \leq σ_{1} \leq 0.0007

and

0 \leq σ_{2} \leq 0.00005

. The stability region for which

Λ \leq 0

and

E = 0

σ_{2} = 0

) and (d) pure 2-simplex (

σ_{1} = 0

) cases.

Figure 6. (a),(c) The phase diagram and (b),(d) the time series of the mHR neuron maps in their synchronous state described in System (21) for

σ_{1} = 10 σ_{2} = 0.0001

(first row) and

σ_{1} = 10 σ_{2} = 0.00055

second row). Other parameters are,

b = 3

c = 1

d = 5

ϵ = 0.1,

and

m = 1.4 .

The initial values are considered randomly around the origin.

Figure 6. (a),(c) The phase diagram and (b),(d) the time series of the mHR neuron maps in their synchronous state described in System (21) for

σ_{1} = 10 σ_{2} = 0.0001

(first row) and

σ_{1} = 10 σ_{2} = 0.00055

second row). Other parameters are,

b = 3

c = 1

d = 5

ϵ = 0.1,

and

m = 1.4 .

The initial values are considered randomly around the origin.

Figure 7. First row: (a) The maximum Lyapunov exponent of System (28) and (b) the synchronization error of Network (25) with

N = 10

for

0 \leq σ_{1} \leq 0.01

and

0 \leq σ_{2} \leq 0.00005

. The stability region for which

Λ \leq 0

and

E = 0

is coded in purples spectra and dark blue in the analytical and numerical approaches, respectively. Second row: The maximum Lyapunov exponent of System (28) (shown in orange) and the synchronization error of Network (25) (shown in navy blue) for the (c)

0 \leq σ_{1} \leq 0.01

and

σ_{2} = 0.00004

and (d)

σ_{1} = 0.008

and

0 \leq σ_{2} \leq 0.00005

Figure 7. First row: (a) The maximum Lyapunov exponent of System (28) and (b) the synchronization error of Network (25) with

N = 10

for

0 \leq σ_{1} \leq 0.01

and

0 \leq σ_{2} \leq 0.00005

. The stability region for which

Λ \leq 0

and

E = 0

0 \leq σ_{1} \leq 0.01

and

σ_{2} = 0.00004

and (d)

σ_{1} = 0.008

and

0 \leq σ_{2} \leq 0.00005

Figure 8. (a),(c) The phase diagram and b,d) the time series of the mHR neuron maps in their synchronous state described in System (26) for

1000 σ_{1} = σ_{2} = 0.00001

(first row) and

40 σ_{1} = σ_{2} = 0.00004

(second row). Other parameters are

a = 1

b = 3

c = 1

d = 5

ϵ = 0.1,

and

m = 1.4 .

The initial values are considered randomly around the origin.

Figure 8. (a),(c) The phase diagram and b,d) the time series of the mHR neuron maps in their synchronous state described in System (26) for

1000 σ_{1} = σ_{2} = 0.00001

(first row) and

40 σ_{1} = σ_{2} = 0.00004

(second row). Other parameters are

a = 1

b = 3

c = 1

d = 5

ϵ = 0.1,

and

m = 1.4 .

The initial values are considered randomly around the origin.

Figure 9. First row: (a) The maximum Lyapunov exponent of System (31) and (b) the synchronization error of Network (29) with

N = 10

for

0 \leq σ_{1} \leq 0.012

and

0 \leq σ_{2} \leq 0.00005

. The stability region for which

Λ \leq 0

and

E = 0

is coded in purples spectra and dark blue in the analytical and numerical approaches, respectively. Second row: The maximum Lyapunov exponent of System (31) (shown in orange) and the synchronization error of Network (29) (shown in navy blue) for the (c)

0 \leq σ_{1} \leq 0.012

and

σ_{2} = 0.000038

and (d).

σ_{1} = 0.0091

and

0 \leq σ_{2} \leq 0.00005

Figure 9. First row: (a) The maximum Lyapunov exponent of System (31) and (b) the synchronization error of Network (29) with

N = 10

for

0 \leq σ_{1} \leq 0.012

and

0 \leq σ_{2} \leq 0.00005

. The stability region for which

Λ \leq 0

and

E = 0

0 \leq σ_{1} \leq 0.012

and

σ_{2} = 0.000038

and (d).

σ_{1} = 0.0091

and

0 \leq σ_{2} \leq 0.00005

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

The Impact of Higher-Order Interactions on the Synchronization of Hindmarsh-Rose Neuron Maps under Different Coupling Functions

Abstract

1. Introduction

2. Higher-Order Network Model

3. Results

3.1. Electrical Pairwise and Electrical Non-Pairwise Interactions

3.2. Inner Linking Pairwise and Inner Linking Non-Pairwise Interactions

3.3. Chemical Pairwise and Chemical Non-Pairwise Interactions

3.4. Electrical Pairwise and Chemical Non-Pairwise Interactions

3.5. Inner Linking Pairwise and Chemical Non-Pairwise Interactions

4. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe