Synchronizability of Multi-Layer-Coupled Star-Composed Networks

Abstract

In this paper, several multi-layer-coupled star-composed networks with similar symmetrical structures are defined by using the theory of graph operation. The supra-Laplacian matrix of the corresponding multi-layer networks is obtained according to the master stability equation (MSF). Two important indexes that reflect the synchronizability of these kinds of networks are derived in the case of bounded and unbounded synchronized regions. The relationships among the synchronizability, the number of layers, the length of the paths, the branchings, and the interlayer and intralayer coupling strengths in the two cases are studied. At the same time, the simulation experiments are carried out with the MATLAB software, and the simulated images of the two symmetrical structure networks’ synchronizability are compared. Finally, the factors affecting the synchronizability of multi-layer-coupled star-composed networks are found. On this basis, optimization schemes are given to improve the synchronizability of multi-layer-coupled star-composed networks and the influences of the number of central nodes on the networks’ synchronizability are further studied.

1. Introduction

Complex networks have appeared in many aspects of our lives; they are used in the prevention of the spread of epidemic diseases, computer networks, social networks, and other fields [1,2,3,4]. These networks are an important interdisciplinary research field that has attracted the interest of many scholars. In recent years, many new advances have been made in the research of complex networks, such as the spread and stability of complex networks [5,6], the synchronization and consensus of multi-layer complex networks [7,8,9,10,11,12,13,14,15], diffusion in multi-layer networks [16,17,18,19,20], and the modeling and robustness of multi-layer complex networks [21,22,23,24,25].

Common network synchronization types are complete synchronization [8] and phase synchronization [26]. This paper studies the complete synchronization of networks according to the master stability equation (MSF) proposed by Pecora and Carroll [27]. Research on the synchronizability of multi-layer networks mostly focuses on structures such as star-like networks, and there has been little research on the synchronizability of multi-layer star-composed networks. On the one hand, the symmetrical structures of multi-layer-coupled star-composed networks are more complex and the calculation of the synchronizability indexes in these symmetrical structures is more difficult than the calculation of the indexes of other types of structure. On the other hand, the characteristic polynomials of many star-composed coupled networks can not be solved analytically.

In recent years, scholars have obtained many classical results on the synchronization of star-like networks. Aguirre et al. [28] studied how the degree of nodes connected between two layers of networks affects the synchronizability of the whole system and proposed an effective synchronization strategy. Zhu et al. [29] studied the relationships among the synchronizability of two k-layer variable-coupling windmill-type networks, the number of nodes, the coupling strength, and other parameters in the case of unbounded and bounded synchronized regions and provided an optimization scheme to improve the synchronizability. Deng et al. [30] gave the characteristic polynomials of multi-layer-coupled chain networks, discussed the relationships between synchronizability and different parameters in the case of unbounded and bounded synchronized regions, and further provided image examples to verify the effectiveness of their theoretical analysis. Li et al. [31] studied the bipartite synchronization of coupled inertial memory neural networks with cooperative and competitive interaction. Two simulation examples were given to verify the effectiveness of the proposed discontinuous control strategy for bipartite synchronization. Inspired by the above works of literature, the main contents of this paper are as follows:

• We provide definitions of star-composed networks and introduce the relevant theoretical knowledge required for the synchronization of multi-layer networks.
• The supra-Laplacian matrix of multi-layer-coupled star-composed networks is given according to the MSF. With the help of graph theory and graph operation, the important indexes reflecting the synchronizability are obtained for the cases of bounded and unbounded synchronized regions.
• According to the important indexes of synchronizability, we obtain the relationships between synchronizability and the parameters of multi-layer-coupled star-composed networks.
• Through simulation experiments, we obtain the relationships between the synchronizability and various parameters, and provide a theoretical basis for optimizing the synchronizability of multi-layer-coupled star-composed networks.

The structure of this paper is as follows. In Section 2, the preliminaries required to study the synchronization of multi-layer-coupled star-composed networks are given. In Section 3, the eigenvalue spectrum and the synchronizability indexes of multi-layer-coupled star-composed networks are studied and extended to a more general case. In Section 4, the synchronizability of multi-layer-coupled star-composed networks is simulated numerically. In Section 5, the main conclusions of this paper are given.

2. Preliminaries

2.1. The Synchronizability of Multi-Layer Networks

According to the dynamic equation of the ith node of the Lth layer in the Q-layer network [32]:

