Topologically protected synchronization in networks

Many physical, biological, chemical, and technological systems can effectively be seen as networks of interacting oscillators, each with its own natural frequency. Two oscillators are considered neighbors if they are coupled by some nonlinear function of their phase difference. Spontaneous synchronization represents one of the most intriguing and ubiquitous aspects of these systems StrogatzC ; Pikovsky ; ArenasReview and has countless applications such as, pacemaker cells in the heart Heart , pacemaker cells in the brain Brain , flashing fireflies Fireflies , arrays of lasers Lasers , and superconducting Josephson junctions Josephson . Put simply, if the coupling constant $J$ is larger than the spread $\Delta\omega$ of the natural frequency distribution in the system, a finite portion of its oscillators tends to synchronize, i.e., they tend to rotate according to a common mean frequency and their phases fall into step with one another (locked in phase). Kuramoto-like models, where the coupling between two oscillators is modulated by the sine of their phase difference, provide the most popular approach for addressing synchronization in complex systems. In its mean-field version, where each oscillator is coupled to all other oscillators, the model allows to be exactly solved Strogatz2000 while, in the general case, despite no exact solution is known, rigorous bounds Strogatz1988 , large scale simulations and theoretical studies ArenasReview ; Munoz ; Kahng have succeeded to provide a general description of the different synchronization scenarios that can take place as a function of the topology of the underlying graph Goltsev , ranging from regular to random, small-world, and even scale-free. For the present manuscript, we emphasize the role played by modularity. Many complex networks are modular, i.e., can be seen as subgraphs (also known as groups or communities) with different internal and external connectivities Barabasi_Hierarchical ; Santo . It has been understood that, in modular networks, the densely connected communities synchronize first, and subsequently, the larger and less densely connected ones also tend to synchronize until full synchronization is achieved ArenasPRL ; ArenasReview .

Most of the above literature analyzes synchronization in the thermodynamic limit, viewed as a phase transition characterized by an emergent phase diagram containing essentially an incoherent region (where most of the oscillators follow their own natural frequencies) and a coherent one (where most of the oscillators get locked within each other). Here, we consider a rather different issue related to a local and protected synchronization: “local” because it concerns a finite group of oscillators embedded in some arbitrary graph containing the group; “protected” because the remainder of the graph cannot affect the synchronization status of the group, even when the former is subject to noise. Formally, our work belongs to what is known as synchronization of chaotic systems via bidirectional coupling Boccaletti , a branch less developed than the unidirectional case (or “slave-master” Note ) Boccaletti2 ; Topological_Non_Hermitian . At any rate, to the best of our knowledge, our analysis and findings — based on topological equivalence combined with modularity — have never been reported. In recent years, topology has been shown to play a major role in many areas of physics, ranging from topological insulators Topological_Insulator to quantum computation Topological_Quantum_Computing ; in this sense, synchronization cannot be an exception. We stress however that, despite a reminiscent terminology, our definition of topologically equivalent (TE) nodes is quite unrelated to these works.

Given a graph with $N$ nodes and symmetric adjacency matrix $\bm{a}$, consider a Kuramoto-like model on the top it, where the phases of $N$ oscillators, $\theta_{1},\ldots,\theta_{N}$, with natural frequencies $\omega_{1},\ldots,\omega_{N}$, evolve according to

$\displaystyle\dot{\theta}_{i}=\omega_{i}+J\sum_{j}a_{i,j}\sin\left(\theta_{j}-\theta_{i}\right),$

(1)

$J>0$ being a coupling constant.

Group of $N^{\prime}=2$ TE nodes. Consider now two specific nodes, say node 1 and node 2 (note that $a_{1,2}=a_{2,1}$). From now on, unless otherwise stated, we assume that the two nodes are connected: $a_{1,2}=1$. For the phase difference variable ${\Theta}=\theta_{2}-\theta_{1}$, Eqs. (1) lead to

$\displaystyle\dot{{\Theta}}=\omega_{2}-\omega_{1}-2J\sin\left({\Theta}\right)+J\left[h_{2}(\theta_{2})-h_{1}(\theta_{1})\right],$

(2)

where we have introduced the two functions

$\displaystyle h_{i}(\theta)=\sum_{j\neq 1,2}a_{i,j}\sin\left(\theta_{j}-\theta\right),\quad i=1,2.$

(3)

