Stability of synchronization in simplicial complexes with multiple interaction layers

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.

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]}$.

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.