# Stability of synchronization in simplicial complexes with multiple
interaction layers

Md Sayeed Anwar  Physics and Applied Mathematics Unit, Indian Statistical
Institute, 203 B. T. Road, Kolkata 700108, India  Dibakar Ghosh
<EMAIL_ADDRESS>Physics and Applied Mathematics Unit, Indian Statistical
Institute, 203 B. T. Road, Kolkata 700108, India

###### Abstract

Understanding how the interplay between higher-order and multilayer structures
of interconnections influences the synchronization behaviors of dynamical
systems is a feasible problem of interest, with possible application in
essential topics such as neuronal dynamics. Here, we provide a comprehensive
approach for analyzing the stability of the complete synchronization state in
simplicial complexes with numerous interaction layers. We show that the
synchronization state exists as an invariant solution and derive the necessary
condition for a stable synchronization state in presence of general coupling
functions. It generalizes the well-known master stability function scheme to
the higher-order structures with multiple interaction layers. We verify our
theoretical results by employing them on networks of paradigmatic Rössler
oscillators and Sherman neuronal models, and demonstrate that the presence of
group interactions considerably improves the synchronization phenomenon in the
multilayer framework.

##  I Introduction

Network theory Boccaletti  _et al._ (2006); Newman (2003) has proven to be an
excellent foundation for modeling a variety of natural and artificial systems.
One of the most appealing aspects of this approach is its ability to recognize
unified characteristics in the framework of connections between fundamental
elements of a system Watts and Strogatz (1998); Dorogovtsev  _et al._ (2008).
The study of dynamical systems on networks has thus piqued scientists’
curiosity, with applications ranging from engineering and physics to ecology
and social science Newman (2018); Latora  _et al._ (2017). These network
models are premised on the hypothesis that the components of a system are
connected only through pairwise links, schematized by graphs. However, this
pairwise representation is merely a first-order estimation Battiston  _et al._
(2020); Majhi  _et al._ (2022) in many real-world systems, such as structural
Sizemore  _et al._ (2018) and functional Petri  _et al._ (2014) brain
networks, ecological networks Grilli  _et al._ (2017); Billick and Case
(1994), collaboration graphs Patania  _et al._ (2017); Vasilyeva  _et al._
(2021), protein interaction networks Estrada and Ross (2018), social systems
Benson  _et al._ (2016); Carletti  _et al._ (2020), and consensus dynamics
Neuhäuser  _et al._ (2021), where many-body interactions between system
components are ubiquitous and of immense interest. To schematize these higher-
order interactions, one then requires a generalization of classical network
theory Bick  _et al._ (2021); Battiston  _et al._ (2021).

Simplicial complexes Torres and Bianconi (2020); Giusti  _et al._ (2016), a
topological structure formed by an ensemble of simplices of different
dimensions, are a natural generalization of graphs to describe these group
interactions. A $d$-simplex is a set of $(d+1)$ nodes
$\sigma=[v_{0},v_{1},\cdots,v_{d}]$, namely the $0$-simplices are the
vertices, the edges are $1$-simplices, while $2$-simplices correspond to full
triangles and so on. A simplicial complex $\mathscr{S}$ is a collection of
simplices with an extra constraint that if simplex $\sigma\in\mathscr{S}$ then
all the subsets of $\sigma$ also belong to $\mathscr{S}$, i.e., a three-body
interaction, for example, necessitates the presence of all pairwise
connections corresponding to the same triangle. A simplicial complex is called
$D$-dimensional if the maximal simplices contained in it are of dimension $D$.
Although the concept of simplicial complexes is not new Aleksandrov (1998);
Berge (1973), the recent advances in real-world data availability Carlsson
(2009) have renewed the interest of researchers in examining the collective
dynamical behavior in simplicial complexes, which shows some interesting
results due to the inclusion of many-body interactions among dynamical
entities Iacopini  _et al._ (2019); Alvarez-Rodriguez  _et al._ (2021);
Skardal and Arenas (2019); Battiston  _et al._ (2020); Petri and Barrat
(2018).

Synchronization Pikovsky  _et al._ (2001); Boccaletti  _et al._ (2002); Arenas
_et al._ (2008), or the formation of the order in ensembles of interacting
individuals, is one such fascinating emergent phenomenon where the connection
between isolated dynamical elements plays an important role. Researchers have
recently expanded the study of synchronization to the framework including
higher-order network structures, with the majority of them opting for
simplicial complexes to simulate group interactions due to their simple
topological representation Courtney and Bianconi (2016); Ghorbanchian  _et
al._ (2021). The presence of many-body interactions has been linked to the
emergence of abrupt synchronization transitions Skardal and Arenas (2019,
2020); Kuehn and Bick (2021); Kachhvah and Jalan (2022a), improvement of
synchronization Skardal  _et al._ (2021); Zhang  _et al._ (2022),
multistability Xu  _et al._ (2020), cluster synchronization Zhang  _et al._
(2021), antiphase synchronization Kachhvah and Jalan (2022b), chimeras Kundu
and Ghosh (2022), etc. These findings have coincided with the development of
analytical paradigms for interpreting coupled oscillators with many-body
interactions, such as Hodge decomposition Millán  _et al._ (2020); Arnaudon
_et al._ (2021), Laplacian operators Gambuzza  _et al._ (2021); Lucas  _et
al._ (2020), and low dimensional descriptions Gong and Pikovsky (2019).

On the other hand, multilayer systems are other generalized network structures
in which interactions of different kinds and meanings can coexist, resulting
in networks of networks Sorrentino (2012); Del Genio  _et al._ (2016). In
these systems, nodes can be linked using various types of connections, and
each of these connection types characterizes an interaction layer Boccaletti
_et al._ (2014). Mobility networks, where different modes of transportation
connect individual elements Cardillo  _et al._ (2013), social networks,
wherein individuals are interconnected and acquainted by various kinds of
relations Szell  _et al._ (2010), and neuronal networks, where the neurons
communicate through chemical and electrical channels Rakshit  _et al._ (2018),
are examples of multilayered systems. Various collective phenomena, such as
synchronization Anwar  _et al._ (2022, 2021a), percolation Bianconi and
Dorogovtsev (2014), diffusion Gomez  _et al._ (2013), epidemic spreading
Granell  _et al._ (2013), and evolutionary games Gómez-Gardenes  _et al._
(2012), have been investigated in multilayer frameworks, all revealing a
phenomenology that differs significantly from that observed in monolayer
structures. In all these studies, the connections between individual nodes of
a layer are considered to be pairwise, schematized by links. However, this
pairwise connectivity representation within layers may not always be able to
illustrate the connection topology accurately among individual nodes. The
neuronal system is an example where, on the one hand, individual neurons are
connected via electrical and chemical synapses Rakshit  _et al._ (2018);
Parastesh  _et al._ (2022), and on the other hand, proof of many-body
interactions between individual neurons have recently been found Battiston
_et al._ (2020); Ince  _et al._ (2009); Tlaie  _et al._ (2019); Amari  _et
al._ (2003). Nevertheless, the interplay between higher-order structures and
multilayer networks Sun and Bianconi (2021), under some aspects has yet to be
investigated, and specifically, the study of synchronization is still in its
early stages Anwar and Ghosh (2022); Jalan and Suman (2022). In the previous
study Anwar and Ghosh (2022) of synchronization on multilayer framework with
higher-order interactions, the interactions between individual elements of a
particular layer are considered to be of a specific functional form (linear
diffusive). These linear interactions factorize the many-body structures into
pairwise interactions due to the superposition property of linear functions,
resulting in a weighted pairwise network. Nevertheless, this weighted pairwise
network representation cannot always capture the group interactions between
individuals, especially when the interactions between unitary elements are
nonlinear Battiston  _et al._ (2020); Neuhäuser  _et al._ (2021) ( e.g., in
neuronal systems). The higher-order interactions cannot be factorized into
weighted pairwise connections in such circumstances, and therefore there is
still room for improvement in the investigation of synchronization on higher-
order multilayered systems.

In this work, we fill this gap by investigating the interplay between
multilayer and higher-order structures and providing a generic multilayer
framework to study the synchronization phenomenon on simplicial complexes. Our
goal is to develop a comprehensive mathematical approach for evaluating the
stability of the synchronization state for nonlinear dynamical systems growing
in simplicial complexes with numerous interaction layers. To do so, we first
derive the condition of invariance for the synchronization state and then
conduct the linear stability analysis of the synchronization solution, which
as a result generalizes the well-known master stability function (MSF)
approach Pecora and Carroll (1998) to the realm of simplicial complexes
Gambuzza  _et al._ (2021) with multiple interaction layers Del Genio  _et al._
(2016). The end result is a system of coupled linear equations ( master
stability equation ) for the evolution of separation of the simplicial
complex, which can be optimally separated for a particular instance. The
proposed approach does not require any restriction on the functional forms of
the coupling schemes between individual elements of the layers. We validate
all our theoretical predictions in a simplicial complex of coupled Sherman
neuronal model and chaotic Rössler oscillators with two layers of interactions
featuring different nonlinear coupling schemes and connection topologies.

The remaining of this article is structured in the following manner. A
mathematical framework for the simplicial complexes with multiple interaction
layers is proposed in Sec. II. Section III shows the theoretical results on
the stability of the synchronization state. The invariance criterion for the
synchronization state is obtained in Sec. III.1. Following that, in Sec.
III.2, we develop the analytical condition for the stable synchronous
solution. To validate our theoretical findings, in Sec. IV, we reported
numerical results on systems of coupled oscillators. Finally, in Sec. V, we
summarize all of our findings and draw a conclusion.

##  II Generalized Mathematical Model of the simplicial complex

We consider a simplicial complex of dimension $D=2$, i.e., the interactions
between dynamical units are not only limited to pairwise but also three-body
interactions are allowed. The simplicial complex is composed of $N$ dynamical
units connected through $M$ distinct interacting layers. It is worth noting
that the nodes communicating on each layer in this setup are essentially the
same nodes. The node $i$ in the layer- $1$ is identical to the node $i$ in all
the other layers. Various sorts of transport between cities or electrical and
chemical synaptic couplings between neurons are examples of complex systems
with multiple interaction layers where essentially the same nodes are
interconnected by different means of connections. This differs from other
works Rakshit  _et al._ (2018); Anwar  _et al._ (2022, 2021b); Anwar and Ghosh
(2022); Anwar  _et al._ (2021a) where nodes in various layers can be identical
or different from each other, and further, the nodes in adjacent layers are
connected through a set of links called the interlayer connections. Now, we
further consider that among all the layers, within $L$ $(1\leq L\leq M)$
number of layers a $2$-simplex $(i,j,k)$ is created by promoting an empty
triangle composed of three $1$-simplices $(i,j)$, $(j,k)$, $(k,i)$ to a full
triangle $(i,j,k)$, i.e., all the cliques of size $3$ are considered as a
$2$-simplex. Therefore, these $L$ layers consider both the two-body and three-
body interactions. However, in the other $(M-L)$ layers, no empty triangle
composed of three $1$-simplices is promoted to a full triangle ($2$\-
simplex), i.e., these $(M-L)$ layers consider only the pairwise interactions
between dynamical units. This implies that the proposed simplicial complex
allows the presence of both empty and full triangles simultaneously, which is
a more realistic model of simplicial complexes as in many real-world systems,
the presence of pairwise interactions between three individuals does not
necessarily imply a three-body interaction between them. It is important to
note that the proposed model can easily be extended to any arbitrary higher-
dimensional simplicial complexes (i.e., $D>2$), but to keep things simple, we
consider only up to $D=2$ dimension. Further, all the derived theoretical
results can be expanded to simplicial complexes of any dimension $D$ without
any difficulty.

