Ginzburg-Landau approximation for self-sustained oscillators weakly coupled on complex directed graphs

TL;DR

This paper derives a Ginzburg-Landau approximation for reaction-diffusion systems on directed graphs near a Hopf bifurcation, enabling analysis of stability and pattern formation in coupled oscillators.

Contribution

It introduces a complex Ginzburg-Landau equation tailored for directed networks, linking network topology with oscillatory stability analysis.

Findings

The CGLE accurately predicts stability and symmetry-breaking instabilities.

Numerical simulations with the Brusselator model validate the theoretical predictions.

Pattern formations from the CGLE match those from the original reaction-diffusion system.

Abstract

A normal form approximation for the evolution of a reaction-diffusion system hosted on a directed graph is derived, in the vicinity of a supercritical Hopf bifurcation. Weak diffusive couplings are assumed to hold between adjacent nodes. Under this working assumption, a Complex Ginzburg-Landau equation (CGLE) is obtained, whose coefficients depend on the parameters of the model and the topological characteristics of the underlying network. The CGLE enables one to probe the stability of the synchronous oscillating solution, as displayed by the reaction-diffusion system above Hopf bifurcation. More specifically, conditions can be worked out for the onset of the symmetry breaking instability that eventually destroys the uniform oscillatory state. Numerical tests performed for the Brusselator model confirm the validity of the proposed theoretical scheme. Patterns recorded for the CGLE resemble closely those recovered upon integration of the original Brussellator dynamics.

Equations

$${\dot{\mathbf{x}}}_j=\mathbf{F}({\mathbf{x}}_j,\mathbf{\mu })+\mathbf{D}\sum _{k=1}^N{\Delta }_{jk}{\mathbf{x}}_k$$

$$\frac{\partial }{dt}{\mathbf{u}}_j={\mathbf{Lu}}_j+\mathbf{D}\sum _{k=1}^N{\Delta }_{jk}{\mathbf{u}}_k+\mathcal{M}{\mathbf{u}}_j{\mathbf{u}}_j+\mathcal{N}{\mathbf{u}}_j{\mathbf{u}}_j{\mathbf{u}}_j+\dots$$

$${\left(\mathcal{M}{\mathbf{u}}_j{\mathbf{u}}_j\right)}_l$$

$${\left(\mathcal{N}{\mathbf{u}}_j{\mathbf{u}}_j{\mathbf{u}}_l\right)}_l$$

$$\frac{d}{dt}⟶\frac{\partial }{\partial t}+{ϵ}^2\frac{\partial }{\partial \tau }.$$

$$\mathbf{L}$$

$$\mathcal{M}$$

$$\mathcal{N}$$

$${\mathbf{u}}_j(t)=ϵ{\mathbf{u}}_j^{(1)}(t,\tau )+{ϵ}^2{\mathbf{u}}_j^{(2)}(t,\tau )+\dots$$

$$\begin{gathered} \left(\frac{\partial }{\partial t}{\mathbb{I}}_2+{ϵ}^2\frac{\partial }{\partial \tau }{\mathbb{I}}_2-{\mathbf{L}}_0-{ϵ}^2{\mathbf{L}}_1-\dots \right)(ϵ{\mathbf{u}}_j^{(1)}+{ϵ}^2{\mathbf{u}}_j^{(2)}+\dots ) \\ -{ϵ}^2\mathbf{D}\sum _{k=1}^N{\Delta }_{jk}(ϵ{\mathbf{u}}_k^{(1)}+{ϵ}^2{\mathbf{u}}_k^{(2)}+\dots )={ϵ}^2{\mathcal{M}}_0{\mathbf{u}}_j^{(1)}{\mathbf{u}}_j^{(1)}+{ϵ}^3(2{\mathcal{M}}_0{\mathbf{u}}_j^{(1)}{\mathbf{u}}_j^{(2)} \\ +{\mathcal{N}}_0{\mathbf{u}}_j^{(1)}{\mathbf{u}}_j^{(1)}{\mathbf{u}}_j^{(1)})+\mathcal{O}({ϵ}^4) \end{gathered}$$

$$\left(\frac{\partial }{\partial t}{\mathbb{I}}_2-{\mathbf{L}}_0\right){\mathbf{u}}_j^{(\nu )}={\mathbf{B}}_j^{(\nu )}$$

$${\mathbf{u}}_j^{(1)}(t,\tau )=W_j(\tau ){\mathbf{U}}_0{e}^{i{\omega }_{0}t}+c.c.$$

$${\mathbf{u}}_j^{(2)}=W_j^2{\mathbf{V}}_2{e}^{2i{\omega }_{0}t}+\gamma {\mathbf{u}}_j^{(1)}+|W_j{|}^2{\mathbf{V}}_0.$$