$$\frac{d{x}_i^L}{dt}=f\left(x_i^L\right)-a_L\sum _{\genfrac{}{}{0pt}{}{j=1}{}}^N{\omega }_{ij}^L\mathsf{\Phi }\left(x_j^L\right)-d\sum _{\genfrac{}{}{0pt}{}{T=1}{}}^Q{d}_i^{LT}\mathrm{Y}\left(x_i^L\right),$$

(1)

$$i=1,2,\cdots ,N,L=1,2,\cdots ,Q.$$

where N is the number of nodes in each layer; Q is the number of layers; $x_i^L$ is the state variable of the ith node of the Lth layer; $f(*)$ is the dynamic function; $a_L$ is the intralayer coupling strength of the Lth layer; d is the interlayer coupling strength; $\mathsf{\Phi }$ is the intralayer coupling function; and $\mathrm{Y}$ is the interlayer coupling function. The property of the Laplacian matrix can ensure the block diagonalization of the Jacobian, which means that the MSF can be used to determine the stability of a synchronized solution by reducing the dimensionality of the synchronization problem [27], $H^L=\left({\omega }_{ij}^L\right)=S^L-W^L$ is the Laplacian matrix of the Lth layer. Additionally, $S^L=\left(s_{ii}^L\right)\in R^{N\times N}$ is the degree matrix of the Lth layer, where $s_{ii}^L=d\left(v_i\right)$, and $d\left(v_i\right)$ is the degree of the vertex $v_i$. $W^L=\left(w_{ij}^L\right)\in R^{N\times N}$ is the adjacency matrix of the Lth layer. If there is an edge between the vertex $v_i$ and $v_j$ of the Lth layer, $w_{ij}^L=1$; otherwise, $w_{ij}^L=0(i,j=1,2,3,\cdots ,N-1,N)$.

In this paper, we consider the coupling strength of multi-layer networks. The weighted Laplacian matrix of the Lth layer is $H^{\left(L\right)}=a_L(S^L-W^L)$. The interlayer-weighted Laplacian matrix is defined as $dD=d\left(d_i^{LT}\right)\in R^{Q\times Q}(L,T=1,2,3,\cdots ,Q)$. If there is an edge between the ith node of the Lth layer and the ith node of the Tth layer, $d_i^{LT}=-1$; otherwise, $d_i^{LT}=0$ and $d_i^{LL}=-{\sum }_{T=1,T\ne L}^Q{d}_i^{LT}(L,T=1,2,3,\cdots ,Q)$.

The supra-Laplacian matrix of the multi-layer network is $H=H_L+H_I$, where $H_I=dD⨂I_N$, $I_N$ is the identity matrix of $N\times N$, and ⨂ is the Kronecker product.

$$H_L=\left(\begin{array}{ccccc}{H}^{\left(1\right)}& 0& \cdots & 0& 0\\ 0& H^{\left(2\right)}& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & H^{(Q-1)}& 0\\ 0& 0& \cdots & 0& H^{\left(Q\right)}\end{array}\right)=\underset{\genfrac{}{}{0pt}{}{L=1}{}}{\overset{Q}{⨁}}{H}^{\left(L\right)}$$

(2)

The eigenvalues of a square matrix A are the roots of its characteristic polynomial $\mathsf{\Psi }\left(\lambda \right)$, $\mathsf{\Psi }\left(\lambda \right)=|\lambda I-A|$. The eigenvalues of matrices have many applications [8,33]. In this paper, the eigenvalues of the supra-Laplacian matrix H are used to study the synchronizability of the networks. The eigenvalues of the supra-Laplacian matrix H are ${\lambda }_1,{\lambda }_2,{\lambda }_3,\cdots ,{\lambda }_{n-1},{\lambda }_n$. ${\lambda }_1=0$ is obtained from the properties of the supra-Laplacian matrix. According to the MSF, when the synchronized region is unbounded, the larger the ${\lambda }_2$ is, the stronger the synchronizability of the network is; when the synchronized region is bounded, the smaller the $r={\lambda }_{max}/{\lambda }_2$ is, the stronger the synchronizability of the network is [29].

2.2. Graph Theory and Star-Composed Networks

In this paper, we consider undirected and connected graph G. The vertex set of G is $V\left(G\right)=\{v_1,v_2,v_3,\cdots ,v_{n-1},v_n\}$, and the edge set of G is $E\left(G\right)=\left\{(v_i,v_j)\right|i,j=1,2,3\cdots ,n;i\ne j\}$. The edge set and the vertex set uniquely determine the graph. Let $K_p$ be a complete graph with p nodes, $P_n$ be a path of length n, and $E_k$ be a graph that has k vertices with an empty edge set.

