Eigenvector-based analysis of cluster synchronization in general complex networks of coupled chaotic oscillators

Abstract

Whereas topological symmetries have been recognized as crucially important to the exploration of synchronization patterns in complex networks of coupled dynamical oscillators, the identification of the symmetries in large-size complex networks remains as a challenge. Additionally, even though the topological symmetries of a complex network are known, it is still not clear how the system dynamics is transited among different synchronization patterns with respect to the coupling strength of the oscillators. We propose here the framework of eigenvector-based analysis to identify the synchronization patterns in the general complex networks and, incorporating the conventional method of eigenvalue-based analysis, investigate the emergence and transition of the cluster synchronization states. We are able to argue and demonstrate that, without a prior knowledge of the network symmetries, the method is able to predict not only all the cluster synchronization states observable in the network, but also the critical couplings where the states become stable and the sequence of these states in the process of synchronization transition. The efficacy and generality of the proposed method are verified by different network models of coupled chaotic oscillators, including artificial networks of perfect symmetries and empirical networks of non-perfect symmetries. The new framework paves a way to the investigation of synchronization patterns in large-size, general complex networks.

Graphical abstract

Keywords

cluster synchronization / complex networks / network symmetry / coupled oscillators

1 Introduction

As a universal phenomenon in complex systems of coupled dynamical units, synchronization has been broadly interested by researchers from different fields in the past decades [1–4]. Briefly, synchronization refers to the coherent motion of coupled oscillators, which is observed normally when the coupling strength between the oscillators exceeds a critical value. Depending on the correlation between the oscillators, synchronization can be realized in different forms, e.g., complete synchronization, phase synchronization and generalized synchronization [3]. In synchronization studies, a key issue is to find, numerically or analytically, the critical coupling for generating synchronization. For systems consisting of linearly coupled identical chaotic oscillators, complete synchronization can be realized and the critical coupling for synchronization can be analyzed by the formalism of master stability function (MSF) [5–7]; for systems composed of nonlinearly coupled non-identical phase oscillators, the critical coupling characterizing the onset of phase synchronization can be estimated by some mean-field approaches [8, 9]. Whereas earlier studies of oscillator synchronization have been focusing on small-size systems of regular coupling structures [10], e.g., the lattices and globally coupled systems, recent studies extend the studies to large-size systems with complex coupling structures, e.g., the synchronization behaviors in complex networks [11–16]. The adoption of complex networks opens up a new window to the study of oscillator synchronization, in which the important roles of network structure on synchronization have been revealed [17–20].

For systems of coupled identical chaotic oscillators, an interesting phenomenon observed in modeling systems and experiments is that under certain circumstances the oscillators can be synchronized in groups, namely the phenomenon of cluster synchronization (CS) [21–31]. In CS, oscillators within each cluster are highly correlated, but not if the oscillators belong to different clusters. Compared to global synchronization, CS is normally observed at weaker couplings in the transition regime from desynchronization to global synchronization. The spatial distribution of the clusters on the network defines a pattern and, by varying the coupling strength, the same system can present different synchronization patterns. In the study of CS, two of the central questions are: (i) how to find the possible synchronization patterns for the given network structure and (ii) how to estimate the critical couplings generating the patterns. For systems possessing regular coupling structures, e.g., the ring-structure network, the patterns are normally generated through the mechanism of symmetry breaking, and the spatial and dynamical properties of the patterns can be analyzed by methods such as eigenvalue-based analysis [32, 33]. Yet challenges arise when the oscillators are coupled on complex networks [34–41]. Different from the regular networks in which the synchronization patterns can be inspected visually, the patterns on complex networks are blurred by the network structures. As such, to identify and analyze the synchronization patterns on complex networks, special techniques and methods different from the conventional ones should be employed [42].

One approach to analyzing the synchronization patterns in complex networks is exploiting the information of network topological symmetries [34–41, 43–50]. Briefly, if a set of oscillators whose permutations on the network do not change the network dynamics, the oscillators are said to be symmetric with each other and the set of symmetric oscillators have the potential to form a synchronization cluster [43]. Whereas the symmetries of small-size networks can be figured out straightforwardly, the finding of all symmetries in large-size complex networks requires some sophisticated techniques [38, 51–54]. Once the network symmetries are identified, the next step is to analyze the stability of the synchronization patterns associated with the symmetries. In doing this, a general approach is to decouple the dynamics of the CS state from the dynamics of the perturbations. This can be done by techniques such as irreducible representations (IRR) and generalized MSF [36–38, 47], by which the critical couplings for generating the CS states can be estimated. Besides the approach of network-symmetry-based analysis, the CS states of complex networks can also be analyzed by methods such as external equitable partition (EEP) and simultaneous block diagonalization (SBD) [52–58]. The EEP method concerns the dynamics of a small-size quotient network, and the contents of each cluster are defined according to the inputting signals [52–54]; the SBD method is based on partitions of the networked nodes, and is featured by a simultaneous block diagonalization of the network coupling matrix [55–58].

The existing methods of CS analysis, however, encounter a major difficulty when dealing with realistic systems: perfect symmetry is hardly observed in real-world complex networks [38, 43, 58]. In grouping network nodes into clusters, the existing methods require a perfect symmetry of the nodes [34–41, 43–50], i.e., the exchanges of the symmetric nodes on the network does not affect the network dynamics. This rigorous requirement makes most of the clusters identified in typical complex networks trivial, in the sense that the clusters contain only several or just a single node [59]. Yet mesoscale synchronization clusters do exist in realistic systems, which are believed as playing crucial roles in realizing the system functions, e.g., the brain networks [60–63]. In specific, empirical studies show that many realistic networks are composed of communities, with nodes inside each community being densely connected and the connections between communities are sparse [64]. When oscillators are coupled on a complex community network, distinct synchronization clusters can be observed and the contents of the clusters are consistent with the partition of the communities [60, 65, 66]. As both the inter- and intra-connections are randomly established, perfect symmetry does not exist in complex community networks in general, which, according to the existing methods of CS analysis, implies that distinct synchronization clusters can not be observed. The contradiction between theory and reality makes it necessary to generalize the current methods of CS analysis to the general complex networks, which is the major objective of our present work. Besides the concern of network symmetry, another question encountered in analyzing CS is the separated analysis of cluster partition and stability. To obtain the critical coupling associated with a specific CS state by the existing methods, one needs to first find the contents of the clusters associated with the state by some node-partition techniques, then construct the quotient network governing the dynamics of this CS state, and finally evaluate the stability of the CS state and obtain the critical coupling. As the quotient network is dependent on the partition of the clusters, the stability of the CS states needs to be evaluated individually and separately, making a complete analysis of the CS states in large-size complex networks time-consuming and inconvenient. The second objective of our present work is to combine the two-step analysis, namely cluster partition and stability analysis, into a one-step analysis, so as to facilitate the exploration of CS in large-size networks.

The mission of our present work is to argue that the above objectives, i.e., extending the CS studies to general complex networks and combining the analysis of cluster identification and stability into a single analysis, can be accomplished within the framework of eigenvector-based analysis. To be more specific, we shall demonstrate in different network models that, without a priori knowledge of the network symmetries, the proposed method is able to predict not only the CS states to be emerged in synchronization transition, but also the range over which each CS state is observable in the parameter space and how the network dynamics is transited among the CS states as the coupling parameter is varying. The rest of the paper is organized as follows. In the following section, we shall present the theoretical framework used for analyzing CS in general complex networks. The applications of the proposed framework to different network models, including small-size networks of perfect symmetries and large-size complex networks of non-perfect symmetries, will be reported in Section 3. Finally, discussion and conclusions will be given in Section 4.

2 Eigenvector-based analysis

The dynamical system we consider here is a complex network of coupled identical chaotic oscillators, with the system dynamics described by the equations:

(1){\dot{\mathit{x}}}_i=\mathit{F}({\mathit{x}}_i)+\epsilon {\displaystyle \sum _{j=1}^N{w}_{ij}\mathit{H}({\mathit{x}}_j).}

Here, i,j=1,\dots ,N are the node (oscillator) indices, {\mathit{x}}_i represents the state vector of the ith oscillator, \mathit{F}(\mathit{x}) denotes the dynamics of the isolated oscillators, \mathit{H}(\mathit{x}) defines the coupling function, and \epsilon is the uniform coupling strength. The coupling relationship of the oscillators is captured by the weighted matrix \mathit{W}=\{w_{ij}\}, with w_{ij}>0 denoting the strength of the coupling that node i is received from node j. If nodes i and j are not connected, we set w_{ij}=w_{ji}=0. The diagonal elements of \mathit{W} are set as w_{ii}=-\sum _{j(j\ne i)}{w}_{ij}, i.e., \mathit{W} is a Laplacian matrix. The model described by Eq. (1) has been widely adopted in the literature for exploring the synchronization behaviors of coupled chaotic oscillators, which is a good approximation to the dynamics of many realistic systems in the vicinity of their functional states [3, 14, 15].

As the oscillators are linearly coupled and \mathit{W} is a Laplacian matrix, the state of global synchronization is a solution to Eq. (1). Denote \mathit{s} as the manifold of global synchronization state, i.e., \mathit{s}={\mathit{x}}_1={\mathit{x}}_2=\dots ={\mathit{x}}_N, the starting point of our analysis is to evaluate the stability of this state in the presence of small, random perturbations. Let \delta {\mathit{x}}_i={\mathit{x}}_i-\mathit{s} be infinitesimal perturbation added onto oscillator i in the global synchronization state, the evolution of \delta {\mathit{x}}_i is governed by the variational equations