$$\frac{d}{d\tau }{W}_j(\tau )=\sigma W_j-g|W_j{|}^2{W}_j+d\sum _{k=1}^N{\Delta }_{jk}{W}_k$$

$$\frac{\tau }{{\sigma }_{Re}}\to \tau \frac{d_{Re}}{{\sigma }_{Re}}{\Delta }_{ij}\to {\Delta }_{ij}\sqrt{\frac{g_{Re}}{{\sigma }_{Re}}}{W}_j\to W_j,$$

$$c_0=\frac{{\sigma }_{Im}}{{\sigma }_{Re}}{c}_1=\frac{d_{Im}}{d_{Re}}{c}_2=\frac{g_{Im}}{g_{Re}}$$

$$\frac{d}{d\tau }{W}_j=W_j-(1+i{c}_2)|W_j{|}^2{W}_j+(1+i{c}_1)\sum _k{\Delta }_{jk}{W}_k$$

$$W_j(\tau )\equiv W^{LC}=e^{-i{c}_2\tau }$$

$$W_j(\tau )=W^{LC}(\tau )(1+{\rho }_j(\tau ))e^{i{\theta }_j(\tau )}$$

$$\left(\begin{array}{c}{\dot{\rho }}_j\\ {\dot{\theta }}_j\end{array}\right)=\left(\begin{array}{cc}-2& 0\\ -2c_2& 0\end{array}\right)\left(\begin{array}{c}{\rho }_j\\ {\theta }_j\end{array}\right)+\sum _{k=1}^N{\Delta }_{jk}\left(\begin{array}{cc}1& -c_1\\ c_1& 1\end{array}\right)\left(\begin{array}{c}{\rho }_k\\ {\theta }_k\end{array}\right).$$

$$\left(\begin{array}{c}{\rho }_j\\ {\theta }_j\end{array}\right)=\sum _{\alpha =1}^N\left(\begin{array}{c}{\rho }^{(\alpha )}\\ {\theta }^{(\alpha )}\end{array}\right)e^{{\lambda }^{(\alpha )}\tau }{\Phi }_j^{(\alpha )}$$

$$\text{det}\left({\mathbf{J}}_{\alpha }-{\lambda }^{(\alpha )}{\mathbb{I}}_2\right)=0$$

$${\mathbf{J}}_{\alpha }=\left(\begin{array}{cc}-2+{\Lambda }^{(\alpha )}& -c_1{\Lambda }^{(\alpha )}\\ -2c_2+c_1{\Lambda }^{(\alpha )}& {\Lambda }^{(\alpha )}\end{array}\right).$$

$$S_2({\Lambda }_{Re}^{(\alpha )})⩽S_1({\Lambda }_{Re}^{(\alpha )}){\left[{\Lambda }_{Im}^{(\alpha )}\right]}^2$$

$$S_2({\Lambda }_{Re}^{\alpha })$$

$$-C_{2,1}{\Lambda }_{Re}^{(\alpha )}$$

$$S_1({\Lambda }_{Re}^{\alpha })$$

$$C_{2,4}$$

$$C_{2,3}$$

$$C_{2,2}$$

$$C_{2,1}$$

$$C_{1,2}$$

$$C_{1,1}$$

$$C_{1,0}$$

