Global Robustness vs. Local Vulnerabilities in Complex Synchronous Networks

TL;DR

This paper analyzes how complex networks of coupled oscillators respond to local perturbations, identifying vulnerable nodes and proposing spectral topological indices to measure and enhance global robustness.

Contribution

It introduces resistance centralities and Kirchhoff indices derived from spectral analysis as novel tools to assess and improve network robustness against perturbations.

Findings

Resistance centralities identify most fragile oscillators.

Kirchhoff indices quantify global robustness.

Inertia has minimal impact on robustness in homogeneous systems.

Abstract

In complex network-coupled dynamical systems, two questions of central importance are how to identify the most vulnerable components and how to devise a network making the overall system more robust to external perturbations. To address these two questions, we investigate the response of complex networks of coupled oscillators to local perturbations. We quantify the magnitude of the resulting excursion away from the unperturbed synchronous state through quadratic performance measures in the angle or frequency deviations. We find that the most fragile oscillators in a given network are identified by centralities constructed from network resistance distances. Further defining the global robustness of the system from the average response over ensembles of homogeneously distributed perturbations, we find that it is given by a family of topological indices known as generalized Kirchhoff indices. Both resistance centralities and Kirchhoff indices are obtained from a spectral decomposition of the stability matrix of the unperturbed dynamics and can be expressed in terms of resistance distances. We investigate the properties of these topological indices in small-world and regular networks. In the case of oscillators with homogeneous inertia and damping coefficients, we find that inertia only has small effects on robustness of coupled oscillators. Numerical results illustrate the validity of the theory.

Equations

$${\Omega }_{ij}^{(1)}={\mathbb{L}}_{ii}^{†}+{\mathbb{L}}_{jj}^{†}-{\mathbb{L}}_{ij}^{†}-{\mathbb{L}}_{ji}^{†},$$

$${\Omega }_{ij}^{(1)}=\sum _{\alpha \ge 2}\frac{(u_{\alpha ,i}-u_{\alpha ,j}{)}^2}{{\lambda }_{\alpha }},$$

$$\begin{array}{r}{C}_1(k)={\left[n^{-1}\sum _j{\Omega }_{kj}^{(1)}\right]}^{-1}={\left[\sum _{\alpha \ge 2}\frac{u_{\alpha ,k}^2}{{\lambda }_{\alpha }}+n^{-2}K{f}_1\right]}^{-1}.\end{array}$$

$$K{f}_1$$

$${\Omega }_{ij}^{(p)}=({\mathbb{L}}^p{)}_{ii}^{†}+({\mathbb{L}}^p{)}_{jj}^{†}-({\mathbb{L}}^p{)}_{ij}^{†}-({\mathbb{L}}^p{)}_{ji}^{†}.$$

$${\Omega }_{ij}^{(p)}=\sum _{\alpha \ge 2}\frac{(u_{\alpha ,i}-u_{\alpha ,j}{)}^2}{{\lambda }_{\alpha }^p}.$$

$$\begin{array}{r}{C}_p(k)={\left[n^{-1}\sum _j{\Omega }_{kj}^{(p)}\right]}^{-1}={\left[\sum _{\alpha \ge 2}\frac{u_{\alpha ,k}^2}{{\lambda }_{\alpha }^p}+n^{-2}K{f}_p\right]}^{-1},\end{array}$$

$$K{f}_p$$

$$K{f}_2=\frac{n}{4}\sum _{i,j}{{\Omega }_{ij}^{(1)}}^2-\frac{1}{2}\sum _{i,j,k}{\Omega }_{ij}^{(1)}{\Omega }_{jk}^{(1)}+\frac{(K{f}_1{)}^2}{n}.$$

$$m_i{\ddot{\theta }}_i+d_i{\dot{\theta }}_i=P_i-\sum _j{b}_{ij}\sin ({\theta }_i-{\theta }_j).$$

$$\mathbf{M}\delta \ddot{\mathbf{\theta }}+\mathbf{D}\delta \dot{\mathbf{\theta }}$$