Let us now suppose that node 1 and node 2 “see” the same remaining graph, i.e., i.e.,

$\displaystyle a_{1,j}=a_{2,j},\quad\forall j\neq 1,2.$

(4)

We shall say that the two nodes are TE. In this case

$\displaystyle h(\theta)=h_{1}(\theta)=h_{2}(\theta),$

(5)

and Eq. (2) becomes

$\displaystyle\dot{{\Theta}}=\omega_{2}-\omega_{1}-2J\sin\left({\Theta}\right)+J\left[h(\theta_{2})-h(\theta_{1})\right].$

(6)

We can rewrite Eq. (6) as

$\displaystyle\dot{{\Theta}}=\omega_{2}-\omega_{1}-2J\sin\left({\Theta}\right)+Jh^{\prime}({\theta}^{*}){\Theta}$

(7)

where ${\theta}^{*}\in(\min\{\theta_{1},\theta_{2}\},\max\{\theta_{1},\theta_{2}\})$ and we have applied the mean-value theorem to the function $h(\theta)$. Note that ${\theta}^{*}$ is itself an unknown function of $\theta_{1}$ and ${\Theta}$ (or, alternatively, $\theta_{2}$ and ${\Theta}$), yet, as we shall see in a moment, it is worth to consider Eq. (7). We will use the following bounds:

$\displaystyle h^{\prime}(\theta)\leq k^{(\mathrm{OUT})}\quad\forall\theta,$

(8)

where $k^{(\mathrm{OUT})}$ is the (common) outgoing degree of the two nodes

$\displaystyle k^{(\mathrm{OUT})}=\sum_{j\neq 1,2}a_{1,j}=\sum_{j\neq 1,2}a_{2,j},$

(9)

and

$\displaystyle-\theta\leq-\sin(\theta)<-\theta(1-\gamma),\quad 0\leq\theta<1,$

(10)

where

$\displaystyle\gamma=\frac{1}{3!}+\frac{1}{7!}+\frac{1}{11!}+\ldots=0.1668651044\ldots.$

(11)

Equal frequencies with $a_{1,2}=1$. Let us assume that the two nodes have also equal natural frequencies, $\omega_{1}=\omega_{2}$. From Eq. (7) we see that, in the manifold $\theta_{1}=\theta_{2}$ ($=\theta^{*}$), the two external fields cancel out, in other words, ${\Theta}=~{}0$ provides a fixed point of the equation. Let us introduce $\lambda=k^{(\mathrm{OUT})}-2(1-\gamma)$ and $\mu=k^{(\mathrm{OUT})}+2$. From Eq. (7), by using the bounds (8) and (10), we see that, if $1>{\Theta}\geq 0$,

$\displaystyle-J\mu{\Theta}\leq\dot{{\Theta}}\leq J\lambda{\Theta},$

(12)

while, if $-1<{\Theta}<0$, hold the opposite inequalities. Let us assume that, for the initial condition $\Theta_{0}=\Theta(t=0)$, we have $0<\Theta_{0}<1$. By continuity, there exists a sufficiently small time $t_{1}$ such that $1>{\Theta}(t)>0$ for $t\in[0,t_{1})$ whereby, from Eq. (12),

$\displaystyle{\Theta}_{0}e^{-J\mu t}\leq{\Theta}\leq{\Theta}_{0}e^{J\lambda t}.$

(13)

On the other hand, if $\lambda<0$, Eq. (13) and the fact that each $\theta_{i}(t)$ (an hence also $\theta_{2}-\theta_{1}$) as well as $\dot{\theta}_{i}(t)$ are continuous, imply that $t_{1}=\infty$, so that Eq. (13) holds for any $t$ SM . For simplicity, in the subsequent cases, on assuming certain constrains on the initial conditions, we shall limit ourselves to check that the upper bounding solution satisfies the same constrains for any $t$. Note however that, in general, simply bounding $|\Theta|$ is not enough for claiming good synchronization; the conservation of the sign of $\Theta$ is crucial.