$$\{\begin{array}{rl}\dot{{\varphi }_j}& =A-(B+1){\varphi }_j+{\varphi }_j^2{\psi }_j+D_{\varphi }\sum _{k=1}^N{\Delta }_{jk}{\varphi }_k\\ \dot{{\psi }_j}& =B{\varphi }_j-{\varphi }_j^2{\psi }_j+D_{\psi }\sum _{k=1}^N{\Delta }_{jk}{\psi }_k\end{array}$$

$$c_1$$

Full text

††journal: CNSNS

1 Introduction

Many real-life phenomena can be ultimately described in terms of mutually interacting entities which can occasionally give rise to collective behaviors [1, 2]. The emerging patterns may play important functional roles, as it is for instance the case for biochemical processes [3] and ecological applications [4, 5, 6]. Synchrony is among the most striking example of self-organized dynamics [7]. It is encountered in a wide gallery of natural systems, think to the beating of the heart [8] and the firing of the firefly [9]. It also plays a role of paramount importance for the correct functioning of man-designed technology, as e.g. in power grids [10, 11].

Synchronization of self-sustained oscillations is a particularly rich field of investigation. Mathematically, self-sustained oscillations correspond to stable limit cycles in the state space of an autonomous continuous-time dynamical system. Oscillators can be embedded in continuum space, being subject to diffusive couplings. The ensuing reaction-diffusion system displays synchronous oscillations, which, under specific conditions, prove robust to external perturbations. Each oscillator can alternatively occupy a node of a complex network [12, 13, 14, 15], a generalization that opens up the perspective to tackle a a large plethora of problems that deal with a discrete and heterogeneous hosting support [15, 16, 17].

For continuous reaction-diffusion systems near the supercritical Hopf bifurcation, one can obtain a simplified, normal form description of the dynamics in terms of an associated Ginzburg-Landau equation (CGLE) [18, 19, 7]. This is a reduced picture which provides an accurate representation of the original model, while allowing for analytical progress to be made. From a more general perspective, it is also important to classify minimal, though effective descriptive frameworks that could guide in the search for universal traits that happen to be shared by multispecies models, besides their intrinsic degree of inherent specificity. Simplified schemes that exemplify an exact underlying dynamics can be effectively employed to shed light on the onset of spatio-temporal chaos and propagation of nonlinear waves. Externally imposed non homogeneous perturbation can, for instance, break the synchrony of the oscillations as displayed by the original reaction-diffusion system, or equivalently its CGLE analogue, so materializing in a colorful density patterns that sustain spatio-temporal propagation. The formal link between continuous reaction diffusion-systems and the CGLE was established in the pioneering work by Kuramoto [18], exploiting a multiple-scale perturbative analysis. More recently the analysis has been extended by Nakao [20] to the relevant setting where the reaction-diffusion system is made to evolve on a symmetric, hence undirected, network. As remarked in [21], directionality matters and can seed the emergence of non trivial collective dynamics which cannot manifest when the scrutinized system is made to evolve on a symmetric discrete support. Motivated by this finding, we considered in [22] the dynamics of a reaction-diffusion system defined on a directed graph which displays a stable fixed point and obtained an effective description for the evolution mode triggered unstable, just above the threshold of criticality. The analysis exploits a multiple time-scale analysis and eventually yields a Stuart-Landau for the amplitude of the unstable mode, whose complex coefficients reflect the topology of the network, the factual drive to the instability.

In this paper we aim at following the similar strategy of [22] to derive an approximate equation for the evolution of a reaction diffusion system on a directed (and balanced) graph, in the vicinity of a supercritical Hopf bifurcation. As a matter of fact we will assume weak the strength of the coupling that links adjacent nodes. In doing so, we will generalize the work of Nakao [20] to the interesting setting where asymmetry in the couplings needs to be accommodated for and, at the same time, reformulate the classical work of Kuramoto [18], on a discrete spatial backing. To anticipate our findings, and at variance with the analysis reported in [22], we will finally obtain a CGLE as a minimal description for the dynamics of the self-sustained oscillators coupled on a complex and asymmetric graph. The obtained CGLE enables one to analytically probe the stability of the synchronous uniform state, as displayed by the reaction-diffusion system. Specifically, it allows to determine the parameters setting that instigates a symmetry breaking transition to non-uniform patterns. Numerical tests made for the Brusselator model, here assumed as a reference model for its pedagogical interests, will confirm the adequacy of the proposed approximate scheme.

2 Diffusive oscillators on networks

We will here consider a generic two dimensional reaction-diffusion system and label with ${\bf x}_{j}(t)=(\phi_{j},\psi_{j})^{T}$ for $j=1,\ldots,N$ the two-dimensional real vector of the concentrations. The index $j$ refers to the node of the network to which the selected components refer to. $N$ stands for the size of the network, i.e. the total number of nodes. The only further assumption that we shall make is the existence of self-sustained oscillations for the system under scrutiny and for this reason we will point to ${\bf x}_{j}(t)$ as to the oscillators’ variables. The dynamics of the system is hence described by the following differential equation

[TABLE]

where the two-dimensional nonlinear function ${\bf F}$ specifies the reaction terms: it depends on the local concentrations of the species ${\bf x}_{j}$ and on $\boldsymbol{\mu}$ which is a vector of arbitrary dimension that gather together the parameters of the model. The second term represents the diffusive coupling: ${\bf D}=K\text{diag}(D_{\phi},D_{\psi})$ denotes the diagonal matrix of the diffusion coefficients. $K$ is a constant parameter that set the strength of the coupling. As we will make clear in the following the perturbative analysis that we shall develop, hold for $K<<1$. In equation (1), ${\bf\Delta}$ is the Laplacian matrix whose elements read $\Delta_{ij}=A_{ij}-\delta_{ij}k_{i}$. Here we focus on directed networks, thus the adjacency matrix ${\bf A}$ is not symmetric. Using standard notations, $A_{ij}=1$ if a link exists that goes from node $i$ to node $j$. Otherwise, $A_{ij}=0$. $k_{i}$ is the number of outgoing edges from node $i$, $\delta_{ij}$ is the Kronecker’s delta. We assume that system (1) admits a homogeneous equilibrium point that we here label ${\bf x}^{*}=(\phi^{*},\psi^{*})^{T}$. This request implies dealing with a balanced network, namely a network where the outgoing and incoming connectivities are equal.

We additionally require that ${\bf x}^{*}$ undergoes a Hopf bifurcation for $\boldsymbol{\mu}=\boldsymbol{\mu}_{0}$. Accordingly, the Jacobian matrix associated to system (1) has a pair of imaginary eigenvalues $\pm i\omega_{0}$. Slightly above the supercritical Hopf bifurcation, ${\bf x}^{*}$ becomes unstable, the reaction-diffusion system admits an time dependent homogeneous solution. This is the uniform state obtained by replicating on each node of the network and in complete synchrony, the limit cycle displayed by the system in its a-spatial limit ($K=0$). The spatial coupling, sensitive to tiny non homogeneities, which configure as injected perturbation, can eventually destabilized the uniform synchronous equilibrium. When diffusion is small ($K=\epsilon^{2}<<1$), the method of multiple timescales [18] constitutes a viable strategy to characterize the nonlinear evolution of the perturbation and hence elaborate on the stability of the time-dependent uniform periodic solution.

To this end, inspired by the analysis carried out in [20], we introduce small inhomogeneous perturbations, $\delta\phi_{j}$ and $\delta\psi_{j}$, to the uniform equilibrium point, namely $(\phi_{j},\psi_{j})=(\phi^{*},\psi^{*})+(\delta\phi_{j},\delta\psi_{j})$ for $j=1,\ldots,N$. We then substitute this ansatz into equations (1). Performing a Taylor expansion of the resulting system and packing $\delta\phi_{j}$ and $\delta\psi_{j}$ into the column vector ${\bf u}_{j}=\left(\delta\phi_{j},\delta\psi_{j}\right)^{T}$, one ends up with the following equation for the time evolution of ${\bf u}_{j}$:

[TABLE]

where ${\bf L}$ is the Jacobian matrix evaluated at the steady state ${\bf x}^{*}$. Adopting the same symbolic notations of [18], $\mathcal{M}{\bf u}_{j}{\bf u}_{j}$ and $\mathcal{N}{\bf u}_{j}{\bf u}_{j}{\bf u}_{j}$ denote two-dimensional vectors whose components respectively read

[TABLE]

for $l=1,2$.

Setting the model above the supercritical Hopf bifurcation, translates into redefining $\boldsymbol{\mu}$ by means of a small parameter $\epsilon$ as $\boldsymbol{\mu}=\boldsymbol{\mu}_{0}+\epsilon^{2}\boldsymbol{\mu}_{1}$, where $\boldsymbol{\mu}_{1}$ is order one. To follow the nonlinear evolution of the amplitude of the most unstable mode, it is appropriate to introduce a slow time variable $\tau=\epsilon^{2}t$ and, accordingly, to rewrite the total derivative with respect to the original time $t$:

[TABLE]

As a key assumption, already alluded to in the preceding discussion, we set $K=\epsilon^{2}$: in other words we suppose that the impact of diffusion is of the same order as the deviation from the bifurcation point. Moreover, near criticality, the matrix ${\bf L}$ and the operators $\mathcal{M}$ and $\mathcal{N}$ may be expanded in powers of $\epsilon^{2}$

[TABLE]

We further assume that ${\bf u}$ can be expressed as a series, function of both $t$ and $\tau$:

[TABLE]

Inserting Eqs. (3), (4) and (5) into Eq. (2) we get:

[TABLE]

Collecting together terms of the same order in $\epsilon$ returns the following set of equations

[TABLE]

for $\nu=1,2,3...$. The existence of a nontrivial solution for the aforementioned linear systems is guaranteed by the Fredholm theorem [23]. As we shall discuss in a moment, the solvability condition (see A) is directly satisfied for $\nu=1$ and $\nu=2$, while it must be explicitly imposed for $\nu=3$.

Let us start by focusing on $\nu=1$. The corresponding right hand side of system (6) is ${\bf B}_{j}^{(1)}=0$. It is easy to see that the analytical solution is

[TABLE]

where ${\bf U}_{0}$ the right eigenvector of ${\bf L}$ corresponding to the eigenvalue $i\omega_{0}$. The complex variable $W_{j}(\tau)$ is the amplitude of the perturbation. Its dynamics is supposed to be slow, i.e. ruled by the rescaled time $\tau$. As we shall see, by imposing the solvability condition for $\nu=3$ will return an equation for the time evolution of $W_{j}(\tau)$.

Moving to $\nu=2$, we find ${\bf B}_{j}^{(2)}=\mathcal{M}_{0}{\bf u}_{j}^{(1)}{\bf u}_{j}^{(1)}$ and we try to cast the solution of the linear system in the form

[TABLE]

Inserting this ansatz into (6) and grouping together terms that do not depend on $t$, one finds ${\bf V}_{0}=-2{\bf L}_{0}^{-1}\mathcal{M}_{0}{\bf U}_{0}\bar{{\bf U}}_{0}$ where the bar stands for the complex conjugate. In the same way, collecting terms proportional to $e^{2i\omega_{0}t}$ yields ${\bf V}_{2}=\left(2i\omega_{0}\mathbb{I}_{2}-L_{0}\right)^{-1}\mathcal{M}_{0}{\bf U}_{0}{\bf U}_{0}$.

The third constant term appearing in Eq. (6) is explicitly given by ${\bf B}_{j}^{(3)}=\left({\bf L}_{1}-\frac{d}{d\tau}\mathbb{I}_{2}\right){\bf u}_{j}^{(1)}$ $+{\bf D}\sum_{k=1}^{N}\Delta_{jk}{\bf u}_{k}^{(1)}+2\mathcal{M}_{0}{\bf u}_{j}^{(1)}{\bf u}_{j}^{(2)}+\mathcal{N}_{0}{\bf u}_{j}^{(1)}{\bf u}_{j}^{(1)}{\bf u}_{j}^{(1)}$. By imposing the solvability condition, one eventually obtains the following CGL equation for $W_{j}(\tau)$

[TABLE]

where $\sigma=\sigma_{Re}+i\sigma_{Im}=({\bf U}_{0}^{*})^{\dagger}{\bf L}_{1}{\bf U}_{0}$, $g=g_{Re}+ig_{Im}=$ $-({\bf U}_{0}^{*})^{\dagger}\left[2\mathcal{M}_{0}{\bf V}_{2}\bar{\bf U}_{0}\right.$ $\left.+2\mathcal{M}_{0}{\bf V}_{0}{\bf U}_{0}+3\mathcal{N}_{0}{\bf U}_{0}{\bf U}_{0}\bar{\bf U}_{0}\right]$ and $d=d_{Re}+id_{Im}=({\bf U}_{0}^{*})^{\dagger}{\bf D}{\bf U}_{0}$ are complex numbers.

Using the following scale transformation

[TABLE]

introducing three new coefficients

[TABLE]

and making the further transformation $W_{j}\rightarrow W_{j}\exp(ic_{o}\tau)$, we find the following reduced CGLE

[TABLE]

for the complex amplitude $W_{j}$ of the $j$-th oscillator.

Equation (9) displays a uniformly synchronized solution given by

[TABLE]

for $j=1,\ldots,N$. In this respect it provides a local approximation for the dynamics of the original reaction-diffusion system near the Hopf bifurcation, namely when self-organized oscillations rise. The dynamics of the ensemble made of $N$ oscillators is sensitive to the coupling imposed, this latter being encoded in the Laplacian operator. Mutual interferences among oscillators, as mediated by the weak diffusive coupling here at play, can disrupt the coherency of the synchronous state, making the system to evolve towards a different attractor. The stability of the synchronized dynamics can be assessed via a linear stability analysis which is carried out for CGLE, and whose conclusion can be readily translated to the original reaction-diffusion framework. Having obtained a universal description of the dynamics of the system in terms of a normal mode representation allows one to obtain general conditions of the onset of the instability, which prescind the specific context of analysis.

To investigate the linear stability of solution (10), we briefly review the method described in [24, 20] and point to the specific aspects that relate to the directed nature of the couplings, as discussed in [25]. Notice that in this latter paper, the CGLE is assumed as the reference model, while here it is obtained as a local approximation of a generic reaction-diffusion system defined on a directed network, under weak diffusive couplings.

To carry out the linear stability calculation, we substitute the following amplitude and phase perturbation [24, 20]

[TABLE]

into Eq. (9) and we linearize around $\rho_{j}=0$ and $\theta_{j}=0$ obtaining

[TABLE]

Using the Laplacian eigenvalues and eigenvectors that, by definition, satisfy ${\bf\Delta}{\bf\Phi}^{(\alpha)}$ $=\Lambda^{(\alpha)}{\bf\Phi}^{(\alpha)}$ for $\alpha=1,\dots,N$, the inhomogeneous perturbations $\rho_{j}$ and $\theta_{j}$ can be expanded as

[TABLE]

where $\rho^{(\alpha)}$ and $\theta^{(\alpha)}$ denote the expansion coefficients, while $\lambda^{(\alpha)}$ controls the exponential growth of the $\alpha$-th mode. If the real part of $\lambda^{(\alpha)}$ ($\lambda^{(\alpha)}_{Re}$) is negative, perturbations shrink to zero, otherwise they grow exponentially. Plugging expansion (12) into equation (11) yields the following condition

[TABLE]

with

[TABLE]

Constrain (13) returns an equation for $\lambda^{(\alpha)}$ as a function of $\Lambda^{(\alpha)}$ known as the dispersion relation. For a symmetric undirected network, as the one explored by Nakao [24, 20], the eigenvalues $\Lambda^{(\alpha)}$ of the Laplacian operators are real and this observation translates in simple condition for the instability to develop. A straightforward calculation, allows one to conclude that the perturbation can grow, and so destroy the synchrony of the system, provided $1+c_{1}c_{2}<0$. When the underlying graph of coupling is instead directed, the spectrum of the Laplacian operator is complex and the instability can set in also when $1+c_{1}c_{2}>0$. For this reason, the generalized class of instability that is encountered when the system is confined on a asymmetric support is referred to as to topology driven. The conditions which result in $\lambda^{(\alpha)}_{Re}>0$ when the network of couplings is directed, or, equivalently, when $\Lambda^{(\alpha)}$ is complex ($\Lambda^{(\alpha)}=\Lambda_{Re}^{(\alpha)}+i\Lambda_{Im}^{(\alpha)}$), have been worked out in [25], building on the recipe outlined in [21]. The conclusion of [25] are here briefly reviewed for for the sake of completeness. Indeed, $\lambda_{Re}^{(\alpha)}>0$ if

[TABLE]

where

[TABLE]

with

[TABLE]

When condition (14) is met for some eigenvalue $\Lambda^{(\alpha)}$, the limit cycle looses its stability and the system evolves towards a different attractors, which is non homogeneous in space.

The next section is aimed at challenging the validity of Eq. (9), as an approximate equation for the non linear evolution of an injected non homogeneous perturbation. To anticipate our findings, we will show that the CGLE allows to correctly delineate the region of parameters that yields stable and synchronized self-sustained oscillators, and so accurately mark the transition to the pattern forming domain. In both cases, i.e. when uniform oscillations or complex patterns in density are displayed, we obtain a satisfying degree of correspondence between the original reaction-diffusion model and its corresponding CGLE. To perform the numerical tests we will make use of the Brusselator model, as a reference case study.

3 Numerical validation of the theory: the Brusselator model

To challenge the analysis performed in the previous section, we here consider a specific reaction-diffusion system, the Brusselator model, which allows for self-sustained oscillations. The governing equations read:

[TABLE]

where $A$, $B$ stand for non negative constants parameters. This system admits a homogeneous steady state given by $(\phi^{*},\psi^{*})=(A,B/A)$. By acting on the parameter $B$, while keeping $A$ fixed, one can induce a Hopf bifurcation, which opens up the route to self-sustained oscillations. The bifurcation occurs for $B_{c}=1+A^{2}$, as explained in the annexed B. Above the bifurcation point, the Brusselator dynamics can be approximated by a CGLE. Following the derivation reported in the B, which builds on the seminal paper by Kuramoto [18], one ends up with the compact expression for $c_{1}$ and $c_{2}$:

[TABLE]

The CGLE is therefore completely characterized and it can be readily employed to asses the stability of the uniform oscillatory state, namely the spatially extended configuration that is obtained as the synchronous replica of $N$ self-sustained oscillators. Fix the values of $D_{\phi}$ and $D_{\psi}$. Then, the stability of the synchronous state is ultimately controlled by $A$, the sole parameter that one can freely tune. We will denote with $A_{c}$ the value of $A$ (if it exists), above which stability is eventually lost. Moreover, we will set the parameters $c_{1}$ and $c_{2}$ so as to enforce stability on a symmetric support, namely $1+c_{1}c_{2}<0$. The instability that materializes for $A>A_{c}$ is hence indirectly reflecting the asymmetry of the imposed couplings.

The results of the analysis is depicted in Fig. 1(a) and (b). The shaded grey area in the plane ($\Lambda_{Re}^{(\alpha)},\Lambda_{Im}^{(\alpha)}$) delimits the region of instability, as it follows inequality (14). Symbols represent the eigenvalues of the Laplacian matrix in the complex plane and follow the specific choice of the network operated. As explained in the caption of Fig. 1, we here work with a balanced Newman-Watts graph, generated by following the recipe discussed in [22]. For the selected value of $A$ (larger than the critical value $A_{c}=1.48$), two eigenvalues of the Laplacian matrix protrude inside the shaded region, thus triggering the instability. This is indeed a topology driven instability, as it cannot realize if the same choice of parameters is operated for an homologous system hosted on a symmetric support. To argue along this line, consider a symmetric network of couplings: the spectrum of the Laplacian matrix is now real ($\Lambda_{Im}^{(\alpha)}=0$) and should hence fall on the real axis of Fig. 1(a). The grey region where the instability materializes, intersects the real axis only in the origin, where the trivial zero eigenvalue (corresponding to the uniform, fully synchronized, limit cycle state) sits. For the chosen values of $A$ and $B$, hence $c_{1}$ and $c_{2}$, it is therefore necessary to activate an imaginary component of $\Lambda_{Im}^{(\alpha)}$, to invade the region of the complex plane deputed to the instability. To state it differently, for the same choice of parameters, the instability sets in only if a suitable degree of asymmetry is enforced in the networks of connections, while it is bound to fade away otherwise. This scenario is confirmed by Fig. 1(b) where the (magenta online) triangles represent the real part of the dispersion relation, $\lambda_{Re}$, as a function of $-\Lambda_{Re}^{(\alpha)}$. Also in this case, the dispersion relation is punctually positive in correspondence of two specific entries $-\Lambda_{Re}^{(\alpha)}$. The solid line curve represents the homologous dispersion relation when the diffusive coupling is instead declined on a symmetric (and continuous) support.

Patterns emerging from the aforementioned instability are displayed in Fig. 1. Panel (c) shows patterns relative to species $\phi$ obtained via a numerical integration of system (15), while panel (d) originates from the CGLE (8). The degree of correspondence is satisfactory and testifies a posteriori on the adequacy of the proposed approximation.

Figure 2 refers to $A<A_{c}$: the synchronized limit cycle solution is predicted to be stable, following the CGL approximation. Direct integration of the original reaction diffusion system confirms this conclusion. The degree of correspondence between original and approximated schemes is, also in this case, satisfying.

4 Conclusions

Synchronization is a fascinating field of investigation that has applications ranging from biology to computer science, just to name few examples. When this phenomenon occurs on complex networks, the stability of the synchronized states is strongly affected by the topology of the graph. When the network is directed, external non homogeneous perturbations can destabilize synchronous oscillations also if the parameters are set to values that fall in the region of stability for the dynamics formulated on a symmetric support.

Starting from these premises, we here considered a reaction diffusion system defined on a directed (and balanced) graph. The system is made to evolve in the vicinity of a supercritical Hopf bifurcation and weak diffusive couplings are assumed to hold between adjacent nodes. By generalizing the work of Nakao [20] to the setting where asymmetry in the couplings is enforced, and adapting the multiple-scales approach discussed in Kuramoto [18] to the case of a discrete, network like support, we derived a normal form equation to approximate the local evolution of the underlying reaction-diffusion system. This is a Complex Ginzburg-Landau equation (CGLE) whose coefficients clearly depend on the parameters of the model, but also on the topological characteristics of the hosting network. The CGLE can be effectively employed to probe the stability of the periodic uniform solution and work out the conditions for the onset of the symmetry breaking instability that eventually destroy the synchronized regime. Numerical tests performed for the Brusselator model, confirm the adequacy of the proposed approximation. Patterns obtained under the CGLE scheme resemble closely those displayed upon integration of the original reaction-diffusion system.

5 Acknowledgments

This work has been supported by the program PRIN 2012 founded by the Italian Ministero dell’Istruzione, della Università e della Ricerca (MIUR). D.F. acknowledges financial support from H2020-MSCA-ITN-2015 project COSMOS 642563. The work of T.C. presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office.

Appendix A The solvability condition

We consider a linear operator $\bf A$, and two complex vectors ${\bf u}(t)$ and ${\bf b}(t)$ of the same length. According to the Fredholm theorem, the linear system ${\bf Au}(t)={\bf b}(t)$ is solvable if $\langle{\bf v}(t),{\bf b}(t)\rangle=0$ for all vectors ${\bf v}(t)$ solution of ${\bf A}^{*}{\bf v}(t)=0$, where ${\bf A}^{*}$ is the adjoint operator satisfying $\langle{\bf A}^{*}\bf y,x\rangle=\langle y,Ax\rangle$ $\forall x,y$. The angular brackets denote the scalar product that we here define as $\langle{\bf v}(t),{\bf b}(t)\rangle$ $=\int_{0}^{2\pi/\omega_{0}}{\bf v}^{\dagger}(t){\bf b}(t)dt$, the symbol $\dagger$ standing for the conjugate transpose.

With reference to equation (6), the first requirement of the Fredholm theorem consists in finding ${\bf v}(t)$ such that $(\partial/\partial t\mathbb{I}_{2}-{\bf L}_{0})^{*}{\bf v}(t)=0$. By partial integration, and recalling that ${\bf L}_{0}$ is a real matrix, we find that $(\partial/\partial t\mathbb{I}_{2}-{\bf L}_{0})^{*}=-(\partial/\partial t\mathbb{I}_{2N}+{\bf L}_{0})^{T}$. As a consequence, the system to be solved is $-(\partial/\partial t\mathbb{I}_{2}+{\bf L}_{0})^{T}{\bf v}(t)=0$. In analogy with equation (7), we search ${\bf v}(t)$ in the form ${\bf v}(t)={\bf U}_{0}^{*}e^{i\omega_{0}t}$ for some vector ${\bf U}_{0}^{*}$. Substituting this ansatz into the previous equation, we find ${\bf L}_{0}^{T}{\bf U}_{0}^{*}=-i\omega_{0}{\bf U}_{0}^{*}$. To simplify calculation, it is convenient to normalize ${\bf U}_{0}^{*}$ so that $({\bf U}_{0}^{*})^{\dagger}{\bf U}_{0}=1$.

Having defined ${\bf U}_{0}^{*}$, one can explicitly write down the solvability condition $\langle{\bf U}_{0}^{*}e^{i\omega_{0}t},{\bf B}_{j}^{(\nu)}(t,\tau)\rangle=0$. Since ${\bf B}_{j}^{(\nu)}(t,\tau)$ turns out to be periodic functions of period $2\pi/\omega_{0}$, it is appropriate to express them in the form ${\bf B}_{j}^{(\nu)}(t,\tau)=\sum_{l=-\infty}^{+\infty}\left({\bf B}_{j}^{(\nu)}(\tau)\right)_{l}e^{il\omega_{0}t}$. If we multiply this series by $({\bf U}_{0}^{*}e^{i\omega_{0}t})^{\dagger}$ we again obtain periodic functions that, when integrated over the period $2\pi/\omega_{0}$ give zero. The only exception holds for $l=1$ which gives $\langle{\bf U}_{0}^{*}e^{i\omega_{0}t},\left({\bf B}_{j}^{(\nu)}(\tau)\right)_{1}e^{i\omega_{0}t}\rangle=\int_{0}^{2\pi/\omega_{0}}({\bf U}_{0}^{*})^{\dagger}\left({\bf B}_{j}^{(\nu)}(\tau)\right)_{1}dt$. The integrand does not depend on time $t$ and therefore the integral is zero only if the integrand itself is identically equal to zero. For this reason the solvability condition reduces to $({\bf U}_{0}^{*})^{\dagger}\left({\bf B}_{j}^{(\nu)}(\tau)\right)_{1}=0$ $\forall j$.

Appendix B Details of the coefficients of the GL for the Brusselator model

In this section we explicitly derive the coefficients of the CGL equation for the Brusselator model.

The Jacobian matrix associate to system (15) evaluated at the equilibrium point reads

[TABLE]

Setting the system at the bifurcation point means requiring ${\bf L}_{0}$ to have a pair of imaginary eigenvalues, which translates into $\text{tr}{\bf L}_{0}=0$. This relation implies a constrain on the critical value of the parameter $B$ which must be chosen as $B_{c}=1+A^{2}$. With this assumption, ${\bf L}_{0}$ can be rewritten as

[TABLE]

and its pair of imaginary eigenvalues are $\lambda_{0}=\pm i\omega_{0}=\pm iA$. The right eigenvectors of $L_{0}$ and $L_{0}^{T}$ corresponding to, respectively, eigenvalues $iA$ and $-iA$ are given by

[TABLE]

To proceed with the multiscale analysis, we perturb $B_{c}$ as $B=B_{c}+\epsilon^{2}B_{c}$ and we immediately find the entries of matrix ${\bf L}_{1}$

[TABLE]

and the value of $\sigma$ which turns out to be $\dfrac{(1+A^{2})}{2}$.

The operators $\mathcal{M}_{0}$ and $\mathcal{N}_{0}$ are given by

[TABLE]

where $(\cdot)_{i}$ defines the $i$-th component of a two-dimensional vector. These preliminary quantities allow us to calculate

[TABLE]

and

[TABLE]

From Eq. (17) we can make explicit ${\bf V}_{0}$ and ${\bf V}_{2}$

[TABLE]

and thus

[TABLE]

Notice that our equations for ${\bf U}_{0}$ and ${\bf V}_{2}$ are slightly different from those reported in [18]. Despite this minor difference, we end up with same values of $d$ and $g$. Finally, we can write down the two last coefficients of the CGL equation:

[TABLE]