Figure 1:  Diagrammatic representation of simplicial complex with two distinct
interaction layers. The two layers (solid red and dashed magenta links,
respectively) are built up of different sorts of interconnections for the same
nodes, where the layer with red links allows three-body interactions between
the nodes. The layers are completely independent; the presence of a link
between two nodes in one layer has no influence on their connection status in
the other.

We suppose that the equation of motion characterizing the simplicial complex
dynamics can be stated in terms of the following equation,

$\begin{array}[]{lcl}\dot{\bf x}_{i}(t)=F({\bf
x}_{i})+\sum\limits_{\beta=1}^{M}\epsilon_{\beta}\sum\limits_{j=1}^{N}{\mathscr{A}_{ij}^{[\beta]}}G_{\beta}({\bf
x}_{i},{\bf x}_{j})\\\
+\sum_{q=1}^{L}\epsilon_{SI}^{q}\sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N}{\mathscr{A}_{ijk}^{[q,SI]}}G_{SI}^{q}({\bf
x}_{i},{\bf x}_{j},{\bf x}_{k}),\hskip 10.0pti=1,2,...,N.\end{array}$ (1)

Here the $p$-dimensional state vector $\mathbf{x}_{i}(t)$ indicates the
dynamics of the unit $i$, $F:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$
illustrates the isolated node dynamics presumed to be identical for all units,
$\epsilon_{\beta}$ $(\beta=1,2,\cdots,M)$ are the real-valued coupling
strengths corresponding to the pairwise interactions within the layers,
whereas $\epsilon_{SI}^{q}$ describes the coupling strength associated with
the three-body ($2$-simplex) interactions among dynamical units of the layer
where higher-order interactions are allowed. $G_{\beta}:\mathbb{R}^{(2\times
p)}\rightarrow\mathbb{R}^{p}$, $G_{SI}^{q}:\mathbb{R}^{(3\times
p)}\rightarrow\mathbb{R}^{p}$ are the real-valued continuously differentiable
functions illustrating the coupling forms corresponding to pairwise and non-
pairwise interactions, respectively. Here we can consider the coupling
functions in arbitrary linear or nonlinear form. Furthermore, the $N\times N$
adjacency matrices $\mathscr{A}^{[\beta]}$ $(\beta=1,2,\cdots,M)$ recount the
network structures of the layers, where $\mathscr{A}^{[\beta]}_{ij}=1$, if the
$i$-th and $j$-th units are interconnected and zero otherwise.
$\mathscr{A}^{[q,SI]}$ describes $N\times N\times N$ adjacency tensors that
account for $2$-simplices, where $\mathscr{A}^{[q,SI]}_{ijk}=1$, if $(i,j,k)$
construct a full triangle and zero otherwise. If one considers $M=1$ and
$L=1$, then our mathematical framework will eventually coincide with the
generalized framework proposed earlier in Gambuzza  _et al._ (2021) for mono-
layer representation of simplicial complexes. Therefore, Eq. (1) represents a
generalized mathematical framework to investigate dynamical processes in
simplicial complexes with multiple interaction layers.

A schematic illustration of the proposed simplicial complex with $N=7$ nodes
and $M=2$ distinct interaction layers is shown in Fig. 1. The connection
mechanisms between the nodes via two different layers are represented in two
different colors; solid red lines correspond to the first layer and dotted
magenta lines for the second layer. In the first layer, both pairwise and
three-body interactions are considered, where all the three cliques (empty
triangles) are promoted to full triangles ($2$-simplices), shaded in blue
color. However, in the second layer, no empty triangles are promoted to full
triangles.

##  III Theoretical approaches

In this section, our main goal is to analytically derive the necessary
conditions for the stable complete synchronization state of the simplicial
complex (1). To do so, we first look into the conditions for invariance of the
synchronization state. The invariance condition gives the necessary
requirement on the network topology by which a stable synchronization state
may observe depending on the critical value of the coupling strengths.

###  III.1 Invariance of synchronization state

A complete synchronization state arises in the simplicial complex (1) when
each node follows in unison with the rest nodes. Mathematically, there exists
a solution $\mathbf{x}_{0}(t)\in\mathbb{R}^{p}$ such that,

$\begin{array}[]{l}\mbox{for}\leavevmode\nobreak\ \leavevmode\nobreak\
t\to\infty,\leavevmode\nobreak\ \leavevmode\nobreak\
\left\lVert\mathbf{x}_{i}(t)-\mathbf{x}_{0}(t)\right\rVert\to
0,\leavevmode\nobreak\ \leavevmode\nobreak\ i=1,2,\dots,N.\end{array}$

Following this, the corresponding synchronized manifold is defined as follows,

$\begin{array}[]{l}\mathcal{S}=\\{\mathbf{x}_{0}(t)\in\mathbb{R}^{p}:\mathbf{x}_{i}(t)=\mathbf{x}_{0}(t),\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{for}\leavevmode\nobreak\
\leavevmode\nobreak\ i=1,2,\dots,N\leavevmode\nobreak\ \leavevmode\nobreak\
\mbox{and}\leavevmode\nobreak\ \leavevmode\nobreak\
t\in\mathbb{R}^{+}\\}.\end{array}$

If the functional forms of coupling schemes are synchronization noninvasive
Gambuzza  _et al._ (2021), i.e.,

$\begin{array}[]{l}G_{\beta}({\bf x}_{0},{\bf
x}_{0})=0,\;\mbox{and}\;G_{SI}^{q}({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})=0,\end{array}$ (2)

then the existence and invariance of the synchronous solution will be assured.
However, we cast away any restriction on functional forms of coupling
mechanisms, consider arbitrary linear or nonlinear coupling functions, and
derive the invariance condition for such board circumstances.

If the simplicial complex starts evolving with the synchronization state at
any instance of time $t=t_{s}$, then
$\mathbf{x}_{i}(t_{s})=\mathbf{x}_{0}{(t_{s})}$, for $i=1,2,\dots,N$. The
evolution of the $i^{\mbox{th}}$ node of the simplicial complex at $t=t_{s}$
can then be expressed from Eq. (1) as follows,

$\begin{array}[]{l}\dot{\bf x}_{i}(t_{s})=F({\bf
x}_{0})+\sum_{\beta=1}^{M}\epsilon_{\beta}k_{(in)_{i}}^{\beta}G_{\beta}({\bf
x}_{0},{\bf x}_{0})\\\ \\\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+2\sum_{q=1}^{L}\epsilon_{SI}^{q}k_{(in)_{i}}^{[q,SI]}G_{SI}^{q}({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0}),\end{array}$ (3)

where
$k_{(in)_{i}}^{\beta}={\sum\limits_{j=1}^{N}{\mathscr{A}_{ij}^{[\beta]}}}$
denotes the number of links incident to node $i$ in each individual layer
$\beta$, i.e., in-degree of node $i$ and
$k_{(in)_{i}}^{[q,SI]}=\frac{1}{2}{\sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N}{\mathscr{A}_{ijk}^{[q,SI]}}}$
indicates the number of $2$-simplices (full triangle) incident to node $i$,
i.e., generalized in-degree of node $i$ Gallo  _et al._ (2022). Now, in order
to maintain the synchronous solution, each node must progress at the same
rate. As a result, any two distinct nodes $i$ and $j$ in the simplicial
complex should have the same velocity, i.e.,

$\begin{array}[]{l}\dot{\mathbf{x}}_{i}(t_{s})=\dot{\mathbf{x}}_{j}(t_{s}),\;\;{\mbox{for
any}\leavevmode\nobreak\ i\neq j,\leavevmode\nobreak\
\mbox{and}}\leavevmode\nobreak\ i,j=1,2,\dots,N.\end{array}$ (4)

Since the functional forms of coupling schemes are arbitrary linear or
nonlinear, Eq. (4) together with Eq. (3) gives the following relation,

$\begin{array}[]{l}k_{(in)_{i}}^{\beta}=k_{(in)_{j}}^{\beta}=k_{(in)}^{\beta}\;\mbox{(say)},\;\beta=1,2,\cdots,M\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\mbox{and}\;k_{(in)_{i}}^{[q,SI]}=k_{(in)_{j}}^{[q,SI]}=k_{(in)}^{[q,SI]}\;\mbox{(say)}.\end{array}$
(5)

Hence for the consistency of complete synchronization state, the same number
of edges must be incident to each node of layer-$\beta$ and for the layers
with three-body interactions, exactly the same number of $2$-simplices must
incident to each node. However, we here stick ourselves with only the
undirected connection structure of the simplicial complex. Therefore, in this
case, $k_{(in)}^{\beta}=k^{\beta}$, the degree of each node and
$k_{(in)}^{[q,SI]}=k^{[q,SI]}$, the generalized degree of each node.

As a result, the synchronous solution evolves according to the following
equation,

$\begin{array}[]{l}\dot{\bf x}_{0}(t_{s})=F({\bf
x}_{0})+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}G_{\beta}({\bf x}_{0},{\bf
x}_{0})\\\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+2\sum_{q=1}^{L}\epsilon_{SI}^{q}k^{[q,SI]}G_{SI}^{q}({\bf x}_{0},{\bf
x}_{0},{\bf x}_{0}).\end{array}$ (6)

At this point, we would like to emphasize that the derived invariance
condition is only necessary for the layers with arbitrary linear or nonlinear
coupling schemes between individual units. If the forms of interaction
functions within any layer are synchronization noninvasive (i.e., satisfies
Eq.(2)), then the invariance condition (Eq. (5)) is not necessary as the
coupling terms vanish at the synchronization solution, and the existence of
synchronous state can be guaranteed.

Therefore, the invariance of the complete synchronization state depends either
upon the synchronization noninvasive form of coupling schemes within the
layers or upon the regular topology of the layers (i.e., equal in-degrees of
the nodes in a particular layer). However, to investigate the stability of the
synchronization solution, we proceed with the case of arbitrary coupling
schemes. The study corresponding to synchronization noninvasive coupling forms
is just a special case of this, as the last two terms of Eq. (6) vanish for
noninvasive couplings.

