Synchronizability in randomized weighted simplicial complexes

Complex systems arise when a large number of dynamical units interact with each other, giving rise to emergent properties distinct from those of the individual subsystems [1]. Such systems are ubiquitous, found naturally in entities like the brain and man-made ones like the internet or financial markets [2]. The study and modelling of these systems have been a longstanding focus of researchers [3, 4, 5]. Complex networks represent the dynamical units, and their interactions as nodes and links are a common way to model such systems. However, these networks typically only account for pairwise interactions, whereas many complex systems involve interactions among three or more units. To gain a more comprehensive understanding of these systems, one must incorporate higher-order networks, such as hypergraphs and simplicial complexes, to capture these higher-order interactions [6].

Simplicial complexes are used to model complex systems with higher-order interactions [4, 7, 8, 9]. A simplicial complex consists of $d$-simplices, where $d$ is the dimension of the simplex. A $0$-dimensional simplex is a node, a $1$-dimensional simplex is a link, a $2$-dimensional simplex is a triangle, and so on. A $d$-dimensional simplicial complex consists of simplices up to dimension $d$. They provide a good representation of complex systems with higher-order interactions [10, 11, 7, 12]. In most studies, these simplicial complexes are unweighted, meaning all the links and triangles have the same weight [13, 14]. However, it has been observed that in most real-world systems, this is not feasible. For instance, in collaboration networks, if there is a three-author paper, not every pair of authors has to produce research papers by themselves. But in unweighted simplicial complexes, if there is a three-author paper, then there is also the presence of all the possible two-author research papers, and this leads to information loss, which one could avoid by introducing proper weights for links [15]. Courtney and Bianconi [16] have developed a model to distribute weights to links in terms of bare affinity weights and topological weights. Every simplex has bare-affinity weights, while the topological weights represent whether a particular link has a contribution other than being part of a triangle. Therefore, weighted simplicial complexes provide a more realistic representation of higher-order networks.

After accurately modelling complex systems, we can better understand their emergent properties, such as synchronization. As a result, we can control the synchronization that occurs in natural or artificial systems. Synchronization occurs when individual dynamical systems adjust their properties to have common dynamics [3]. Neuronal synchronization can cause epilepsy, while synchronization from a healthy to an infected state can cause an epidemic outbreak, and so on. Synchronization can be both constructive and destructive. However, it is possible to promote synchronization in desired systems and inhibit it where it leads to destruction. This can be achieved by determining under what circumstances the system goes into synchronization and desynchronization and whether synchronization is possible. Two main parameters are used to calculate the synchronizability of the system: eigenratio and the cost of the connections of nodes in the system. The lower the eigenratio and cost, the system has more synchronizability [17].

To calculate the eigenratio and costs, we need to construct the connectivity matrices, such as adjacency matrices and Laplacian matrices, and then find the eigenvalues of these matrices. Writing connectivity matrices for large, complex systems, such as connection networks, collaboration networks, and neuronal networks, can be challenging, and collecting the necessary data is not an easy task. It becomes even more difficult when we include higher-order interactions, as we need to construct the adjacency tensors for these interactions. Even for three-body interactions, the adjacency tensor contains $N\times N\times N$ elements. In 2006, Zhou et al. [17] developed an approach to calculate eigenratios and costs based on the heterogeneity of the system intensities. This approach significantly reduced the computational cost of constructing connectivity matrices and evaluating their eigenvalues. We extend this approach to randomized weighted simplicial complexes. As a result, we can determine the synchronizability parameters solely from the heterogeneity of the intensities, specifically the maximum and minimum intensities of a node as part of links ($d=1$), triangles ($d=2$), and so on and from their coupling strengths. In this paper, we discuss the mathematical formulation, derive general formulas for determining the synchronization without explicitly computing the eigenvalues of huge matrices and corroborate with numerical results. We validate the results by applying the general formula to Rössler oscillators and real-world connection networks. We check the emergence of complete synchronization in diffusively coupled identical Rössler oscillators, and we show that the results are independent of the number of oscillators in the networks, the components through which they are connected and also the distributions of the degrees and intensities.

In order to derive the general formulas for the eigenratio and the cost, we first write the dynamical equations for the weighted simplicial complex. Then, we deduce the corresponding variational equations and modify these equations in terms of the effective matrix $M$. The eigenvalues of the effective matrix determine the stability of the synchronized state. Since the effective matrix is a zero row sum matrix, the first eigenvalue will be zero, corresponding to the mode along the synchronization manifold. We can find the synchronizability of the system from the ratio between $\lambda_{N}$ and $\lambda_{2}$. For simplicity, we consider the simplicial complex of dimension $2$, which we can extend to any dimension.