In conclusion, within a basin of attraction for the initial conditions contained in the region $|{\Theta}_{0}|<1$, a sufficient condition for the fixed point ${\Theta}=0$ to be stable is $k^{(\mathrm{OUT})}<2(1-\gamma)$. Of course, since $1<2(1-\gamma)<2$ and $k^{(\mathrm{OUT})}$ is integer, it follows that $k^{(\mathrm{OUT})}$ can be either 0 or 1, but, in view of generalizations to subsystems with $N^{\prime}>2$ TE oscillators, it is useful to keep in mind the inequality. Note that this synchronization between node 1 and node 2 occurs regardless of the dynamics of all the other oscillators, which in particular do not need to be synchronized. Note also that we have imposed only the equality of the frequencies of the two nodes, but the other nodes can have arbitrary frequencies and the result does not change. We stress that this is a consequence of the topological equivalence (4); when (4) does not hold, $h_{1}$ and $h_{2}$ remain two different functions so that, in the manifold $\theta_{1}=\theta_{2}$, they do not cancel out in Eq. (2). We can even imaging to modify the inner part of the remainder of the graph by dynamically removing, adding, or rewiring some of its links, as well as by allowing for the presence of any site-dependent noise: as far as such links are not those arriving at nodes 1 and 2, and as far as such noise applies only to the other nodes, the sufficient condition $k^{(\mathrm{OUT})}\leq 2(1-\gamma)$ remains satisfied, i.e., the subsystem fixed point $\theta_{1}=\theta_{2}$ keeps being stable.

Equal frequencies with $a_{1,2}=0$. Before analyzing the most general case, it is worth also considering the sub-case in which the two nodes are disconnected,i.e., $a_{1,2}=0$ and have equal frequencies. In this case, the evolution equation for ${\Theta}$ is of no help because $\theta_{1}=\theta_{2}$ is no longer a stable attractor, in other words, in general, they do not get synchronized. However, from Eq. (1) we see that, if nodes 1 and 2 are TE, $\theta_{1}$ and $\theta_{2}$ obey the same identical equation: we simply have that one oscillator follows the other along the same identical trajectory and, in particular, if they started with the same initial condition, they will remain identical at any instant. In this specific case, ${\Theta}=0$ is no longer an attractor but rather a constant of the motion.

Different frequencies with $a_{1,2}=1$. Let us now allow for the two frequencies to be different. If the pair is isolated, i.e., $k^{(\mathrm{OUT})}=0$, we have $h^{\prime}(\theta^{*})=0$ and Eq. (7) returns the known critical condition for the synchronization of $N^{\prime}=2$ oscillators: $2J>|\omega_{2}-\omega_{1}|$, the stability condition of the fixed point ${\Theta}$ being $\cos({\Theta})>0$. In the general case, $h^{\prime}(\theta^{*})\neq 0$, we must observe that, although Eq. (7) admits (under suitable conditions) a formal fixed point, due to the fact that $\theta^{*}$ depends on $\theta_{1}$ (besides ${\Theta}$), such a fixed point is conditioned on the value of $\theta_{1}$, which in turn is a function of time, rendering this formal fixed point useless. Our approach here is different: we assume that both $|\omega|=|\omega_{2}-\omega_{1}|/J$ and the initial condition $|{\Theta}_{0}|=|\theta_{2}(t=0)-\theta_{1}(t=0)|$ are sufficiently small and look for a bounding solution by exploiting again the bounds (8) and (10). Let us suppose ${\Omega}=\omega_{2}-\omega_{1}>0$ and $0<{\Theta}_{0}<1$. By continuity, at small enough times, we also have $0<{\Theta}<1$ and from Eq. (7) we get

$\displaystyle\dot{{\Theta}}\leq{\Omega}+J\lambda{\Theta}.$

(14)

Equation (14) implies that, if $\lambda=k^{(\mathrm{OUT})}-2(1-\gamma)<0$, then ${\Theta}$, with initial condition $0<{\Theta}_{0}<1$, will remain bounded as follows

$\displaystyle{\Theta}\leq\frac{{\Omega}}{J|\lambda|}+\left[{\Theta}_{0}-\frac{{\Omega}}{J|\lambda|}\right]\exp\left[-J|\lambda|t\right].$

(15)

Note that, for $0<\Omega/J|\lambda|<\Theta_{0}$, all the above procedure turns out to be consistent with the required bound $0<{\Theta}<1$ for any $t$.

