Cluster and group synchronization in delay-coupled networks
We investigate the stability of synchronized states in delay-coupled networks where synchronization takes place in groups of different local dynamics or in cluster states in networks with identical local dynamics. Using a master stability approach, we find that the master stability function shows a discrete rotational symmetry depending on the number of groups. The coupling matrices that permit solutions on group or cluster synchronization manifolds show a very similar symmetry in their eigenvalue spectrum, which helps to simplify the evaluation of the master stability function. Our theory allows for the characterization of stability of different patterns of synchronized dynamics in networks with multiple delay times, multiple coupling functions, but also with multiple kinds of local dynamics in the networks’ nodes. We illustrate our results by calculating stability in the example of delay-coupled semiconductor lasers and in a model for neuronal spiking dynamics.
The scientific field of synchronization in coupled systems has evolved rapidly in the last decades . Complete or isochronous synchronization of coupled chaotic units  as well as of time-periodic systems has been extensively studied . In general, more complicated synchronization patterns may be observed including cluster, group, and sublattice synchronization . Cluster synchronization, where certain clusters inside the network show isochronous synchronization, will be investigated in this paper. Additionally, we describe the stability of group synchronization, i.e., a generalization of cluster synchronization where the local dynamics of the nodes in each group differ.
The characterization of stability of isochronous synchronization has been widely studied, and the ground-breaking work by Pecora and Carroll  which allows for a separation of network topology and local dynamics of the nodes was recently also applied to networks with delays in the links . Such delay times can greatly change the synchronization properties and appear in many natural coupled systems. For example, in optical applications delay times arise from the finite speed of light and in neuronal networks delays play a role due to finite distances between interacting neurons, but also due to processing lags in the neurons.
For group and cluster synchronization, attempts have been made to treat stability within a master stability approach. Sorrentino and Ott  considered two groups of nodes governed by different local dynamics. In the present paper, we show how this can be generalized to a higher number of groups and what restrictions for the topology of the network arise. Moreover, our framework allows us to have multiple delay times in the network. Making use of a separation of the topologies into multiple coupling matrices we can lift the restriction that no coupling may exist inside groups or clusters, i.e., a restriction to multipartite topologies. This makes our theory accessible for a wide range of topologies.
After introducing the notion of cluster and group dynamics in Section 2, we derive the master stability function and show the restrictions that arise upon the topology in Section 3. In Section 4 we investigate the symmetries that group and cluster synchronization impose on the master stability function. In Section 5, we demonstrate this symmetry for networks of delay-coupled lasers. Multiple coupling matrices are introduced in Section 6, where we use a hierarchical network structure as an example. The effect of different delay times is shown for the example of neuronal networks in Section 7. Finally, we conclude with Section 8.
2Cluster and group dynamics
In a network consisting of identical nodes, we refer to cluster synchronization as a state where clusters of nodes exist that show isochronous synchronization internally, but synchronization between these cluster does not occur, or is of non-isochronous type, i.e., there may be a phase lag between clusters .
Group synchronization describes a similar state of synchrony, but the node dynamics – determined by the functional form of the local dynamics – differs from cluster to cluster. We refer to these clusters as groups. As cluster synchronization is a special case of group synchronization, we use the more general notion of groups in the following.
Assume the number of groups to be where numbers the individual groups. The dynamical variables of the nodes in each group are then given by with , where denotes the number of nodes in the -th group. The dimension of the is given by the particular node model, e.g., the complex Hopf normal-form (Stuart-Landau) oscillator , the two-dimensional FitzHugh-Nagumo model , or the three-dimensional Lang-Kobayashi equations .
In general the dimension of the nodes may be different for each group . Consequently, also the local dynamics can be different for each group, but must be identical for all nodes in a given group . For example, consider a network of neurons, where one group contains inhibitory neurons and another group contains excitatory ones. The local dynamics will be different for each group, and depending on the model used to describe both types of neurons also the dimension of the node dynamics may be different.
Let be the coupling strength for the coupling from the -th to the -th group. In the same sense, let be an coupling matrix, such that its entries represent the coupling of node (which is in the -th group) to node (which is in the -th group). By this construction we obtain a multipartite topology in which one cluster has incoming links from only one neighbor while having outgoing links to another one. The stability analysis performed in this Section works for these topologies; but we will lift this restriction by allowing multiple coupling matrices in Section 6. Without loss of generality we assume the row sums of the coupling matrices to be unity, which corresponds to the condition of unity or constant row sum needed in the special case of complete isochronous synchronization . If a coupling matrix has arbitrary non-zero but constant row sum, unity row sum can easily be obtained by rescaling the corresponding coupling strength .
As coupling schemes we introduce matrices, given that and are the dimensions of and , i.e., the dimensions of the local dynamics in the -th and -th group, respectively. Note that, as a generalization, nonlinear coupling functions may also be used instead of matrices .
Finally, we allow the coupling delays to be different for any pair of groups being connected. A schematic diagram of the variables and matrices is shown in Figure 1(a). At this point we consider only multipartite topologies, i.e., only the dashed arrows in the Figure.
The dynamics of any single node in the network can then be described by the differential equation
for , . This type of coupling is applicable for optical systems  and electronic circuits. In other cases, for instance neural dynamics, a diffusive-like coupling term of the form is used. Both forms are equivalent since the local dynamics can be transformed by . In the following, we will use the form of Eq. .
The group synchronization manifold is then given by
which follows by inserting into Eq. (, ). For the example of two groups, Figure 1(b) illustrates the synchronization manifold, where Eq. corresponds to the dashed arrows only.
Note that each group may exhibit different synchronous dynamics. Even if the functions , the coupling matrices , and the delay times are identical for each group, different initial conditions can lead to different dynamics.
3Stability of group synchronization
In order to investigate the stability of the synchronous state, we linearize Eq. (Equation 1) around the group synchronization manifold :
Now assume that for each group each of the solutions of Eq. can be written in the form
with time-independent scalars . We show that the vectors formed from the possible combinations of the span a space of dimension , thus the form yields all solutions of Eq. as linear combinations. Using the form , Eq. becomes
Eq. can be rewritten as
is independent of . This is the case if a set of linearly independent vectors , , can be found, which we show in the following. Using these vectors , the conditions can be written as
One particular solution of Eq. (and equivalently of Eq. ) is obtained when setting :
Introducing rescaling factors and a fixed , this can be rewritten as
Substituting , Eq. becomes
Setting , it follows that Eq. yields all possible solutions of Eq. assuming that are free parameters and .
The scaling factors change only the magnitude of the variational vectors, thus their particular choice is not important for the stability of synchronization. Therefore setting in Eq. yields
which, in conclusion, qualifies as a master stability equation for this network topology. Here, is chosen from the set of eigenvalues of the block matrix
because Eq. is equivalent to the eigenvalue problem .
The largest Lyapunov exponent calculated from Eq. (Equation 12) as a function of the parameter is called the master stability function (MSF). It determines the stability of group synchronization if evaluated at the eigenvalues of .
4Symmetry of the master stability function
Note that the master stability equation () is of dimension and thus independent of the sizes of the individual groups and the particular coupling topologies . Because of the structure of , there always exist eigenvalues corresponding to dynamics inside the group synchronization manifold. We will refer to these as longitudinal eigenvalues.
Besides these longitudinal eigenvalues, the spectrum of shows a more general symmetry: For a given eigenvalue of , is also an eigenvalue of for any . See Appendix A for a detailed survey on the spectrum of the coupling matrix .
Looking closely at the master stability equation , we find another symmetry. The equation is invariant with respect to the transformation :
With the basis transformation , which leaves the Lyapunov spectrum unchanged, the original equation is regained Consequently, the master stability equation is invariant with respect to rotations .
Combining both results – the invariance of the MSF and the spectrum of against rotations of – we can conclude that it is sufficient to evaluate the MSF in an angular sector given by .
In the next Section, we demonstrate this symmetry and calculate the MSF for the example of delay-coupled laser networks.
5Example: laser networks
For semiconductor lasers subjected to optical feedback, the Lang-Kobayashi (LK) model  is a paradigmatic model. This model is based on simple rate equations and includes as variables the carrier inversion and the complex electric field , which is reduced to its slowly varying envelope. The LK model in its dimensionless form includes the local dynamics
where denotes the excess carrier density and the complex electric field . denotes the ratio of carrier and photon lifetimes, is the normalized pump current in excess of the laser threshold, and is the linewidth enhancement factor. The dynamics of a solitary laser – without any feedback or coupling – is described by . Coupling groups of lasers in a network of the form Eq. , we consider identical local dynamics and focus on all-optical coupling, thus
The cluster synchronization manifold and thus the master stability equation are -dimensional.
Figure 2 shows the MSF for one, two, three, and four clusters in panel (a), (b), (c), and (d), respectively. The black asterisks mark the position of the longitudinal eigenvalues , , of . Clearly visible in panels (a)-(d) is the symmetry with respect to discrete rotations of as discussed in the last Section. In particular, the Lyapunov exponent is identical at all longitudinal eigenvalues , . Note that this result is independent of a particular topology. Choosing any topology that has the structure , its eigenvalue will always show the discrete rotational symmetry discussed in Section 4.
For large delay, as shown in the right part of the Figure, the MSF has a circular shape for one cluster (panel (e)). This was recently shown to be a universal feature of networks where the coupling delay is large compared to the time scale of the local dynamics . Due to the discrete rotational symmetry discussed above, the circular shape cannot change when increasing the number of clusters, hence the shape of the MSF is independent of the number of clusters for large coupling delay (see Figure 2(f-h)), but we observe that the size of the disc of stability is shrinking with increasing number of clusters. This shrinking can be explained as follows: The dimension of the synchronization manifold Eq. is proportional to the number of clusters. Since the blocks of the matrix are arranged in a unidirectional ring, the dynamics inside the synchronization manifold lives inside such a unidirectional ring. Hence, the time that a signal takes to travel through this ring scales linearly with the number of groups . This signal traveling-time can be seen as an effective time-delay governing the degree of chaos, i.e., the longitudinal Lyapunov exponent. As was shown in Refs. , a larger longitudinal Lyapunov exponent yields a smaller radius of the stable region.
The above example used identical local dynamics in all of the groups, which corresponds to the case of cluster synchronization. In order to illustrate our theory for group synchronization, we now consider two groups of lasers, where the pump current is in the first group and in the second group. Figure 3 shows the resulting master stability function for delay times and in panels (a) and (b), respectively. Compared to Figs. Figure 2(b,f), only the pump current in one of the groups is increased. In the case of a small delay time (Fig. Figure 3(a)) this does not change the master stability function, because both groups still lock to the same dynamics. In the case of a large delay time (Fig. Figure 3(b)), the stable region shrinks compared to Figure 2(f) due to the higher pump current in one of the groups leading to a larger longitudinal Lyapunov exponent . Note that the discrete symmetry of the master stability function is also present for group synchronization.
6Beyond multipartite topologies
So far we have developed a master stability formalism to determine the stability of group and cluster synchronization. In order to utilize the master stability framework, one major restriction has been made: Each group must receive input from one and only one other group, i.e., the network topology has to be multipartite. More complex topologies beyond multipartite structures like, for instance, lattices  could not be dealt with. In the following we will derive the master stability equation for group synchronization with multiple coupling matrices. Thereby some of the former stringent restrictions can be dropped.
The network dynamics for group synchronization with one coupling matrix has been written in the form of Eq. , and the synchronization manifold and the master stability equation were given by Eqs. and , respectively. Generalizing this to two coupling matrices yields for the network dynamics of groups:
where the matrix describes the coupling from the -th to the -th group as before and describes the coupling from the -th to the -th group. That is, the -th group now receives input from two groups, and . The row sums of all and must be unity. Any constant non-zero row sum can be rescaled by means of the coupling strengths.
For the sake of simplicity and readability, we use identical coupling schemes and identical time delays for both coupling terms. In general, our framework works for different time delays and coupling schemes. The sum of and must yield the overall coupling strength used before in order to arrive at the same dynamical regime: . Figure 1(a) shows schematically the coupling parameters and matrices that are present in Eq. . In the case of two groups shown here, and represent the coupling within the groups, depicted by solid arrows.
From the above, the synchronization manifold is obtained as
for . See Figure 1(b) for a schematic diagram of the synchronization manifold for the example of two groups. The coupling inside a group translates into a self-feedback loop, depicted by solid arrows. Let be the matrix containing the blocks at positions and the matrix containing the blocks at positions . If and commute, i.e., , it is possible to obtain a master stability equation
for , where and are chosen from the eigenvalue spectrum of the matrices matrices and , respectively. These eigenvalues have to be evaluated in pairs corresponding to one eigenvector. Since and commute they always have a set of identical eigenvectors.
6.1Example: two groups
Let us first consider the simplest example, namely only two groups. Above, we have shown results for synchronization in two groups using a single coupling matrix of the form
where the matrices and describe the coupling from the second to the first group and vice versa, respectively, see also Ref. . We will now elaborate what happens when we introduce a second coupling matrix
i.e., and . Figure 4 shows the master stability function for the structure given by these matrices and and for laser parameters in the regime of low-frequency fluctuations with , , . The coupling strengths are chosen as and with . This resembles strong coupling in the clusters, but weak coupling between clusters. We use only one time delay for simplicity. For the same reason, we investigate cluster synchronization, i.e., the local dynamics and the coupling scheme are identical for both groups in this example. We consider matrices with real eigenspectrum only and set . The eigenvalue pairs depicted by the black (blue) dots correspond to a particular network topology that will be discussed below.
Using the forms and , the commutation relation is equivalent to
These conditions are fulfilled for certain classes of matrices only. We will give an example of hierarchical coupling that yields matrices which fulfill these conditions.
6.2Towards hierarchical networks
A hierarchical network usually consists of topological clusters that are densely coupled inside, while links to other such topological clusters are sparse. The hierarchy is then built by larger topological clusters that contain the smaller ones . This procedure can be continued over many levels of hierarchy. It is important to distinguish these topological clusters from the dynamical cluster states that are investigated in this paper.
The simplest hierarchical structure consists of just two topological clusters. Figure 5 illustrates this in a schematic sketch of a graph of nodes with two topological clusters, . Solid (blue) arrows correspond to links inside one cluster while dashed (red) arrows denote links between both clusters.
In this Section we will show that each cluster can exhibit isochronous synchronization under certain conditions. In this sense, the notions of topological cluster and of dynamical cluster coincide at this point.
The graph in Figure 5 is modeled by the coupling matrices
where is the identity matrix and is an undirected Erdős-Rényi random graph with a certain link probability . The undirectedness is necessary to obtain a real-valued eigenvalue spectrum. Then it is sufficient to calculate the master stability function in the plane as done in Figure 4.
In order to comply with the link density of a hierarchical network, we choose the coupling strength for the two matrices and to be different as used for the calculation of the master stability function in Figure 4. The coupling strengths and are chosen as , where is the overall coupling strength corresponding to the regime of low-frequency fluctuations of the laser dynamics. The coupling strengths corresponding to are chosen as . For a high link probability in the random matrix , the matrix contains comparatively more links than , which has only one link per row. Given that both matrices are renormalized to unity row sum, it is therefore a reasonable choice that and are significantly larger.
The black (blue) dots in Figure 4 show the eigenvalue pairs of the hierarchical example given by Eqs. and using a link probability of in the random matrix for . It can be seen that this network shows stable synchronization in this 2-cluster state. That is, each topological cluster exhibits synchronization internally. All eigenvalue pairs transversal to the synchronization manifold are inside the stable region, while the longitudinal eigenvalue pairs and do not affect the stability of synchronization. For any choice of the eigenvalues will always be lined up on the dotted vertical lines, which are determined by the matrix being constructed from identity-matrix blocks.
The link probability is just above the threshold of stable synchronization. Using lower values, some eigenvalues will cross the boundary of the stable region of the master stability function, leading to desynchronization.
Since the eigenvalues are always aligned along the lines in this example, the stability for other choices of the coupling strength can easily be obtained by evaluating the master stability function at a fixed value of as a function of and or . The other value yields identical results and therefore need not be considered, which is a result of the symmetry discussed in Section 4 and observable also in Figure 4.
Figure 6 shows the master stability function in the plane, where . The other coupling strength is set using the relation , where the overall coupling strength is chosen as corresponding to the laser regime of low-frequency fluctuations. This relation ensures that the overall coupling strength leads to an operation in this regime.
The dashed black (blue) lines enclose the region where no 2-cluster state can exist. This can be seen by considering the two-node network motif described by
, where the coupling matrix
describes the behavior on the 2-cluster synchronization manifold. For coupling strengths between the black (blue) lines, this motif shows stable synchronization for the chosen laser parameters and thus the dynamics in both clusters will be identical. The stability of the 2-cluster state is therefore only meaningful below the lower and above the upper dashed black (blue) line.
The boundaries of stability for the 2-cluster state are nearly independent of in the lower range of , which corresponds to a high coupling strength inside the clusters, but a low coupling strength between clusters. The upper range of , which corresponds to low coupling strength inside the clusters, but high coupling strength between them, also allows for the existence of the 2-cluster state, but this state cannot be stable for any topology. In conclusion, the coupling strength must be comparatively large inside the clusters to allow for a stable 2-cluster state.
7Example: neural networks
Synchronization in the brain can be related to cognitive capacities  as well as to pathological conditions, e.g., epilepsy . Therefore, there has been tremendous interest in the study of synchronization in neural networks . The master stability approach has been applied to the study of synchronization patterns independently of a specific network topology . The brain is organized in different brain areas leading to different delay times between neurons of different areas and neurons within the same area. Furthermore, different types of neurons exist, corresponding to different local dynamics. Therefore we propose that the master stability function for group synchronization introduced here will be especially useful for investigating complex neural synchronization phenomena.
Here we apply our method to a neural network where the nodes are modeled as FitzHugh-Nagumo (FHN) systems. As in the last Section we consider a network of two groups coupled via two coupling matrices (intergroup coupling) and (intragroup coupling). We use a diffusive-like coupling. As discussed above this can be transformed to the coupling used in the previous sections by transforming the local dynamics of the -th node in the -th cluster as follows:
with and . Here and denote the activator and inhibitor variables, respectively. The parameter determines the threshold of excitability. A single FHN oscillator is excitable for and exhibits self-sustained periodic firing beyond the Hopf bifurcation at . We will focus on the excitable regime with . The time-scale parameter is chosen as . The synchronized dynamics and the master stability equation are then given by Eq. and Eq. , respectively. We assume the coupling scheme .
The cluster synchronized dynamics is equivalent to a system of two coupled nodes with self-feedback. In Ref.  it was shown that depending on the delay times, the coupling strength, and the strength of the self-feedback different dynamical scenarios, i.e., in-phase synchronization, anti-phase synchronization, or bursting can arise. Figure 7 shows the master stability function in panels (a)-(c) for in-phase synchronization, anti-phase synchronization and for synchronization in two bursting groups, respectively. The right hand panels of Figure 7 depict the corresponding time series: In panel (d), (f), and (h) for the activator variables and in panel (e), (g), and (i) for the inhibitor for in-phase, anti-phase, and bursting dynamics, respectively. Because the different dynamical scenarios yield distinctively different stable regions, topologies might arise which show stable synchronization for one of the patterns but not for the others. However, for all scenarios the stable region contains the unity square, i.e., . With Gershgorin’s circle theorem  it can easily be shown that the eigenvalues of symmetrical matrices with positive entries and unity row sum are always contained in the interval . Thus, if and have only positive entries, i.e., if the coupling is excitatory, synchronization is stable for the dynamics and parameters shown here. As a consequence, only the introduction of inhibitory links can lead to desynchronization. In Ref.  it has been shown that for and this is the case for all and for which the synchronized dynamics is periodic. A detailed study of these phenomena for the eight-dimensional parameter space of () is beyond the scope of this paper. Note that the symmetry discussed in Section 4 does not show up in Figure 7(a), because both clusters synchronize to and the invariance of Eq. does not hold in this case of in-phase synchronized spiking.
As an example of a network with inhibitory links which will exhibit stable synchronization only in one of the patterns discussed above, but not in the other ones, we choose and as
where with is an all-to-all coupling matrix with self-coupling, and is an undirected random matrix with both excitatory (positive entries) and inhibitory links (negative entries). The matrix describes a fixed node degree with 12 excitatory and 9 inhibitory links for each node. The number of nodes is chosen as . The black dots in Figure 7 denote the corresponding eigenvalue pairs. In panels (a) and (b) some eigenvalues are located outside the stable region, while in panel (c) they are all inside, which means that the zero-lag and anti-phase synchronized solutions will be unstable in such a network, while synchronization in the bursting state will be stable.
Based on a master stability approach, we have studied patterns of cluster and group synchronization in delay-coupled networks and determined their stability. We have shown that the master stability function applied to cluster and group synchronization exhibits a discrete -fold rotational symmetry for dynamical clusters. This reduces the numerical effort, such that for a larger number of clusters the master stability function must be evaluated only on a smaller angular sector in the complex plane. Within our approach we can treat a wide range of multipartite network topologies. Using multiple commuting coupling matrices, we have generalized our stability analysis beyond multipartite topologies, for instance towards hierarchical network structures. As concrete examples we have focused on delay-coupled lasers and neural networks. The interplay of complex topologies, multiple delay times, and possibly different local dynamics and different coupling functions extends the scope of the master stability framework and is a step towards understanding complex patterns of synchronization in real-world networks.
ASpectrum of the coupling matrix
We investigate the eigenvalue spectrum of the coupling matrix
has at least zero eigenvalues, where
arises solely due to the block structure of : if all have maximum rank , there are exactly those zeros. In general, the exact number of zeros is given by , which may be larger than due to the particular structure of the .
Consider the matrix , which is of block diagonal structure
Note that each block on the diagonal is a product of all , only the order differs.
Assume that the groups are arranged such that (), which can always be achieved by an index permutation, and that each has maximum rank
The non-zero eigenvalues of are then given by the -th roots of the non-zero eigenvalues of , and the whole spectrum of reads
Note, in particular, that the eigenvalue of corresponds to the longitudinal eigenvalues of , which are related to directions longitudinal to the group synchronization manifold. Their existence can already be seen solely by looking at itself, because its eigenvectors
where each corresponds to the longitudinal eigenvalue , do not depend on the inner structure of the blocks .
Given that the MSF is invariant with respect to rotations and that each of the multiple roots of are also rotations by multiples of with respect to the roots ,, , we can restrict ourselves to evaluating the master stability function at the location of the eigenvalues
which lie all inside the angular sector . Note that the zero eigenvalue is only added here if , i.e., if at least one block differs from the others in size.
- The latter assumption simplifies the argument regarding the zero eigenvalues, but the final result is valid for arbitrary ranks of the block matrices
- (, , )bibitemNoStop
- (, , )bibitemNoStop
- (, , )bibitemNoStop
- , bibitemNoStop