###  III.2 Linear stability analysis

To derive the conditions for stable synchronization state, one considers a
small perturbation $\delta\mathbf{x}_{i}=\mathbf{x}_{i}-\mathbf{x}_{0}$ around
the synchronous solution $\mathbf{x}_{i}=\mathbf{x}_{0}\leavevmode\nobreak\
(i=1,2,\cdots,N)$ and performs linear stability analysis of Eq. (1). This
gives

$\begin{array}[]{l}\delta\dot{{\bf x}}_{i}=JF({\bf x}_{0})\delta{\bf
x}_{i}+\sum_{\beta=1}^{M}\epsilon_{\beta}\sum_{j=1}^{N}\mathscr{A}^{[\beta]}_{ij}\bigg{[}J_{1}G_{\beta}{({\bf
x}_{0},{\bf x}_{0})}\delta{\bf x}_{i}+J_{2}G_{\beta}{({\bf x}_{0},{\bf
x}_{0})}\delta{\bf x}_{j}\bigg{]}\\\
\;\;\;\;\;\;\;\;\;+\sum_{q=1}^{L}\epsilon_{SI}\sum_{j=1}^{N}\sum_{k=1}^{N}\mathscr{A}^{[q,SI]}_{ijk}\bigg{[}J_{1}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}\delta{\bf x}_{i}+J_{2}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}\delta{\bf x}_{j}+J_{3}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}\delta{\bf x}_{k}\bigg{]},\end{array}$ (7)

where $JF({\bf x}_{0})$ indicates the Jacobian of $F$ calculated at the
synchronization solution $\mathbf{x}_{0}$. $J_{1}$, $J_{2}$, and $J_{3}$ are
the Jacobian operators with respect to the first, second and third variables,
respectively. $\mathscr{L}^{[\beta]}$, $\beta=1,2,\cdots,M$ are standard zero
row-sum graph Laplacian matrices, defined by

$\begin{array}[]{l}\mathscr{L}^{[\beta]}_{ij}=\begin{cases}-\mathscr{A}^{[\beta]}_{ij},&i\neq
j\\\ k^{\beta}_{i},&i=j.\end{cases}\end{array}$

Now one has that $\sum_{j=1}^{N}\mathscr{A}^{[\beta]}_{ij}=k_{i}^{\beta}$,
$\sum_{j=1}^{N}\sum_{k=1}^{N}\mathscr{A}^{[q,SI]}_{ijk}=2k_{i}^{[q,SI]}$ and
the invariance condition (5). Hence plugging back the terms in the Eq. (7),
one eventually obtains,

$\begin{array}[]{l}\delta\dot{{\bf x}}_{i}=JF({\bf x}_{0})\delta{\bf
x}_{i}+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}\bigg{[}J_{1}G_{\beta}{({\bf
x}_{0},{\bf x}_{0})}+J_{2}G_{\beta}{({\bf x}_{0},{\bf
x}_{0})}\bigg{]}\delta{\bf
x}_{i}-\sum_{\beta=1}^{M}\epsilon_{\beta}\sum_{j=1}^{N}\mathscr{L}^{[\beta]}_{ij}J_{2}G_{\beta}{({\bf
x}_{0},{\bf x}_{0})}\delta{\bf x}_{j}\\\
\;\;\;\;\;\;+2\sum_{q=1}^{L}\epsilon^{q}_{SI}k^{[q,SI]}J_{1}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}\delta{\bf
x}_{i}+\sum_{q=1}^{L}\epsilon^{q}_{SI}\sum_{j=1}^{N}\sum_{k=1}^{N}\mathscr{A}^{[q,SI]}_{ijk}\bigg{[}J_{2}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}\delta{\bf x}_{j}+J_{3}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}\delta{\bf x}_{k}\bigg{]}.\end{array}$ (8)

At this stage of derivation, it is important to note that our approach does
not require any specific functional form of coupling functions, for instance,
diffusive or synchronization noninvasive Gambuzza  _et al._ (2021) functional
forms to derive the stability condition of the synchronization state, and as a
result, we are literally encompassing arbitrary sorts of coupling functions.
The derivation for diffusive or synchronization noninvasive functional forms
can be considered as some special cases of our approach.

Let us now use the fact that the Jacobians $J_{2}G^{q}_{SI}{({\bf x}_{0},{\bf
x}_{0},{\bf x}_{0})}$, $J_{3}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})}$ are independent of $k$ and $j$, respectively. This simplifies the
Eq. (8) as,

$\begin{array}[]{l}\delta\dot{{\bf x}}_{i}=JF({\bf x}_{0})\delta{\bf
x}_{i}+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}\bigg{[}J_{1}G_{\beta}{({\bf
x}_{0},{\bf x}_{0})}+J_{2}G_{\beta}{({\bf x}_{0},{\bf
x}_{0})}\bigg{]}\delta{\bf
x}_{i}-\sum_{\beta=1}^{M}\epsilon_{\beta}\sum_{j=1}^{N}\mathscr{L}^{[\beta]}_{ij}J_{2}G_{\beta}{({\bf
x}_{0},{\bf x}_{0})}\delta{\bf x}_{j}\\\
\;\;\;\;\;\;\;+2\sum_{q=1}^{L}\epsilon^{q}_{SI}k^{[q,SI]}J_{1}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}\delta{\bf
x}_{i}+\sum_{q=1}^{L}\epsilon^{q}_{SI}\bigg{[}\sum_{j=1}^{N}J_{2}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}\delta{\bf
x}_{j}\sum_{k=1}^{N}\mathscr{A}^{[q,SI]}_{ijk}\\\ \hskip
200.0pt+\sum_{k=1}^{N}J_{3}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})}\delta{\bf
x}_{k}\sum_{j=1}^{N}\mathscr{A}^{[q,SI]}_{ijk}\bigg{]}.\end{array}$ (9)

Then the symmetric property of the adjacency tensors $\mathscr{A}^{[q,SI]}$,
i.e.,
$\sum_{j=1}^{N}\mathscr{A}^{[q,SI]}_{ijk}=\sum_{j=1}^{N}\mathscr{A}^{[q,SI]}_{ikj}$
and the fact that $\sum_{k=1}^{N}\mathscr{A}^{[q,SI]}_{ijk}=K^{[q,SI]}_{ij}$,
where $K^{[q,SI]}_{ij}$ counts the number of $2$-simplices to which the link
$(i,j)$ participates, modifies Eq. (9) as follows,

$\begin{array}[]{l}\delta\dot{{\bf x}}_{i}=JF({\bf x}_{0})\delta{\bf
x}_{i}+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}\bigg{[}J_{1}G_{\beta}{({\bf
x}_{0},{\bf x}_{0})}+J_{2}G_{\beta}{({\bf x}_{0},{\bf
x}_{0})}\bigg{]}\delta{\bf
x}_{i}-\sum_{\beta=1}^{M}\epsilon_{\beta}\sum_{j=1}^{N}\mathscr{L}^{[\beta]}_{ij}J_{2}G_{\beta}{({\bf
x}_{0},{\bf x}_{0})}\delta{\bf x}_{j}\\\ \\\
\;\;\;\;\;\;\;\;\;\;\;\;\;+2\sum_{q=1}^{L}\epsilon^{q}_{SI}k^{[q,SI]}\bigg{[}J_{1}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}+J_{2}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})}+J_{3}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})}\bigg{]}\delta{\bf x}_{i}\\\ \\\
\;\;\;\;\;\;\;\;\;\;\;\;\;-\sum_{q=1}^{M}\epsilon^{q}_{SI}\sum_{j=1}^{N}\mathscr{L}^{[q,SI]}_{ij}\bigg{[}J_{2}G^{q}_{SI}{({\bf
x}_{0},{\bf x}_{0},{\bf x}_{0})}+J_{3}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})}\bigg{]}\delta{\bf x}_{j},\end{array}$ (10)

where $\mathscr{L}^{[q,SI]}$ illustrates the generalized zero row-sum
Laplacian matrix associated with three-body interactions, defined by

$\begin{array}[]{l}\mathscr{L}^{[q,SI]}_{ij}=\begin{cases}-{K}^{[q,SI]}_{ij},&i\neq
j\\\ 2k^{[q,SI]}_{i},&i=j.\end{cases}\end{array}$

Let us now rewrite the Eq. (10) in vector form by introducing
$\delta\mathbf{X}=[\delta\mathbf{x}_{1}^{tr},\delta\mathbf{x}_{2}^{tr},\cdots,\delta\mathbf{x}_{N}^{tr}]^{tr}$,
where $[\;\;]^{tr}$ denotes the matrix transpose and denoting the terms as
$h=JF({\bf x}_{0})$, $H_{\beta}=J_{1}G_{\beta}{({\bf x}_{0},{\bf
x}_{0})}+J_{2}G_{\beta}{({\bf x}_{0},{\bf x}_{0})}$,
$\Phi_{\beta}=J_{2}G_{\beta}{({\bf x}_{0},{\bf x}_{0})}$,
$\Psi_{q}=J_{1}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})}+J_{2}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})}+J_{3}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf x}_{0})}$,
$\chi_{q}=J_{2}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf
x}_{0})}+J_{3}G^{q}_{SI}{({\bf x}_{0},{\bf x}_{0},{\bf x}_{0})}$. Eventually,
the equation in terms of stake vectors becomes as follows,

$\begin{array}[]{l}\delta\dot{{\bf X}}=\bigg{[}(I_{N}\otimes
h)+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}(I_{N}\otimes
H_{\beta})-\sum_{\beta=1}^{M}\epsilon_{\beta}(\mathscr{L}^{[\beta]}\otimes\Phi_{\beta})+2\sum_{q=1}^{L}\epsilon^{q}_{SI}k^{[q,SI]}(I_{N}\otimes\Psi_{q})\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-\sum_{q=1}^{L}\epsilon^{q}_{SI}(\mathscr{L}^{[q,SI]}\otimes\chi_{q})\bigg{]}\delta{\bf
X},\end{array}$ (11)