We consider a simplicial complex of identically coupled oscillators, described by the following equations,

where $\mathbf{x}_{i}$ is the state vector of $i$-th oscillator of dimension $m$, and $\mathbf{f}(\mathbf{x}_{i})$ represents the dynamics of the uncoupled oscillators. ${\omega_{ij}}^{(1)}$ and $\omega_{ijk}^{(2)}$ are the topological weights of links and triangles, respectively. $a_{ij}^{(1)}$ and $a_{ijk}^{(2)}$ are the elements of the adjacency matrix $A^{(1)}$ and adjacency tensor $A^{(2)}$. $a_{ij}^{(1)}=1$ if the nodes $i$ and $j$ form a link, or $a_{ij}^{(1)}=0$ otherwise. Also $a_{ijk}^{(2)}=1$ if the nodes $i$, $j$, and $k$ form a triangle, or $a_{ijk}^{(2)}=0$ otherwise. Here $\sigma_{1}$ and $\sigma_{2}$ are the coupling strengths of the links (pairwise) and triangles (non-pairwise), respectively. Here, $\mathbf{g}^{(1)}(\mathbf{x}_{i},\mathbf{x}_{j})$ and $\mathbf{g}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k})$ are the synchronization noninvasive functions [7] and are chosen in the forms $\mathbf{g}^{(1)}(\mathbf{x}_{i},\mathbf{x}_{j})=\mathbf{h}^{(1)}(\mathbf{x}_{j})-\mathbf{h}^{(1)}(\mathbf{x}_{i})$ and $\mathbf{g}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k})=\mathbf{h}^{(2)}(\mathbf{x}_{j},\mathbf{x}_{k})-\mathbf{h}^{(2)}(\mathbf{x}_{i},\mathbf{x}_{i})$. At the synchronization state, they tend to zero, making the dynamics of the system resemble that of the uncoupled oscillators.

For the above diffusive coupling functions, we can rewrite the above equation as

where $\mathbf{h}^{(1)}(\mathbf{x}_{j})$ and $\mathbf{h}^{(2)}(\mathbf{x}_{j},\mathbf{x}_{k})$ are coupling functions that couples nodes in links and triangles, respectively. The choice of coupling functions affects the synchronizability of the system [7]. So we have to choose the coupling functions in a way that the oscillators tend to synchronize and we can write,

Here $S_{i}^{(1)}$ and $S_{i}^{(2)}$ are the intensities of node $i$ to form links and triangles, i.e., the number of weighted links and weighted triangles incident on node $i$. For randomized simplicial complexes with $K^{(1)}_{\mbox{min}}\gg 1$, we can use mean-field approximation and this equation can be rewritten as,

Here, $K_{i}^{(1)}$ is the degree of pairwise interactions, i.e., the total number of links that a node participates in. $K_{i}^{(2)}$ is the total number of triangles that a node participates in, also called the degree for three-body interactions. Now, we consider pairwise and non-pairwise local mean fields as

due to the interaction between the nodes through links and

due to the interaction between nodes through triangles. Then, Eq. (3) can be written in terms of mean field as,

So, the interaction between nodes in a large randomized weighted simplicial complex can be approximated as the interaction of a node with a mean field. This mean field is not only generated by the interaction of links, as in the case of complex networks, but also by the interaction of nodes as triangles.

According to the condition for natural coupling, close to the synchronized state, the interaction between nodes in links and triangles will be similar since, upon synchronization, the two nodes that are in the same state can be considered as a single node, and the three-body interactions reduce to the two-body interactions [7]. The above assumptions enable us to write $\mathbf{h}^{(2)}(\mathbf{x},\mathbf{x})=\mathbf{h}^{(1)}(\mathbf{x})$ and $\overline{\mathbf{H}}^{(2)}{(\mathbf{x},\mathbf{x})}=\overline{\mathbf{H}}^{(1)}{(\mathbf{x})}=\mathbf{H(x)}$. Equation (4) can be rewritten as,

The variational equation of the above Eq. (5) is of the form

where $\mathbf{x}^{s}$ is the synchronization state. The eigenvalues are approximately equal to $\sigma_{1}S_{i}^{(1)}+\sigma_{2}S_{i}^{(2)}$, $i=1,2,\ldots,N$. From the above equation, we can write the eigenratio [17] as,

