Engineering mesoscale structures with distinct dynamical implications

Over the past decade networks have become the principal tool for the analysis of complex systems [1]. By describing a given complex system as a network of discrete nodes that interact via discrete links, a considerable simplification is achieved, but the complexity of the topology of interactions (and thus many of the emergent properties of the real world system) is retained.

In some of the biggest success stories of graph and network theory, the simplification that is achieved by treating a system as a network enables a subsequent analytical investigation that links an emergent phenomenon to local properties of the nodes. Examples include for instance Euler's solution of the Königsberg bridge problem and the giant-component transition in random graphs [2].

It is apparent that not all emergent phenomena can be traced back to node properties alone. Instead, progress has been made by linking phenomena to well-studied global properties such as the spectra of matrices. In particular, the spectra of the adjacency matrix and the different graph Laplacians contain information about the net's diameter [3], the degree of modularity [4] and the isoperimetry (i.e. how many edges have to be removed to cut it in pieces of a certain size) [3]. Also in studies of the dynamics on networks clever mapping to a spectral problem has often proved instrumental. Prominent examples include for instance the computation of percolation thresholds of networks [5] and the master stability function approach to synchronization[6].

A direct connection between dynamics and spectral properties can be made for dynamical systems close to steady states [7]. Here the relevant matrix is the Jacobian, which consists of the coefficients of a linearized system at the steady state. In a given dynamical system a steady state is stable if all eigenvalues of the Jacobian matrix have negative real parts. Notably, this connection between the stability of steady states and the spectrum of the Jacobian holds in all dynamical systems including strongly heterogeneous ones.

For instance the seminal paper by May [8] used spectral reasoning to link network dynamics to global structural properties. The paper showed that large densely connected networks are unlikely to be stable, a result that was recently confirmed for a broad class of dynamics [9, 10].

While linking emergent phenomena to spectral properties has led to considerable progress, it does not provide a satisfactory solution for the investigation of many real world networks, because the size of the network can be such that computation of the spectrum is infeasible and, more importantly, information on the network structure is often incomplete. Perhaps inspired by Alon [11], this has triggered a recent search for mesoscale subgraphs that permit conclusions on the spectra of important matrices and by extension emergent phenomena. In the dynamical context one thus seeks so-called dynamical motifs, whose presence in a given network has well-defined dynamical implications [12, 13].

A direct search for such motifs was undertaken for instance in [14], which studied the impact of three-node subgraphs on the stability of ecological food webs. Furthermore, it was observed that, compared to statistical null-models, the spectra of real world networks [15–20] display characteristic peaks, which can in part be linked to the statistical over-representation of certain subgraphs [18–21]. Thus, for instance, the multiplicity of the eigenvalue 0 of a graph's Laplacian matrix has been found to increase with the number of star-like subgraphs with two or more leaf nodes [15, 22, 23].

Mathematically speaking, an eigenvalue of a network can be linked to a specific subgraph if the corresponding eigenvector v is localized on the subgraph. In this case the eigenvector v and the corresponding localized eigenvalue ${\mathrm { Ev}}({\boldsymbol{v}})$ are independent of the topology outside the subgraph [20]. Thus, the localized eigenvector and the corresponding eigenvalue are preserved in all networks, in which the subgraph occurs.

For a long time, the only well-studied examples for subgraphs with localized eigenvalues were star-like subgraphs, and complete subgraphs (cliques) [18]. Starting from the observation that spectral peaks are often associated with network symmetries [24], MacArthur and co-workers determined generic symmetry properties that lead to the localization of eigenvectors [25, 26]. By describing these properties in terms of a net's automorphism group, they provided means to construct and classify arbitrarily large subgraphs with localized eigenvalues.

In the present paper we analyse the dynamical implications of subgraphs with localized eigenvectors in a system of delay-coupled delay oscillators, which has been previously proposed as a model for the ecological metapopulations and gene-regulatory nets [28]. In contrast to the food web model, for which some dynamical implications of symmetric structures are discussed in [27], the metapopulation model lends itself to a deeper analysis as it is a so-called factorizing system [28]. In these systems the effect of topology can be 'factored-out' such that eigenvalues of the Jacobian, which characterize the dynamics, can be written as functions of eigenvalues of the adjacency matrix, which characterizes the structure. While we exploit this property extensively to provide analytical results that illustrate our findings, we note that these findings are not limited to factorizing systems. As the main result of the present paper, we show that it is possible to engineer subgraphs that have distinct dynamical consequences, irrespective of the embedding network.