(2)\delta {\dot{\mathit{x}}}_i=\mathit{D}\mathit{F}(\mathit{s})\delta {\mathit{x}}_i+\epsilon {\displaystyle \sum _{j=1}^N{w}_{ij}\mathit{D}\mathit{H}(\mathit{s})\delta {\mathit{x}}_j,}

with \mathit{D}\mathit{F}(\mathit{s}) and \mathit{D}\mathit{H}(\mathit{s}) the Jacobian matrices of \mathit{F} and \mathit{H} evaluated on the synchronization manifold \mathit{s}. Transforming Eq. (2) into the space spanned by the eigenvectors of \mathit{W}, we have

(3)\delta {\dot{\mathit{y}}}_i=[\mathit{D}\mathit{F}(\mathit{s})+\epsilon {\lambda }_i\mathit{D}\mathit{H}(\mathit{s})]\delta {\mathit{y}}_i.

Here, \mathrm{\Delta }\mathit{Y}=[\delta {\mathit{y}}_1,\delta {\mathit{y}}_2,\dots ,\delta {\mathit{y}}_N{]}^{\mathrm{T}}={\mathit{V}}^{-1}[\delta {\mathit{x}}_1,\delta {\mathit{x}}_2,\dots ,\delta {\mathit{x}}_N{]}^{\mathrm{T}} are the perturbation modes in the new space, \mathit{V} is the transformation matrix composed by the eigenvectors of \mathit{W}, and 0={\lambda }_1>{\lambda }_2\ge \dots \ge {\lambda }_N are the eigenvalues of \mathit{W}. Among the N modes, the mode associated with {\lambda }_1 describes the motion parallel to the synchronization manifold, and the modes associated with {\lambda }_{2,\dots ,N} describe the motions transverse to the synchronization manifold. Following the MSF formalism [5–7], we introduce the generic coupling strength \sigma \equiv -\epsilon \lambda and rewrite Eq. (3) as

(4)\delta {\dot{\mathit{y}}}_i=[\mathit{D}\mathit{F}(\mathit{s})-{\sigma }_i\mathit{D}\mathit{H}(\mathit{s})]\delta {\mathit{y}}_i.

For the global synchronization state to be stable, the necessary condition is that all the perturbation modes transverse to the synchronization state should be damping with time. That is, the largest conditional Lyapunov exponent, {\mathrm{\Lambda }}_i, of Eq. (4) should be negative for modes i=2,\dots ,N. As Eq. (4) applies to all the perturbation modes, it is named the master equation of the perturbation dynamics [5]. By solving the master equation numerically (or analytically in some special cases), we can obtain the variation of \mathrm{\Lambda } with respect to \sigma, which defines the MSF curve. In the MSF curve, the regions with constitute the stable domain in the parameter space, and the regions with \mathrm{\Lambda }>0 constitute the unstable domain. Depending on the nodal dynamics and the coupling function, the stable domain may have different forms, e.g., bounded or unbounded [7]. For simplicity, we consider here only the case of unbounded stable domain, in which \mathrm{\Lambda } becomes negative when \sigma is larger than a critical value {\sigma }_c. In this case, global synchronization will be achieved when \epsilon >{\epsilon }_c=-{\sigma }_c/{\lambda }_2, with {\lambda }_2 the 2nd largest eigenvalue of \mathit{W}. So far, our analysis is identical to the standard MSF formalism [5–7], in which the only role of the 2nd largest eigenvector of the network coupling matrix \mathit{W} is to construct the transformation matrix \mathit{V}, while the stability of the global synchronization state is dependent on only the eigenvalues of the network coupling matrix. In what follows, we are going to argue that the eigenvectors play an important role in exploring the CS states in general complex networks of coupled chaotic oscillators.

We move on to analyze the desynchronization state of the network when . We consider first the situation when a single mode in the transverse space is unstable. As we are focusing on MSF of unbounded stable domain, the most unstable mode is associated with {\lambda }_2. That is, by decreasing \epsilon from {\epsilon }_c, the mode \delta {\mathit{y}}_2 is first emerged from the uniform background of global synchronization. Since the mode associated with {\lambda }_1=0 is always unstable, we thus have for this situation two unstable modes in fact, \delta {\mathit{y}}_1 and \delta {\mathit{y}}_2. Denote \delta {\mathit{y}}_1(t) and \delta {\mathit{y}}_2(t) as the modes of {\lambda }_1 and {\lambda }_2 at time t respectively, we transfer them into the node space and obtain

(5)\delta {\mathit{x}}_i(t)=v_{1,i}\delta {\mathit{y}}_1(t)+v_{2,i}\delta {\mathit{y}}_2(t).

Here, {\mathit{v}}_1=\{v_{1,i}{\}}_{i=1,\dots ,N} is the eigenvector associated with {\lambda }_1, and {\mathit{v}}_2=\{v_{2,i}{\}}_{i=1,\dots ,N} is the eigenvector associated with {\lambda }_2. As {\lambda }_1=0, we have v_{1,i}=1/\sqrt{N} for all the elements of {\mathit{v}}_1. That is, the first term on the right-hand-side (RHS) of Eq. (5) is identical for all the oscillators, and the perturbation that oscillator i is away from the global synchronization state is solely determined by the second term on the RHS of the equation. More specifically, we have

(6)\delta {\mathit{x}}_i(t)=c(t)+v_{2,i}\delta {\mathit{y}}_2(t),

with c(t) a variable independent of the oscillator index. Therefore, by checking the values of v_{2,i}, the stability of the oscillators can be evaluated individually: the larger is |v_{2,i}|, the more unstable is oscillator i. In particular, if v_{2,i}=v_{2,j}, the states of the pair of oscillators i and j will be identical during the process of system evolution, i.e., they are completely synchronized. The set of oscillators with the identical eigenvector element thus forms a synchronization cluster. That is, by checking just the elements of the eigenvector {\mathit{v}}_2, we are able to identify all the synchronization clusters when the mode of {\lambda }_2 is unstable.

We consider next the situation when several perturbation modes are unstable. This occurs when \epsilon is much less than {\epsilon }_c and means that the network is deeply desynchronized. We assume that the statistical properties of the oscillators, e.g., the measure of the oscillator trajectory, can be approximated by that of the isolated oscillator. (The validity of this approximation will be verified by numerical simulations later.) With this approximation, the stability of the perturbation modes is still governed by Eq. (4). We assume further that \tilde{m} modes are unstable, and the eigenvectors of the unstable modes are {\mathit{v}}_k, with k=2,\dots ,\tilde{m}+1. Note that for the case of unbounded MSF curve, the perturbation modes are destabilized in sequence by the descending order of the eigenvalues, with the mode of {\lambda }_2 being the first one and the mode of {\lambda }_N being the last one. Specifically, the kth perturbation mode becomes unstable when . When \tilde{m} modes are unstable in the network, the perturbation of the ith oscillator can be written as

(7)\delta {\mathit{x}}_i(t)=c(t)+{\displaystyle \sum _{k=2}^{\tilde{m}+1}{v}_{k,i}\delta {\mathit{y}}_k(t).}

Still, c(t) represents the perturbation in parallel to the synchronization manifold, which is independent of the oscillator index. Now, for oscillators i and j to be completely synchronized, the necessary condition is that v_{k,i}=v_{k,j} for all the unstable modes k=2,\dots ,\tilde{m}+1. Hence, by checking only the eigenvectors of the \tilde{m} unstable modes, we are able to find all the synchronization clusters on the desynchronized network.

The above analysis applies to only networks of perfect symmetries. In specific, for two oscillators to be synchronized, their elements in the eigenvectors of all the unstable modes should be exactly the same. From the point of view of network topology, this means that the two oscillators are identical and indistinguishable. This restricts seriously the application of the above analysis to the general complex networks, in which links are normally weighted and perfect symmetry does not exist in general. To cope with this problem, we loosen the requirement for cluster synchronization, and regard oscillators i and j as synchronized if the time-averaged error \delta {\tilde{x}}_{i,j}=⟨|{\mathit{x}}_i-{\mathit{x}}_j|⟩ between them is smaller to a small threshold. To investigate theoretically the CS behaviors in the general complex networks, now the question becomes: for a moderate coupling strength by which global network synchronization is impossible, can we predict the contents of the loosely defined synchronization clusters based on the information of the network coupling matrix and, furthermore, the ranges over which the synchronization clusters are observed in the parameter space of the uniform coupling strength? Noticing that {\mathit{x}}_i-{\mathit{x}}_j=\delta {\mathit{x}}_i-\delta {\mathit{x}}_j, we have from Eq. (7) the relation

(8)\delta {\tilde{x}}_{i,j}=⟨\left|{\displaystyle \sum _{k=2}^{\tilde{m}+1}(v_{k,i}-v_{k,j})\delta {\mathit{y}}_k(t)}\right|⟩.

For a deeply desynchronized network, the amplitudes of the unstable modes are approximately the same, i.e., |\delta \hat{\mathit{y}}(t)|\approx |\delta {\mathit{y}}_k(t)|. With this approximation, Eq. (8) can be simplified as

(9)\delta {\tilde{x}}_{i,j}=⟨\left|\delta {\tilde{e}}_{i,j}\delta \hat{\mathit{y}}(t)\right|⟩,

with

(10)\delta {\tilde{e}}_{i,j}={\displaystyle \sum _{k=2}^{\tilde{m}+1}|v_{k,i}-v_{k,j}|}