The star-composed networks studied in this paper are introduced below. We first give two graph operations.

Definition 1.

([34]) Let $G_1,G_2$ represent two simple graphs with $n_1,n_2$ vertices and $m_1,m_2$ edges, respectively. The corona operation $G_1\circ G_2$ is the graph obtained from one copy of $G_1$ and $n_1$ copies of $G_2$, before joining the ith vertex of $G_1$ to every vertex in the ith copy of $G_2$$(i=1,2,\cdots ,n_1)$.

Definition 2.

([35]) Let $G_1,G_2$ represent two simple graphs. The conjunction operation $G_1▽G_2$ is defined as the graph, which is obtained from the disjoint union of $G_1$ and $G_2$ by adding the edges $\{ab:a\in V\left(G_1\right),b\in V\left(G_2\right)\}$.

$G_1(1,k,n)=(K_1▽E_k)\circ P_n$ is composed of one center node, k branchings, and paths of length n. $G_2(2,k-1,n)=(K_2▽E_{k-1})\circ P_n$ is composed of two center nodes, $k-1$ branchings, and paths of length n. Its connection relationships are shown in Figure 1. The red nodes represent the central nodes, the yellow nodes represent the branchings nodes, the green nodes represent the nodes of the paths, and the black solid lines represent the interlayer connection relationships between the nodes.

$G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ are composed of Q layers $G_1(1,k,n)$ and $G_2(2,k-1,n)$, respectively. The total number of nodes is $Q(nk+n+k+1)$, and the connection relationships are shown in Figure 2. The colored dotted lines represent the connection relationships between each layer.

Lemma 1.

([30]) Let $A,B$ be two $M\times M$ square matrices, and K be an integer; then,

$${\left|\begin{array}{cccc}A& B& \cdots & B\\ B& A& \cdots & B\\ ⋮& ⋮& \ddots & ⋮\\ B& B& \cdots & A\end{array}\right|}_{K\times K}={|A+(K-1)B|·|A-B|}^{(K-1)}.$$

(3)

3. The Eigenvalue Spectrum and Synchronizability Indexes of Multi-Layer Coupled Star-Composed Networks

3.1. The Eigenvalue Spectrum and Synchronizability Indexes of $G_1^Q(1,k,n)$

Similar to the network structure of $G_1^2(1,8,4)$ in Figure 2, each layer has the same topology. We assume that a is the intralayer coupling strength and d is the interlayer coupling strength. According to the concept of supra-Laplacian matrix, it is not difficult to obtain the supra-Laplacian matrix of $G_1^Q(1,k,n)$:

$$H_1={\left(\begin{array}{ccccc}(k+n)a& -a& -a& \cdots & -a\\ -a& {\mathsf{\Theta }}_1& 0& \cdots & 0\\ -a& 0& {\mathsf{\Theta }}_1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -a& 0& 0& \cdots & {\mathsf{\Theta }}_1\end{array}\right)}_{(k+1)\times (k+1)},$$

where ${\mathsf{\Theta }}_1=(n+1)a$.

$$H_2={\left(\begin{array}{ccccc}-a& 0& 0& \cdots & 0\\ 0& -a& 0& \cdots & 0\\ 0& 0& -a& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & -a\end{array}\right)}_{(k+1)\times (k+1)},$$

$$H_3={\left(\begin{array}{ccccc}2a& 0& 0& \cdots & 0\\ 0& 2a& 0& \cdots & 0\\ 0& 0& 2a& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 2a\end{array}\right)}_{(k+1)\times (k+1)},$$

$$H_4={\left(\begin{array}{ccccc}3a& 0& 0& \cdots & 0\\ 0& 3a& 0& \cdots & 0\\ 0& 0& 3a& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 3a\end{array}\right)}_{(k+1)\times (k+1)},$$

The supra-Laplacian matrix of Lth layer is:

$$H_1^{\left(L\right)}=\left(\begin{array}{ccccccc}{H}_1& H_2& H_2& H_2& \cdots & H_2& H_2\\ H_2& H_3& H_2& 0& \cdots & 0& 0\\ H_2& H_2& H_4& H_2& \cdots & 0& 0\\ H_2& 0& H_2& H_4& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ H_2& 0& 0& 0& \cdots & H_4& H_2\\ H_2& 0& 0& 0& \cdots & H_2& H_3\end{array}\right),$$

Combined with the interlayer supra-Laplacian matrix:

$$H=\left(\begin{array}{ccccc}{\mathsf{\Theta }}_2& -dI& -dI& \cdots & -dI\\ -dI& {\mathsf{\Theta }}_2& -dI& \cdots & -dI\\ -dI& -dI& {\mathsf{\Theta }}_2& \cdots & -dI\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -dI& -dI& -dI& \cdots & {\mathsf{\Theta }}_2\end{array}\right),$$

where ${\mathsf{\Theta }}_2=H_1^{\left(L\right)}+(Q-1)dI$.

Then, the characteristic polynomial is:

$${\mathsf{\Psi }}_I\left(\lambda \right)=det(\lambda I-H)=\left(\begin{array}{ccccc}{\mathsf{\Theta }}_3& dI& dI& \cdots & dI\\ dI& {\mathsf{\Theta }}_3& dI& \cdots & dI\\ dI& dI& {\mathsf{\Theta }}_3& \cdots & dI\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ dI& dI& dI& \cdots & {\mathsf{\Theta }}_3\end{array}\right),$$

where ${\mathsf{\Theta }}_3=\lambda I-H_1^{\left(L\right)}-(Q-1)dI$.

According to Lemma 1, we have:

$${\mathsf{\Psi }}_I\left(\lambda \right)=|\lambda I-H_1^{\left(L\right)}\left|\right|\lambda I-H_1^{\left(L\right)}{-QdI|}^{Q-1}.$$

(4)

The eigenvalue spectrum of H is:

$$0,(n+1)a,\frac{k+n+2\pm \sqrt{{(k+n+2)}^2-4(k+1)}}{2}a,\underset{k-1}{\underbrace{\frac{n+2\pm \sqrt{{(n+2)}^2-4}}{2}a}},$$

$$\underset{Q-1}{\underbrace{Qd,(n+1)a+Qd}},\underset{k+1}{\underbrace{4asi{n}^2\left(\frac{\alpha \pi }{2n}\right)+a}}(\alpha =1,2,\cdots ,n-1),$$