In conclusion, we have proven that, under the three conditions $k^{(\mathrm{OUT})}<2(1-\gamma)$ (equivalent to $k^{(\mathrm{OUT})}\leq 1$), $0<{\Theta}_{0}<1$, and $0<\Omega/J|\lambda|<\Theta_{0}$, we have $0\leq\Theta\leq\Omega/J|\lambda|$ when $t\to\infty$. In other words, the phases of the two oscillators get asymptotically close to each other and the larger is $J$ the closer they stay. Moreover, as in the case of equal frequencies, we see that this bound holds regardless of the status of all the other oscillators located on the remainder of the graph and, again, any change in it, cannot affect the synchronization of the pair.

It is important to note that the nature of the condition $k^{(\mathrm{OUT})}\leq 1$, unlike the others related to the initial value of ${\Theta}$ and to the ratio $|\Omega|/J$, is strictly topological. We shall call it the topological condition (TC).

Generalization to $N^{\prime}>2$ TE oscillators. Let us now consider a subsystem of $N^{\prime}=3$ TE nodes, say nodes 1, 2 and 3. Besides being TE with respect to the remainder of the graph, we want them to be TE among each other; in particular, each of them must have the same number of links pointing to the other two TE nodes. For $N^{\prime}=3$ there exists only one possibility, the one where $k^{(\mathrm{IN})}=2$, i.e., the three nodes form a triangle. Here we have introduced $k^{(\mathrm{IN})}$ as the internal (common) connectivity of the group, i.e., $k^{(\mathrm{IN})}$ is the number of links emanating from a node of the group and pointing to other nodes of the same group. Note that, formally, both the subsystem $N^{\prime}=3$ and the one already seen case $N^{\prime}=2$ (where it was $k^{(\mathrm{IN})}=1$), are fully connected (FC) graphs with $N^{\prime}$ nodes. However, as Fig. 1 shows, when $N^{\prime}>3$, there exist more configurations in which the $N^{\prime}$ TE nodes can be arranged. In fact, the number of possible arrangements tends to grow exponentially with $N^{\prime}$, but with a smaller rate for $N^{\prime}$ odd (see also SM ).

We introduce some notation and point out a few crucial points. Given $N^{\prime}$ TE nodes with indices $1,2,\ldots,N^{\prime}$, we indicate their corresponding phase differences by ${\Theta}_{i,j}=\theta_{i}-\theta_{j}$ and frequency differences by ${\Omega}_{i,j}=\omega_{i}-\omega_{j}$. Observe that ${\Theta}_{j,i}=-{\Theta}_{i,j}$ and that the variables $\{{\Theta}_{i,j}\}$ with $j>i$, are not all independent. For example, for $N^{\prime}=3$ we have the constrain ${\Theta}_{3,1}={\Theta}_{3,2}+{\Theta}_{3,1}$. In general, given $N^{\prime}$, we can always write a system of $N^{\prime}-1$ independent equations involving $N^{\prime}-1$ independent variables. We shall also assume $0<{\Theta}_{i,j}(0)<1/(N^{\prime}-1)$ so that, at small enough times, we also have $0<{\Theta}_{i,j}(t)<1/(N^{\prime}-1)$. Our general strategy is to use the bounds (8), (10); if the found bounding solution keeps satisfying the above constrains for any time $t$, the procedure is consistent. The resulting bounding system can be written vectorially:

$\displaystyle\dot{\bm{{\Theta}}}\leq\bm{{\Omega}}+J\bm{B}\cdot\bm{{\Theta}},$

(16)

where $\bm{{\Theta}}=({\Theta}_{2,1},{\Theta}_{3,2},\ldots,{\Theta}_{N^{\prime},N^{\prime}-1})^{T}$, $\bm{{\Omega}}=({\Omega}_{2,1},{\Omega}_{3,2},\ldots,{\Omega}_{N^{\prime},N^{\prime}-1})^{T}$, and $\bm{B}$ is a $(N^{\prime}-1)\times(N^{\prime}-1)$ matrix that depends on the parameters $N^{\prime}$, $k^{(\mathrm{OUT})}$, $k^{(\mathrm{IN})}$, and $\gamma$. Let $\lambda_{i}$ and $\bm{u}_{i}$ be the eigenvalues and normalized eigenvectors of $\bm{B}$, respectively. The coefficients $c_{i}(t)=\bm{{\Theta}}\cdot\bm{u}_{i}$ satisfy the decoupled system $dc_{i}(t)/dt\leq\bm{{\Omega}}\cdot\bm{u}_{i}+J\lambda_{i}c_{i}(t)$. As a consequence, each component of the vector $\bm{{\Theta}}$ remains bounded if all the eigenvalues of $\bm{B}$ are negative, leading to the TC (note that the eigenvalues do not depend on $J$ or $\bm{{\Omega}}$). For example, in the FC case $N^{\prime}=3$ we have