the distance between nodes i and j in the eigenvector space. To facilitate the analysis of CS in general complex networks (so that the thresholds used for finding the clusters are independent of the coupling strength), here we adopt the normalized distance \delta e_{i,j}=\delta {\tilde{e}}_{i,j}/\tilde{m} to characterize the topological difference between nodes i and j. As \delta \hat{\mathit{y}} is independent of the node index, we finally have

(11)\delta {\tilde{x}}_{i,j}\propto \delta e_{i,j}.

Eq. (11) is our main theoretical result, which tells how the synchronization degree of two oscillators (characterized by \delta {\tilde{x}}_{i,j}) can be inferred from the network topology (characterized by \delta e_{i,j}). We note that both \delta {\tilde{x}}_{i,j} and \delta e_{i,j} are dependent on the coupling strength, but the relation described by Eq. (11) holds for any coupling strength.

In applications, the above framework of eigenvector-based analysis is implemented as follows. The first step is to cacluate the eigenvalues \{{\lambda }_i{\}}_{i=1,\dots ,N} and the corresponding eigenvectors \{{\mathit{v}}_i{\}}_{i=1,\dots ,N} of the network coupling matrix \mathit{W}. The second step is to find the boundary of the stable regime of the MSF curve, i.e., the value of {\sigma }_c, which can be obtained by solving Eq. (4) numerically. The third step is to find the set of unstable modes, \{{\lambda }_k{\}}_{k=2,\dots ,\tilde{m}+1}, for a specific coupling strength of interest, and then calculate the matrix of eigenvector distance, \{\delta e_{i,j}{\}}_{N\times N}, according to Eq. (11). Finally, we find from the distance matrix the set of oscillators forming a cluster. For complex networks of perfect symmetries, the contents of each cluster are identified by the requirement \delta e_{i,j}=0 (i.e., the eigenvector distance between the oscillators within each cluster is 0); while for the general complex networks, the contents of each cluster are identified by requiring \delta e_{i,j} to be smaller than a predefined threshold.

3 Applications

We next apply the theoretical framework to investigate the CS behaviors in different network models, including small-size networks of perfect symmetries, large-size complex networks of community structures, and two empirical neural networks.

3.1 Small-size networks

We start with a toy network model of perfect symmetries. The network structure is shown in Fig.1(a1), which contains N=6 nodes and 9 links. Following Ref. [37], we set the weight of the link between nodes 2 and 6 as 1.5, and the same weight is arranged for the link between nodes 3 and 5. The weights of the other links are all set as unity. The coupling matrix of the network reads

(12)\mathit{W}=\left(\begin{array}{cccccc}-3& 1& 0& 1& 0& 1\\ 1& -3.5& 1& 0& 0& 1.5\\ 0& 1& -3.5& 1& 1.5& 0\\ 1& 0& 1& -3& 1& 0\\ 0& 0& 1.5& 1& -3.5& 1\\ 1& 1.5& 0& 0& 1& -3.5\end{array}\right).

The eigenvalues of \mathit{W} are