$$\underset{Q-1}{\underbrace{\frac{k+n+2\pm \sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd}},$$

$$\underset{(Q-1)(k-1)}{\underbrace{\frac{n+2\pm \sqrt{{(n+2)}^2-4}}{2}a+Qd}},\underset{(Q-1)(k+1)}{\underbrace{4asi{n}^2\left(\frac{\alpha \pi }{2n}\right)+a+Qd}}(\alpha =1,2,\cdots ,n-1).$$

Deduced from the preliminary knowledge, the minimum nonzero eigenvalue is:

$${\lambda }_2=min\{\frac{n+2-\sqrt{{(n+2)}^2-4}}{2}a,Qd\}.$$

The maximum eigenvalue is:

$${\lambda }_{max}=\frac{k+n+2+\sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd.$$

Then, we have:

$$r=(\frac{k+n+2+\sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd)/{\lambda }_2.$$

(5)

3.2. The Eigenvalue Spectrum and Synchronizability Indexes of $G_2^Q(2,k-1,n)$

Similar to the network structure of $G_2^2(2,7,4)$ in Figure 2, we examine the eigenvalue spectrum of $G_2^Q(2,k-1,n)$ using a similar method:

$$H_5={\left(\begin{array}{ccccc}(k+n)a& -a& -a& \cdots & -a\\ -a& (k+n)a& -a& \cdots & -a\\ -a& -a& {\mathsf{\Theta }}_4& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -a& -a& 0& \cdots & {\mathsf{\Theta }}_4\end{array}\right)}_{(k+1)\times (k+1)},$$

where ${\mathsf{\Theta }}_4=(n+2)a$.

The supra-Laplacian matrix of the Lth layer is:

$$H_2^{\left(L\right)}=\left(\begin{array}{ccccccc}{H}_5& H_2& H_2& H_2& \cdots & H_2& H_2\\ H_2& H_3& H_2& 0& \cdots & 0& 0\\ H_2& H_2& H_4& H_2& \cdots & 0& 0\\ H_2& 0& H_2& H_4& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ H_2& 0& 0& 0& \cdots & H_4& H_2\\ H_2& 0& 0& 0& \cdots & H_2& H_3\end{array}\right),$$

Combined with the interlayer supra-Laplacian matrix:

$$\tilde{H}=\left(\begin{array}{ccccc}{\mathsf{\Theta }}_5& -dI& -dI& \cdots & -dI\\ -dI& {\mathsf{\Theta }}_5& -dI& \cdots & -dI\\ -dI& -dI& {\mathsf{\Theta }}_5& \cdots & -dI\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -dI& -dI& -dI& \cdots & {\mathsf{\Theta }}_5\end{array}\right),$$

where ${\mathsf{\Theta }}_5=H_2^{\left(L\right)}+(Q-1)dI$.

The characteristic polynomial is:

$${\mathsf{\Psi }}_{II}\left(\lambda \right)=|\lambda I-H_2^{\left(L\right)}\left|\right|\lambda I-H_2^{\left(L\right)}{-QdI|}^{Q-1}.$$

(6)

The eigenvalue spectrum of $\tilde{H}$ is: $0,(n+1)a,\frac{k+n+2\pm \sqrt{{(k+n+2)}^2-4(k+1)}}{2}a,\underset{k-2}{\underbrace{\frac{n+3\pm \sqrt{{(n+3)}^2-8}}{2}a}},$

$$\underset{Q-1}{\underbrace{Qd,(n+1)a+Qd}},\underset{k+1}{\underbrace{4asi{n}^2\left(\frac{\beta \pi }{2n}\right)+a}}(\beta =1,2,\cdots ,n-1),$$

$$\underset{(Q-1)}{\underbrace{\frac{k+n+2\pm \sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd}},$$

$$\underset{(k-2)(Q-1)}{\underbrace{\frac{n+3\pm \sqrt{{(n+3)}^2-8}}{2}a+Qd}},\underset{(k+1)(Q-1)}{\underbrace{4asi{n}^2\left(\frac{\beta \pi }{2n}\right)+a+Qd}}(\beta =1,2,\cdots ,n-1).$$

By comparison, the minimum nonzero eigenvalue is:

$${\lambda }_2=min\{\frac{n+3-\sqrt{{(n+3)}^2-8}}{2}a,Qd\}.$$

The maximum eigenvalue is:

$${\lambda }_{max}=\frac{k+n+2+\sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd.$$

Then, we have:

$$r=(\frac{k+n+2+\sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd)/{\lambda }_2.$$

(7)

3.3. The Eigenvalue Spectrum and Synchronizability Indexes of $G_p^Q(p,k+1-p,n)$

The multi-layer-coupled star-composed networks with p central nodes are considered, $(3\le p\le k)$, and the multi-center-coupled star-composed networks can be written as $G_p^Q(p,k+1-p,n)=(K_p\nabla E_{k+1-p})\circ P_n$.

The eigenvalue spectrum of the supra-Laplace matrix obtained by the same method is: $0,(n+1)a,\underset{p}{\underbrace{\frac{k+n+2\pm \sqrt{{(k+n+2)}^2-4(k+1)}}{2}a}},\underset{k-p}{\underbrace{\frac{n+p+1\pm \sqrt{{(n+p+1)}^2-4p}}{2}a}},$

$$\underset{Q-1}{\underbrace{Qd,(n+1)a+Qd}},\underset{k+1}{\underbrace{4asi{n}^2\left(\frac{\gamma \pi }{2n}\right)+a}}(\gamma =1,2,\cdots ,n-1),$$

$$\underset{p(Q-1)}{\underbrace{\frac{k+n+2\pm \sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd}},$$

$$\underset{(k-p)(Q-1)}{\underbrace{\frac{n+p+1\pm \sqrt{{(n+p+1)}^2-4p}}{2}a+Qd}},\underset{(k+1)(Q-1)}{\underbrace{4asi{n}^2\left(\frac{\gamma \pi }{2n}\right)+a+Qd}}(\gamma =1,2,\cdots ,n-1).$$

We can easily obtain the required indicators. The minimum nonzero eigenvalue is as follows:

$${\lambda }_2=min\{\frac{n+p+1-\sqrt{{(n+p+1)}^2-4p}}{2}a,Qd\}.$$

The maximum eigenvalue is:

$${\lambda }_{max}=\frac{k+n+2+\sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd.$$

Then, we have:

$$r=(\frac{k+n+2+\sqrt{{(k+n+2)}^2-4(k+1)}}{2}a+Qd)/{\lambda }_2.$$

(8)

4. Numerical Simulation Experiment and Analysis

In the analysis of the above sections, we can obtain the relationships between the synchronizability and various parameters of multi-layer-coupled star-composed networks. This section tests the accuracy of the above mentioned theoretical results through a large number of numerical simulation experiments. Firstly, the state trajectories of star-composed networks are described to illustrate that their synchronization can be realized by MATLAB. Secondly, considering the actual initial value, the change images of the synchronizability of the star-composed networks are obtained by MATLAB. Finally, we analyze the relationships between the synchronizability and parameters of multi-layer-coupled star-composed networks.

We consider an example of a coupling system composed of two layers. The network topology is shown in Figure 3. The number of layers is $Q=2$, $N=9$ is the number of nodes in each layer, and the dynamics of each node are given as follows [36]:

$$\frac{d{x}_i^L}{dt}=\frac{|x_i^L+1|-|x_i^L-1|}{2}+\sum _{j\in N_i}{a}_{ij}(x_j^L-x_i^L)+\sum _{j\in N_i}{b}_{ij}(x_j^{3-L}-x_i^L)$$

where $x_i^L$ is the state of a node, $\frac{|x+1|-|x-1|}{2}$ represents the inherent nonlinear dynamics of the node, $i=1,2,3,\cdots ,8,9$, $L=1,2$; $a_{ij}$ represents the intralayer connection, and the connection weight is 1. $b_{ij}$ denotes the interlayer connection, and the connection weight is $2.5$. The initial value of the system is randomly selected. Through the simulation, the state trajectory of the system is shown in Figure 4 and Figure 5. It can be seen that the state trajectory of nodes achieves synchronization under the coupling effect.

In the following analysis, in order to distinguish the synchronizability indexes of multi-layer-coupled star-composed networks, when the synchronized region is unbounded, we use ${\lambda }_{2a}$ and ${\lambda }_{2b}$ to represent the synchronizability indexes of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$. When the synchronized region is bounded, the synchronizability indexes of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ are represented by $r_a$ and $r_b$, respectively.

• When the synchronized region is unbounded, we set $n=5,a=2,d=0.008,Q=40$, with the increase in k as in Figure 6. ${\lambda }_{2a}$ remains unchanged at 0.2918 and ${\lambda }_{2b}$ remains unchanged at 0.32. The synchronizability is not affected by k. We set $k=10,a=2,$$d=0.008,Q=40$. With the increase in n in Figure 7, ${\lambda }_{2a}$ remains unchanged at 0.32 and then decreases monotonically to 0.0196, while ${\lambda }_{2b}$ remains unchanged at 0.32 and then decreases monotonically to 0.0486. The synchronizability of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ reached the maximum at $n=\lfloor a/Qd+Qd/a-2\rfloor$ and $n=\lfloor 2a/Qd+Qd/a-3\rfloor$, respectively, and the length of the paths reaches the maximum. We set $k=10,n=5,d=0.008,Q=40$, with the increase in a as in Figure 8. ${\lambda }_{2a}$ monotonically increases from 0.0729 to 0.32 and then remains unchanged, while ${\lambda }_{2b}$ monotonically increases from 0.1629 to 0.32 and then remains unchanged. The synchronizability of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ reaches the maximum at $a=(n+2+\sqrt{{(n+2)}^2-4})Qd/2$ and $a=(n+3+\sqrt{{(n+3)}^2-8})Qd/4$, respectively, and the intralayer coupling strength is the minimum. We set $k=10,$$n=5,$$a=2,$$Q=40$. With the increase in d in Figure 9, ${\lambda }_{2a}$ monotonically increases from 0.04 to 0.2918 and then remains unchanged, while ${\lambda }_{2b}$ monotonically increases from 0.04 to 0.6515 and then remains unchanged. The synchronizability of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ reaches the maximum at $d=(n+2-\sqrt{{(n+2)}^2-4})a/2Q$ and $d=(n+3-\sqrt{{(n+3)}^2-8})a/2Q$, respectively, and the interlayer coupling strength is the minimum. We set $k=10,n=5,a=2,d=0.008$, with the increase in as Q in Figure 10. ${\lambda }_{2a}$ monotonically increases from 0.016 to 0.2918 and then remains unchanged, while ${\lambda }_{2b}$ monotonically increases from 0.016 to 0.6515 and then remains unchanged. The synchronizability of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ reaches the maximum at $Q=(n+2-\sqrt{{(n+2)}^2-4})a/2d$ and $Q=(n+3-\sqrt{{(n+3)}^2-8})a/2d$, respectively, and the number of layers is the minimum.
• When the synchronized region is bounded, we set $n=5,a=2,d=0.008,Q=40$, with the increase in k as in Figure 6. $r_a$ monotonically increases from 91.2508 to 1094.7 and $r_b$ monotonically increases from 83.2082 to 998.1725. When k is the minimum, the synchronizability is the maximum. We set $k=10,a=2,d=0.008,Q=40$. With the increase in n in Figure 7, $r_a$ rapidly increases from 127.0276 to 17141 and $r_b$ rapidly increases from 127.0276 to 6922.5. When n is the minimum, the synchronizability is the maximum. We set $k=10,n=5,d=0.008,Q=40$. With the increase in a in Figure 8, $r_a$ decreases slowly from 174.5482 to 171.1677 and then increases rapidly to 233.7460, while $r_b$ decreases slowly from 78.1736 to 77.2113 and then increases rapidly to 233.7460. The synchronizability of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ reaches the maximum at $a=(n+2+\sqrt{{(n+2)}^2-4})Qd/2$ and $a=(n+3+\sqrt{{(n+3)}^2-8})Qd/4$, respectively. We set $k=10,n=5,a=2,Q=40$. With the increase in d in Figure 9, $r_a$ decreases rapidly from 1242.3 to 171.2582 and then increases monotonically to 581.4077, while $r_b$ decreases slowly from 1242.3 to 77.2526 and then increases rapidly to 260.3906. The synchronizability of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ reaches the maximum at $d=(n+2-\sqrt{{(n+2)}^2-4})a/2Q$ and $d=(n+3-\sqrt{{(n+3)}^2-8})a/2Q$, respectively. We set $k=10,n=5,a=2,d=0.008$. With the increase in Q in Figure 10, $r_a$ rapidly decreases from 3104.3 to 171.1760 and then slowly increases to 172.9032, while $r_b$ decreases rapidly from 3104.3 to 77.2158 and then increases slowly to 77.4368. The synchronizability of $G_1^Q(1,k,n)$ and $G_2^Q(2,k-1,n)$ reaches the maximum at $Q=(n+2-\sqrt{{(n+2)}^2-4})a/2d$ and $Q=(n+3-\sqrt{{(n+3)}^2-8})a/2d$, respectively.

5. Main Results

5.1. The Parameters That Affect the Synchronizability of Star-Composed Networks

Based on the above situation, for $G_1^Q(1,k,n)$ networks when the synchronized region is unbounded, the synchronizability of the networks is only related to the number of layers, the length of paths, the interlayer coupling strength, and the intralayer coupling strength. When the synchronized region is bounded, the synchronizability of the networks is related to the interlayer coupling strength, the intralayer coupling strength, the length of paths, the branchings, and the number of layers. See

Table 1. Relationships among all parameters affecting $G_1^Q(1,k,n)$ networks’ synchronizability.

		$\mathit{k}\uparrow$	$\mathit{n}\uparrow$	$\mathit{a}\uparrow$	$\mathit{d}\uparrow$	$\mathit{Q}\uparrow$
${\lambda }_2=min\{Qd,{\mathrm{\Delta }}_1\}$	${\mathrm{\Delta }}_1<Qd$	−	↓	↑	−	−
	${\mathrm{\Delta }}_1>Qd$	−	−	−	↑	↑
$r={\lambda }_{max}/{\lambda }_2$	${\mathrm{\Delta }}_1<Qd$	↑	↑	↓	↑	↑
	${\mathrm{\Delta }}_1>Qd$	↑	↑	↑	↓	↓

${\mathsf{\Delta }}_1=\frac{n+2-\sqrt{{(n+2)}^2-4}}{2}a$, ↑ strengthen, ↓ weaken, − unchanged.

for the specific relationships. For $G_2^Q(2,k-1,n)$ networks, whether the synchronized region is bounded or unbounded, the factors affecting the networks’ synchronizability are the same as $G_1^Q(1,k,n)$; see

Table 2. Relationships among all parameters affecting $G_2^Q(2,k-1,n)$ networks’ synchronizability.

		$\mathit{k}\uparrow$	$\mathit{n}\uparrow$	$\mathit{a}\uparrow$	$\mathit{d}\uparrow$	$\mathit{Q}\uparrow$
${\lambda }_2=min\{Qd,{\mathrm{\Delta }}_2\}$	${\mathrm{\Delta }}_2<Qd$	−	↓	↑	−	−
	${\mathrm{\Delta }}_2>Qd$	−	−	−	↑	↑
$r={\lambda }_{max}/{\lambda }_2$	${\mathrm{\Delta }}_2<Qd$	↑	↑	↓	↑	↑
	${\mathrm{\Delta }}_2>Qd$	↑	↑	↑	↓	↓

${\mathsf{\Delta }}_2=\frac{n+3-\sqrt{{(n+3)}^2-8}}{2}a$, ↑ strengthen, ↓ weaken, − unchanged.

for the specific relationships.

5.2. Optimizing Synchronizability of $G_p^Q(p,k+1-p,n)$ Networks

We further discuss the synchronizability of the multi-layer-coupled star-composed networks with the central nodes p. Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 show the relationships between the synchronizability of $G_p^Q(p,k+1-p,n)(p=1,2)$ and the affecting factors. $G_p^Q(p,k+1-p,n)$ have similar figures; the specific relationship is shown in

Table 3. Relationships between all parameters affecting $G_p^Q(p,k+1-p,n)(3\le p\le k)$ networks’ synchronizability.

		$\mathit{k}\uparrow$	$\mathit{n}\uparrow$	$\mathit{a}\uparrow$	$\mathit{d}\uparrow$	$\mathit{Q}\uparrow$
${\lambda }_2=min\{Qd,{\mathrm{\Delta }}_3\}$	${\mathrm{\Delta }}_3<Qd$	−	↓	↑	−	−
	${\mathrm{\Delta }}_3>Qd$	−	−	−	↑	↑
$r={\lambda }_{max}/{\lambda }_2$	${\mathrm{\Delta }}_3<Qd$	↑	↑	↓	↑	↑
	${\mathrm{\Delta }}_3>Qd$	↑	↑	↑	↓	↓

${\mathsf{\Delta }}_3=\frac{n+p+1-\sqrt{{(n+p+1)}^2-4p}}{2}a$, ↑ strengthen, ↓ weaken, − unchanged.

.

• When the synchronized region is unbounded, the synchronizability is not affected by k. With the increase in the length of paths n, the synchronizability of the networks remains unchanged and then decreases. The synchronizability of $G_p^Q(p,k+1-p,n)$ reaches the maximum at $n=\lfloor ap/Qd+Qd/a-(p+1)\rfloor$ and the length of paths reaches the maximum. With the increase in intralayer coupling strength a, the synchronizability of the networks increases first and then remains unchanged. The synchronizability of $G_p^Q(p,k+1-p,n)$ reaches the maximum at $a=(n+p+1+\sqrt{{(n+p+1)}^2-4p})Qd/2p$, and the intralayer coupling strength is the minimum. With the increase in interlayer coupling strength d, the synchronizability of the networks increases first and then remains unchanged. The synchronizability of $G_p^Q(p,k+1-p,n)$ reaches the maximum at $d=(n+p+1-\sqrt{{(n+p+1)}^2-4p})a/2Q$ and the interlayer coupling strength is the minimum. With the increase in the number of layers Q, the synchronizability of the networks increases first and then remains unchanged. The synchronizability of $G_p^Q(p,k+1-p,n)$ reaches the maximum at $Q=(n+p+1-\sqrt{{(n+p+1)}^2-4p})a/2d$ and the number of layers is the minimum. It can be seen that in order to improve the synchronizability of the networks, the length of paths n should be appropriately reduced, and the intralayer coupling strength a, interlayer coupling strength d, and number of layers Q should be appropriately increased.
• When the synchronized region is bounded, with the increase in the branchings k and the length of paths n, the synchronizability of the networks continues to weaken. With the increase in the intralayer coupling strength a, the synchronizability of the networks increases first and then decreases. The synchronizability of $G_p^Q(p,k+1-p,n)$ reaches the maximum at $a=(n+p+1+\sqrt{{(n+p+1)}^2-4p})Qd/2p$. With the increase in the interlayer coupling strength d, the synchronizability of the networks increases first and then decreases. The synchronizability of $G_p^Q(p,k+1-p,n)$ reaches the maximum at $d=(n+p+1-\sqrt{{(n+p+1)}^2-4p})a/2Q$. With the increase in the number of layers Q, the synchronizability of the networks increases first and then decreases. The synchronizability of $G_p^Q(p,k+1-p,n)$ reaches the maximum at $Q=(n+p+1-\sqrt{{(n+p+1)}^2-4p})a/2d$. From the above analysis, in order to improve the synchronizability of the networks, we reduce the branchings k and the length of paths n and appropriately increase the intralayer coupling strength a, interlayer coupling strength d, and number of layers Q to the corresponding inflection points.

5.3. Impact of the Number of Central Nodes on the Synchronizability

In the numerical simulation, we determined that the synchronizability of the networks with two central nodes is better than that of the networks with one central node. Whether the synchronized region is bounded or unbounded, with the increase in the central nodes p, the synchronizability of the networks will not be weaker than the synchronizability of the networks with fewer central nodes.

There are still many problems that need to be solved in multi-layer-coupled star-composed networks. For example, what will happen to the synchronizability after giving different coupling strengths to multi-layer-coupled star-composed networks? For multi-layer-coupled star-composed networks, we need to further study how to keep the synchronizability unchanged when changing multiple parameters.