$\displaystyle\bm{B}=\begin{bmatrix}k^{(\mathrm{OUT})}-3(1-\gamma)&\gamma\\ \gamma&k^{(\mathrm{OUT})}-3(1-\gamma),\end{bmatrix}$

(17)

and the eigenvalues of $B$ are $\lambda_{1}=k^{(\mathrm{OUT})}-3+4\gamma$ and $\lambda_{2}=k^{(\mathrm{OUT})}-3+2\gamma$. We conclude that the group of $N^{\prime}=3$ TE FC nodes synchronize if $k^{(\mathrm{OUT})}-3+4\gamma<0$, which, on observing that $0<4\gamma<1$, amounts to $k^{(\mathrm{OUT})}\leq 2$, or $k^{(\mathrm{OUT})}\leq k^{(\mathrm{IN})}$, where we have made use of the fact that, here, $k^{(\mathrm{IN})}=2$. Note that also for the previously seen case $N^{\prime}=2$, where it was $k^{(\mathrm{IN})}=1$, the TC for the synchronization can be written as $k^{(\mathrm{OUT})}\leq k^{(\mathrm{IN})}$. In general, for a group of $N^{\prime}$ TE nodes FC, the diagonal elements of the matrix $\bm{B}$ are all equal to $k^{(\mathrm{OUT})}-N^{\prime}(1-\gamma)$, while the off-diagonal ones are in the form $p\gamma$ where $p\in\{1,\ldots,N^{\prime}-2\}$. We have verified that for $N^{\prime}=4$ the TC remains $k^{(\mathrm{OUT})}\leq k^{(\mathrm{IN})}$ but for $N^{\prime}=5$ the TC becomes $k^{(\mathrm{OUT})}\leq k^{(\mathrm{IN})}-1$. For details see SM where we also analyze cases where the group forms regular polygons. We stress that these are sufficient conditions; they might be not necessary.

We illustrate our analysis by solving numerically Numerical_Recipes model (1) in a system of moderate size but involving most of the ideas so far discussed. The graph associated to this system is depicted in Fig. 1. It contains five groups of TE nodes of various kind for which the TC is satisfied. We performed several simulations with different choices of the parameters $\omega_{1},\ldots,\omega_{N}$ and $\theta_{1}(0),\ldots,\theta_{N}(0)$ ranging from situations in which the sufficient conditions for synchronization discussed previously are strictly satisfied, to situations in which only the TC is satisfied. Remarkably, even in these latter situations, the synchronization scenario predicted analytically (via stricter conditions) holds. Figs. 2 and 3 show the behavior of $\theta_{i+1}(t)-\theta_{i}(t)$ in the two situations. Inset of Fig. 3 shows also the rotating numbers $\theta_{i}(t)/t$. For more details see SM .

In conclusion, simulations confirm the local synchronization scenario predicted analytically and extend its region of validity for a quite wider range of initial conditions and spread of natural frequencies: groups of TE nodes with $k^{(\mathrm{OUT})}$ sufficiently smaller than $k^{(\mathrm{IN})}$ (depending on the specific group), synchronize and remain protected from the remainder of the system, which could be even noisy. Moreover, the plots of the rotating numbers strengthen the idea that these groups of TE nodes play the role of independent pacemakers of the system. In fact, in these plots, the number of lines of convergence seems to scale as the number of TE groups. As the coupling $J$ increases, these lines approach each other and eventually converge to a unique line establishing global synchronization, consistently with ArenasPRL ; Note2 . However, for any finite value of $J$, each group keeps its own synchronization status, regardless of the others; this fact seems to hinder global synchronization while allowing for the local one. Of course, this hypothesis needs to be tested in larger systems. Many other questions arise, particularly the following: i) Do real-world networks contain groups of TE nodes? ii) Do real-world networks exist where the number of different groups, say $Q$, scales with the system size $N$? iii) Is it possible to probe efficiently very large systems to quantify $Q$ as well as to identify the nodes of these groups? We anticipate that these three questions meet a positive answer and will be the subject of a subsequent publication.