(13)({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5,{\lambda }_6)=(0,-2,-3,-4,-5,-6),

and the eigenvector matrix is

(14)\begin{array}{rl}\mathit{V}=& ({\mathit{v}}_1,{\mathit{v}}_2,{\mathit{v}}_3,{\mathit{v}}_4,{\mathit{v}}_5,{\mathit{v}}_6)\\ =& \left(\begin{array}{cccccc}0.41& -0.41& 0.58& 0& 0.58& 0\\ 0.41& -0.41& -0.29& 0.5& -0.29& 0.5\\ 0.41& 0.41& -0.29& 0.5& 0.29& -0.5\\ 0.41& 0.41& 0.58& 0& -0.58& 0\\ 0.41& 0.41& -0.29& -0.5& 0.29& 0.5\\ 0.41& -0.41& -0.29& -0.5& -0.29& -0.5\end{array}\right).\end{array}

We adopt the chaotic Lorenz oscillator to describe the nodal dyamimcs [67]. The dynamics of the network reads

(15)\{\begin{array}{l}\dot{x_i}=\alpha (y_i-x_i),\\ \dot{y_i}=x_i(r-z_i)-y_i+\epsilon \sum w_{ij}{x}_j,\\ \dot{z_i}=x_i{y}_i-b{z}_i.\end{array}

Here, i,j=1,\dots ,N are the oscillator indices, and the coupling function is chosen as \mathit{H}(\mathit{x})=[0,x,0{]}^{\mathrm{T}}, i.e., the x variable is coupled to the y variable. The parameters of the oscillators are all set as (\alpha ,r,b)=(10,35,8/3), with which the oscillators present chaotic motion (the largest Lyapunov exponent is about 1.05). To characterize the synchronization relationship among the oscillators, we set oscillator 2 as the reference oscillator and calculate the synchronization error \delta x_i=⟨|x_i-x_2|{⟩}_{\mathrm{T}}, where ⟨\cdot {⟩}_{\mathrm{T}} denotes the time average function. In simulations, the system dynamics is evolved numerically by the fourth-order Runge−Kutta method with the time step \delta t=1\times {10}^{-2}. In calculating \delta x_i, the system is firstly evolved for a transient period of T_0=1\times {10}^4 so as to remove the impact of the initial conditions, and the results are averaged over a period of T=1\times {10}^4.

Following the MSF formalism, we first calculate numerically the variation of the largest Lyapunov exponent, \mathrm{\Lambda }, of the master equation [i.e., Eq. (4)] with respect to the generic coupling strength, \sigma. The results show that when \sigma >{\sigma }_c\approx 8.3. The critical couplings where the modes of {\lambda }_2, {\lambda }_3, {\lambda }_4, {\lambda }_5 and {\lambda }_6 become unstable thus are {\epsilon }_2={\epsilon }_c\approx 4.2, {\epsilon }_3\approx 2.8, {\epsilon }_4\approx 2.1, {\epsilon }_5\approx 1.7 and {\epsilon }_6\approx 1.4, respectively. By the framework of eigenvector-based analysis, we next predict the possible CS states to be observed in the process of network desynchronization, starting from the global synchronization state generated at a larger coupling strength \epsilon >{\epsilon }_c={\epsilon }_2. When \epsilon crosses {\epsilon }_2 from above, the mode of {\lambda }_2 will be unstable. According to the analysis presented in Section 2, the specific form of the CS state to be generated from the uniform background of global synchronization is determined by the eigenvector {\mathit{v}}_2. As shown in Eq. (14) (the 2nd column of the eigenvector matrix), we have v_{2,1}=v_{2,2}=v_{2,6}=-0.41 and v_{2,3}=v_{2,4}=v_{2,5}=0.41. The eigenvector distances therefore are \delta e_{i,j}=0 for i,j\in \{1,2,6\} or i,j\in \{3,4,5\}. That is, the theory predicts that when only the mode of {\lambda }_2 is destabilized, the network nodes are synchronized into two clusters, {\mathit{C}}_1=\{1,2,6\} and {\mathit{C}}_2=\{3,4,5\}. Decreasing further \epsilon, the mode of {\lambda }_3 will be destabilized when \epsilon crosses {\epsilon }_3. By checking the eigenvector {\mathit{v}}_3 Eq. (14) (the 3rd column of the eigenvector matrix), we see that v_{3,1}=v_{3,4}=0.58, v_{3,2}=v_{3,3}=v_{3,5}=v_{3,6}=-0.29. As both the modes of {\lambda }_2 and {\lambda }_3 are unstable [i.e., \tilde{m}=2 in Eq. (10)], the eigenvector distances between the nodes thus are \delta e_{2,6}=0, \delta e_{3,5}=0, and \delta e_{i,j}>0 for other node pairs. Hence, the theory predicts that in this case the network contains two synchronization pairs, {\mathit{C}}_1=\{2,6\} and {\mathit{C}}_2=\{3,5\}, and two desynchronized nodes, {\mathit{C}}_3=\{1\} and {\mathit{C}}_4=\{4\}. Finally, as \epsilon falls below {\epsilon }_4, the mode associated with {\lambda }_4 will be unstable. As v_{4,1}=v_{4,4}=0, v_{4,2}=v_{4,3}=0.5 and v_{4,5}=v_{4,6}=-0.5 [see Eq. (14), the 4th column in the eigenvector matrix], the eigenvector distance between any pair of nodes is larger than 0. That is, the theory predicts that when the mode of {\lambda }_4 is unstable (), no synchronization is observed between the oscillators, and the network is completely desynchronized.

The above predictions are well verified by numerical simulations. Shown in Fig.1(a2) is the variation of the synchronization errors, \delta x_i=⟨|x_i-x_2|{⟩}_{\mathrm{T}}, with respect to the uniform coupling strength, \epsilon. We see that the network is globally synchronized when \epsilon >{\epsilon }_c^{\prime }\approx 4.2. In the range of \epsilon \in ({\epsilon }_3^{\prime }\approx 2.7,{\epsilon }_c^{\prime }), the oscillators are organized into two synchronizastion clusters, {\mathit{C}}_1=\{1,2,6\} and {\mathit{C}}_2=\{3,4,5\}. In the range of \epsilon \in ({\epsilon }_4^{\prime }\approx 1.8,{\epsilon }_3^{\prime }), the network contains two synchronization pairs, (2,6) and (3,5), and two desynchronized nodes (oscillators 1 and 4). Finally, when , no synchronization is observed between the oscillators, and the network is completely desynchronized. The critical couplings of the CS states, as well as the order of their emergences in the process of network desynchronization, are in good agreement with the theoretical predictions.

The second network model of perfect symmetries employed in our studies is the Nepal power-grid [38]. The network structure of the Nepal power-grid is plotted in Fig.1(b1), which consists of N=15 nodes and L=62 unweighted, non-directed links. Still, we describe the nodal dynamics by the chaotic Lorenz oscillator and couple the oscillators through the function \mathit{H}(\mathit{x})=[0,x,0{]}^{\mathrm{T}}. The eigenvalues of the network coupling matrix are

(16)\begin{array}{rl}\{{\lambda }_i{\}}_{i=1,\dots ,15}=& (0,-0.94,-3.10,-8,-8,-8,-8,-9,\\ & -9,-9,-9,-9.94,-14,-14,-14.03),\end{array}

and the eigenvector matrix, \mathit{V}=({\mathit{v}}_1,\dots ,{\mathit{v}}_{15}), reads

(17)\left(\begin{array}{ccccccccccccccc}0.26& -0.11& -0.32& 0.71& -0.12& 0& -0.53& 0& 0& 0& 0& 0.01& 0& 0& 0.14\\ 0.26& -0.11& -0.32& -0.71& -0.12& 0& -0.53& 0& 0& 0& 0& 0.01& 0& 0& 0.14\\ 0.26& -0.11& -0.32& 0& 0.48& 0.71& 0.27& 0& 0& 0& 0& 0.01& 0& 0& 0.14\\ 0.26& -0.11& -0.32& 0& 0.48& -0.71& 0.27& 0& 0& 0& 0& 0.01& 0& 0& 0.14\\ 0.26& -0.11& -0.32& 0& -0.72& 0& 0.53& 0& 0& 0& 0& 0.01& 0& 0& 0.14\\ 0.26& -0.07& 0.01& 0& 0& 0& 0& 0& 0& 0& 0& -0.03& 0.72& -0.39& -0.51\\ 0.26& -0.07& 0.01& 0& 0& 0& 0& 0& 0& 0& 0& -0.03& -0.02& 0.82& -0.51\\ 0.26& -0.07& 0.01& 0& 0& 0& 0& 0& 0& 0& 0& -0.03& -0.70& -0.43& -0.51\\ 0.26& -0.05& 0.29& 0& 0& 0& 0& 0.74& 0.21& 0.43& 0.14& 0.17& 0& 0& 0.14\\ 0.26& -0.05& 0.29& 0& 0& 0& 0& -0.15& -0.19& -0.38& 0.77& 0.17& 0& 0& 0.14\\ 0.26& -0.05& 0.29& 0& 0& 0& 0& 0.20& 0.10& -0.69& -0.53& 0.17& 0& 0& 0.14\\ 0.26& -0.05& 0.29& 0& 0& 0& 0& -0.20& -0.74& -0.36& -0.30& 0.17& 0& 0& 0.14\\ 0.26& -0.05& 0.29& 0& 0& 0& 0& -0.59& 0.61& 0.28& -0.09& 0.17& 0& 0& 0.14\\ 0.26& 0.06& 0.23& 0& 0& 0& 0& 0& 0& 0& 0& -0.92& 0& 0& 0.17\\ 0.26& 0.95& -0.11& 0& 0& 0& 0& 0& 0& 0& 0& 0.10& 0& 0& -0.01\end{array}\right).

While CS in the Nepal power-grid network has been studied in the literature [38, 44], the existing studies rely on a priori knowledge of the network symmetries. In what follows, we are going to show that, without knowing the network symmetries, the CS states can be well analyzed by the framework of eigenvector-based analysis.

We now analyze the CS states emerged in the process of desynchronization transition for the Nepal power-grid network, starting from the global synchronization state. The critical coupling for the whole network to be synchronized is {\epsilon }_c={\epsilon }_2=-{\sigma }_c/{\lambda }_2\approx 8.8. When \epsilon crosses {\epsilon }_c from above, the mode of {\lambda }_2 will be unstable and, according to the method of eigenvector analysis, the synchronization clusters can be identified from the elements of {\mathit{v}}_2 (the 2nd column of the matrix \mathit{V}). As v_{2,1}=\dots =v_{2,5}=-0.11, v_{2,6}=v_{2,7}=v_{2,8}=-0.07 and v_{2,9}=\dots =v_{2,13}=-0.05, we have three clusters in this case, {\mathit{C}}_1=(1,\dots ,5), {\mathit{C}}_2=(9,\dots ,13) and {\mathit{C}}_3=(6,7,8), with the eigenvector distance being 0 for nodes within the same cluster. The contents of the clusters are kept unchanged when the mode of {\lambda }_3 becomes unstable, as the partition of the elements of {\mathit{v}}_3 is identical to that of {\mathit{v}}_2. However, when the mode of {\lambda }_4 is destabilized (), the elements of {\mathit{v}}_4 (the 4th column of the matrix \mathit{V}) suggest that the 1st cluster will be broken, while the 2nd and 3rd clusters hold still. That is, the theory predicts that in the range of \epsilon \in ({\epsilon }_4,{\epsilon }_c), three synchronization clusters will be observed on the network, with the contents of the 1st, 2nd and 3rd clusters being defined by {\mathit{C}}_1, {\mathit{C}}_2 and {\mathit{C}}_3, respectively. By decreasing \epsilon further, the eigenvector distances of the nodes are kept unchanged till the critical coupling {\epsilon }_8=0.92 is met, where the mode of {\lambda }_8 becomes unstable. By checking the elements of {\mathit{v}}_8 in Eq. (17) (the 8th column), we see that elements of the nodes in cluster 2 are non-identical, suggesting the breaking of the 2nd cluster at {\epsilon }_8. As the 1st cluster has already been broken at {\epsilon }_4, only cluster 3 survives in this case. Finally, at the critical coupling {\epsilon }_{13}\approx 0.6, the mode of {\lambda }_{13} is destabilized and, according to the elements of {\mathit{v}}_{13}, the eigenvector distances between nodes in the 3rd cluster will be all larger than 0. As such, the 3rd cluster will be desynchronized at {\epsilon }_{13}. Summarizing up the above results, the theory predicts that: (i) the network contains three synchronization clusters, {\mathit{C}}_{1,2,3}, when \epsilon \in ({\epsilon }_4,{\epsilon }_c); (ii) two synchronization clusters, {\mathit{C}}_{2,3}, are formed when \epsilon \in ({\epsilon }_8,{\epsilon }_4); (iii) a single synchronization cluster, {\mathit{C}}_3, is left when \epsilon \in ({\epsilon }_{13},{\epsilon }_8); and (iv) no synchronization cluster is observed when .

By the approach of numerical simulations, we plot in Fig.1(b2) the variation of the cluster-based synchronization error, \delta x^c, with respect to the coupling strength, \epsilon. Here, cluster-based synchronization error is defined as \delta x_m^c=\sum _{i,j\in C_m}⟨|x_i-x_j|{⟩}_{\mathrm{T}}/n_m(n_m-1), with m=1,2,3 the cluster index and n_m the number of nodes in cluster m. Clearly, cluster m is synchronized when \delta x_m^c=0. Fig.1(b2) shows that: (i) three synchronization clusters are generated when \epsilon >{\epsilon }_4^{\prime }\approx 0.9, with the contents of the clusters identical to the ones predicted by the theory; (ii) clusters 2 and 3 are synchronized when \epsilon \in ({\epsilon }_8^{\prime }\approx 0.7,{\epsilon }_4^{\prime }); (iii) only cluster 3 is synchronized when \epsilon \in ({\epsilon }_{13}^{\prime }\approx 0.32,{\epsilon }_8^{\prime }); and (iv) no cluster is synchronized when . We see that, though the framework of eigenvector-based analysis cannot predict precisely the critical couplings where the clusters become unstable, it does predict accurately the patterns of the CS states and, additionally, the sequence of these states in the process of network desynchronization.

3.2 Complex community networks

So far the framework of eigenvector-based analysis has been utilized to investigate the CS states in only small-size networks of perfect symmetries. It remains as not clear whether the framework can be applied to general complex networks in which perfect symmetries do not exist. We next check the efficacy of the proposed framework in analyzing the CS behaviors of complex community networks – the type of network that is representative to many realistic systems [64]. In constructing the complex community network, we first divide the N nodes into M communities of equal size, and then connect nodes within each community by the probability p_1 and nodes from different communities by the probability p_2, with . To demonstrate the generality of the theoretical framework, here we adopt the chaotic Hindmarsh−Rose (HR) oscillator as the nodal dynamics [68]. The dynamics of isolated HR oscillators is governed by equations [\dot{x},\dot{y},\dot{z}]=[y-x^3+3x^2-z+I,1-5x^2-y,r[s(x+1.6)-z]{]}^{\mathrm{T}}. The parameters of the oscillators are set as (r,s,I)=(6\times {10}^{-3},4,3.2), by which the oscillators present the chaotic motion (the largest Lyapunov exponent is about 1.3\times {10}^{-2}). The oscillators are coupled by gap junctions, with the coupling function being described as \mathit{H}(\mathit{x})=[x,0,0{]}^{\mathrm{T}}. To make the network dynamic stable, here we adopt the strategy of normalized couplings [69, 70], i.e., w_{ij}=a_{ij}/k_i, with \{a_{ij}{\}}_{N\times N} the binary adjacency matrix capturing the network structure and k_i=\sum a_{ij} the degree of node i. The diagonal elements of \mathit{W} are all set as −1. Though the network links are weighted and directed, the eigenvalues of the coupling matrix are all real [69, 70].

We start by analyzing the CS behaviors in a trivial complex community network consisting of N=90 nodes and M=3 communities. To distinguish the community feature, here we set the connecting probabilities as p_1=0.9 (for the intra-community links) and p_2=0.1 (for the inter-community links). The network structure is plotted in Fig.2(a1), in which nodes are ordered according to the communities, {\text{embassy}\text{C}}_1=(1,\dots ,30), {\text{embassy}\text{C}}_2=(31,\dots ,60) and {\text{embassy}\text{C}}_3=(61,\dots ,90). As the links are randomly added among the nodes, there is no perfect symmetry in the network structure, which has been confirmed by the algorithm developed from computational group theory [38]. While previous studies have shown that the synchronization clusters in community networks can be inferred from the eigenvector associated with the leading mode [65, 66], we are going to demonstrate in the following that the framework of eigenvector-based analysis is capable of offering more information about the CS behaviors. In applying the theoretical framework, we first calculate from the network coupling matrix the eigenvalues and the associated eigenvectors. The three largest non-trivial eigenvalues are {\lambda }_2=-0.12, {\lambda }_3=-0.14 and {\lambda }_4=-0.88. By simulating the master equation of the HR oscillator [see Eq. (4)], we have for \sigma >{\sigma }_c\approx 0.94. The critical coupling for the modes of {\lambda }_2, {\lambda }_3 and {\lambda }_4 therefore are {\epsilon }_2\approx 7.8 (which is also the critical coupling for global synchronization), {\epsilon }_3\approx 6.7 and {\epsilon }_4\approx 1.1. Shown in Fig.2(a2) are the elements of the eigenvectors {\mathit{v}}_2, {\mathit{v}}_3 and {\mathit{v}}_4. We see that the elements of {\mathit{v}}_2 are clearly seperated into two groups, v_{2,i}\approx -0.15 for i=1,\dots ,30 and v_{2,j}\approx 0.15 for j=31,\dots ,90. In terms of the eigenvector distance defined by Eq. (10), this means that \delta e_{ij}\approx 0 only for i,j\in {\text{embassy}\text{C}}_1 or i,j\in {\text{embassy}\text{C}}_2\cup {\text{embassy}\text{C}}_3. To have a global picture on the distribution of eigenvector distances, we plot in Fig.2(b1) the matrix \{\delta e_{ij}{\}}_{i,j=1,\dots ,N}, for the coupling strength \epsilon =7.6 (by which the mode of {\lambda }_2 is unstable). Clearly, the oscillators are divided into two groups. According to the method of eigenvector-based analysis, we thus predict that as \epsilon decreases and crosses {\epsilon }_2, the oscillators are synchronized into two clusters, {\mathit{C}}_1=(1,\dots ,30), {\mathit{C}}_2=(31,\dots ,90). Decreasing \epsilon to {\epsilon }_3, the mode of {\lambda }_3 will be unstable. By checking the elements of {\mathit{v}}_3 in Fig.2(a2), it is found that the elements are divided into three distinct groups: v_{3,i}\approx 0 for i=1,\dots ,30, v_{3,j}\approx -0.14 for j=31,\dots ,60 and v_{3,l}\approx 0.14 for l=61,\dots ,90. The matrix of eigenvector distance for this case is plotted in Fig.2(c1), which shows that the nodes are divided into three distinct groups, {\mathit{C}}_1=(1,\dots ,30), {\mathit{C}}_3=(31,\dots ,60) and {\mathit{C}}_3=(61,\dots ,90), with \delta e\approx 0 for nodes within the same group but is larger than zero for nodes from different groups. The theory thus predicts that as \epsilon decreases from {\epsilon }_3, the giant cluster formed in the previous stage [\epsilon \in ({\epsilon }_3,{\epsilon }_2)] will be broken into two small-size clusters, and now there are three synchronization clusters on the network. Finally, as \epsilon crosses {\epsilon }_4, the mode of {\lambda }_4 will be unstable. As depicted in Fig.2(a2), the elements of {\mathit{v}}_4 are randomly distributed and, as a result of this, no eigenvector distance is close to 0 and no clear structure is observed in the matrix [see Fig.2(d1)], indicating that the network is fully desynchronized in the case.

The above predictions are well verified by numerical simulations. Setting \epsilon =7.6, we plot in Fig.2(b2) the matrix of the synchronization error \delta x_{i,j}={⟨|x_i-x_j|⟩}_{\mathrm{T}}. We see that, in consistent with the theoretical prediction shown in Fig.2(b1), the oscillators are synchronized into two clusters, {\mathit{C}}_1={\text{embassy}\text{C}}_1 and {\mathit{C}}_2={\text{embassy}\text{C}}_2\cup {\text{embassy}\text{C}}_3. Plotted in Fig.2(c2) is the synchronization-error matrix for the coupling strength \epsilon =6, which corresponds to the situation when the modes of {\lambda }_2 and {\lambda }_3 are unstable. We see that, in agreement with the predictions shown in Fig.2(c1), the oscillators are synchronized three clusters, {\mathit{C}}_1={\text{embassy}\text{C}}_1, {\mathit{C}}_2={\text{embassy}\text{C}}_2 and {\mathit{C}}_3={\text{embassy}\text{C}}_3. Finally, we plot in Fig.2(d2) the synchronization-error matrix for the coupling strength \epsilon =0.5. In consistent with the predictions in Fig.2(d1), we see that there is no clear synchronization cluster and the network is fully desynchronized.

We next investigate the CS behaviors in a large-size community network of hierarchical structures. The network structure is plotted in Fig.3(a1), which consists of M=6 communities and each community contains n=30 nodes. The network is generated by replacing each node in the 6-node network plotted in Fig.1(a) with a community. That is, the inter-connections are established only between the paired communities. To capture the weighted feature of the 6-node network, the inter-connecting probability between communities 2 and 6 (and also between communities 3 and 5) is set as p_2=0.15, while the inter-connecting probability for other paired communities is set as p_2=0.1. The connecting probability between nodes within each community is set uniformly as p_1=0.9. Still, the couplings are normalized and the nodal dynamics is described by the chaotic HR oscillator. The network nodes are ordered according to the partition of the communities, {\text{embassy}\text{C}}_m=((m-1)n+1,\dots ,mn), with m=1\dots ,M the community index. The three non-trivial largest eigenvalues of the network coupling matrix are {\lambda }_2=-8.2\times {10}^{-2}, {\lambda }_3=-0.136 and {\lambda }_4=-0.167, and the corresponding critical couplings are {\epsilon }_2\approx 11.46, {\epsilon }_3\approx 6.9 and {\epsilon }_4\approx 5.6. Shown in Fig.3(a2) are the elements of the eigenvectors {\mathit{v}}_2, {\mathit{v}}_3 and {\mathit{v}}_4. As \epsilon decreases and crosses {\epsilon }_2, the mode of {\lambda }_2 will be unstable. A check of the elements of {\mathit{v}}_2 in Fig.3(a2) shows that the nodes are divided into two groups, {\mathit{C}}_1={\text{embassy}\text{C}}_1\cup {\text{embassy}\text{C}}_2\cup {\text{embassy}\text{C}}_6 (in which v_{2,i}\approx 0.075) and {\mathit{C}}_2={\text{embassy}\text{C}}_3\cup {\text{embassy}\text{C}}_4\cup {\text{embassy}\text{C}}_5 (in which v_{2,i}\approx -0.075). The corresponding matrix of eigenvector distance is plotted in Fig.3(b1), which, according to the framework of eigenvector-based analysis, implies that the oscillators will be synchronized into two clusters, {\mathit{C}}_1 and {\mathit{C}}_2, when the mode of {\lambda }_2 is destabilized. As \epsilon crosses {\epsilon }_3, the mode of {\lambda }_3 will be also unstable. By checking the elements of both {\mathit{v}}_2 and {\mathit{v}}_3 in Fig.3(a2), we see that the nodes are divided into four distinct groups, {\mathit{C}}_1={\text{embassy}\text{C}}_1, {\mathit{C}}_2={\text{embassy}\text{C}}_2\cup {\text{embassy}\text{C}}_6, {\mathit{C}}_3={\text{embassy}\text{C}}_4 and {\mathit{C}}_4={\text{embassy}\text{C}}_3\cup {\text{embassy}\text{C}}_5. The matrix of eigenvector distance is plotted in Fig.3(c1), which suggests that the oscillators will be synchronized into four clusters when the modes of {\lambda }_2 and {\lambda }_3 are unstable. Finally, when the mode of {\lambda }_4 is unstable (as \epsilon crosses {\epsilon }_4), the elements of {\mathit{v}}_2, {\mathit{v}}_3 and {\mathit{v}}_4 in Fig.3(a2) show that the nodes are divided into six distinct groups, {\mathit{C}}_m={\text{embassy}\text{C}}_m, with m=1,\dots ,6. The corresponding matrix of eigenvector distance is plotted in Fig.3(d1).

The above analysis thus predicts the following scenario of network desynchronization. With the decrease of the coupling strength, the whole network is firstly separated into two larger synchronization clusters (at {\epsilon }_2), with one cluster formed by communities 1, 2 and 6 and the other one formed by communities 3, 4 and 5. Then, as \epsilon crosses {\epsilon }_3, each cluster is broken into two small-size clusters, and there are in total 4 synchronization clusters on the network. In this stage, communities 2 and 6 are still synchronized, and also communities 3 and 5. Finally, as \epsilon falls below {\epsilon }_4, community 2 (3) is desynchronized from community 6 (5), and there are 6 synchronization clusters on the network in total, with each cluster corresponding to one community. These predictions are well verified by numerical simulations, as depicted in Fig.3(b2) [corresponding to the prediction in Fig.3(b1)], Fig.3(c2) [corresponding to the prediction in Fig.3(c1)] and Fig.3(d2) [corresponding to the prediction in Fig.3(d1)].

3.3 Empirical neural networks

We finally employ the framework of eigenvector-based analysis to investigate the CS behaviors in two empirical neural networks. The first example is the cortical network of the cat brain [60, 71]. The network structure is plotted in Fig.4(a1), which are constructed by N=53 nodes (cortex areas) and L=830 links (fiber connections). According to their functions, the cortex areas are divided into four divisions of variant size: 16 areas in the visual division [{\text{embassy}\text{C}}_V=(1,\dots ,16)], 7 areas in the auditory division [{\text{embassy}\text{C}}_A=(17,\dots ,23)], 16 areas in the somatomotor division [{\text{embassy}\text{C}}_{SM}=(24,\dots ,39)], and 14 areas in the frontolimbic division [{\text{embassy}\text{C}}_{FL}=(40,\dots ,53)]. (We note that the order of the neurons does not affect the synchronization behaviors of the network.) Still, the strategy of normalized couplings is adopted in constructing the network coupling matrix, and the chaotic HR oscillator is adopted to describe the nodal dynamics. The two largest non-trivial eigenvalues of the network coupling matrix are {\lambda }_2=-0.395 and {\lambda }_3=-0.432. The critical couplings of the two modes are {\epsilon }_2\approx 2.38 (which is also the critical coupling for global synchronization) and {\epsilon }_3\approx 2.18. The elements of the eigenvenctors {\mathit{v}}_2 and {\mathit{v}}_3 are plotted in Fig.4(a2). We see that, compared with the eigenvectors of the artificial networks [e.g., the results in Fig.2(a2) and Fig.3(a2)], the elements of {\mathit{v}}_2 and {\mathit{v}}_3 are not clearly grouped.

Shown in Fig.4(b1) is the matrix of eigenvector distance calculated for the coupling strength \epsilon =2.2, by which only the mode of {\lambda }_2 is unstable. We see that the nodes are partitioned into three distinct clusters, {\mathit{C}}_1={\text{embassy}\text{C}}_V\cup {\text{embassy}\text{C}}_{FL}, {\mathit{C}}_2={\text{embassy}\text{C}}_A and {\mathit{C}}_3={\text{embassy}\text{C}}_{SM}. Shown in Fig.4(c1) is the matrix of eigenvector distance for \epsilon =1.95, by which both the modes of {\lambda }_2 and {\lambda }_3 are unstable. We see that in this case, in agreement with the neural divisions, the nodes are partitioned into four clusters, {\mathit{C}}_1={\text{embassy}\text{C}}_V, {\mathit{C}}_2={\text{embassy}\text{C}}_A, {\mathit{C}}_3={\text{embassy}\text{C}}_{SM} and {\mathit{C}}_4={\text{embassy}\text{C}}_{FL}. To check the accuracy of the predictions, we plot in Fig.4(b2) and Fig.4(c2) the matrix of synchronization error, \{\delta x_{i,j}{\}}_{i,j=1,\dots ,N}, for the coupling strengths \epsilon =2.2 and 1.95, respectively. We see that the synchronization patterns predicted by the framework of eigenvector-based analysis [Fig.4(b1) and (c1)] are in good agreement with the ones obtained by model simulations [Fig.4(b2) and (c2)].

The second empirical neural network we consider is the cerebral cortex of the human brain [63, 72, 73], which contains N=989 nodes (496 nodes in the right hemisphere and 493 nodes in the left hemisphere) and L=17865 links. According to the cytoarchitecture and functional parcellation, the nodes are partitioned into M=64 cortical regions (communities). Still, we adopt the normalized coupling strategy and use the chaotic HR oscillator to describe the nodal dynamics. The leading non-trivial eigenvalues of the network coupling matrix are ({\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5,{\lambda }_6)=(-3.7\times {10}^{-2},-7.7\times {10}^{-2},-0.1,-0.12,-0.16,-0.18), and the corresponding critical couplings are about ({\epsilon }_2,{\epsilon }_3,{\epsilon }_4,{\epsilon }_5,{\epsilon }_6)=(25,12,9.4,7.8,5.9,5.2). Plotted in Fig.5(a1) are the synchronization patterns predicted by the framework of eigenvector-based analysis for the coupling strength \epsilon =15 (by which only the mode of {\lambda }_2 is unstable). In plotting Fig.5(a1), we first calculate for each cortical region the averaged eigenvector distance {⟨\delta e⟩}_m=\sum _{i,j\in {\text{embassy}\text{C}}_m}\delta e_{i,j}/[n_m(n_m-1)], with {\text{embassy}\text{C}}_m denoting the set of nodes in the mth region, n_m is the size of the mth region, and \delta e_{i,j} is the eigenvector distance between nodes i and j. The nodes within the mth region are regarded as synchronizable if {⟨\delta e⟩}_m is smaller than the threshold {⟨\delta e⟩}_c=0.01. By doing so, the M=64 regions are divided into two groups, the synchronizable and non-synchronizable groups. The nodes in the non-synchronizable regions are shown in grey, which are excluded from the further analysis. Having identified the set of synchronizable regions, we next partition the regions into clusters based on the eigenvector distance between them. Here the eigenvector distance between regions m and m^{\prime } is defined as {⟨\delta e⟩}_{m,m^{\prime }}=\sum _{i\in {\text{embassy}\text{C}}_m;j\in {\text{embassy}\text{C}}_{m^{\prime }}}\delta e_{i,j}/(n_m{n}_{m^{\prime }}), with {\text{embassy}\text{C}}_m and {\text{embassy}\text{C}}_{m^{\prime }} denoting the sets of nodes in the mth and m^{\prime }th regions, respectively, and \delta e_{i,j} is the eigenvector distance betwteen nodes i and j. Finally, regions with are considered as synchronizable and nodes in these regions are marked by the same color. We see in Fig.5(a1) that most of the nodes are colored, indicating a higher synchronization degree of the whole network. Meanwhile, the colored nodes are organized into clusters of different sizes and specific spatial distributions. By the same procedure, we plotted in Fig.5(b1) and (c1) the synchronization patterns predicted by the theoretical framework for the coupling strengths \epsilon =8.5 (by which the modes of {\lambda }_2, {\lambda }_3 and {\lambda }_4 are unstable) and \epsilon =5.0 (by which the modes of {\lambda }_{2,\dots ,6} are unstable). We see that, as more modes become unstable, more nodes are shown in grey, indicating the degraded synchronization performance of the whole network. In the meantime, with the decrease of the coupling strength, the sizes of the clusters are shrunk gradually.

The above predictions are well verified by simulations. Setting \epsilon =15, we plot in Fig.5(a2) the synchronization pattern obtained by model simulations. In plotting Fig.5(a2), we first calculate the node-averaged synchronization error of each region, {⟨\delta x⟩}_m=\sum _{i,j\in {\text{embassy}\text{C}}_m}\delta x_{i,j}/[n_m(n_m-1)], with \delta x_{i,j}={⟨|x_i-x_j|⟩}_{\mathrm{T}} is the time-averaged synchronization error between nodes i and j. The nodes within the mth region are regarded as synchronizable if {⟨\delta x⟩}_m is smaller than the threshold value {⟨\delta x⟩}_c=0.035. Still, nodes in the non-synchronizable regions are shown in grey. We then partition the synchronizable regions into clusters based on the synchronization error {⟨\delta x⟩}_{m,m^{\prime }}=\sum _{i\in {\text{embassy}\text{C}}_m;j\in {\text{embassy}\text{C}}_{m^{\prime }}}\delta x_{i,j}/[n_m{n}_{m^{\prime }}]. The mth and m^{\prime }th regions are regarded as synchronized if , and nodes in the synchronized regions are marked by the same color. We see that the synchronization pattern obtained by simulations [Fig.5(a2)] is in good agreement with the pattern predicted by the theory [Fig.5(a1)]. The good agreement between simulation and prediction is also observed for the coupling strengths \epsilon =8.5 and \epsilon =5.0, as depicted in Fig.5(b2) [corresponding to the prediction in Fig.5(b1)] and Fig.5(c2) [corresponding to the prediction in Fig.5(c1)]. We note that by changing the thresholds {⟨\delta e⟩}_c and {⟨\delta x⟩}_c, the sizes of the clusters will be changed, but the consistency between the theoretically predicted and numerically obtained patterns is kept unchanged, which has been checked by additional simulations (not shown).

4 Discussion and conclusions

Whereas the significant impacts of network structure on synchronization have been well recognized and demonstrated in the literature, attention has been mainly focused on the roles of the extreme modes of the network coupling matrix. Specifically, for complex networks of linearly coupled identical chaotic oscillators, the MSF formalism suggests that the synchronizability of the networks can be well characterized by the extreme eigenvalues of the network coupling matrix [5], i.e., the largest non-trivial eigenvalue ({\lambda }_2) or the ratio between the largest and smallest non-trivial eigenvalues (\rho ={\lambda }_N/{\lambda }_2). For this reason, a large body of studies on network synchronization has been focusing on the dependence of the extreme eigenvalues on the network structure or the coupling strategy [14, 15], with the impacts and roles of the other modes being largely overlooked. While analysis based on the extreme eigenvalues is efficient and convenient in applications, it applies to only the special case of global synchronization state but not the CS states appeared in the transition to global synchronization. Our studies in the present work show that, to fully characterize the synchronization behaviors on complex networks and the CS states in synchronization transition in particular, it is equally important to consider the impacts and roles of the modes with moderate eigenvalues.

Our present work also highlights the significance of eigenvectors in exploring the CS behaviors of complex networks. Previously, eigenvectors of the network coupling matrix have been exploited to investigate the synchronization dynamics of networked systems, but the studies are limited to slightly desynchronized networks or small-size network motifs [74–76]. In Ref. [74], the authors studied the synchronization behavior of chaotic oscillators in large-size complex networks at the boundary of global synchronization, and found that the stabilities of the oscillators can be well predicted from the elements of the eigenvector associated with the unstable mode. In Ref. [75], the authors investigated the dynamical patterns on small-size networks of up to 5 nodes, and found that, when the eigenvectors of the network coupling matrix satisfy some special properties, the dynamical patterns can be well inferred from the network topology and the local dynamics. In Ref. [76], the authors demonstrated that, by the information of the eigenvector of the leading mode, all synchronization clusters in a complex network can be identified without knowing the network symmetries and, in addition, their sequence in the process of synchronization transition can be predicted. Compared to the existing studies, the framework proposed in our current study is featured by the capability of analyzing the synchronization patterns (CS states) formed in deeply desynchronized networks (with the number of unstable modes unlimited) and, more importantly, in large-size complex networks of non-perfect symmetries. The latter makes the two different approaches currently employed in CS studies, namely the symmetry-based (in which perfect network symmetry is required and the clusters are defined on complete synchronization) [34–50] and community-based approaches (in which no perfect network symmetry is required and the synchronization clusters are loosely defined) [60, 61, 63, 65, 66], unified within the same framework, thereby paving a way to the exploration of the synchronization patterns in real-world systems.

The key result we have obtained in the theoretical analysis is given by Eq. (11), which reveals how the synchronization patterns are connected to the eigenvectors of the network coupling matrix. In obtaining Eq. (11), we have made the assumption that the statistical properties of the coupled oscillators in CS are close to that of the isolated oscillator. The purpose of this assumption (approximation) is to make the same master equation, i.e., Eq. (4), applicable to all the perturbation modes, so that the critical couplings where the modes become unstable can be analytically estimated ({\epsilon }_i=-{\sigma }_c/{\lambda }_i). While the validity of this approximation has been verified by numerical simulations in our current work, there are cases where the approximation is invalid [78]. In these invalid cases, we are still able to identify the synchronization patterns (the CS states) according to the eigenvector matrix, but are unable to predict analytically the critical couplings where the patterns are generated. (In cases like this, the stability of each CS state should be evaluated individually, e.g., by calculating the largest conditional Lyapunov exponent [38, 49].)

An alternative approach to analyzing the CS behaviors in networked chaotic oscillators could be the energy-balance-based method [79]. Like the method of Lyapunov function, the CS states stand as the local minima of the energy landscape. Our study suggests that these local minima are determined by the network symmetries and can be inferred from the eigenvectors of the network coupling matrix. A systematic analysis of the application of the energy-balance-based method to CS behaviors is out of the scope of our current study.

To summarize, exploiting the eigenvectors of the network coupling matrix, we have proposed a new theoretical framework to explore the CS behaviors in general complex networks. The new framework, which requires no prior knowledge of the network symmetries, is able to predict not only all the CS states generated in the process of synchronization transition, but also the critical couplings where the CS states are generated and the sequence of the CS states in the transition. The efficacy and efficiency of the framework have been verified by a variety of network models of identical chaotic oscillators, including small-size artificial networks of perfect symmetries, large-size complex networks of distinct community structures, and large-size empirical neural networks. Our studies highlight the importance of eigenvectors in exploring CS in complex networks, and the framework we have proposed provides a powerful tool to the exploration of synchronization patterns in real-world complex systems.

References

[1] Y.Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer, Berlin, 1984

[2] A.T. Winfree, Timing of Biological Clocks, W H Freeman & Co, 1987

[3] A.S. PikovskyM.G. RosenblumJ.Kurths, Synchronization: A Universal Concept in Nonlinear Science, Cambridge University Press, Cambridge, 2001

[4] S.Strogatz, Sync: The Emerging Science of Spontaneous Order, Hyperion, New York, 2003

[5] L. M. Pecora , T. L. Carroll . Master stability functions for synchronized coupled systems. Phys. Rev. Lett., 1998, 80(10): 2109

[6] G. Hu , J. Z. Yang , W. Liu . Instability and controllability of linearly coupled oscillators: Eigenvalue analysis. Phys. Rev. E, 1998, 58(4): 4440

[7] L. Huang , Q. Chen , Y. C. Lai , L. M. Pecora . Generic behavior of master-stability functions in coupled nonlinear dynamical systems. Phys. Rev. E, 2009, 80(3): 036204

[8] J. A. Acebrón , L. L. Bonilla , C. J. Pérez Vicente , F. Ritort , R. Spigler . The Kuramoto model: A simple paradigm for synchronization phenomena. Rev. Mod. Phys., 2005, 77(1): 137

[9] E. Ott , T. M. Antonsen . Low dimensional behavior of large systems of globally coupled oscillators. Chaos, 2008, 18(3): 037113

[10] K.Kaneko, Theory and Application of Coupled Map Lattice, Wiley, Chichester, 1993

[11] D. J. Watts , S. H. Strogatz . Collective dynamics of “small-world” networks. Nature, 1998, 393(6684): 440

[12] A. L. Barabási , R. Albert . Emergence of scaling in random networks. Science, 1999, 286(5439): 509

[13] M.E. J. Newman, Networks: An Introduction, Oxford University Press, 2010

[14] S. Boccaletti , V. Latora , Y. Moreno , M. Chavez , D. U. Hwang . Complex networks: Structure and dynamics. Phys. Rep., 2006, 424(4−5): 175

[15] A. Arenas , A. Diaz-Guilera , J. Kurths , Y. Moreno , C. S. Zhou . Synchronization in complex networks. Phys. Rep., 2008, 469(3): 93

[16] T. Wu , X. Zhang , Z. Liu . Understanding the mechanisms of brain functions from the angle of synchronization and complex network. Front. Phys., 2022, 17(3): 31504

[17] X. Wang , G. Chen . Synchronization in small-world dynamical networks. Int. J. Bifurcat. Chaos, 2002, 12(1): 187

[18] M. Barahona , L. M. Pecora . Synchronization in small-world systems. Phys. Rev. Lett., 2002, 89(5): 054101

[19] T. Nishikawa , A. E. Motter , Y. C. Lai , F. C. Hoppensteadt . Heterogeneity in oscillator networks: Are smaller worlds easier to synchronize. Phys. Rev. Lett., 2003, 91(1): 014101

[20] A. Arenas , A. Díaz-Guilera , C. J. Pérez-Vicente . Synchronization reveals topological scales in complex networks. Phys. Rev. Lett., 2006, 96(11): 114102

[21] D. Hansel , G. Mato , C. Meunier . Clustering and slow switching in globally coupled phase oscillators. Phys. Rev. E, 1993, 48(5): 3470

[22] M. Hasler , Yu. Maistrenko , O. Popovych . Simple example of partial synchronization of chaotic systems. Phys. Rev. E, 1998, 58(5): 6843

[23] Y. Zhang , G. Hu , H. A. Cerdeira , S. Chen , T. Braun , Y. Yao . Partial synchronization and spontaneous spatial ordering in coupled chaotic systems. Phys. Rev. E, 2001, 63(2): 026211

[24] A. Pikovsky , O. Popovych , Yu. Maistrenko . Resolving clusters in chaotic ensembles of globally coupled identical oscillators. Phys. Rev. Lett., 2001, 87(4): 044102

[25] I. A. Heisler , T. Braun , Y. Zhang , G. Hu , H. A. Cerdeira . Experimental investigation of partial synchronization in coupled chaotic oscillators. Chaos, 2003, 13(1): 185

[26] C. R. S. Williams , T. E. Murphy , R. Roy , F. Sorrentino , T. Dahms , E. Schöll . Experimental observations of group synchrony in a system of chaotic optoelectronic oscillators. Phys. Rev. Lett., 2013, 110(6): 064104

[27] J. Zhang , Y. Z. Yu , X. G. Wang . Synchronization of coupled metronomes on two layers. Front. Phys., 2017, 12(6): 120508

[28] M. M. Norton , N. Tompkins , B. Blanc , M. C. Cambria , J. Held , S. Fraden . Dynamics of reaction-diffusion oscillators in star and other networks with cyclic symmetries exhibiting multiple clusters. Phys. Rev. Lett., 2019, 123(14): 148301

[29] H. Fan , L. W. Kong , X. G. Wang , A. Hastings , Y. C. Lai . Synchronization within synchronization: Transients and intermittency in ecological networks. Natl. Sci. Rev., 2021, 8(10): nwaa269

[30] E. Rodriguez , N. George , J. P. Lachaux , J. Martinerie , B. Renault , F. J. Varela . Perception’s shadow: Long-distance synchronization of human brain activity. Nature, 1999, 397(6718): 430

[31] S. Kitsunai , W. Cho , C. Sano , S. Saetia , Z. Qin , Y. Koike , M. Frasca , N. Yoshimura , L. Minati . Generation of diverse insect-like gait patterns using networks of coupled Rössler systems. Chaos, 2020, 30(12): 123132

[32] J. F. Heagy , L. M. Pecora , T. L. Carroll . Short wavelength bifurcations and size instabilities in coupled oscillator systems. Phys. Rev. Lett., 1995, 774(21): 4185

[33] L. M. Pecora . Synchronization conditions and desynchronizing patterns in coupled limit-cycle and chaotic systems. Phys. Rev. E, 1998, 58(1): 347

[34] B. Ao , Z. G. Zheng . Partial synchronization on complex networks. Europhys. Lett., 2006, 74(2): 229

[35] F. Sorrentino , E. Ott . Network synchronization of groups. Phys. Rev. E, 2007, 76(5): 056114

[36] C. Fu , Z. Deng , L. Huang , X. G. Wang . Topological control of synchronous patterns in systems of networked chaotic oscillators. Phys. Rev. E, 2013, 87(3): 032909

[37] C. Fu , W. Lin , L. Huang , X. G. Wang . Synchronization transition in networked chaotic oscillators: The viewpoint from partial synchronization. Phys. Rev. E, 2014, 89(5): 052908

[38] L. M. Pecora , F. Sorrentino , A. M. Hagerstrom , T. E. Murphy , R. Roy . Cluster synchronization and isolated desynchronization in complex networks with symmetries. Nat. Commun., 2014, 5(1): 4079

[39] M. T. Schaub , N. O’Clery , Y. N. Billeh , J. C. Delvenne , R. Lambiotte , M. Barahona . Graph partitions and cluster synchronization in networks of oscillators. Chaos, 2016, 26(9): 094821

[40] F. Sorrentino , L. M. Pecora , A. M. Hagerstrom , T. E. Murphy , R. Roy . Complete characterization of the stability of cluster synchronization in complex dynamical networks. Sci. Adv., 2016, 2(4): e1501737

[41] J. D. Hart , Y. Zhang , R. Roy , A. E. Motter . Topological control of synchronization pattern: Trading symmetry for stability. Phys. Rev. Lett., 2019, 122(5): 058301

[42] D. M. Abrams , L. M. Pecora , A. E. Motter . Introduction to focus issue: Patterns of network synchronization. Chaos, 2016, 26(9): 094601

[43] M. Golubitsky , I. Stewart . Recent advances in symmetric and network dynamics. Chaos, 2015, 25(9): 097612

[44] W. Lin , H. Fan , Y. Wang , H. Ying , X. G. Wang . Controlling synchronous patterns in complex networks. Phys. Rev. E, 2016, 93(4): 042209

[45] W. Lin , H. Li , H. Ying , X. G. Wang . Inducing isolated-desynchronization states in complex network of coupled chaotic oscillators. Phys. Rev. E, 2016, 94(6): 062303

[46] T. Nishikawa , A. E. Motter . Network-complement transitions, symmetries, and cluster synchronization. Chaos, 2016, 26(9): 094818

[47] Y. Cho , T. Nishikawa , A. E. Motter . Stable chimeras and independently synchronizable clusters. Phys. Rev. Lett., 2017, 119(8): 084101

[48] B. Cao , Y. F. Wang , L. Wang , Y. Z. Yu , X. G. Wang . Cluster synchronization in complex network of coupled chaotic circuits: An experimental study. Front. Phys., 2018, 13(5): 130505

[49] Y. F. Wang , L. Wang , H. Fan , X. G. Wang . Cluster synchronization in networked nonidentical chaotic oscillators. Chaos, 2019, 29(9): 093118

[50] L. Wang , Y. Guo , Y. Wang , H. Fan , X. G. Wang . Pinning control of cluster synchronization in regular networks. Phys. Rev. Res., 2020, 2(2): 023084

[51] Y. Long , Z. Zhai , M. Tang , Y. Liu , Y. C. Lai . Structural position vectors and symmetries in complex networks. Chaos, 2022, 32(9): 093132

[52] D. M. Cardoso , C. Delorme , P. Rama . Laplacian eigenvectors and eigenvalues and almost equitable partitions. Eur. J. Combin., 2007, 28(3): 665

[53] M. A. D. Aguiar , A. P. S. Dias , M. Golubitsky , M. C. A. Leite . Bifurcations from regular quotient networks: A first insight. Physica D, 2009, 238(2): 137

[54] N. O’Clery , Y. Yuan , G. B. Stan , M. Barahona . Observability and coarse graining of consensus dynamics through the external equitable partition. Phys. Rev. E, 2013, 88(4): 042805

[55] D. Irving , F. Sorrentino . Synchronization of dynamical hypernetworks: Dimensionality reduction through simultaneous block-diagonalization of matrices. Phys. Rev. E, 2012, 86(5): 056102

[56] Y. Zhang , A. E. Motter . Symmetry-independent stability analysis of synchronization patterns. SIAM Rev., 2020, 86: 056102

[57] Y. Zhang , A. E. Motter . Unified treatment of synchronization patterns in generalized networks with higher-order, multilayer, and temporal interactions. Commun. Phys., 2021, 4(1): 195

[58] S. Panahi , N. Amaya , I. Klickstein , G. Novello , F. Sorrentino . Failure of the simultaneous block diagonalization technique applied to complete and cluster synchronization of random networks. Phys. Rev. E, 2022, 105(1): 014313

[59] M.GolubitskyI.StewartD.G. Schaeffer, Singularities and Groups in Bifurcation Theory, Springer-Verlag, 1985

[60] C. Zhou , L. Zemanová , G. Zamora , C. C. Hilgetag , J. Kurths . Hierarchical organization unveiled by functional connectivity in complex brain networks. Phys. Rev. Lett., 2006, 97(23): 238103

[61] C. Zhou , L. Zemanová , G. Zamora-López , C. C. Hilgetag , J. Kurths . Structure-function relationship in complex brain networks expressed by hierarchical synchronization. New J. Phys., 2007, 9(6): 178

[62] R. Wang , P. Lin , M. Liu , Y. Wu , T. Zhou , C. Zhou . Hierarchical connectome modes and critical state jointly maximize human brain functional diversity. Phys. Rev. Lett., 2019, 123(3): 038301

[63] S. Huo , C. Tian , M. Zheng , S. Guan , C. Zhou , Z. Liu . Spatial multi-scaled chimera states of cerebral cortex network and its inherent structure-dynamics relationship in human brain. Natl. Sci. Rev., 2020, 8(1): nwaa125

[64] M. E. J. Newman . Modularity and community structure in networks. Proc. Natl. Acad. Sci. USA, 2006, 103(23): 8577

[65] L. Huang , K. Park , Y. C. Lai , L. Yang , K. Yang . Abnormal synchronization in complex clustered networks. Phys. Rev. Lett., 2006, 97(16): 164101

[66] X. G. Wang , L. Huang , Y. C. Lai , C. H. Lai . Optimization of synchronization in gradient clustered networks. Phys. Rev. E, 2007, 76(5): 056113

[67] E. N. Lorenz . Deterministic nonperiodic flow. J. Atmos. Sci., 1963, 20(2): 130

[68] J. L. Hindmarsh , R. M. Rose . A model of neuronal bursting using three coupled first order differential equations. Proc. R. Soc. Lond. B, 1984, 221(1222): 87

[69] A. E. Motter , C. S. Zhou , J. Kurths . Enhancing complex-network synchronization. Europhys. Lett., 2005, 69(3): 334

[70] X. G. Wang , Y. C. Lai , C. H. Lai . Enhancing synchronization based on complex gradient networks. Phys. Rev. E, 2007, 75(5): 056205

[71] J. W. Scannell , G. A. P. C. Burns , C. C. Hilgetag , M. A. O’Neil , M. P. Young . The connectional organization of the cortico-thalamic system of the cat. Cereb. Cortex, 1999, 9(3): 277

[72] P. Hagmann , L. Cammoun , X. Gigandet , R. Meuli , C. J. Honey , V. J. Wedeen , O. Sporns . Mapping the structural core of human cerebral cortex. PLoS Biol., 2008, 6(7): e157

[73] C. J. Honey , O. Sporns , L. Cammoun , X. Gigandet , J. P. Thiran , R. Meuli , P. Hagmann . Predicting human resting-state functional connectivity from structural connectivity. Proc. Natl. Acad. Sci. USA, 2009, 106(6): 2035

[74] C. Fu , H. Zhang , M. Zhan , X. Wang . Synchronous patterns in complex systems. Phys. Rev. E, 2012, 85(6): 066208

[75] W. Poel , A. Zakharova , E. Schöll . Partial synchronization and partial amplitude death in mesoscale network motifs. Phys. Rev. E, 2015, 91(2): 022915

[76] P. Khanra , S. Ghosh , K. Alfaro-Bittner , P. Kundu , S. Boccaletti , C. Hens , P. Pal . Identifying symmetries and predicting cluster synchronization in complex networks. Chaos Solitons Fractals, 2022, 155: 111703

[77] F. B. Denton , S. J. Parke , T. Tao , X. Zhang . Eigenvectors from eigenvalues: A survey of a basic identity in linear algebra. Bull. Am. Math. Soc., 2022, 59(1): 31

[78] Y. Wang , D. Zhang , L. Wang , Q. Li , H. Cao , X. G. Wang . Cluster synchronization induced by manifold deformation. Chaos, 2022, 32(9): 093139

[79] J. Ma . Biophysical neurons, energy, and synapse controllability: A review. J. Zhejiang Univ. – Sci. A, 2023, 24: 109

RIGHTS & PERMISSIONS

Higher Education Press