where $\sigma_{1}S_{\mbox{max}}^{(1)}+\sigma_{2}S_{\mbox{max}}^{(2)}\geq\mbox{max}\{\sigma_{1}S_{1}^{(1)}+\sigma_{2}S_{1}^{(2)},\sigma_{1}S_{2}^{(1)}+\sigma_{2}S_{2}^{(2)},\cdots,\sigma_{1}S_{N}^{(1)}+\sigma_{2}S_{N}^{(2)}\}$ will give the maximum bound for eigenvalues and similar for minimum $\sigma_{1}S_{\mbox{min}}^{(1)}+\sigma_{2}S_{\mbox{min}}^{(2)}\leq\mbox{min}\{\sigma_{1}S_{1}^{(1)}+\sigma_{2}S_{1}^{(2)},\sigma_{1}S_{2}^{(1)}+\sigma_{2}S_{2}^{(2)},\cdots,\sigma_{1}S_{N}^{(1)}+\sigma_{2}S_{N}^{(2)}\}$. These maximum and minimum bounds may not necessarily belong to any particular node index $i$. It further simplifies the problem as it is sufficient to consider only the maximum and minimum values of the intensities to calculate the eigenratio and it eliminates the necessity to calculate all the eigenvalues. It would be of great advantage when $N$ is very large.

For the case of pairwise interaction only ($\sigma_{2}=0.0$), the eigenratio $R$ is very similar to the previous result [17]. For non-zero pairwise coupling strength ($\sigma_{1}\neq 0.0,$), we can rewrite the above equation as,

In the similar way, we can derive the eigenratio for the $d$-dimensional simplical complex as,

Thus, the eigenratio of a randomized simplicial complex depends on the coupling strengths (pairwise and non-pairwise), and maximum and minimum intensities of the nodes. If we assume all the non-pairwise coupling strengths are identical to the pairwise coupling strength (i.e., $\sigma_{2}=\sigma_{3}=\cdots=\sigma_{d}=\sigma_{1}$), then the eigenratio $R$ is independent on the coupling strengths and depends on the maximum and minimum intensities as observed for pairwise weighted random networks [17].

To find the tight bounds of the above mean-field approximation, we can write Eq. (5) as,

where $G_{ij}^{(1)}$ and $G_{ij}^{(2)}$ are the elements of matrices $G^{(1)}=S^{(1)}D^{(1)-1}L^{(1)}$ and $G^{(2)}=S^{(2)}D^{(2)-1}L^{(2)}$, respectively. Here $L^{(1)}$ and $L^{(2)}$ are generalized Laplacian matrices, while $S$ and $D$ are diagonal matrices of strengths and generalized degrees, respectively. It is a good point that we have separated weights from topology. We can write normalized Laplacian matrices as $\overline{L}^{(1)}=D^{-1}L^{(1)}$ and $\overline{L}^{(2)}=D^{-1}L^{(2)}$, then $G^{(1)}=S^{(1)}\overline{L}^{(1)}$ and $G^{(2)}=S^{(2)}\overline{L}^{(2)}$. Thus, the largest and smallest eigenvalues of the matrices $G^{(1)}$ and $G^{(2)}$ are bounded by the eigenvalues of $\overline{L}^{(1)}$ and $\overline{L}^{(2)}$. We can write the linearized variational Eq. (10) by following the procedure outlined in [18] and is given by

We can define an effective matrix $M$ as,

Equation (11) can be rewritten, in terms of effective matrix $M$, as

The eigenvalues of the matrix $M$ depend on the ratio of coupling strengths, degrees and intensities, and the eigenratio of this matrix can be expressed as $R=\lambda_{N}/\lambda_{2},$ where $\lambda_{2}$ and $\lambda_{N}$ are the 2nd and $N$-th eigenvalues of $M$ such that $0=\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\leq\cdots\leq\lambda_{N}$.

Another measure of synchronizability is the cost $C$ involved in the coupling of nodes in a simplicial complex, and it is the total strength of connections of all nodes. It can be written as,

By using the eigenvalues of $M$ in Eq. (13), the normalized cost can be written as,

where $\alpha_{1}=\sigma_{1}\lambda_{2}(M)$ and mean intensity

From Eq. (6), the normalized cost will be,

where $r={\sigma_{2}}/{\sigma_{1}}$. We can verify Eq. (8) and Eq. (17) numerically by computing the eigenvalues of the effective matrix $M$, and we can obtain an explicit bound from the eigenvalues of $\overline{L}^{(1)}$ and $\overline{L}^{(2)}$. Using this bound, we can derive formulas for $R$ and $C_{0}$ as a function of $S_{\mbox{max}}^{(1)}$, $S_{\mbox{min}}^{(1)}$ and the generalized degrees.

We can demonstrate that the upper and lower bounds of the nonzero eigenvalues of the effective matrix $M$ are given by the eigenvalues $\mu^{(1)}_{l}$ and $\mu^{(2)}_{l}$ of the matrices $G^{(1)}$ and $G^{(2)}$, respectively, as follows,