This paper is organized as follows. We start in section 2 by recalling the conditions under which a mesoscale structure has localized eigenvalues. In section 3, we introduce a class of models of delay-coupled delay oscillators and derive bifurcation conditions that explicitly depend on the topological eigenvalues. In section 4, we study model realizations that include symmetric structures and show that the localized eigenvalues give rise to distinct instabilities of the system.

The topology of a network of N nodes can be described by an N × N matrix, the so-called adjacency matrix

$$\begin{equation} {\rm \bf A}: A_{ij}=\left\{\begin{array}{@{}ll@{}} 1 & {\rm if~there~is~a~link~pointing~from}~{\it j}~{\rm to}~{\it i}, \\ 0 & {\rm otherwise}, \end{array}\right. \end{equation}\noindent \tag{ 1 }$$

where A_ij = A_ji, for undirected networks.

We say that a vector ${\boldsymbol{v}}\in \mathbb {C}^N$ is a topological eigenvector of a network if it is an eigenvector of the network's adjacency matrix A. In the remainder of this section, we discuss the conditions under which topological eigenvectors can be attributed to a specific subgraph. The discussion presented here aims to provide an intuitive understanding. For a rigorous analysis we refer the reader to [25, 26].

Let us consider a small subgraph that is part of a bigger network. We group the nodes of the network into three sets: the set I of inner nodes, which are in the subgraph and adjacent only to nodes of the subgraph; the set C of connecting nodes, which are in the subgraph but not in I and thus constitute the interface between the subgraph and the residual graph; and the set R of residual nodes, which are not in the subgraph (cf figure 1).

Zoom In Zoom Out Reset image size

Below we say that a topological eigenvector is localized on a given subgraph, if the components of the eigenvector vanish in all nodes except those belonging to the set I, i.e. v_j = 0 ∀j∈R,C.⁴ The defining condition for a localized eigenvector is thus

$$\begin{equation} \sum_{k}A_{jk}v_k=c\cdot v_j \left\{\begin{array}{@{}ll@{}} = 0 & \mathrm{for \ all} \ j\in R,C,\\ \in \mathbb{C} & \mathrm{for \ all} \ j\in I. \end{array}\right. \end{equation} \tag{ 2 }$$

Let us discuss the implications of equation (2). First, consider nodes j∈R. By definition, these nodes can only have neighbours k∈R,C. As v_k = 0 ∀k∈R,C, we can conclude that (i) equation (2) is trivially fulfilled for j∈R, and that (ii) the values of all coefficients A_jk, j∈R and k∈R∪C can be changed without changing the right hand side of any of equations (2). In other words, changes of the residual graph's topology, its size and/or its coupling to connecting nodes do not affect the eigenvector property of v.

Next, consider nodes j∈C, i.e. the connecting nodes. These nodes have neighbours ∈R,C,I, which contribute to the sum in equation (2). As v_k = 0 for k∈R,C, contributions from all R-neighbours and C-neighbours of a C-node vanish. Thus, $\sum _{k}A_{jk}v_k=0$ requires that contributions from all I-neighbours of a C-node cancel each other⁵.

Finally, consider nodes j∈I. The respective entries v_j can be non-zero but are subject to two conditions imposed by equation (2): firstly, as shown above, the first line of equation (2) requires that there are entries v_j j∈I that sum up to zero. Secondly, for the general case of complex $v_{j}=:p_{j}{\rm e}^{\imath q_{j}}$, the second line of equation (2) requires that $\sum _{k}\frac { A_{jk}p_k}{p_j}\;\mathrm {e}^{\imath (q_k-q_j)}=c$ for all j∈I, i.e. that the accumulated phase difference of a node with respect to its neighbours is identical for all nodes j∈I. Both conditions can be fulfilled if the v_j are equally spaced on the unit circle, which is the case if the nodes j∈I are subject to symmetry relations [25].

Let us briefly summarize. The reasoning above showed that eigenvectors localized on a subgraph are retained in any network, in which an arbitrary residual graph is attached in an arbitrary manner to the connecting node(s). Further, we motivated that localized eigenvectors occur on all subgraphs whose inner nodes are subject to symmetry relations. In the following we call those subgraphs basic symmetric subgraphs (BSS). Each BSS thus has one or more localized eigenvectors, whose corresponding eigenvalues appear in the spectrum of every network that contains the BSS.

A classification of BSS with only undirected links and a recipe for their construction can be found in [25, 26]. Because a generalization of the scheme to directed BSS is straightforward, the present paper focuses on dynamical implications, rather than providing a wealth of examples for topological localization. For illustration, we thus restrict ourselves to three generic BSS with both, directed and undirected links, which are shown in table 1.

Table 1. Topological eigenvalues and eigenvectors of three BSS. Shown are a star-like subgraph, a directed triangle and a directed hexagon with one or two connecting nodes respectively. Inner nodes are denoted by black circles and enumerated counter-clockwise beginning from the upper left corner. Connecting nodes are denoted by open circles and enumerated last. Lines with arrows denote directed links, lines without arrows undirected links. We use the following abbreviations: $\phi _n:=\exp (2\pi n \mathrm {i}/3)$, $\xi _\pm :=(1\pm \sqrt {13})/2$, and $\zeta _\pm :=(1\pm \sqrt{7})$; + and − denote +1,−1 respectively. Localized eigenvectors are distinguished by vanishing values of the eigenvector on the connecting nodes (after the '|'). Additionally we have marked the localized eigenvectors by a star (*).

In the following, we present a generalized model of delay-coupled delay oscillators, which we use to study the dynamical implications of localized topological eigenvalues. The class of models has previously been proposed in [28], but is here extended to the case of directed links which has not previously been considered.

Delay networks constitute a powerful, and intensively studied framework for modelling biological systems [29–35]. Modelling aspects of biology often requires us to reduce complex processes to simple steps within a mechanism. Delay terms in models can account for time requirements of complex processes, such as transcription, gestation, etc, without resolving the process in detail [36]. In the context of the present paper, delay networks also offer another advantage: even small delay systems can exhibit complex dynamics, which allows us to study nontrivial dynamics in relatively small networks [37]. In systems without delay, very similar results could be produced, but would require larger examples and a considerably more cumbersome notation.

The analytical treatment of delay networks is generally demanding as even simple delay differential equations constitute infinite dimensional dynamical systems and lead to transcendental equations [38]. We therefore focus on a class of delay-coupled delay networks on degree homogeneous networks, which are known to fall into a class of factorizing systems [28] and thus can be treated analytically.

3.1. The metapopulation model

We consider a directed network of N nodes, corresponding to different patches that can be inhabited by a given species. Each node i has an internal dynamical variable X_i representing the abundance of the population within the patch. In time the variable X_i changes due to internal dynamics and due to coupling to other variables X_j according to

$$\begin{equation} \dot{X}_i=G(X_i^{\tau})-L(X_i)+\sum_{j} \left[A_{ij}F(X_j^{\delta})-A_{ji}F(X_i)\right], \end{equation} \tag{ 3 }$$

where X^δ_i and X^τ_i denote values of the variable at an earlier time given by the travel time delay δ and the growth delay τ. Moreover, A_ij and A_ji denote entries of the adjacency matrix, and G, L, and F are positive functions describing growth, loss, and migration, respectively. In the following we do not restrict the functions G, L, and F to specific functional forms, but consider formally the whole class of models, using the approach of generalized modelling [39–41]. The class of models under consideration thus includes several well-studied examples such as the Mackey–Glass [29] and the Ikeda model [30].

3.2. Generalized stability analysis

For simplicity we study equation (3) on a network, in which every node has as many incoming links dⁱⁿ_i as outgoing links d^out_i, such that dⁱⁿ_i = d^out_i = d_i. In this class of networks, which includes all undirected nets, flows tend to balance unequal node loads. Thus, if the functions G, L and F are identical for all nodes, there is a homogeneous steady state, in which all variables X*_i = X* are identical.

We analyse the stability of the steady state at X* by introducing normalized variables x_i = X_i/X* and normalized functions f(x_i) = F(x_iX*)/F(X*), g(x_i) = G(x_iX*)/G(X*), and l(x_i) = L(x_iX*)/L(X*). Using X* = X^τ* = X^δ*, we can rewrite equation (3) as

$$\begin{equation} \dot{x_{i}}=\alpha(g(x_{i}^{\tau})-l(x_i))+\beta\sum_{j}\left[A_{ij}f(x_j^{\delta})-A_{ji}f(x_i)\right], \end{equation} \tag{ 4 }$$

where α = G(X*)/X* = L(X*)/X* and β = F(X*)/X* are parameters that can be interpreted as characteristic turnover and transport rates⁶. We set α = 1 by timescale normalization.

The Jacobian matrix that governs the stability of steady states in the system can be expressed as a function of the quantities g', l', and f' that are defined as derivatives of the normalized functions with respect to the normalized variable s_i, or equivalently the logarithmic derivatives of the functions before normalization

$$\begin{equation} {g'}=\left.{\frac{{\rm \partial} g}{{\rm \partial x}}}\right|_{x=1}=\left.\frac{X^*}{G^*}{\frac{{\rm \partial} G}{{\rm \partial X}}}\right|_{X=X^*}=\left.{\frac{{\rm \partial} \log{G}}{{\rm \partial \log{X}}}}\right|_{X=X^*}. \end{equation}\noindent \tag{ 5 }$$

As we did not restrict the functional forms in the model, the three quantities g', l' and f' are unknown constants that can be treated as parameters of the system. These parameters are generally easily interpretable in the context of the application [39].

Following [28] we write the Jacobian as

$$\begin{equation} {\rm \bf J}: J_{ij}(\lambda)=J_{i}^{\mathrm{d}}\delta_{ij}+J^{\mathrm{o}}A_{ij}, \end{equation}\noindent \tag{ 6 }$$

where δ_ij is the Kronecker delta,

$$\begin{eqnarray} &&J_i^{\mathrm{d}}=(g'\exp(-\lambda\tau)-l')-d_i\beta f', \end{eqnarray}\noindent \tag{ 7 }$$

$$\begin{eqnarray} &&J^{\mathrm{o}}=\beta f' \exp(-\lambda \delta) \end{eqnarray} \tag{ 8 }$$

and λ is a self-consistent eigenvalue of J(λ).

3.3. Bifurcation conditions

In a continuous-time dynamical system, a steady state is locally asymptotically stable if all eigenvalues, Ev(J), have negative real-parts and it is unstable if at least one eigenvalue has a positive real part [7, 42]. Because the Jacobian of a delay system contains exponentials of its own eigenvalue λ, the characteristic polynomial of the Jacobian is typically a transcendental equation that can admit an infinite number of eigenvalues.

To provide some analytical intuition we now focus on the case of degree homogeneous networks, where d_i = d* for all i. For simplicity of notation we will for the moment pretend that degree homogeneity is satisfied for the entire network. We emphasize however that for the study of localized eigenvalues and eigenvectors, we only need to assume degree homogeneity within the I-nodes, which is often true because of symmetry. As the localized eigenvectors are insensitive to the network structure outside the BSS, they are also insensitive to a violation of the degree homogeneity assumption if the violation occurs outside the BSS. It can be shown that essentially the same is true for the C-nodes.

For degree-homogeneous networks all diagonal elements of the Jacobian are identical such that J^d_i = :J^d. We can then write

$$\begin{equation} {\rm \bf J} = J^{\mathrm{d}}{\rm \bf I} + J^{\mathrm{o}}{\rm \bf A}, \end{equation} \tag{ 9 }$$

where I is the identity matrix. The eigenvalue computation then becomes

$$\begin{equation} {\rm Ev}({\rm \bf J}) = {\rm Ev}(J^{\mathrm{d}}{\rm \bf I} + J^{\mathrm{o}}{\rm \bf A}) = J^{\mathrm{d}} + J^{\mathrm{o}} {\rm Ev}({\rm \bf A}). \end{equation} \tag{ 10 }$$

This equation establishes a relationship between the topological and dynamical eigenvalues of the system. However, we emphasize that both J^d and J^o depend on λ. The relationship (10) thus constitutes a self-consistency condition that typically admits an infinite set of dynamical eigenvalues for each topological eigenvalue. For a given topological eigenvalue c this condition reads

$$\begin{equation} \lambda=J^{\mathrm{d}}(\lambda)+c J^{\mathrm{o}}(\lambda). \end{equation} \tag{ 11 }$$

Instead of trying to compute the dynamical eigenvalues, we focus on the computation of bifurcation points where the stability changes and eigenvalues cross the imaginary axis in the complex plane, so that λ = iω. Separating the real and imaginary part of equation (11) yields

$$\begin{eqnarray} &&0={g'}\cos(\phi)-{l'}-d{k}+|c|{k}\cos(\psi), \end{eqnarray}\noindent \tag{ 12 }$$

$$\begin{eqnarray} &&\omega=-{g'}\sin(\phi)-|c|{k}\sin(\psi), \end{eqnarray} \tag{ 13 }$$

with k = βf', ϕ = ωτ and ψ: = ωδ − ψ^c, where ψ^c is the complex phase of c.

The bifurcation condition is thus given by an equation system with three unknown variables, ϕ, ψ and ω. For every solution triplet (ϕ,ψ,ω), we find another solution (−ϕ,−ψ,−ω), providing the same bifurcation lines. We therefore only consider solutions with ω > 0. Further, we are only interested in solutions for positive delays, τ > 0, and thus only consider solutions with ϕ > 0. We now choose ϕ as a free parameter to find parametric representations of the bifurcation lines. This is done separately for three cases |c| = 0, |c| = d, and 0 < |c| < d.

For c = 0, equations (12) and (13) are independent of ψ and hence δ. Thus, the bifurcations are vertical lines in the (k,δ)-plane, which are solely characterized by

$$\begin{equation} {k}=\frac{1}{\tau}\frac{h}{d}, \end{equation} \tag{ 14 }$$

with h: = (g'cos(ϕ) − l')τ, where ϕ has to satisfy q(ϕ) = 0, with q: = ϕ + g'τsin(ϕ). For the parameter values used throughout this paper, the condition q(ϕ) = 0 has a unique solution.

For |c| = d, we find

$$\begin{equation} {k}=\frac{1}{\tau}\frac{q^2+h^2}{2dh}. \end{equation} \tag{ 15 }$$

In order to obtain positive solutions, h(ϕ) needs to be positive. Therefore, we have to restrict ϕ to the intervals $I^{k}_r=\left (2\pi r-\phi ^{k}, 2\pi (r+1) + \phi ^{k}\right )$, where r is a non-negative integer and ϕ^k = cos⁻¹(l'/g'). The corresponding δ values can be calculated using equations (12) and (13), which yield

$$\begin{eqnarray} &&\delta^{\mathrm{L,R}}=\frac{\psi^{\mathrm{L,R}}+\psi_{\rm c}+2\pi s}{\phi}\tau, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&\psi^{\mathrm{L,R}}=\pm\cos^{-1}\left(1-\frac{h}{\tau |c|{k}}\right) \end{eqnarray} \tag{ 17 }$$

for non-negative integers s. Thus, after calculating k with equation (15), we are able to calculate δ. In order to obtain valid solutions, we need to apply the L-branch for all q(ϕ) < 0 and the R-branch otherwise.

Finally, for 0 < |c| < d, we find

$$\begin{equation} {k}_{\mathrm{1,2}}=\frac{1}{\tau}\frac{dh\pm\sqrt{a}}{d^2-{\left| c \right|}^2}, \end{equation} \tag{ 18 }$$

with

$$\begin{equation} a=d^2h^2-(d^2-|c|^2)(q^2+h^2). \end{equation} \tag{ 19 }$$

In order to obtain real and positive solutions, not only h(ϕ) needs to be larger than 0 but also a(ϕ). This sets additional restrictions to ϕ, which can be calculated numerically. In contrast to the case |c| = d, valid solutions can only be found inside a finite number of intervals I^k_r. The corresponding values for δ are calculated as for |c| = d.

In the latter two cases, two integers r and s enumerate different solution branches. From equation (16), we see that increasing s shifts the bifurcation lines towards larger values of δ. By contrast, the influence of r is less obvious. While different values of r correspond to distinct bifurcation lines in the parameter space, solutions that differ only by s often connect to each other such that they correspond to different segments of the same bifurcation line.

The results from sections 2 and 3 enable us to calculate the bifurcation points of instabilities that are created in BSS in which all C-nodes have an identical degree d. Even in large and complex networks, instabilities that are generated entirely in certain symmetric subgraphs can thus be understood analytically.

The results of the previous sections imply that a given BSS has a characteristic set of bifurcations that will appear in every network that contains the BSS. In the present section we support this conclusion with numerical evidence.

Because the computation of the potentially infinite set of eigenvalues is difficult, we use the approach proposed in [28, 43], which uses Cauchy's argument principle to compute the number of eigenvalues with positive real parts. In the following we apply this approach to analyse networks with a given topology. For each of these networks we consider a large ensemble of (k, δ) parameter pairs that are drawn independently from a uniform distribution. We show this ensemble as a scatter plot in the (k, δ)-plane, in which we colour-code the number of eigenvalues with positive real parts. As a result we obtain plots such as the ones shown in figure 2, where changes in colour indicate bifurcation lines.

Zoom In Zoom Out Reset image size

As a first example, we consider two networks containing the triangular BSS from table 1. For the directed triangle d = 2 as every inner node of the BSS has two incoming and two outgoing links. The bifurcation lines induced by the localized topological eigenvalues, $c=\exp (2\pi \mathrm {i}/3)$ and $c=\exp (4\pi \mathrm {i}/3)$, can thus be calculated using equation (18). The scatter plot for the isolated BSS (figure 2, centre left) shows a number of bifurcation lines, some of which arise from delocalized instabilities. However, two of the bifurcation lines correspond to characteristic instabilities of the triangle. These lines agree exactly with the analytical results.

Let us now consider the more complicated example shown on the right hand side of figure 2, which also contains the directed triangle BSS. For the larger network the corresponding scatter plot is significantly more complex. Note that the delocalized instabilities of the isolated triangle BSS do not carry over to the larger network, while the bifurcation lines corresponding to the localized instabilities reappear exactly.

The effect of the characteristic instabilities on the system can be visualized by simulations (figure 2, bottom). In contrast to the results above, which are valid for the general class of models, simulations require restricting the functional form of the equations of motion. We use coupled Mackey–Glass systems

$$\begin{equation} \dot{X}_i=\frac{aX_i^\tau}{1+(X_i^\tau)^b}-c X_i+\epsilon \sum_j\left(A_{ij}X_j-A_{ji}X_i\right), \end{equation}\noindent \tag{ 20 }$$

with a = 2, b = 5, c = 1,  = 10, τ = 5 corresponding to the generalized parameters g' = −1.5, l' = 1, k = τ = 1. The simulation results show that the two characteristic instabilities correspond to clockwise and counter-clockwise oscillations on the triangle. One can therefore think of the directed triangle BSS as a network motor that will induce a characteristic forward or reverse oscillation if specific conditions are met, independently of the embedding topology. Let us emphasize that these oscillations do not generally remain confined to the triangle subgraph. After the onset of oscillations the propagation of the oscillatory solutions across the network is not described by the same Jacobian and the delicate balance which isolates the instability can be destroyed in the non-stationary state.

Two more complicated examples are shown in figure 3. The figure shows the scatter plots obtained for model realization on two different networks, which both contain the hexagonal BSS with two connecting nodes. The network on the right-hand side moreover contains the star-like BSS and the directed-triangle BSS. In both scatter plots, we find the characteristic bifurcation lines of the hexagonal BSS, corresponding to the topological eigenvalues c = 0, $2\exp (2\pi {\mathrm { i}}/3)$ and $2\exp (4\pi {\mathrm { i}} /3)$. The two-fold degeneracy of the eigenvalue c = 0 is reflected in the number of positive eigenvalues changing by 4 instead of 2 when crossing the corresponding bifurcation line (solid vertical line). In the scatter plot on the right, we further find the bifurcation lines of the triangular BSS (see figure 2) and the bifurcation line originating from the localized eigenvalue c = 0 of the star-like BSS. This confirms that bifurcations induced by different BSS can be superimposed without interfering.

Zoom In Zoom Out Reset image size

In the present paper we showed that symmetric subgraphs induce characteristic localized instabilities. The bifurcations arising from these instabilities occur independently of the structure of the embedding network. Symmetric structures thus provide a wealth of examples for mesoscale structures that have precise and distinct dynamical implications. Because these implications are independent of the embedding structure, they can be analysed in small, structurally simple systems. Therefore, analytical solutions can be obtained even for relatively complicated dynamics such as the delay-coupled delay systems studied here.

The beauty of the results lies perhaps in their simplicity. Essentially a symmetry in topology translates into a symmetry of topological eigenvalues in the complex plane, which then translates into a symmetry of dynamic eigenvalues, and onward to a corresponding dynamical mode. Here we have mainly considered the case of cyclic structures such as the directed triangle, which cause corresponding cyclic modes. It can be easily extrapolated that for instance bipartite structures cause characteristic instabilities that correspond to an imbalance between the partitions. For instance in a directed four-cycle, one would expect four cyclic instabilities and an additional instability where two nodes start to suppress the other two.

Although the connection between mesoscale structure and dynamics may be intuitive, we note that it is only exactly true in the special case of exact symmetry that is not generic in the mathematical sense. In a completely random system, one would not expect that any dynamical instability localizes exactly on a finite set of nodes. It has been pointed out by MacArthur and Sánchez-García [25] that symmetric structures are relatively common in real world networks, with up to 80% of the nodes belonging to symmetric structures in some examples. However, the symmetry of the Jacobian matrix on which we build here, requires additionally that the symmetric nodes also follow the same dynamical rules. In many biological systems this requirement is only met exactly in special situations. Nevertheless, we expect the present results to remain valid to good approximation if the dynamical differences between topologically symmetric nodes are reasonably small (e.g. in a metapopulation where topologically symmetric nodes correspond to patches with similar but not identical environmental conditions). Furthermore, it is easily conceivable that almost exact symmetries exist in many technical systems. Thus, the exploration of both technological and biological implications of the present work, appear to be promising targets for future work, including the analysis of real-world data.

Despite the caveat above, the connection between mesoscale structure and dynamics can be expected to remain valid in a wide variety of different systems. In the present paper we illustrated this connection by the example of delay-coupled delay oscillators. For simplicity we focused on topologies where the number of a node's outgoing links is identical to the number of incoming links. Furthermore, we only considered symmetric subgraphs in which all inner nodes had the same number of links, such that the system was factorizing. While these assumptions were essential for a clean illustration of the results, very similar results could have been obtained in other systems. Firstly, we have focused here on delay-coupled delay oscillators, because they provide non-trivial examples, where already very simple structures can sustain complex dynamics. The same reasoning could have been presented for ordinary differential equations or discrete time maps, although perhaps slightly more complex symmetric subgraphs would have been needed to provide interesting dynamics. Secondly, the balance condition for links allowed us to avoid the computation of steady states which is impossible with the chosen degree of generality, but does not otherwise affect the subsequent reasoning. Thirdly, the factorization enabled us to establish an explicit relationship between the topological and dynamical eigenvalues. We expect that the same will be possible for many symmetric structures and systems. However, even if this is not the case, a symmetry in the structure will generally induce a corresponding symmetry in the Jacobian if the symmetric nodes follow the same dynamical rules. Thus, even when the connection between structure and dynamics cannot be made explicit, a structural symmetry in a subgraph will generally induce characteristic instabilities that are independent of the embedding network. Because these instabilities originate in small sub-systems their numerical investigation should in general be possible with relative ease (see [27] for an example).

One difficulty that may be encountered in future work is that the Jacobian can be sensitive to the location of the steady state. While an embedding network cannot affect the characteristic bifurcation lines in the bifurcation diagrams shown here, it can, by shifting the steady state, alter the parameter values at which the system is located within the bifurcation diagram. Our results would thus have been less clear if we did not use the generalized modelling formalism that conceptually separates the effects of shifting steady state concentrations from the effects of changes in the nonlinearities [39]. In practice this means that subgraphs, such as the directed-triangle, can be constructed to sense some information from the surrounding network if this information is communicated via a change in the steady state values of the variables associated with the symmetric nodes. Conversely, one can also construct networks that are completely insensitive to the surrounding network. This can be done for instance by choosing the functional forms of the processes in the model to follow power laws. For such functions the generalized parameters that appear in the Jacobian are independent of the steady state under consideration (see equation (5)).

In conclusion, the work presented here points to a general connection between the structure of certain mesoscale subgraphs and system-level dynamics. It is conceivable that, in future works, this relationship may prove instrumental for understanding or engineering certain dynamics in complex nonlinear systems. For instance our results seem to suggest that networks can be constructed such that some function that is sustained on a certain part of the network can be isolated from functions that are localized on other parts of the net. Such localization is indeed observed in many biological networks, but significantly more work will be necessary to determine whether the symmetries described here contribute to this localization. Furthermore, attempts could be made to engineer simple information processing systems based on communicating symmetric subgraphs that are embedded in a larger network. While much work remains to be done to explore these possibilities we believe that, if successful, these efforts might lead to design principles for information processing in applications such as synthetic biology or chemical computing.