$${\mathbb{L}}_{ij}=\left\{\begin{array}[]{cc}-b_{ij}\cos(\theta_{i}^{(0)}-\theta_{j}^{(0)})\,,&i\neq j\,,\\ \sum_{k}b_{ik}\cos(\theta_{i}^{(0)}-\theta_{k}^{(0)})\,,&i=j\,.\end{array}\right.$$

$$m{\ddot{c}}_{\alpha }(t)+d{\dot{c}}_{\alpha }(t)=\delta \mathbf{P}(t)\cdot {\mathbf{u}}_{\alpha }-{\lambda }_{\alpha }{c}_{\alpha }(t),$$

$$c_{\alpha }(t)=$$

$$\times {\int }_0^{t_1}\delta \mathbf{P}(t_2)\cdot {\mathbf{u}}_{\alpha }{e}^{(\gamma -{\Gamma }_{\alpha })t_2/2}\mathrm{d}{t}_2\mathrm{d}{t}_1,$$

$${\mathcal{P}}_1(T)$$

$${\mathcal{P}}_2(T)$$

$${\mathcal{P}}_1(T)$$

$${\mathcal{P}}_2(T)$$

$${\mathcal{P}}_1^{\infty }$$

$${\mathcal{P}}_2^{\infty }$$

$$\begin{array}{rl}{\mathcal{P}}_1^{\infty }& =\frac{{\tau }_0^2}{2d}\sum _{\alpha \ge 2}\frac{(\delta {\mathbf{P}}_0\cdot {\mathbf{u}}_{\alpha }{)}^2}{{\lambda }_{\alpha }},\end{array}$$

$$\begin{array}{rl}{\mathcal{P}}_2^{\infty }& =\frac{{\tau }_0^2}{2md}\sum _{\alpha \ge 2}(\delta {\mathbf{P}}_0\cdot {\mathbf{u}}_{\alpha }{)}^2.\end{array}$$

$$\begin{array}{rl}{\mathcal{P}}_1^{\infty }& ={\tau }_0\sum _{\alpha \ge 2}\frac{(\delta {\mathbf{P}}_0\cdot {\mathbf{u}}_{\alpha }{)}^2}{{\lambda }_{\alpha }^2},\end{array}$$

$$\begin{array}{rl}{\mathcal{P}}_2^{\infty }& =d^{-1}\sum _{\alpha \ge 2}\frac{(\delta {\mathbf{P}}_0\cdot {\mathbf{u}}_{\alpha }{)}^2}{{\lambda }_{\alpha }}\end{array}.$$

$$\begin{array}{rl}{\mathcal{P}}_1^{\infty }(k)& =\frac{\delta P_0^2{\tau }_0^2}{2d}\sum _{\alpha \ge 2}\frac{u_{\alpha ,k}^2}{{\lambda }_{\alpha }}=\frac{\delta P_0^2{\tau }_0^2}{2d}[C_1^{-1}(k)-n^{-2}K{f}_1],\end{array}$$

$${\mathcal{P}}_2^{\infty }(k)$$

$$\begin{array}{rl}{\mathcal{P}}_1^{\infty }(k)& =\delta P_0^2{\tau }_0\sum _{\alpha \ge 2}\frac{u_{\alpha ,k}^2}{{\lambda }_{\alpha }^2}=\delta P_0^2{\tau }_0[C_2^{-1}(k)-n^{-2}K{f}_2],\end{array}$$

$$\begin{array}{rl}{\mathcal{P}}_2^{\infty }(k)& =\frac{\delta P_0^2}{d}\sum _{\alpha \ge 2}\frac{u_{\alpha ,k}^2}{{\lambda }_{\alpha }}=\frac{\delta P_0^2}{d}[C_1^{-1}(k)-n^{-2}K{f}_1]\end{array}.$$

$$⟨{\mathcal{P}}_1^{\infty }⟩$$

$$⟨{\mathcal{P}}_2^{\infty }⟩$$

$$⟨{\mathcal{P}}_1^{\infty }⟩$$

$$⟨{\mathcal{P}}_2^{\infty }⟩$$

$${\lambda }_{\alpha }=4-2\cos (k_{\alpha })-\cos (k_{\alpha }q),\alpha =1,...,n,$$

$$K{f}_p=n\sum _{\alpha \ge 2}[4-2\cos (k_{\alpha })-\cos (k_{\alpha }q){]}^{-p}.$$

Full text

I Introduction

Complex networks are widely used to model nature- as well as man-made coupled dynamical systems Rodrigues et al. (2016). Physical realizations of such systems range from microscopic Josephson junction arrays Wiesenfeld et al. (1998) and interacting molecules in chemical reactions Kuramoto (1984, 1975) to macroscopic high voltage electric power grids Bergen and Hill (1981) and communication or social networks Stankovski et al. (2017); Barabási (2016). Individual elements are represented by nodes in a complex network, which have internal parameters and degrees of freedom. The latter are governed by differential equations that depend on both the internal dynamics of the individual elements and the coupling to the adjacent nodes. Two central questions are (i) how to identify nodes, which, once attacked, perturbed or removed, have the most dramatic effect on the overall dynamics of the coupled system and (ii) how to devise a coupling network guaranteeing robustness of the system against random external perturbations. Attempts to answer such questions are often based on complex network theory, numerically relating dynamical effects to graph-theoretic metrics. This approach has been often criticized, e.g. in Ref. Borgatti (2006); Boldi and Vigna (2014); Hines et al. (2010), because (i) it gives no a priori criterion for which metric should be considered in which situation and (ii) it does not directly indorporate the intrinsic dynamics of the network-coupled system.

Here we propose an altogether different analytical approach. First, we use robustness performance measures that quantify the excursion during the transient dynamics following a perturbation. Second, we spectrally decompose the coupling matrix to calculate the response of the system to some external perturbations.

Third, by direct calculation, we relate the obtained analytical expressions for performance measures (i) to local centralities when analyzing local vulnerabilities, and (ii) to global topological indices when assessing global robustness of the networked system. Following these steps, we identify a new class of local and global topological indices that characterize robustness of synchrony of complex network-coupled oscillators. Our method builds up on investigations of consensus algorithms Bamieh et al. (2012); Grunberg and Gayme (2018), electric power systems Poolla et al. (2017); Paganini and Mallada (2017); Siami and Motee (2016) and coupled oscillators systems Tyloo et al. (2018a, b). Already implicitly present in Refs. Grunberg and Gayme (2018); Paganini and Mallada (2017); Siami and Motee (2016), the Kirchhoff index was first identified as a global robustness quantifier in our earlier work, Ref. Tyloo et al. (2018a). Local vulnerabilities have been more recently connected to centralities related to the resistance distance Grunberg and Gayme (2018); Tyloo et al. (2018b).

In this manuscript, we investigate vulnerabilities and global robustness of synchronous network-coupled oscillators. Frequency synchronization often occurs in such systems when the coupling between individual oscillators is strong enough that they start to oscillate at the same frequency, even when their natural frequency is not homogeneous Strogatz (2004); Pikovsky et al. (2003). Frequency synchronization has attracted a large interest, in particular, the robustness of the synchronous state has been studied from a variety of points of view. One may for instance consider the linear stability of the synchronous state Pecora and Carroll (1998), the range of network parameters where synchrony occurs Barahona and Pecora (2002); Chavez et al. (2005); Zhou et al. (2006), the volume of the basin of attraction of the synchronous state Wiley et al. (2006); Menck et al. (2013); Delabays et al. (2017), the influence of noise on the synchronous state, in particular how it can lead to desynchronization or drive the system to another synchronous state DeVille (2012); Hindes and Schwartz (2016); Schäfer et al. (2017); Hindes and Schwartz (2018); Tyloo et al. (2018c); Hindes et al. (2019), how disturbances spread across the network Kettemann (2016); Wolter et al. (2018); Pagnier and Jacquod (2019), or even how topological changes affect synchrony Coletta et al. (2016); Soltan et al. (2017); Coletta and Jacquod (2019). Here, we investigate the robustness of the synchronous state against external perturbations. For both local and ensemble-averaged perturbations on oscillators with identical dynamical parameters, we find that the robustness of the synchronous state is given by a new family of topological indices based on the resistance distance Klein and Randić (1993); Tyloo et al. (2018a, b). This is illustrated in Fig. 1 which shows the ratio between performance measures defined in Eq. (15a), numerically obtained by perturbing each node of graphs (a) and (e) shown in the inset. Even if the average value of the performance measure is lower for graph (e), the latter can be more strongly sensitive to certain local perturbations. Below we show that specific local vulnerabilities and global averaged robustness are determined by nodal centralities and global topological indices (orange).

The manuscript is organized as follows. Section II recalls the definition of the resistance distance and generalizes it to graphs corresponding to powers of the Laplacian matrix. Section III describes our model of coupled oscillators, briefly discusses synchronized states and evaluates how they respond to external perturbations. Performance measures quantifying this response are also introduced and calculated for quench perturbations. Sections IV numerically illustrates the theory on different graphs for local and global vulnerabilities. An analysis of Kirchhoff indices in both small-world and regular networks is also done. We conclude in Section V.

II Resistance Distances, Centralities and Kirchhoff Indices

The resistance distance $\Omega_{ij}^{(1)}$ is a graph-theoretic metric with an intuitive physical interpretation Klein and Randić (1993). To any graph, one associates an electrical network of resistors whose capacities are given by the inverse of the edge weights. In this case, $\Omega_{ij}^{(1)}$ is the effective resistance between $i$ and $j$, i.e. the voltage that develops between $i$ and $j$ when a unit current is injected at $i$ and collected at $j$ with no injection nor collection at any other node. The resistance distance can be expressed with the network Laplacian matrix $\mathbb{L}$ as

[TABLE]

where $\mathbb{L}^{\dagger}$ is the Moore-Penrose pseudo inverse of $\mathbb{L}$. The resistance distance can be formulated in a convenient way using eigenvectors ${\bf u}_{\alpha}$ and eigenvalues $\lambda_{\alpha}$ of $\mathbb{L}$. It is given by Xiao and Gutman (2003); Coletta and Jacquod (2019),

[TABLE]

where the zero-eigenvector of $\mathbb{L}$ corresponding to $\lambda_{1}=0$ is omitted in the sum. The resistance distance is a graph-theoretic distance metric because (i) $\Omega_{ii}^{(1)}=0$, $\forall i$, (ii) $\Omega_{ij}\geq 0$ $\forall i,j$, and (iii) $\Omega_{ij}^{(1)}+\Omega_{jk}^{(1)}\geq\Omega_{ik}^{(1)}$, $\forall i,j,k$ (triangle inequality) Klein and Randić (1993).

A measure of nodal centrality is given by the inverse of the average resistance distance from any node $k$ to all other network nodes,

[TABLE]

It is a closeness centrality in the usual sense Boldi and Vigna (2014), meaning in particular that large values of $C_{1}(k)$ indicate nodes $k$ that are central in the network according to the resistance distance $\Omega^{(1)}_{ij}$. The second term in bracket on the right-hand-side of Eq. (3) is a graph topological index known as the Kirchhoff index of the network and defined by Klein and Randić (1993),

[TABLE]

where the second equality follows from Eq. (2) Tyloo et al. (2018a).

Until now we have introduced global topological indices and local centralities expressed through resistance distances of the original coupling network. In the upcoming sections, we show how resistance distances naturally come out when quantifying robustness of network-coupled oscillators, but that new distance metrics related to powers of the Laplacian matrix also emerge. We therefore generalize Eqs. (1)–(4) to quantities corresponding to the $p^{\rm th}$ power $\mathbb{L}^{p}$ of the Laplacian matrix ($p\in\mathbb{N}$). This matrix is still a Laplacian matrix, and the associated resistance distance is defined as

[TABLE]

Still using the eigenvectors and eigenvalues of $\mathbb{L}$ we have,

[TABLE]

One can easily check that $\Omega_{ij}^{(p)}$ is still a graph-theoretic distance metric satisfying the properties mentioned between Eqs. (2) and (3). We finally have generalized resistance centralities Tyloo et al. (2018b)

[TABLE]

and generalized Kirchhoff indices Tyloo et al. (2018a)

[TABLE]

We note that generalized resistance distances can in principle be expressed as function of $\Omega_{ij}^{(1)}$. For instance one has

[TABLE]

Below we show how global robustness and local vulnerabilities quantified with performances measures can be expressed in terms of the resistance distance-based centralities and the generalized Kirchhoff indices just introduced.

III Synchronized oscillators under external perturbations

III.1 The Kuramoto model with inertia and its linearization

We consider a set of network-coupled oscillators defined by the following set of coupled differential equations,

[TABLE]

Oscillators labeled $i=1,...,n$ sit on the $n$ nodes of a weighted graph defined by the adjacency matrix with elements $b_{ij}\geq 0$. They have compact angle coordinates $\theta_{i}\in(-\pi,\pi]$, natural frequencies $P_{i}/d_{i}$ rem and inertia as well as damping parameters $m_{i}$ and $d_{i}$. For $m_{i}=0$, Eq. (10) gives the celebrated Kuramoto model on a complex network, for which it is known that when the couplings are sufficiently strong, a finite fraction of, or all oscillators synchronize, i.e. with $\dot{\theta}_{i}-\dot{\theta}_{j}=0$, depending on the distribution of the natural frequencies Kuramoto (1975); Pikovsky et al. (2003); Jadbabaie et al. (2004); Acebrón et al. (2005). Here, we consider $P_{i}$ defined on a bounded, real interval and set $\sum_{i}P_{i}=0$ without loss of generality, so that synchronous states have $\dot{\theta}_{i}=0$, $\forall i$.

Eq. (10) is governed by three sets of time scales. The first one consists of the inverse natural frequencies $d_{i}/P_{i}$. The second one is given by ratios $m_{i}/d_{i}$ and corresponds to the relaxation time of individual oscillators. Finally, the third one is given by the network relaxation times $d_{i}/\lambda_{\alpha}$ defined by the damping parameters and the eigenvalues $\lambda_{\alpha}$ of the weighted Laplacian matrix defined in Eq. (12) below. The first of these sets essentially determines the synchronous state, together with the coupling network. Depending on the other two sets of time scales, perturbations are locally damped or they propagate across the network Pagnier and Jacquod (2019).

We consider a stable fixed-point solution ${\bm{\theta}}^{(0)}=(\theta_{1}^{(0)},\ldots,\theta_{n}^{(0)})$ to Eq. (10) with unperturbed natural frequencies $\bm{P}^{(0)}$. We subject this state to a time-dependent perturbation ${\bm{P}}(t)={\bm{P}}^{(0)}+\delta{\bm{P}}(t)$, which renders angles time-dependent, ${\bm{\theta}}(t)={\bm{\theta}}^{(0)}+\delta{\bm{\theta}}(t)$. Linearizing the dynamics of Eq. (10) about $\bm{\theta}^{(0)}$, one obtains

[TABLE]

where we introduced inertia and damping matrices, ${\bm{M}}={\rm diag}\{m_{i}\}$ and ${\bm{D}}={\rm diag}\{d_{i}\}$, respectively, and the weighted Laplacian matrix ${\mathbb{L}}(\{\theta_{i}^{(0)}\})$ with matrix elements

[TABLE]

This Laplacian is minus the stability matrix of the linearized dynamics, and since we consider a stable synchronous state, it is positive semidefinite, with a single eigenvalue $\lambda_{1}=0$ with eigenvector ${\bf u}_{1}=(1,1,1,...1)/\sqrt{n}$, and $\lambda_{\alpha}>0$, $\alpha=2,3,...n$. From here on, we order the Lyapunov exponents $\lambda_{\alpha}$ in increasing order, i.e. $\lambda_{1}=0<\lambda_{2}<\ldots<\lambda_{n}$.

Eq. (11) can be solved analytically through a spectral expansion if (i) both $\bm{M}$ and $\bm{D}$ commute with $\mathbb{L}$ or (ii) if ${{\bm{M}}^{-1}{\bm{D}}}=\gamma\mathbb{I}$. In case (i), the spectral expansion is over the eigenmodes of $\mathbb{L}$, while in case (ii) it is over the eigenmodes of ${\bm{M}^{-1/2}}\mathbb{L}{\bm{M}^{-1/2}}$ Paganini and Mallada (2017); Coletta and Jacquod (2019). Here, we focus on case (i) with $m_{i}=m$, $d_{i}=d$ $\forall i$.

Expanding the angle deviations over the eigenmodes of $\mathbb{L}$ as $\delta{\bm{\theta}}(t)=\sum_{\alpha}c_{\alpha}(t){\bf u}_{\alpha}$, Eq. (11) leads to a Langevin equation,

[TABLE]

whose general solution reads

[TABLE]

with $\Gamma_{\alpha}=\sqrt{\gamma^{2}-4\lambda_{\alpha}/m}$ and $\gamma=d/m$. Similar expressions have been derived using the transfer function formalism Paganini and Mallada (2017); Guo et al. (2018) or within linear response Manik et al. (2017); Tyloo et al. (2018a); Kettemann (2016). When $\gamma^{2}<4\lambda_{\alpha}/m$, $\Gamma_{\alpha}\in i\mathbb{R}$ and accordingly, $|\Gamma_{\alpha}|$ corresponds to the angular frequency of oscillations along the eigenmode ${\bf u}_{\alpha}$ of $\mathbb{L}$. When on the other hand $\gamma^{2}>4\lambda_{\alpha}/m$, $\Gamma_{\alpha}\in\mathbb{R}$ and gives an additional damping beyond $\gamma$. From Eq. (14), angle and frequency deviations can be calculated as $\delta{\bm{\theta}}(t)=\sum_{\alpha}c_{\alpha}(t){\bf u}_{\alpha}$.

III.2 Performance Measures

The perturbation $\delta{\bm{P}}(t)$ moves the oscillators angles and frequencies away from their value at synchrony and renders them time dependent. For not too strong, finite-time perturbations, they eventually relax to their synchronous values and to assess the magnitude of the excursion away from synchrony, we introduce the following quadratic performance measures

[TABLE]

Similar measures have been discussed in the context of consensus algorithms Bamieh et al. (2012); Grunberg and Gayme (2018), electric power systems Poolla et al. (2017); Paganini and Mallada (2017); Siami and Motee (2016) and coupled oscillators systems Tyloo et al. (2018a, b). The results we are about to present directly connect these performance measures to resistance-distance based centralities and Kirchhoff indices introduced in Section II. While similar connections may have been inferred from some of these works (in particular Refs. Paganini and Mallada (2017); Siami and Motee (2016); Tyloo et al. (2018a)), to the best of our knowledge, it was first unambiguously stated in Ref. Tyloo et al. (2018b).

Because synchronous states are defined modulo any homogeneous angle shift, they are unaffected by the transformation $\theta_{i}^{(0)}\rightarrow\theta_{i}^{(0)}+C$. Accordingly, only angle shifts with $\sum_{i}\delta\theta_{i}(t)=0$ matter, which is incorporated in the definitions of ${\mathcal{P}}_{1,2}$ by subtracting averages $\Delta(t)=n^{-1}\sum_{j}\delta\theta_{j}(t)$ and $\dot{\Delta}(t)=n^{-1}\sum_{j}\delta\dot{\theta}_{j}(t)$. If the perturbation is not too strong and finite in time, both $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ are finite even for $T\rightarrow\infty$. Low values for ${\mathcal{P}}_{1,2}^{\infty}\equiv{\mathcal{P}}_{1,2}(T\rightarrow\infty)$ indicate then that the system absorbs the perturbation with little fluctuations, while large values indicate a temporary fragmentation of the system into independent pieces – qualitatively speaking, ${\mathcal{P}}_{1,2}^{\infty}$ measures the coherence of the synchronous state Bamieh et al. (2012).

Using the spectral expansion with coefficients given in Eq. (14), the performance measures of Eqs. (15) read, in our case of homogeneous inertia and damping coefficients

[TABLE]

Performance measures depend on the perturbation vector $\delta\bm{P}(t)=\delta\bm{P}_{0}\,f(t)$, which may have different time dependences $f(t)$ – such as, for instance, noisy fluctuations or instantaneous, Dirac-delta perturbations – and different geographical dependences encoded in $\delta\bm{P}_{0}$. In this manuscript we consider quenches where $f(t)$ vanishes outside some time interval, inside which it is constant but nonzero. In the next section we calculate performance measures for general perturbation vectors $\delta\bm{P}_{0}$ for such quenches. As for geographical dependences, we then consider two cases of (i) nodal vulnerabilities, with local perturbations $\delta\bm{P}_{0}=(0,...,\delta P_{0,k},...,0)$ and (ii) global robustness, where performance measures are averaged over all possible locations $k$ for the perturbation.

III.3 Quench Perturbation

We compute both performance measures ${\mathcal{P}}_{1,2}$ for a quench perturbation $\delta\bm{P}(t)=\delta\bm{P}_{0}\,\Theta(t)\,\Theta(\tau_{0}-t)$ with the Heaviside function $\Theta(t)$ and a perturbation vector $\delta{\bm{P}}_{0}$ encoding the geographical distribution of the perturbation. The duration $\tau_{0}$ of the quench allows to explore the different time scales of the system and we show below that ${\mathcal{P}}_{1,2}$ varies significantly depending on $\tau_{0}$. Using Eq. (14), Eqs.(15) give

[TABLE]

It is easily checked that ${\mathcal{P}}_{1,2}^{\infty}\in\mathbb{R}$ in both cases $\gamma^{2}>4\lambda_{\alpha}/m$ (with $\Gamma_{\alpha}\in\mathbb{R}$) and $\gamma^{2}<4\lambda_{\alpha}/m$ (with $\Gamma_{\alpha}\in i\mathbb{R}$).

Both performance measures are given by a spectral sum of terms corresponding to the eigenmodes of the network Laplacian matrix $\mathbb{L}$. Each term in this sum depends on the scalar product of the perturbation vector $\delta\bm{P}_{0}$ with the eigenmodes ${\bf u}_{\alpha}$ of $\mathbb{L}$ times a mode-dependent factor. The latter is an almost always decreasing function of the eigenvalues $\lambda_{\alpha}$. Therefore, Eqs. (17) suggest that the largest excursion can be obtained by overlapping $\delta{\bm{P}}_{0}$ with few of the lowest-lying eigenmodes of $\mathbb{L}$, in particular ${\bf u}_{2}$, the so-called Fiedler mode of the network Fiedler (1973).

To get more insight into Eqs. (17), we compute their two asymptotic limits of long and short $\tau_{0}$. For perturbations with very short duration i.e. $\tau_{0}\ll m/d$, $(\gamma\pm\Gamma_{\alpha})^{-1}$, we have,

[TABLE]

Each term in the sum over modes depends on $\lambda_{\alpha}$ for $\mathcal{P}_{1}^{\infty}$ but not for $\mathcal{P}_{2}^{\infty}$. Consequently, $\mathcal{P}_{1}^{\infty}$ depends explicitly on the location of the perturbation, while there is no such dependence for $\mathcal{P}_{2}^{\infty}$, which depends only on the squared norm of the perturbation vector $\delta{\bm{P}}_{0}$ orthogonal to ${\bf u}_{1}$. This reflects the fact that in the regime of short $\tau_{0}$, the perturbation does not act long enough to change the kinetic energy of inertiafull oscillators, which $\mathcal{P}_{2}^{\infty}$ essentially measures. Consequently, the perturbation is quickly damped locally, with little dependence on its location in the situation we consider of homogeneously distributed inertia. We note that similar topology-independent results were obtained for instantaneous, Dirac-delta perturbations Bamieh and Gayme (2013).

In the other limit $\tau_{0}\gg m/d$, $(\gamma\pm\Gamma_{\alpha})^{-1}$, the performance measures read

[TABLE]

In this case of a long-lasting perturbation, both $\mathcal{P}_{1}^{\infty}$ and $\mathcal{P}_{2}^{\infty}$ depend on the location of the perturbation. Furthermore, and perhaps more importantly, the inertia affects neither $\mathcal{P}_{1}^{\infty}$ nor $\mathcal{P}_{2}^{\infty}$. This is so since, for long quenches, oscillators have the time to synchronize at a new frequency with zero angular acceleration before the perturbation is over. We further note that $\tau_{0}$ no longer appears in $\mathcal{P}_{2}^{\infty}$, since the latter considers deviations orthogonal to ${\bf u}_{1}$. Consequently, the whole time spent by the oscillators at the new frequency does not contribute to $\mathcal{P}_{2}^{\infty}$. Most importantly, Eqs. (18) and (19) suggest that in both asymptotic limits of short and long perturbations, $\mathcal{P}_{1,2}^{\infty}\propto\sum_{\alpha\geq 2}(\delta{\bm{P}}_{0}\cdot{\bf u}_{\alpha})^{2}/\lambda_{\alpha}^{p}$ with $p=0,1,2$. That result was already hinted at in Ref. Tyloo et al. (2018a) for inertialess oscillators and various types of perturbations. Below we show how this dependence leads to performance measures depending on the resistance distances, centralities and Kirchhoff indices introduced in Section II.

Eqs. (17) and their asymptotic limits, Eqs. (18) and (19), give the performance measures $\mathcal{P}_{1,2}^{\infty}$ for any perturbation vector $\delta{\bm{P}}_{0}$. We next discuss two important cases of (i) a single-node perturbation, $\delta P_{0,i}=\delta P_{0}\delta_{ik}$, where large values of the node-dependent performance measures $\mathcal{P}_{1,2}^{\infty}\rightarrow\mathcal{P}_{1,2}^{\infty}(k)$ identify local vulnerabilities and (ii) averaged perturbation over ensemble of homogeneously distributed perturbation vectors $\delta{\bm{P}}_{0}$, where large values of $\mathcal{P}_{1,2}^{\infty}\rightarrow\langle\mathcal{P}_{1,2}^{\infty}\rangle$ indicate globally fragile networks.

III.4 Specific Local Vulnerabilities

To assess local vulnerabilities of the coupled oscillators, we apply a quench perturbation on a single node. The vulnerability of node $k$ is then given by Eqs. (17) with the components of the perturbation vector given by $\delta P_{0,i}=\delta P_{0}\delta_{ik}$. In the limit of short duration of perturbation, $\tau_{0}\ll m/d$, $(\gamma\pm\Gamma_{\alpha})^{-1}$, one obtains

[TABLE]

where the right-hand side of Eq. (20a) directly follows from Eq. (3). For a perturbation on node $k$, $\mathcal{P}_{1}^{\infty}(k)$ is expressed in terms of the centrality, $C_{1}(k)$, a local nodal descriptor, and the Kirchhoff index $Kf_{1}$, a global network descriptor. Consequently, the most vulnerable nodes in a given network, according to $\mathcal{P}_{1}^{\infty}(k)$, are identified by their resistance-distance based centrality $C_{1}(k)$.

In the other limit of long perturbations, $\tau_{0}\gg m/d$, $(\gamma\pm\Gamma_{\alpha})^{-1}$, Eqs. (19) give

[TABLE]

This time $\mathcal{P}_{1}^{\infty}$ is given by the higher order centrality $C_{2}(k)$ and Kirchhoff index $Kf_{2}$.

When considering a given, fixed network, Eqs. (20) and (21) show that perturbations on the most central nodes – as measured by either centrality $C_{1}$ or $C_{2}$ – give the smallest overall responses, except when considering $\mathcal{P}^{\infty}_{2}(k)$ for a short-time perturbation. In that latter case, the response is homogeneous and perturbing any node leads to the same performance measure $\mathcal{P}^{\infty}_{2}(k)$. When comparing two nodes with similar centrality on two different networks, on the other hand, Eqs. (20) and (21) indicate that the largest response occurs on the network with smallest generalized Kirchhoff index – except again for $\mathcal{P}^{\infty}_{2}(k)$ and a short-time perturbation. We show below that the overall network robustness is actually given by these generalized Kirchhoff indices, which makes this observation quite counterintuitive : when perturbing two nodes of equal centrality on two different networks, the largest response is actually recorded on the overall more robust network ! We will come back to this point below.

III.5 Averaged Global Robustness

We next assess the global robustness of synchrony in a given network, by averaging Eqs. (17) over an homogeneously distributed ensemble of perturbation vectors defined by $\langle\delta P_{0,i}\delta P_{0,j}\rangle=\delta_{ij}\langle\delta P_{0}^{2}\rangle$ Tyloo et al. (2018a). Averaging Eqs. (17) gives, in the limit of short perturbations, $\tau_{0}\ll m/d$, $(\gamma\pm\Gamma_{\alpha})^{-1}$

[TABLE]

We see that $\langle\mathcal{P}_{1}^{\infty}\rangle$ is given by the Kirchhoff index $Kf_{1}$ which is proportional to the network’s average resistance distance $\Omega_{ij}^{(1)}$ [see Eq. (4)]. Similarly to the local vulnerability in this limit, $\langle\mathcal{P}_{2}^{\infty}\rangle$ depends on the network only marginally through the number of nodes.

In the other limit $\tau_{0}\gg m/d$, $(\gamma\pm\Gamma_{\alpha})^{-1}$, Eqs. (19) give

[TABLE]

Both performance measures depend on generalized Kirchhoff indices. Quite remarkably and as for local vulnerabilities, the only average performance measure that depends on inertia is $\langle\mathcal{P}_{2}^{\infty}\rangle$ in the short $\tau_{0}$ limit. In the next Section, we numerically confirm the validity of the analytical theory presented in this Section.

IV Numerical Results

IV.1 Local Vulnerabilities and Resistance Centralities

We numerically investigate local vulnerabilities by perturbing individual nodes with the quench perturbation discussed above. Our theory applies to network of any geometry with any number $n$ of nodes. However in order to better visualize the agreement between analytical predictions and numerical results we restrict ourselves to relatively small graphs with $n=20$ nodes of the kind shown in Fig. 2.

We check Eqs. (17) for the model defined in Eq. (10) with $b_{ij}=1$ on the edge of the graph considered and $b_{ij}=0$ otherwise, $m_{i}\equiv m=1$ and $d_{i}\equiv d=1$. We numerically time-evolve Eq. (10) with a fourth-order Runge-Kutta method, following a perturbation $\delta{P_{i}}(t)=\delta{P}_{0}\delta_{ik}\,\Theta(t)\,\Theta(\tau_{0}-t)$ away from ${\bm{P}}^{(0)}=0$ and starting from the corresponding synchronous state ${\bm{\theta}}^{(0)}=0$. Fig. 3 shows that the theory of Eqs. (17) is in perfect agreement with numerical results. In particular, one clearly sees the crossover from $[C_{1}^{-1}(k)-n^{-2}Kf_{1}]$ to $[C_{2}^{-1}(k)-n^{-2}Kf_{2}]$ for $\mathcal{P}_{1}^{\infty}$ (dotted to dashed lines on the left panels) and from a constant to $[C_{1}^{-1}(k)-n^{-2}Kf_{1}]$ (dotted to dashed line on the right panels) for $\mathcal{P}_{2}^{\infty}$, as $\tau_{0}$ increases. This fully confirms our theoretical predictions, Eqs. (20)-(21). We conclude that, generally speaking (i.e. except for $\mathcal{P}_{2}^{\infty}$ and short perturbations), the most central nodes are the most robust. They are connected by multiple paths to the rest of the network, and when they are perturbed, the disturbance quickly diffuses across the network with small angle differences. In contrast, the most peripheral nodes such as dead ends have only few paths connecting them to the bulk of the network and the disturbance diffuses across the network with large angle differences. It has been numerically found that dead ends undermine grid stability Menck et al. (2014), and our results shed some analytical light on that observation.

We further illustrate this strong connection between resistance centralities and response of the system. We show in Fig. 4 resistance centralities $C_{1}(k)$ and $C_{2}(k)$ for the six graphs of Fig. 2. One sees that $C_{1}(k)$ and $C_{2}(k)$ tend to become higher while going from graph (a) to (f) indicating that graphs with more rewired edges (and thus with more long-range couplings) have shorter distances between nodes and thus lower Kirchhoff indices. Interestingly, several nodes with a high centrality $C_{1}(k)$ do not necessarily have a high centrality $C_{2}(k)$, and vice-versa. We then show in Fig. 5 the time-evolution of angles and frequencies following a local quench perturbation on two different nodes of graph (f) with very different resistance centralities. One clearly sees that for a perturbed node with low resistance centrality (Fig. 5, top), angles and frequencies spread more during the perturbation than for a node with higher centrality (Fig. 5, bottom).

Generally speaking, networks with higher rewiring probabilities have smaller global topological indices $Kf_{1}$ and $Kf_{2}$ and thus smaller $\langle\mathcal{P}_{1,2}^{\infty}\rangle$ according to our theory. This is confirmed numerically in the four left panels in Fig. 6, where we apply the same quench perturbation on nodes with resistance centralities $C_{1}(k)$ close to their median value in the corresponding graph. One observes that angles and frequencies spread more and take more time to return to the initial fixed point in the graph with higher $Kf_{1}$ and $Kf_{2}$ (top) compared to the one with more rewired edges (bottom).

While this is a rather general rule, it does not forbid exceptions. As a matter of fact, specific perturbations can lead to higher response in a network with lower Kirchhoff index than in a network with higher Kirchhoff index. Such an exception is illustrated in the four right panels in Fig. 6, where the same quench perturbation is applied on nodes with similar resistance centralities $C_{1}(k)$ but belonging to graphs with very different Kirchhoff indices (see insets of Fig. 6). As expected from Eqs. (21), if two nodes on different networks have the same centralities, then, a perturbation applied on the one in the network with lower Kirchhoff index produces the largest response. Another illustration of this effect is given in Fig. 1, where graph (e) is more robust than graph (a) on average (dashed lines). But if we compare the response to specific local perturbations, some nodes of graph (a) are more robust than those of graph (e) (red crosses). Both the generic and the exceptional behaviors are accurately captured by our theory.

IV.2 Global Robustness and Generalized Kirchhoff Indices

We next investigate global robustness by averaging performance measures over an ensemble of perturbation vectors located on a single node, $\delta{\bm{P}}_{0}=(0,..,\delta P_{0,k},0,...)$ with $k=1,...,n$. Fig. 7 compares the resulting numerical averages $\langle\mathcal{P}_{1,2}^{\infty}\rangle$ with the average of the theoretical prediction of Eqs. (17). Numerics and theory agree well. In particular the left panel confirms nicely the crossover between $Kf_{1}$ and $Kf_{2}$ predicted by Eqs. (22a) and (23a). A similar behavior is visible in the right panel, where $\langle\mathcal{P}_{2}^{\infty}\rangle$ does not depend on the network topology for short duration of perturbation (black and blue lines and symbols) but crosses over to $Kf_{1}$ as $\tau_{0}$ increases, as predicted by Eqs. (22b) and (23b).

We finally note that networks with high $Kf_{1}$ do not necessarily have a high $Kf_{2}$, and vice-versa. This is illustrated in Fig. 7 where the chosen network with $p=0.15$ has a higher $Kf_{2}$ but a lower $Kf_{1}$ than the chosen network with $p=0.3$. Below we analyze in more details $Kf_{1,2}$ in randomly rewired networks.

IV.3 Generalized Kirchhoff Indices in Small-World Networks

The results obtained above relate local vulnerabilities to nodal centralities and global network robustness to generalized Kirchhoff indices. This connection is powerful : it gives a vulnerability ranking of nodes and provides robustness assessment based on well-defined, easily calculated network descriptors. To gain qualitative insight on what favors robustness in a graph, we investigate the behavior of the Kirchhoff indices for Watts-Strogatz, randomly rewired networks. Following Ref. Watts and Strogatz (1998), we consider initially regular, circular graphs where nodes are coupled to their nearest, second-nearest aso. up to their 10th neighbors. Each edge in the corresponding coupling network is then rewired with probability $p$. Fig. 8 compares the standard measures of ”nearest-neighborness” of geodesic centrality $l$ and clustering coefficient $Cl$ with the generalized Kirchhoff indices $Kf_{1}$ and $Kf_{2}$, as a function of $p$.

Both Kirchhoff indices drop, roughly following $l$, as $p$ is increased, with $Kf_{2}$ decreasing significantly faster than $Kf_{1}$ and $l$. Traditionnally, the ”small-world” behavior occurs around $p=0.01$, where $l$ is significantly smaller than its initial value, while $Cl$ has not yet changed much. In that region, $Kf_{1}$ has been reduced to $\sim 40$% of its initial value, while $Kf_{2}$ reaches only few percents of its initial value. Accordingly, small-world networks are significantly more robust to external perturbations than regular networks, particularly when considering $\mathcal{P}_{1}^{\infty}$ for long quenches. Only a fraction of edges need to be rewired to achieve a level of robustness comparable to that of random networks. As a side-remark, we note that the ratio of Kirchhoff indices provides for a clear identification of small-world networks, which correspond to values of $p$ where $Kf_{1}(p)/Kf_{2}(p)$ is fast increasing with $p$.

IV.4 Regular Networks

We finally comment on regular networks. In such networks, all the nodes are equivalent and therefore global robustness is equivalent to local vulnerability, $\mathcal{P}_{1,2}^{\infty}(k)=\langle\mathcal{P}_{1,2}^{\infty}\rangle$, $\forall k$, furthermore, Kirchhoff indices can be calculated analytically. The Laplacian matrix can be diagonalized with a Fourier transform, and its spectrum is given by

[TABLE]

with $k_{\alpha}=2\pi(\alpha-1)/n$. Kirchhoff indices Eq. (8) are then given by,

[TABLE]

Fig. 9 shows $Kf_{1}$ and $Kf_{2}$ for such regular networks with $n=50$ nodes. When extending the coupling range $q$, Kirchhoff indices are generally decreasing, indicating the standard trend that longer-range couplings reduce centralities. However, for some values $q=10,17,24$ equal or close to integer divisors of $n$, Kirchhoff indices suddenly become larger. This is so, since then, paths made of few long-range interactions form either closed or almost closed loops (see the inset of Fig. 9 for $q=17,24$), which do not reduce the geodesic distance between many pairs of nodes, compared to long range coupling with $n/q$ not integer (e.g. $q=19$ in Fig. 9). Consequently, graphs that may appear similar, such as those with $q=17$ and $q=19$ or with $q=23$ and $q=24$ may exhibit Kirchhoff indices differing by factors of 2-4 or even more. This illustrates how assessing global robustness is hard to do from a network’s general appearance and/or from arguments solely based on the existence of long-range couplings.

V Conclusion

Building up on earlier works Bamieh et al. (2012); Grunberg and Gayme (2018); Poolla et al. (2017); Paganini and Mallada (2017); Siami and Motee (2016); Tyloo et al. (2018a, b), we have investigated the response under external perturbations of network-coupled dynamical systems initially in a stable synchronous state. We proposed to assess network robustness and identify nodal vulnerabilities through quadratic performance metrics quantifying the magnitude of the perturbation-induced transient excursion away from the synchronous state. As we reported earlier for first-order oscillators Tyloo et al. (2018a), we found that the response of inertiaful, second-order oscillators depends on the overlap between the perturbation vector and the eigenmodes of the weighted Laplacian. In particular, the set of nodes located on the slowest eigenmode corresponding to the smallest eigenvalue produces the largest excursions when perturbed. Considering disturbances localized on a single node we found that, oscillators which, once perturbed, induce the largest transient excursion are the ones with smallest resistance centralities. Extending the results of Ref. Tyloo et al. (2018a) to second-order oscillators, we found that global robustness, assessed by averaging performance measures over ergodic ensembles of perturbation vectors, is also given by generalized Kirchhoff indices. A network can then be made more robust to perturbations by minimizing its average resistance distances, for instance by introducing long-range edges. Quite remarkably, except for $\mathcal{P}_{2}$ and short time perturbation, asymptotic behaviors of performance measures in either limit of long or short perturbations do not depend on the inertia of the oscillators.

Our findings are rather general. Together with Refs. Tyloo et al. (2018a, b), they make it clear that, almost regardless of the presence of inertia, and of the type of perturbation chosen, quadratic performance measures are given by the generalized resistance distance-based centralities or, once averaged over ergodic ensembles of perturbations, by the generalized Kirchhoff indices that we introduced in Section II. These local and global network characteristics therefore provide well-defined, numerically easy to calculate robustness descriptors and local vulnerability indicators.

Further studies could consider the effect of spatially correlated perturbations and go beyond the assumption of homogeneous inertia and damping.

VI Acknowledgments

This work has been supported by the Swiss National Science Foundation under Grant PYAPP2_154275.