where $\otimes$ and $I_{N}$ symbolize the Kronecker product and
$N$-dimensional identity matrix, respectively. This linearized set of
variational Eq. (11) has two parts: one for motion along the synchronization
manifold, termed parallel modes, and another for motion across the manifold,
called transverse modes. All transverse modes must converge to zero in time
for the synchronization state to be stable. The linear stability analysis is
then carried out by decoupling the variational equation into parallel and
transverse modes and determining whether or not the latter is extinct. Now all
the Laplacian matrices $\mathscr{L}^{[\beta]}$ $(\beta=1,2,\cdots,M)$ as well
$\mathscr{L}^{[q,SI]}$ $(q=1,2,\cdots,L)$ are zero row-sum real symmetric
matrices. Hence they are diagonalizable, and all their eigenvalues are non-
negative real numbers with one common smallest eigenvalue $\gamma_{1}=0$, and
the set of eigenvectors form an orthonormal basis of $\mathbb{R}^{N}$. To
separate transverse and parallel modes from Eq. (11), we project the stack
variable $\delta\mathbf{X}$ onto the basis of eigenvectors
$V^{[1]}=\\{v_{1},v_{2},\cdots,v_{N}\\}$ corresponding to $\mathscr{L}^{[1]}$,
the Laplacian associated with layer-$1$ by defining new variable
$\eta=(V^{[1]}\otimes I_{p})^{-1}\delta\mathbf{X}$, where
$v_{1}=\big{(}\frac{1}{\sqrt{N}},\frac{1}{\sqrt{N}},\cdots,\frac{1}{\sqrt{N}}\big{)}^{tr}$
is the eigenvector corresponding to the smallest eigenvalue $\gamma_{1}=0$.
The set of eigenvectors chosen is fully arbitrary; any other basis can be
chosen, and all other sets of eigenvectors will ultimately change to it via
unitary matrix transformation. Then in terms of the new variable, the
variational Eq. (11) transforms as follows,

$\begin{array}[]{l}\dot{\eta}=\bigg{[}(I_{N}\otimes
h)+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}(I_{N}\otimes
H_{\beta})-\epsilon_{1}(\Gamma\otimes\Phi_{1})-\sum_{\beta=2}^{M}\epsilon_{\beta}(\tilde{\mathscr{L}}^{[\beta]}\otimes\Phi_{\beta})\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+2\sum_{q=1}^{L}\epsilon^{q}_{SI}k^{[q,SI]}(I_{N}\otimes\Psi_{q})-\sum_{q=1}^{L}\epsilon^{q}_{SI}(\tilde{\mathscr{L}}^{[q,SI]}\otimes\chi_{q})\bigg{]}\eta,\end{array}$
(12)

where
${V^{[1]}}^{-1}\mathscr{L}^{[1]}V^{[1]}=\Gamma=diag\\{\gamma_{1}=0,\gamma_{2},\cdots,\gamma_{N}\\}$,
${V^{[1]}}^{-1}\mathscr{L}^{[\beta]}V^{[1]}=\tilde{\mathscr{L}}^{[\beta]}$ and
$\tilde{\mathscr{L}}^{[q,SI]}={V^{[1]}}^{-1}\mathscr{L}^{[q,SI]}V^{[1]}$. Now,
as all the Laplacians are zero row-sum matrices, the transformed Eq. (12) can
be decoupled into two parts as follows

$\displaystyle\dot{\eta}_{1}=(h+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}H_{\beta}+2\sum_{q=1}^{L}\epsilon^{q}_{SI}k^{[q,SI]}\Psi_{q})\eta_{1},$
(13a)
$\displaystyle\dot{\eta}_{i}=(h+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}H_{\beta}+2\sum_{q=1}^{L}\epsilon^{q}_{SI}k^{[q,SI]}\Psi_{q}-\epsilon_{1}\gamma_{i}\Phi_{1})\eta_{i}-\sum_{\beta=2}^{M}\epsilon_{\beta}\sum_{j=2}^{N}\tilde{\mathscr{L}}^{[\beta]}_{ij}\Phi_{\beta}\eta_{j}-\sum_{q=1}^{L}\epsilon^{q}_{SI}\sum_{j=2}^{N}\tilde{\mathscr{L}}^{[q,SI]}_{ij}\chi_{q}\eta_{j},$
(13b)

where $\eta_{1}$ indicates the parallel mode and the transverse modes are
symbolized by $\eta_{i}$, $i=2,3,\cdots,N$. Therefore, the problem of
synchronization stability is reduced in solving the coupled linear
differential Eq. (13b) associated with transverse mode to calculate the
maximum Lyapunov exponent (MLE). The Eq. (13b) is called the master stability
equation. For the synchronization state to be stable, the MLE must be
negative.

Due to more complexity of multiple layers and higher-order structures, the
transverse variational Eq. (13b) is a $(N-1)p$-dimensional linear coupled
differential equation, and in general, it is difficult to decouple it to lower
dimensional form. However, there is a suitable instance for which the coupled
transverse modes can be optimally separated, and as a result, the
$(N-1)p$-dimensional coupled equation transforms into $(N-1)$ numbers of
$p$-dimensional linear differential equations. The relevant instance is as
follows-

When all the pairwise Laplacians $\mathscr{L}^{[\beta]}$ and the generalized
Laplacians $\mathscr{L}^{[q,SI]}$ commute with each other then the transverse
Eq. (13b) can be decoupled into $(N-1)$ numbers of $p$-dimensional equations
as,

$\begin{array}[]{l}\dot{\eta}_{i}=\bigg{(}h+\sum_{\beta=1}^{M}\epsilon_{\beta}k^{\beta}H_{\beta}+2\sum_{q=1}^{L}\epsilon^{q}_{SI}k^{[q,SI]}\Psi_{q}\\\
\\\
-\sum_{\beta=1}^{M}\epsilon_{\beta}\gamma^{\beta}_{i}\Phi_{\beta}-\sum_{q=1}^{L}\epsilon^{q}_{SI}\gamma_{i}^{[q,SI]}\chi_{q}\bigg{)}\eta_{i},\leavevmode\nobreak\
i=2,3,\cdots,N,\end{array}$ (14)

where
${V^{[1]}}^{-1}\mathscr{L}^{[\beta]}V^{[1]}=\Gamma^{\beta}=diag\\{0=\gamma^{\beta}_{1}<\gamma^{\beta}_{2}\cdots\leq\gamma^{\beta}_{N}\\}$
and
${V^{[1]}}^{-1}\mathscr{L}^{[q,SI]}V^{[1]}=\Gamma^{[q,SI]}=diag\\{0=\gamma^{[q,SI]}_{1}<\gamma^{[q,SI]}_{2}\cdots\leq\gamma^{[q,SI]}_{N}\\}$,
as in this scenario each of the Laplacians are diagonalizable with respect to
the eigenvector matrix of the Laplacian of 1st layer $\mathscr{L}^{[1]}$.

##  IV Numerical validation

A series of outcomes is presented below, all of which support the validity and
broad applicability of our method. Here we consider two three-dimensional
chaotic oscillators, namely the Sherman neuronal model Jalil  _et al._ (2012);
Sherman (1994) associated with pancreatic $\beta$-cells and the Rössler
oscillators Rössler (1976), coupled through various nonlinear pairwise and
non-pairwise coupling functions interacting via two $(M=2)$ distinct
connection layers. For numerical simulations, we investigate the complete
synchronization of the simplicial complex composed of $N$ nodes. We employ the
synchronization error to analyze complete synchronization as,

$\begin{array}[]{l}E=\lim\limits_{T\to\infty}\dfrac{1}{T}\bigintss_{t_{trans}}^{t_{trans}+T}\sum\limits_{i,j=1}^{N}\dfrac{\|\mathbf{x}_{j}-\mathbf{x}_{i}\|}{N(N-1)},\end{array}$
(15)

where $\|\cdot\|$ symbolizes the Euclidean norm, $t_{trans}$ determines the
transient of the numerical simulation, and $T$ is a sufficiently large
positive number. The zero value of the synchronization error $E$ indicates the
emergence of complete synchrony. The simplicial complex (1) is integrated
using the $4$th order Runge-Kutta scheme subject to integration time step
$dt=0.001$ for neuron dynamics and $dt=0.01$ for Rössler system, over a period
of $t=50000$ time units with a transient of $40000$ time units. Each dynamical
node’s initial condition is randomly selected from the solitary node dynamics
phase space.

Throughout this section, we study the impact of higher-order interactions on
complete synchrony with arbitrary coupling mechanisms and verify all the
results with our analytical conjectures. In a chaotic bursting regime, we
select the system parameters of individual nodes and focus on variation of
coupling strengths along with different coupling schemes and network
topologies.

###  IV.1 Application to neuron dynamics

Recent studies in neuroscience have revealed that neurons may communicate
through many-body interactions Battiston  _et al._ (2020); Parastesh  _et al._
(2022); Ince  _et al._ (2009); Amari  _et al._ (2003). Specifically,
astrocytes and other glial cells are considered a possible source of many-body
interactions Fellin  _et al._ (2004); Tlaie  _et al._ (2019), as they make
contact with thousands of synapses and actively modulate their functions Allen
and Barres (2009). Although it is still unclear how to account for these many-
body interactions in nonlinear neuronal models, we think our simplicial
complex framework may suitably fit the neuronal network systems to study
synchronization in the presence of many-body interactions.

One neuron generally transfers information to other neurons through gap
junctions and chemical synapses. A signal is sent chemically through a
chemical synapse by neurochemical compounds like acetylcholine, gamma-
aminobutyric acid, dopamine, and serotonin, which are contained inside
microscopic synaptic vesicles. Exocytoses release neurotransmitters
probabilistically from a presynaptic neuron into the synaptic gap. Following
that, these compounds attach to certain receptor molecules in nearby
postsynaptic neuronal cells. The gap between the pre-and postsynaptic ends
could be as much as 20–40 nm Hormuzdi  _et al._ (2004) in this scenario. A gap
junction connects the cytoplasm of neighboring cells in electrical synapses.
Electric current, calcium, cyclic AMP, and inositol-1,4,5 triphosphate pass in
both directions between the presynaptic end and the postsynaptic neuron. The
pre-and postsynaptic neurons’ membranes are extraordinarily close to each
other, about 3.5 nm Kandel  _et al._ (2000), and signal transmission is
significantly faster than chemical transmission. In most neurological systems,
both types of synapses exist. In interneuronal transmission, both types of
interaction may not always be present at all times; rather, they operate
independently Pereda (2014) throughout time. A mixed synapse occurs when two
neurons communicate through both chemical and electrical synapses. A
heterosynaptic connection, on the other hand, occurs when a neuron is coupled
to two different neurons, one by a chemical synapse and the other through an
electrical synapse.

Here we analyze a group of pancreatic $\beta$ cells, represented by the
paradigmatic Sherman model Jalil  _et al._ (2012); Sherman (1994), subjected
to both pairwise and three-body interactions. The velocity profile of a single
$\beta$ cell is given by

$\begin{array}[]{l}F(\mathbf{x})=F\begin{bmatrix}V\\\ n\\\ s\end{bmatrix}\\\
=\begin{bmatrix}f(V,n,s)=\frac{1}{\tau}\\{-I_{Ca}(V)-I_{K}(V,n)-I_{S}(V,s)\\}\\\
g(V,n,s)=\frac{1}{\tau}\\{\mu[n^{\infty}(V)-n]\\}\\\
h(V,n,s)=\frac{1}{\tau_{s}}\\{s^{\infty}(V)-s\\}\end{bmatrix},\end{array}$