For sufficiently random simplicial complexes, the spectra of the matrices $G^{(1)}$ and $G^{(2)}$ tend to follow the semicircle law [17] and we can write $\mu_{2}^{(1)}\approx 1-2/\sqrt{K^{(1)}}$, $\mu_{N}^{(1)}\approx 1+2/\sqrt{K^{(1)}}$, $\mu_{2}^{(2)}\approx 1-2/\sqrt{K^{(2)}}$ and $\mu_{N}^{(2)}\approx 1+2/\sqrt{K^{(2)}}$, provided that $K_{\mbox{min}}^{(d)}\gg\sqrt{K^{(d)}}$. Here $K^{(d)}$ is the mean degree of the $d$-simplex. By introducing these in Eq. (18) and Eq. (19), we get the following approximations for the bounds of $R$ and $C_{0}$ as

For a given value of $K$, the synchronizability of randomly weighted simplicial complexes with a large value of $K_{\mbox{min}}^{(1)}$ is expected to follow a general formula as follows,

and

where $A_{R}=\frac{1+2/\sqrt{K^{(1)}}}{1-2/\sqrt{K^{(1)}}}$ and $A_{C}=\frac{1}{1-2/\sqrt{K^{(1)}}}$. Here, $A_{R}\rightarrow 1$ and $A_{C}\rightarrow 1$ in the limit $K^{(1)}\rightarrow\infty$. Also, $B_{R_{1}}=r\frac{1+2/\sqrt{K^{(2)}}}{1+2/\sqrt{K^{(1)}}}$ and $B_{R_{2}}=r\frac{1-2/\sqrt{K^{(2)}}}{1-2/\sqrt{K^{(1)}}}$. $B_{R_{1}}\rightarrow r$ and $B_{R_{2}}\rightarrow r$ in the limit $K^{(d)}\rightarrow\infty(d=1,2)$. They are well agreed with Eqs. (8) and (17).

For $N$-dimensional simplicial complex,

and

where

In this paper, we have analyzed the nature of synchronization in randomized weighted simplicial complexes of identical oscillators with natural coupling. We derived general formulas for the eigenratio and the cost and verified these theoretical results for various network topologies with diffusive couplings and intensity distributions. Our observations indicated that synchronization is independent of the structure of simplicial complexes, the way the oscillators interact with one another and the total number of oscillators present in the system. We observed complete synchronization for fixed values of $\sigma_{1}$, $\sigma_{2}$, $S^{(1)}_{\mbox{max}}$, $S^{(2)}_{\mbox{max}}$, and $S^{(2)}_{\mbox{min}}$, when $S^{(1)}_{\mbox{min}}=0.2\times S^{(1)}_{\mbox{max}}$ (i.e., $20$ percentage of $S^{(1)}_{\mbox{max}}$). Notably, the synchronization error values for $N=150$ and $N=200$ are nearly identical, indicating that synchronisation is independent of the number of oscillators when $N$ is large. Our general formula provides a good approximation for complete synchronization in these cases.

Furthermore, our observations highlighted that the effect of coupling strengths varies depending on the largest intensity. When $S^{(1)}_{\mbox{max}}>S^{(2)}_{\mbox{max}}$, an increase in $\sigma_{2}$ enhances synchronizability. Conversely, when $S^{(2)}_{\mbox{max}}>S^{(1)}_{\mbox{max}}$, synchronization decreases with an increase in $\sigma_{2}$. The behaviours of constants $A_{R}$, $A_{C}$, $B_{R_{1}}$, and $B_{R_{2}}$ suggest that the approximated formula, containing only coupling strengths and intensities, can be used.

We have validated our results by applying the general formula to determine the synchronizability of weighted simplicial complexes of diffusively coupled identical Rössler oscillators in Sec. III.4. We further verified our results using real-world connectivity structure where the Rössler oscillators are placed on each node in Sec. III.5. From numerical simulations in Sec. III.4 and III.5, we see that the results are robust on the number of nodes, components through which they are coupled and the distributions of degree and intensities.

In conclusion, we demonstrated that the synchronizability of large randomized weighted simplicial complexes with natural coupling can be determined using mean degrees, coupling strengths, and intensities without the need to construct connectivity matrices or explicitly calculate their eigenvalues. These formulas substantially reduce computation time and cost when assessing a system’s tendency to synchronize. In essence, these formulas enable the prediction and control of synchronizability in complex networks with higher-order interactions, achieved by manipulating the degrees, weights, and coupling strengths.