where two fast currents: calcium $I_{Ca}$, persistent potassium $I_{K}$, and a
slow potassium current $I_{S}$ are given by $I_{S}(V,s)=g_{S}s(V-E_{K})$,
$I_{K}(V,n)=g_{K}n(V-E_{K})$, and $I_{Ca}(V)=g_{Ca}m^{\infty}(V-E_{Ca})$,
respectively. $V$ is the membrane potential corresponding to the reversal
potential $E_{K}=-0.075$, $E_{Ca}=0.025$. $m$, $n$, and $s$ are the voltage-
dependent gating variables. The maximum conductance and time constants are
$g_{Ca}=3.6$, $g_{K}=10$, $g_{S}=4$ and $\tau=0.02$, $\tau_{s}=5$. $\mu=1$, an
auxiliary scaling factor, manages the time scale of the persistent potassium
channels. The values of the gating variables at steady state are

$\begin{array}[]{l}m^{\infty}(V)=\\{1+\exp[-83.34(V+0.02)]\\}^{-1},\\\
n^{\infty}(V)=\\{1+\exp[-178.57(V+0.016)]\\}^{-1},\;\mbox{and}\\\
s^{\infty}(V)=\\{1+\exp[-100(V+0.035345)]\\}^{-1}.\end{array}$

The parameter values are taken in such a way that the isolated node dynamics
exhibits chaotic behavior (spiking bursting).

We assume that the neurons are interconnected concurrently through two (M = 2)
distinct interacting layers, subject to gap junction coupling via electrical
synapses and chemical ion transit via chemical synapses. Therefore, the
pairwise connections in the layers are described by the coupling functions
$G_{1}(\mathbf{x}_{i},\mathbf{x}_{j})=[V_{j}-V_{i},0,0]^{tr}$ to describe gap
junctions, and
$G_{2}(\mathbf{x}_{i},\mathbf{x}_{j})=[(E_{S}-V_{i})\Gamma(V_{j}),0,0]^{tr}$
to describe chemical ion transportation through chemical synapses. For three-
body interactions, we consider two different instances: (1) only the layer
with electrical synapse coupling allows three-body connection, described by
the diffusive coupling function
$G^{e}_{SI}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k})=[V_{j}+V_{k}-2V_{i},0,0]^{tr}$
and (2) three-body interactions are permitted only in the layer with chemical
synapse coupling, represented by the nonlinear coupling form
$G^{c}_{SI}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k})=[(E_{S}-V_{i})\\{\Gamma(V_{j})+\Gamma(V_{k})\\},0,0]^{tr}$.
Here, the sigmoidal input-output function is,

$\displaystyle\Gamma(V)=\frac{1}{1+\exp(\lambda_{s}(V-\Theta_{s}))},$

describes the procedure for activation and deactivation of nonlinear chemical
synapse. The synaptic reversal potential is $E_{s}=-0.02$. Slope of the
sigmoidal function, and synaptic firing threshold are determined by the real-
valued constants $\lambda_{s}$ and $\Theta_{s}$, which are fixed at
$\Theta_{s}=-0.045$ and $\lambda_{s}=-1000$. Therefore, the equation of
motions describing the dynamics of the neuronal simplicial complex for two
instances are given by,

$\begin{array}[]{l}\tau\dot{V}_{i}=-I_{Ca}(V_{i})-I_{K}(V_{i},n_{i})-I_{S}(V_{i},s_{i})\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\epsilon\sum\limits_{j=1}^{N}{\mathscr{A}_{ij}^{[e]}}(V_{j}-V_{i})+g_{c}(E_{S}-V_{i})\sum\limits_{j=1}^{N}{\mathscr{A}_{ij}^{[c]}}\Gamma(V_{j})\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\epsilon^{e}_{SI}\sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N}{\mathscr{A}_{ijk}^{[e,SI]}}(V_{j}+V_{k}-2V_{i}),\\\
\tau\dot{n}_{i}=\mu[n^{\infty}(V_{i})-n_{i}],\\\
\tau_{s}\dot{s}_{i}=s^{\infty}(V_{i})-s_{i},\end{array}$ (16)

and

$\begin{array}[]{l}\tau\dot{V}_{i}=-I_{Ca}(V_{i})-I_{K}(V_{i},n_{i})-I_{S}(V_{i},s_{i})\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\epsilon\sum\limits_{j=1}^{N}{\mathscr{A}_{ij}^{[e]}}(V_{j}-V_{i})+g_{c}(E_{S}-V_{i})\sum\limits_{j=1}^{N}{\mathscr{A}_{ij}^{[c]}}\Gamma(V_{j})\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+\epsilon^{c}_{SI}(E_{S}-V_{i})\sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N}{\mathscr{A}_{ijk}^{[c,SI]}}[\Gamma({V}_{j})+\Gamma({V}_{k})],\\\
\tau\dot{n}_{i}=\mu[n^{\infty}(V_{i})-n_{i}],\\\
\tau_{s}\dot{s}_{i}=s^{\infty}(V_{i})-s_{i},\end{array}$ (17)

where $\epsilon$, $g_{c}$ symbolize electrical and chemical synaptic coupling
strengths, and $\epsilon^{e}_{SI}$, $\epsilon^{c}_{SI}$ are the corresponding
higher-order coupling strengths. The adjacency matrices $\mathscr{A}^{[e]}$
and $\mathscr{A}^{[c]}$ describes the pairwise network connectivity of
electrical and chemical synapses, and $\mathscr{A}^{[e,SI]}$,
$\mathscr{A}^{[c,SI]}$ represent adjacency tensors corresponding to three-body
interactions. In general, the connectivity network of electrical and chemical
synapses are taken to be bidirectional and unidirectional, respectively. But
for our case, we consider both the connectivity network to be bidirectional.
Furthermore, since the coupling forms associated with chemical ion
transportation are nonlinear synchronization invasive (i.e.,
$G_{2}(\mathbf{x},\mathbf{x})\neq 0$ and
$G^{c}_{SI}(\mathbf{x},\mathbf{x},\mathbf{x})\neq 0$), we consider identical
node degree $k_{i}^{c}=k^{c}$ corresponding to $1$-simplex and
$k_{i}^{[c,SI]}=k^{[c,SI]}$, the number of $2$-simplices adjacent to a node
for synchronization invariance condition to be satisfied. However, since the
coupling forms corresponding to gap junctions are linear diffusive, no further
restriction on connection topology is needed in this scenario.

Figure 2: Synchronization in coupled Sherman model with higher-order gap
junction coupling. (a) Synchronization error $E$ as a function of pairwise gap
junction coupling strength $\epsilon$ for a fixed value of pairwise chemical
synaptic coupling $g_{c}=0.01$. The synchronization error in the absent
$(\epsilon^{e}_{SI}=0)$ and present $(\epsilon^{e}_{SI}=0.01)$ of higher-order
gap junction coupling is displayed by black and red circle curves,
respectively. The dashed black and red lines in (b) depict the corresponding
maximum Lyapunov exponent curves. The variation of synchronization error $E$
in $(\epsilon,g_{c})$ parameter planes for (c) $\epsilon^{e}_{SI}=0$ and (d)
$\epsilon^{e}_{SI}=0.01$ are represented. The solid white curves in (c) and
(d) denote the theoretically derived synchronization threshold curves obtained
from Eq. (18). Here the other parameter values are $k_{c}=4$, $k_{sw}=4$ and
$p_{sw}=0.15$.

####  IV.1.1 Higher-order gap junction coupling

We first consider the layer with a connection through electrical synapses to
allow three-body interactions among the neurons. Therefore, the equations of
motion governing the dynamics of our simplicial complex satisfy Eq. (16). We
consider the connectivity topology of electrical synapses to be small-world
Muldoon  _et al._ (2016), given by the adjacency matrix $\mathscr{A}^{[e]}$.
Particularly, we use the Watts–Strogatz (WS) graph model Watts and Strogatz
(1998), which starts with a non-local ring with $N=200$ nodes, $k_{sw}=4$
nearest neighbor on each side, and the edge rewiring probability
$p_{sw}=0.15$. The $2$-simplices are constructed by promoting all the empty
triangles of the small-world network to full triangles, represented by the
adjacency tensor $\mathscr{A}^{[SI]}$. On the other hand, the chemical
synaptic network is considered to be a random network, described by the
adjacency matrix $\mathscr{A}^{[c]}$ with constant node degree $k_{c}=4$.

We begin by examining the development of synchronization error $E$ (given by
(15)) for varying coupling strengths to understand the emergence of neuronal
synchrony in our simplicial complex (16). Figure 2(a) shows the variation in
$E$ as a function of the pairwise electrical synaptic coupling strength
$\epsilon$ for various values of higher-order coupling strength
$\epsilon^{e}_{SI}$, with fixed chemical synaptic coupling $g_{c}=0.01$. The
figure delineates two different curves corresponding to two different values
of $\epsilon^{e}_{SI}$. Solid black circle corresponds synchronization error
for $\epsilon^{e}_{SI}=0$, i.e., when the higher-order interactions are not
considered. The critical point to achieve synchronization is $\epsilon\approx
0.58$ in this scenario. While for a non-zero value of higher-order coupling
$\epsilon^{e}_{SI}=0.01$ ( shown in the red circles ), the synchrony is
achieved at a much lower value of the coupling $\epsilon=0.35$. Therefore, our
observation indicates that synchronization can be improved in the neuronal
simplicial complex by introducing three-body interactions through electrical
synaptic coupling.

To verify our numerical result, we then proceed through the linear stability
analysis of the synchronization state based on the MSF approach. Followed by a
series of calculations detailed in Sec. III.2, the variational equation
transverse to synchronization manifold can be written analogous to Eq. (13b)
as,

$\begin{array}[]{l}\tau\delta\dot{V}_{i}=Jf(\mathbf{x}_{0})\delta\mathbf{x}_{i}-\epsilon\gamma_{i}\delta
V_{i}-2\epsilon^{e}_{SI}\sum\limits_{j=2}^{N}\tilde{\mathscr{L}}_{ij}^{[e,SI]}\delta
V_{j}\\\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+g_{c}k_{c}[-\Gamma(V_{0})+(E_{S}-V_{0})\Gamma^{\prime}(V_{0})]\delta V_{i}\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-g_{c}(E_{S}-V_{0})\sum\limits_{j=2}^{N}{\tilde{\mathscr{L}}_{ij}^{[c]}}\Gamma^{\prime}(V_{0})\delta
V_{i},\\\ \tau\delta\dot{n}_{i}=Jg(\mathbf{x}_{0})\delta\mathbf{x}_{i},\\\
\tau_{s}\delta\dot{s}_{i}=Jh(\mathbf{x}_{0})\delta\mathbf{x}_{i},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ i=2,3,\cdots,N,\end{array}$ (18)

where $J$ denotes the Jacobian operator, ′ represents derivative with respect
to $V$, $\gamma_{1}(=0)<\gamma_{2}\leq\cdots\leq\gamma_{N}$ are the
eigenvalues of the Laplacian matrix $\mathscr{L}^{[e]}$.
$\tilde{\mathscr{L}}^{[c]}$ and $\tilde{\mathscr{L}}^{[e,SI]}$ are the
$(N-1)\times(N-1)$ real matrices obtained through transformation the
Laplacians $\mathscr{L}^{[c]}$ and $\mathscr{L}^{[SI]}$ by the matrix of
orthonormal eigenvectors of $\mathscr{L}^{[e]}$. Here
$\mathbf{x}_{0}=(V_{0},n_{0},s_{0})^{tr}$ represents the states of the
synchronized solution given by,

$\begin{array}[]{l}\tau\dot{V}_{0}=-I_{Ca}(V_{0})-I_{K}(V_{0},n_{0})-I_{S}(V_{0},s_{0})\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ +g_{c}k_{c}(E_{S}-V_{0})\Gamma(V_{0}),\\\
\tau\dot{n}_{0}=\mu[n^{\infty}(V_{0})-n_{0}],\\\
\tau_{s}\dot{s}_{0}=s^{\infty}(V_{0})-s_{0}.\end{array}$ (19)

Solving $3(N-1)$-dimensional linearized Eq. (18) along with the
$3$-dimensional nonlinear Eq. (19) for the calculation of maximum Lyapunov
exponent gives the condition for stable synchronization state. Let $\Lambda$
be the maximum Lyapunov exponent transverse to the synchronization manifold.
Then the stable synchronous solution is obtained if $\Lambda$ becomes negative
with varying coupling strengths. Figure 2(b) portrays the variation of
$\Lambda$ with respect to $\epsilon$ for same set of other coupling strengths
$g_{c}=0.01$ and $\epsilon^{e}_{SI}=0,0.01$ used in Fig. 2(a). The dashed
black and red curves display the variation of $\Lambda$ for
$\epsilon^{e}_{SI}=0$ and $\epsilon^{e}_{SI}=0.01$, respectively. It is
clearly observable that both the curves cross the zero line and become
negative exactly at the same critical values of $\epsilon$ where the
corresponding curves of synchronization error $E$ achieved zero value as
reported in Fig. 2(a). Hence it affirms our obtained result about the
enhancement of synchrony in the neuronal simplicial complex by the inclusion
of higher-order electrical synaptic coupling.

Thereafter, we concentrate on the effect of higher-order interaction in the
$(\epsilon,g_{c})$ parameter plane by varying pairwise electrical and chemical
synaptic coupling strengths for fixed values of non-pairwise electrical
synaptic coupling parameter $\epsilon^{e}_{SI}$. We plot the corresponding
synchronization error $E(\epsilon,g_{c})$ in Figs. 2(c) and 2(d), where
colorbar indicates the variation of $E$. Figure 2(c) depicts how
synchronization emerges for concurrent variation in $\epsilon\in[0,0.7]$ and
$g_{c}\in[0,0.03]$ whenever there is no higher-order coupling, i.e.,
$\epsilon^{e}_{SI}=0$. As we can notice that as $g_{c}$ is increased, lower
values of $\epsilon$ are found to be sufficient for the achievement of
complete synchrony. A similar phenomenon is observed when we increase
$\epsilon^{e}_{SI}$ to $0.01$, depicted in Fig. 2(d). But in this case, the
critical values of $\epsilon$ for increasing $g_{c}$ are sufficiently smaller
than in the previous scenario. The threshold for synchrony lowers down from
$\epsilon\approx 0.66$ to $\epsilon\approx 0.35$ even with no chemical
synaptic effect, i.e., when $g_{c}=0$. However, the threshold value of $g_{c}$
remains the same in the absence of pairwise electrical synaptic coupling,
i.e., $\epsilon=0$. The solid white curves in Figs. 2(c) and 2(d) are
theoretical predictions of the synchronization thresholds by means of maximum
Lyapunov exponents obtained from Eq. (18), which indeed fully matches with the
numerical results obtained by solving Eq. (16) for synchronization error $E$.

Figure 3: Synchronization in coupled Sherman model with higher-order gap
junction coupling and all-to-all pairwise chemical synaptic connections.
Synchronization error $E$ in $(\epsilon,g_{c})$ plane for (a)
$\epsilon^{e}_{SI}=0$, and (b) $\epsilon^{e}_{SI}=0.01$. The other parameter
values are $p_{sw}=0.15$ and $k_{sw}=4$. The solid white curves depict the
critical curves for which $\Lambda=0$ obtained from Eq. (18), while the
regions below and above it denote the unstable and stable synchronization
states, respectively.

Now we investigate the effect of higher-order electrical synaptic coupling on
neuronal synchronization by changing the network topology of the layer
interacting via chemical synapses. Therefore, we consider that each neuron is
connected to all the other neurons through chemical synapses, i.e., the
connection topology of the layer interacting through chemical synapses is
considered to be an all-to-all network. Figure 3 represents the results for
variation of synchronization error $E$ in $(\epsilon,g_{c})$ parameter plane
for two values of higher-order coupling strength- $\epsilon^{e}_{SI}=0$, i.e.,
when no higher-order effect is considered (Fig. 3(a)) and a non-zero value of
higher-order electrical synaptic coupling $\epsilon^{e}_{SI}=0.01$ (Fig.
3(b)). It is clearly observable that the phenomena remain almost the same as
discovered previously in Fig. 2, i.e., the critical value of $\epsilon$ for
synchrony decreases with increasing $g_{c}$. Further, the critical values for
$\epsilon$ remain almost the same, only the critical values of $g_{c}$
decrease as the number of chemical synaptic connections increases in this
scenario.

Figure 4:  Synchronization in coupled Sherman model with higher-order chemical
synaptic coupling. Synchronization error $E$ in $(\epsilon,g_{c})$ parameter
plane for fixed value of higher-order chemical synaptic coupling strength
$\epsilon^{c}_{SI}=10^{-6}$. The other parameter values are $k_{sw}=4$ and
$p_{sw}=0.15$. The solid white curve depicts the critical curve for which
$\Lambda=0$ obtained from Eq. (20), while the regions below and above it
denote the unstable and stable synchronization states, respectively.

####  IV.1.2 Higher-order chemical synaptic coupling

We now consider that the layer interacting through chemical synapses allows
the three-body interactions among neurons, and the layer interacting via gap
junction allows only pairwise connections. Then the equation of motion
governing the dynamics of the simplicial complex is given by Eq. (17). Now one
can notice that both the pairwise and non-pairwise coupling functions
corresponding to chemical ion transportation are nonlinear in form and further
$G_{2}(\mathbf{x},\mathbf{x})\neq 0$,
$G^{c}_{SI}(\mathbf{x},\mathbf{x},\mathbf{x})\neq 0$. So to achieve complete
neuronal synchrony according to the invariance condition (5), the topology of
the layer interacting through chemical synapses should qualify further
restriction about the same node degree $k_{i}^{c}=k^{c}$ and
$k_{i}^{[c,SI]}=k^{[c,SI]}$. We, therefore, consider the simplest
configuration that the neurons are all-to-all connected. We consider this
scenario so that we can compare the effect of higher-order chemical synaptic
coupling with the higher-order electrical synaptic coupling discussed in the
previous section (Fig. 3). One can choose other network topologies that
satisfy the synchronization invariance condition, for example, the non-local
network also fulfills the demand of the invariance condition.

Figure 4 delineates the synchrony and desynchrony regions in terms of
synchronization error $E$ in $(\epsilon,g_{c})$ parameter plane for fixed non-
pairwise chemical synaptic coupling strength $\epsilon^{c}_{SI}=10^{-6}$. It
is easily discoverable that the synchronization region is greatly enhanced in
this case. Furthermore, in the presence of higher-order chemical synaptic
coupling, not only the critical value of $g_{c}$ lowers down when $\epsilon=0$
but also the critical value of $\epsilon$ decreases when $g_{c}=0$, which is
not the case with only higher-order electrical synaptic coupling as elucidated
in Fig. 3(b). In that scenario, critical value of $g_{c}$ when $\epsilon=0$ is
unaffected with increasing non-pair electrical synaptic coupling strength
$\epsilon^{e}_{SI}$ (Figs. 3(a), 3(b)).

Now, since the neurons are connected globally through the chemical synapses,
the pairwise and generalized Laplacians $\mathscr{L}^{[c]}$ and
$\mathscr{L}^{[c,SI]}$ must satisfy the relation,
$\mathscr{L}^{[c,SI]}=(N-2)\mathscr{L}^{[c]}$. Furthermore, the Laplacians
$\mathscr{L}^{[c]}$ and $\mathscr{L}^{[e]}$ commute with each other. So the
master stability Eq. (13b) decouples optimally in this instance, and as a
result, the transverse variational equation can be written analogous to Eq.
(14) as,

$\begin{array}[]{l}\tau\delta\dot{V}_{i}=Jf(\mathbf{x}_{0})\delta\mathbf{x}_{i}+\bigg{\\{}g_{c}(N-1)[-\Gamma(V_{0})+\Gamma^{\prime}(V_{0})(E_{S}-V_{0})]\\\
\\\
+2(N-1)(N-2)\epsilon^{c}_{SI}[-\Gamma(V_{0})+\Gamma^{\prime}(V_{0})(E_{S}-V_{0})]-\epsilon\gamma^{e}_{i}\\\
\\\
-g_{c}\gamma^{c}_{i}\Gamma^{\prime}(V_{0})(E_{S}-V_{0})-2\epsilon^{c}_{SI}\gamma^{[c,SI]}_{i}\Gamma^{\prime}(V_{0})(E_{S}-V_{0})\bigg{\\}}\delta
V_{i},\\\ \\\ \tau\delta\dot{n}_{i}=Jg(\mathbf{x}_{0})\delta\mathbf{x}_{i},\\\
\\\
\tau_{s}\delta\dot{s}_{i}=Jh(\mathbf{x}_{0})\delta\mathbf{x}_{i},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ i=2,3,\cdots,N,\end{array}$ (20)

where the eigenvalues of the Laplacian $\mathscr{L}^{[e]}$ are symbolized by
$\gamma^{e}_{1}(=0)<\gamma^{e}_{2}\leq\cdots\leq\gamma^{e}_{N}$.
$\gamma^{c}_{i}=N$ and $\gamma^{[c,SI]}_{i}=N(N-2)$, $i=2,3,\cdots,N$ are the
eigenvalues of the Laplacians $\mathscr{L}^{[c]}$ and $\mathscr{L}^{[c,SI]}$,
respectively. In Fig. 4, the solid white line depicts the theoretical
conjecture obtained from Eq. (20) solving for maximum Lyapunov exponents,
which confirms that the numerical result obtained by solving Eq. (17) is in
good agreement with our analytical derivation.

Figure 5: Synchronization in coupled Rössler oscillators with higher-order
nonlinear coupling function. Synchronization error $E$ in
$(\epsilon_{1},\epsilon_{2})$ parameter plane for (a) $\epsilon_{SI}=0$ and
(b) $\epsilon_{SI}=0.00005$. The contour plots of synchronization error $E$ in
the parameter space (c) $(\epsilon_{1},\epsilon_{SI})$ and (d)
$(\epsilon_{2},\epsilon_{SI})$ for fixed values of $\epsilon_{2}=0.005$ and
$\epsilon_{1}=0.05$ are depicted, respectively. The solid white curves
correspond to the theoretical predictions obtained from Eq. (22).

###  IV.2 Application to simplicial complexes of coupled Rössler oscillators

To further illustrate how higher-order interactions play a crucial role in the
emergence of the synchronization state, we investigate the simplicial complex
Eq. (1) composed of $N=100$ nodes having individual node dynamics as Rössler
oscillators, interacting with nonlinear coupling functions through two $(M=2)$
distinct interaction layers, where three-body interactions are allowed in only
one layer. Then the equation of motion describing the dynamics of the
simplicial complex can be written as follows,

$\begin{array}[]{l}\dot{x}_{i}=-y_{i}-z_{i}+\epsilon_{1}\sum_{j=1}^{N}\mathscr{A}^{[1]}_{ij}(\alpha-
x_{i})(x_{j}-\beta)^{2},\\\ \\\
\dot{y}_{i}=x_{i}+ay_{i}+\epsilon_{2}\sum_{j=1}^{N}\mathscr{A}^{[2]}_{ij}(\alpha-
y_{i})(y_{j}-\beta)^{2}\\\ \\\
+\epsilon_{SI}\sum_{j=1}^{N}\sum_{k=1}^{N}\mathscr{A}^{[SI]}_{ijk}(\alpha-
y_{i})[(y_{j}-\beta)^{2}+(y_{k}-\beta)^{2}],\\\ \\\
\dot{z}_{i}=b+z_{i}(x_{i}-c),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ i=1,2,\cdots,N,\end{array}$ (21)

where $G_{1}=[(\alpha-x_{i})(x_{j}-\beta)^{2},0,0]^{tr}$, $G_{2}=[0,(\alpha-
y_{i})(y_{j}-\beta)^{2},0]^{tr}$ are the pairwise coupling forms in the layers
with corresponding coupling strengths $\epsilon_{1}$, $\epsilon_{2}$ and
$G_{SI}=[0,(\alpha-y_{i})\\{(y_{j}-\beta)^{2}+(y_{k}-\beta)^{2}\\},0]^{tr}$ is
the functional form of the non-pairwise coupling acting in the second layer
with coupling strength $\epsilon_{SI}$. The connection topology of the first
layer is considered to be the random network with constant node degree
$k_{1}=5$ to fulfill the synchronization invariance condition, and for the
second layer, we consider the nodes are all-to-all connected. Here, the system
parameters are chosen as $a=b=0.2$, $c=5.7$ for the isolated node dynamics to
be in a chaotic state and the coupling parameters are taken as $\alpha=0.37$,
$\beta=-0.37$.

Figure 5 portrays the results corresponding to the complete synchronization by
varying the coupling strengths in two-dimensional parameter planes in terms of
synchronization error $E$. In the absence of higher-order coupling
$(\epsilon_{SI}=0)$, the synchronization and de-synchronization region is
depicted in $(\epsilon_{1},\epsilon_{2})$ parameter plane in Fig. 5(a). As
perceived, increasing $\epsilon_{2}$ reduces the critical values of
$\epsilon_{1}$ to achieve the synchrony. By introducing
$\epsilon_{SI}=0.00005$, a significant enhancement in synchronization region
is obtained in Fig. 5(b). Thereafter, we investigate the complete synchrony in
$(\epsilon_{1},\epsilon_{SI})$ and $(\epsilon_{2},\epsilon_{SI})$ parameter
planes with fixed values of $\epsilon_{2}=0.005$ and $\epsilon_{1}=0.05$,
respectively. Figures 5(c) and 5(d) depict the corresponding results in terms
of $E$. It is certainly observable that the critical values of $\epsilon_{1}$
and $\epsilon_{2}$ decrease with increasing higher-order coupling strength
$\epsilon_{SI}$. Therefore, our results indicate the enhancement in complete
synchrony in the simplicial complex of Rössler oscillators with the inclusion
of higher-order nonlinear coupling mechanism. We validate the numerical
results with the theoretical conjecture obtained from Eq. (14) in terms of the
maximum Lyapunov exponent, similarly to the previous example of the neuronal
system. Thus, the decoupled master stability equation can be written as,

$\begin{array}[]{l}\delta\dot{x}_{i}=-\delta y_{i}-\delta
z_{i}+\epsilon_{1}\bigg{[}k_{1}\big{\\{}2(\alpha-
x_{0})(x_{0}-\beta)-(x_{0}-\beta)^{2}\big{\\}}\\\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-2\gamma^{[1]}_{i}(\alpha-x_{0})(x_{0}-\beta)\bigg{]}\delta x_{i},\\\
\delta\dot{y}_{i}=\delta x_{i}+a\delta
y_{i}+\epsilon_{2}\bigg{[}k_{2}\big{\\{}2(\alpha-
y_{0})(y_{0}-\beta)-(y_{0}-\beta)^{2}\big{\\}}\\\ -2\gamma^{[2]}_{i}(\alpha-
y_{0})(y_{0}-\beta)\bigg{]}\delta
y_{i}+\epsilon_{SI}\bigg{[}2k^{[SI]}\\{2(\alpha-y_{0})(y_{0}-\beta)\\\
-(y_{0}-\beta)^{2}\\}-4\gamma^{[SI]}_{i}(\alpha-
y_{0})(y_{0}-\beta)\bigg{]}\delta y_{i},\\\ \\\ \delta\dot{z}_{i}=z_{0}\delta
x_{i}+(x_{0}-c)\delta z_{i},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ i=2,3,\cdots,N,\end{array}$ (22)

where the eigenvalues of the Laplacian $\mathscr{L}^{[1]}$ are symbolized by
$\gamma^{[1]}_{1}(=0)<\gamma^{[1]}_{2}\leq\cdots\leq\gamma^{[1]}_{N}$.
$\gamma^{[2]}_{i}=N$ and $\gamma^{[SI]}_{i}=N(N-2)$, $i=2,3,\cdots,N$ are the
eigenvalues of the Laplacians $\mathscr{L}^{[2]}$ and $\mathscr{L}^{[SI]}$,
respectively. The number of links adjacent to each node in second layer is
given by $k_{2}=(N-1)$, and each node is involved in
$k^{[SI]}=\begin{pmatrix}N-1\\\ 2\end{pmatrix}$ numbers of $2$-simplices. The
solid black curves imposed in Figs. 5(a)-5(d) represent analytical
synchronization thresholds obtained from Eq. (22), solved for zero value of
maximum Lyapunov exponent, which affirms that the numerical findings accord
well with our theoretical predictions.

##  V Conclusion

To conclude, here we have developed the most generic model that accounts for
higher-order interactions between dynamical units in simplicial complexes with
multiple interaction layers and studied the complete synchronization
phenomena. In the presence of arbitrary sort of coupling schemes, we derive
the invariance condition for the synchronization state. Under the requirements
of invariance condition, we determine the necessary condition for the stable
synchronization state, which certainly generalizes the well-known MSF approach
Pecora and Carroll (1998) to multilayer structures of simplicial complexes.
The intricacy of multiple layers and arbitrary higher-order coupling in the
considered system is reflected by the fact that our formalism yields a set of
coupled linear differential equations (Eq. (13b)) instead of a single
parametric variational equation. However, we have shown that, for a certain
instance, our formalism gives a set of uncoupled parametric variational
equations (Eq. (14)) with dimensions equal to the dimension of a single
dynamical unit. Finally, a set of numerical results have been added to our
theoretical derivations, confirming the legitimacy and applicability of the
methodology.

Our investigation is based on the most prevalent and well-studied
synchronization behavior, namely, the complete synchronization phenomenon.
However, many additional types of synchrony occurs in the systems with
multiple interaction layers such as chimeras Majhi  _et al._ (2016), cluster
synchronization Della Rossa  _et al._ (2020), relay synchronization Anwar  _et
al._ (2021b), etc. All of these synchronization states have been investigated
in multilayer systems with solely pairwise interactions. The formation of such
states, or indeed fresh ones, in multilayered systems with higher-order
interactions, along with their stability, are both fascinating challenges that
provide way to future investigation. Apart from that, our work can be extended
to higher-order structures with multiple interaction layers having directed
connection topology, which has very recently been investigated for mono-layer
systems with higher-order interactions Gallo  _et al._ (2022). Furthermore, we
have focused only on the situation of clique complexes, it is still unclear
how our numerical conclusions (e.g., enhancement of synchronous region) would
alter if we took into account other higher-order structures, such as
hypergraphs, and more generalized simplicial complexes.

## References

  * Boccaletti  _et al._ (2006) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Physics Reports 424, 175 (2006).
  * Newman (2003) M. E. Newman, SIAM Review 45, 167 (2003).
  * Watts and Strogatz (1998) D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998).
  * Dorogovtsev  _et al._ (2008) S. N. Dorogovtsev, A. V. Goltsev, and J. F. Mendes, Reviews of Modern Physics 80, 1275 (2008).
  * Newman (2018) M. Newman,  _Networks_ (Oxford University Press, 2018).
  * Latora  _et al._ (2017) V. Latora, V. Nicosia, and G. Russo,  _Complex Networks: Principles, Methods and Applications_ (Cambridge University Press, 2017).
  * Battiston  _et al._ (2020) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Physics Reports 874, 1 (2020).
  * Majhi  _et al._ (2022) S. Majhi, M. Perc, and D. Ghosh, Journal of the Royal Society Interface 19, 20220043 (2022).
  * Sizemore  _et al._ (2018) A. E. Sizemore, C. Giusti, A. Kahn, J. M. Vettel, R. F. Betzel, and D. S. Bassett, Journal of Computational Neuroscience 44, 115 (2018).
  * Petri  _et al._ (2014) G. Petri, P. Expert, F. Turkheimer, R. Carhart-Harris, D. Nutt, P. J. Hellyer, and F. Vaccarino, Journal of The Royal Society Interface 11, 20140873 (2014).
  * Grilli  _et al._ (2017) J. Grilli, G. Barabás, M. J. Michalska-Smith, and S. Allesina, Nature 548, 210 (2017).
  * Billick and Case (1994) I. Billick and T. J. Case, Ecology 75, 1529 (1994).
  * Patania  _et al._ (2017) A. Patania, G. Petri, and F. Vaccarino, EPJ Data Science 6, 1 (2017).
  * Vasilyeva  _et al._ (2021) E. Vasilyeva, A. Kozlov, K. Alfaro-Bittner, D. Musatov, A. Raigorodskii, M. Perc, and S. Boccaletti, Scientific Reports 11, 1 (2021).
  * Estrada and Ross (2018) E. Estrada and G. J. Ross, Journal of Theoretical Biology 438, 46 (2018).
  * Benson  _et al._ (2016) A. R. Benson, D. F. Gleich, and J. Leskovec, Science 353, 163 (2016).
  * Carletti  _et al._ (2020) T. Carletti, F. Battiston, G. Cencetti, and D. Fanelli, Physical Review E 101, 022308 (2020).
  * Neuhäuser  _et al._ (2021) L. Neuhäuser, R. Lambiotte, and M. T. Schaub, Physical Review E 104, 064305 (2021).
  * Bick  _et al._ (2021) C. Bick, E. Gross, H. A. Harrington, and M. T. Schaub, arXiv preprint arXiv:2104.11329 (2021).
  * Battiston  _et al._ (2021) F. Battiston, E. Amico, A. Barrat, G. Bianconi, G. Ferraz de Arruda, B. Franceschiello, I. Iacopini, S. Kéfi, V. Latora, Y. Moreno, _et al._ , Nature Physics 17, 1093 (2021).
  * Torres and Bianconi (2020) J. J. Torres and G. Bianconi, Journal of Physics: Complexity 1, 015002 (2020).
  * Giusti  _et al._ (2016) C. Giusti, R. Ghrist, and D. S. Bassett, Journal of Computational Neuroscience 41, 1 (2016).
  * Aleksandrov (1998) P. S. Aleksandrov,  _Combinatorial Topology_ , Vol. 1 (Courier Corporation, 1998).
  * Berge (1973) C. Berge,  _Graphs and Hypergraphs_ (North-Holland Pub. Co., 1973).
  * Carlsson (2009) G. Carlsson, Bulletin of the American Mathematical Society 46, 255 (2009).
  * Iacopini  _et al._ (2019) I. Iacopini, G. Petri, A. Barrat, and V. Latora, Nature Communications 10, 1 (2019).
  * Alvarez-Rodriguez  _et al._ (2021) U. Alvarez-Rodriguez, F. Battiston, G. F. de Arruda, Y. Moreno, M. Perc, and V. Latora, Nature Human Behaviour 5, 586 (2021).
  * Skardal and Arenas (2019) P. S. Skardal and A. Arenas, Physical Review Letters 122, 248301 (2019).
  * Petri and Barrat (2018) G. Petri and A. Barrat, Physical Review Letters 121, 228301 (2018).
  * Pikovsky  _et al._ (2001) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A Universal Concept in Nonlinear Sciences_ , Cambridge Nonlinear Science Series (Cambridge University Press, 2001).
  * Boccaletti  _et al._ (2002) S. Boccaletti, J. Kurths, G. Osipov, D. Valladares, and C. Zhou, Physics Reports 366, 1 (2002).
  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Physics Reports 469, 93 (2008).
  * Courtney and Bianconi (2016) O. T. Courtney and G. Bianconi, Physical Review E 93, 062311 (2016).
  * Ghorbanchian  _et al._ (2021) R. Ghorbanchian, J. G. Restrepo, J. J. Torres, and G. Bianconi, Communications Physics 4, 1 (2021).
  * Skardal and Arenas (2020) P. S. Skardal and A. Arenas, Communications Physics 3, 1 (2020).
  * Kuehn and Bick (2021) C. Kuehn and C. Bick, Science Advances 7, eabe3824 (2021).
  * Kachhvah and Jalan (2022a) A. D. Kachhvah and S. Jalan, New Journal of Physics 24, 052002 (2022a).
  * Skardal  _et al._ (2021) P. S. Skardal, L. Arola-Fernández, D. Taylor, and A. Arenas, Physical Review Research 3, 043193 (2021).
  * Zhang  _et al._ (2022) Y. Zhang, M. Lucas, and F. Battiston, arXiv preprint arXiv:2203.03060 (2022).
  * Xu  _et al._ (2020) C. Xu, X. Wang, and P. S. Skardal, Physical Review Research 2, 023281 (2020).
  * Zhang  _et al._ (2021) Y. Zhang, V. Latora, and A. E. Motter, Communications Physics 4, 1 (2021).
  * Kachhvah and Jalan (2022b) A. D. Kachhvah and S. Jalan, Physical Review E 105, L062203 (2022b).
  * Kundu and Ghosh (2022) S. Kundu and D. Ghosh, Phys. Rev. E 105, L042202 (2022).
  * Millán  _et al._ (2020) A. P. Millán, J. J. Torres, and G. Bianconi, Physical Review Letters 124, 218301 (2020).
  * Arnaudon  _et al._ (2021) A. Arnaudon, R. L. Peach, G. Petri, and P. Expert, arXiv preprint arXiv:2111.11073 (2021).
  * Gambuzza  _et al._ (2021) L. Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora, and S. Boccaletti, Nature Communications 12, 1 (2021).
  * Lucas  _et al._ (2020) M. Lucas, G. Cencetti, and F. Battiston, Physical Review Research 2, 033410 (2020).
  * Gong and Pikovsky (2019) C. C. Gong and A. Pikovsky, Physical Review E 100, 062210 (2019).
  * Sorrentino (2012) F. Sorrentino, New Journal of Physics 14, 033035 (2012).
  * Del Genio  _et al._ (2016) C. I. Del Genio, J. Gómez-Gardeñes, I. Bonamassa, and S. Boccaletti, Science Advances 2, e1601679 (2016).
  * Boccaletti  _et al._ (2014) S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardenes, M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin, Physics Reports 544, 1 (2014).
  * Cardillo  _et al._ (2013) A. Cardillo, J. Gómez-Gardenes, M. Zanin, M. Romance, D. Papo, F. Del Pozo, and S. Boccaletti, Scientific Reports 3, 1 (2013).
  * Szell  _et al._ (2010) M. Szell, R. Lambiotte, and S. Thurner, Proceedings of the National Academy of Sciences 107, 13636 (2010).
  * Rakshit  _et al._ (2018) S. Rakshit, B. K. Bera, and D. Ghosh, Physical Review E 98, 032305 (2018).
  * Anwar  _et al._ (2022) M. S. Anwar, S. Rakshit, D. Ghosh, and E. M. Bollt, Phys. Rev. E 105, 024303 (2022).
  * Anwar  _et al._ (2021a) M. S. Anwar, S. Kundu, and D. Ghosh, Chaos, Solitons & Fractals 142, 110476 (2021a).
  * Bianconi and Dorogovtsev (2014) G. Bianconi and S. N. Dorogovtsev, Physical Review E 89, 062814 (2014).
  * Gomez  _et al._ (2013) S. Gomez, A. Diaz-Guilera, J. Gomez-Gardenes, C. J. Perez-Vicente, Y. Moreno, and A. Arenas, Physical Review Letters 110, 028701 (2013).
  * Granell  _et al._ (2013) C. Granell, S. Gómez, and A. Arenas, Physical Review Letters 111, 128701 (2013).
  * Gómez-Gardenes  _et al._ (2012) J. Gómez-Gardenes, I. Reinares, A. Arenas, and L. M. Floría, Scientific Reports 2, 1 (2012).
  * Parastesh  _et al._ (2022) F. Parastesh, M. Mehrabbeik, K. Rajagopal, S. Jafari, and M. Perc, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 013125 (2022).
  * Ince  _et al._ (2009) R. A. Ince, F. Montani, E. Arabzadeh, M. E. Diamond, and S. Panzeri, in  _Journal of Physics: Conference Series_ , Vol. 197 (IOP Publishing, 2009) p. 012013.
  * Tlaie  _et al._ (2019) A. Tlaie, I. Leyva, and I. Sendiña-Nadal, Physical Review E 100, 052305 (2019).
  * Amari  _et al._ (2003) S.-i. Amari, H. Nakahara, S. Wu, and Y. Sakai, Neural Computation 15, 127 (2003).
  * Sun and Bianconi (2021) H. Sun and G. Bianconi, Physical Review E 104, 034306 (2021).
  * Anwar and Ghosh (2022) M. S. Anwar and D. Ghosh, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 033125 (2022).
  * Jalan and Suman (2022) S. Jalan and A. Suman, arXiv preprint arXiv:2206.10852 (2022).
  * Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, Physical Review Letters 80, 2109 (1998).
  * Anwar  _et al._ (2021b) M. S. Anwar, D. Ghosh, and N. Frolov, Mathematics 9, 2135 (2021b).
  * Gallo  _et al._ (2022) L. Gallo, R. Muolo, L. V. Gambuzza, V. Latora, M. Frasca, and T. Carletti, arXiv preprint arXiv:2202.08707 (2022).
  * Jalil  _et al._ (2012) S. Jalil, I. Belykh, and A. Shilnikov, Physical Review E 85, 036214 (2012).
  * Sherman (1994) A. Sherman, Bulletin of Mathematical Biology 56, 811 (1994).
  * Rössler (1976) O. E. Rössler, Physics Letters A 57, 397 (1976).
  * Fellin  _et al._ (2004) T. Fellin, O. Pascual, S. Gobbo, T. Pozzan, P. G. Haydon, and G. Carmignoto, Neuron 43, 729 (2004).
  * Allen and Barres (2009) N. J. Allen and B. A. Barres, Nature 457, 675 (2009).
  * Hormuzdi  _et al._ (2004) S. G. Hormuzdi, M. A. Filippov, G. Mitropoulou, H. Monyer, and R. Bruzzone, Biochimica et Biophysica Acta (BBA)-Biomembranes 1662, 113 (2004).
  * Kandel  _et al._ (2000) E. R. Kandel, J. H. Schwartz, T. M. Jessell, S. Siegelbaum, A. J. Hudspeth, S. Mack, _et al._ ,  _Principles of Neural Science_ , Vol. 4 (McGraw-hill New York, 2000).
  * Pereda (2014) A. E. Pereda, Nature Reviews Neuroscience 15, 250 (2014).
  * Muldoon  _et al._ (2016) S. F. Muldoon, E. W. Bridgeford, and D. S. Bassett, Scientific Reports 6, 1 (2016).
  * Majhi  _et al._ (2016) S. Majhi, M. Perc, and D. Ghosh, Scientific reports 6, 1 (2016).
  * Della Rossa  _et al._ (2020) F. Della Rossa, L. Pecora, K. Blaha, A. Shirin, I. Klickstein, and F. Sorrentino, Nature communications 11, 1 (2020).