The theory of Turing patterns on time varying networks

TL;DR

This paper investigates how Turing-like pattern formation occurs on time-varying networks, revealing a novel instability mechanism influenced by network dynamics and providing analytical predictions for pattern onset.

Contribution

It introduces a new analytical framework for understanding Turing pattern formation on dynamic networks, extending classical models to time-dependent structures.

Findings

Network dynamics can induce Turing-like instabilities.

Analytical predictions match numerical simulations.

Pattern onset depends on network evolution frequency.

Abstract

The process of pattern formation for a multi-species model anchored on a time varying network is studied. A non homogeneous perturbation superposed to an homogeneous stable fixed point can amplify, as follows a novel mechanism of instability, reminiscent of the Turing type, instigated by the network dynamics. By properly tuning the frequency of the imposed network evolution, one can make the examined system behave as its averaged counterpart, over a finite time window. This is the key observation to derive a closed analytical prediction for the onset of the instability in the time dependent framework. Continuously and piecewise constant periodic time varying networks will be analysed, to set the ground for the proposed approach. The extension to non periodic settings will also be discussed.

Equations

$${\dot{u}}_i(t)$$

$${\dot{v}}_i(t)$$

$${\dot{u}}_i(t)$$

$${\dot{v}}_i(t)$$

$$\vec{x}(\tau )-\vec{y}(\tau )=\mathcal{O}(ϵ),$$

$$\vec{y}{}^{\prime }(\tau )=ϵ\left[F(\vec{y})+⟨\mathcal{L}⟩\vec{y}\right],\vec{y}(0)={\vec{x}}_0\in {\mathbb{R}}^{2N},$$

$$\mathbf{A}(t)=\{\begin{array}{ll}{\mathbf{A}}_1& \text{if {t/T}∈[0,γ)}\\ {\mathbf{A}}_2& \text{if {t/T}∈[γ,1)},\end{array}$$

$$\frac{d\delta X}{dt}=\mathbf{M}(t)\delta X,$$

$${ϵ}^{\ast }=\min \{ϵ>0:\forall s\ge ϵ:\rho ({\mathbf{Q}}_s)\le 1\}.$$

$$\delta X_{ϵ}(ϵT)=e^{{\mathbf{\Gamma }}_{ϵ}(ϵT)}{e}^{ϵT{\sum }_{n\ge 1}{ϵ}^{n-1}{\mathbf{H}}_n}\delta X_{ϵ}(0).$$

$${\mathbf{H}}_1=\frac{1}{T}{\int }_0^T\mathbf{M}(t)dt={\partial }_{X}F(\bar{u},\bar{v})+⟨\mathcal{L}⟩,$$

$${\dot{u}}_i$$

$${\dot{v}}_i$$

$$\dot{\vec{x}}(t)=F(\vec{x})+\mathcal{L}(t/ϵ)\vec{x},$$

$$\mathcal{L}(t)\vec{x}:=\left(\begin{array}{cc}{D}_u\mathbf{L}(t)& 0\\ 0& D_v\mathbf{L}(t)\end{array}\right)\vec{x}.$$

$$\dot{x}=ϵf(t,x),x(0)=x_0\in D\subset {\mathbb{R}}^n,$$

$$⟨f⟩(y)=\frac{1}{T}{\int }_0^{T}f(t,y)dt,$$

$$\dot{y}=ϵ⟨f⟩(y),y(0)=x_0\in D\subset {\mathbb{R}}^n.$$

$$⟨\mathbf{L}⟩=\frac{1}{T}{\int }_0^T\mathbf{L}(t)dt,$$

$$⟨\mathbf{L}⟩=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\mathbf{L}(t)dt.$$

$$\dot{\vec{x}}(t)=F(\vec{x})+⟨\mathcal{L}⟩\vec{x},\vec{x}(0)={\vec{x}}_0\in {\mathbb{R}}^{2N},$$

$$\dot{\vec{x}}(t)=F(\vec{x})+\mathcal{L}(t/ϵ)\vec{x},$$

$${\vec{w}}^{\prime }(\tau )=ϵ\left[F(\vec{w})+\mathcal{L}(\tau )\vec{w}\right].$$

$$\vec{w}(\tau )-\vec{y}(\tau )=\mathcal{O}(ϵ),$$

$$\vec{y}{}^{\prime }(\tau )=ϵ\left[F(\vec{y})+⟨\mathcal{L}⟩\vec{y}\right],\vec{y}(0)=\vec{w}(0)\in {\mathbb{R}}^{2N}.$$

$$\vec{x}(t)-\vec{z}(t)=\mathcal{O}(ϵ)\text{ for all 0<ϵ<ϵ∗ and t=O(1).}$$

$$\dot{\vec{x}}(t)=F(\vec{x})+\mathcal{L}(t/ϵ)\vec{x},$$

$$\dot{\delta \vec{x}}={\mathbf{M}}_{ϵ}(t)\delta \vec{x}\text{mbox}with{\mathbf{M}}_{ϵ}(t)={\partial }_{x}F(\bar{x})+\mathcal{L}(t/ϵ),$$

$$\mathcal{L}(t)={\mathcal{L}}^{[\sigma (t)]}.$$

$$\dot{\delta \vec{x}}={\mathbf{M}}_{ϵ}^{[\sigma (t)]}\delta \vec{x}\text{mbox}with{\mathbf{M}}_{ϵ}^{[\sigma (t)]}={\partial }_{x}F(\bar{x})+{\mathcal{L}}^{[\sigma (t/ϵ)]}.$$

$$\delta \vec{x}(ϵt_{k+1})=\exp (ϵ{\tau }_{k+1}{\mathbf{M}}_{ϵ}^{[\sigma (ϵt_k)]})\delta \vec{x}(ϵt_k),k\in \mathbb{N}.$$

$$\rho \left\{\exp (ϵ{\tau }_k{\mathbf{M}}_{ϵ}^{[\sigma (ϵt_k)]}),k\in \mathbb{N})\right\}<1.$$

Full text

Spatially extended systems can spontaneously yield a multitude of patterns, archetypal imprint of an inherent drive to coordinated self-organisation pikovsky ; murray2 ; cross . In many cases of interest, the subtle interplay between nonlinearities and diffusion seeds a symmetry breaking instability (discovered by Turing in his pioneering work on morphogenesis), which paves the way to the emergence of a rich and colourful gallery of patchy motifs turing ; mimura ; baurmann . Distinct populations of homologous constituents often interact via an intricate architecture of nested couplings, which can be adequately represented as complex heterogeneous networks. Elaborating on the mechanisms that instigate pattern formation for reaction-diffusion systems hosted on discrete network-like supports is hence central to a wide gallery of phenomena of broad applied and fundamental impact nakao .

Topology is known to play, in this respect, a role of paramount importance. Surprisingly, patterns can rise on a directed support directed , even if they are formally impeded on a regular, continuum or discrete, spatial medium. Further, self-organisation may proceed across multiple, interlinked networks, by exploiting the multifaceted nature of resources and organisational skills. Multiplex networks in layers have been therefore introduced as a necessary leap forward in the modelling effort: these are particularly relevant to transportation systems, the learning organisation in the brain and to understanding the emergent dynamics in social communities. Interestingly, the interaction between adjacent layers can yield self-organised patterns which are instead impeded in the limit of decoupled layers AsllaniPRE2014 . Patterns on individual layers can also fade away due to cross-talking between contiguous slabs. All these examples, point to the pivotal role exerted by the topology of the underlying network in shaping the system response and driving the onset of non-trivial instabilities.

In several realms of application, the networks that specify the routing of the spatial or physical interactions are not static, but, instead, do evolve in time holme2012 ; holme ; delvenne ; Masuda . This is an important additional ingredient that should be eventually accommodated for, so to improve on our current understanding of the mechanisms that are ultimately responsible for the appearance of structured patterns across heterogenous networks. In the framework of temporal networks, diffusion processes perra as well as random walks starnini2012 ; MasudaKlemmEguiluz2013 have been studied in recent years, mostly under the assumption of piecewise constant time varying networks Vankeerberghen2014 ; Liberzon1999 ; Stilwell2006 : the topology of the assigned connections is fixed over a finite window in time, whose duration can either depend on the dynamics of the system or exogenously controlled. The couplings are then instantaneously modified (created/destroyed/rewired) and the ensuing network frozen for the subsequent time interval. Synchronisation phenomena have also been considered for nonlinear oscillators Barreto hosted on time varying networks, as well as for systems displaying generic reactions terms Kohar ; Stilwell2006 ; Boccaletti2006 .

In this Letter we consider for the first time the process of pattern formation for reaction-diffusion systems anchored on networks that evolve over time. To this end, we will initially inspect the case of a network that is periodically rearranged in time. We will in particular prove that a symmetry breaking instability which anticipates the pattern derive can be incited, by properly tuning the frequency of the imposed network dynamics. The proposed framework includes the case of a continuously time varying network (links weights change as a smooth function of time) and the one where sudden switches between distinct discrete networks configurations are to be accounted for (links are suddenly created/destroyed/rewired). Surprisingly enough, patterns can emerge for a reaction-diffusion system hosted on a piecewise constant time varying network, also when they are formally impeded on each network snapshot. The extension to non periodic settings will be also analysed, allowing us to draw a general interpretative scenario.

Imagine two different species living on a network that evolves over time and denote by $u_{i}$ and $v_{i}$ their respective concentrations, as seen on node $i$. Links can change their (positive) weights, they can fade away (the corresponding weight goes to [math]), or be created (a null weight is turned positive) or even rewired (a combination of both preceding moves). For the sake of simplicity, and without losing generality, we deal here with networks made of a constant number of nodes, that totals in $N$. The network structure is stored in a time varying $N\times N$ weighted adjacency matrix, whose elements are labelled $A_{ij}(t)$. The (weighted) connectivity, $s_{i}(t)$, of node $i$ is computed as $s_{i}(t)=\sum_{j}A_{ij}(t)$. Species can relocate in space, as follows a standard diffusive mechanism, ruled by the (time dependent) Laplacian operator $L_{ij}(t)=A_{ij}(t)-s_{i}(t)\delta_{ij}$. Remark that, for all $t$ and all $i$, $\sum_{j}L_{ij}(t)=0$, namely $\mathbf{L}(t)$ is a row stochastic matrix. When species happen to share the same node, they mutually interact via non linear terms that reflect the specificity of the problem at hand. We will deal at first with networks that are periodically reset in time, and show how adjusting the period of the modulation can stimulate the onset of the patterns.

The model can be mathematically cast in the form:

[TABLE]

where $f$ and $g$ stand for the generic nonlinear reactions terms and $D_{u}$ and $D_{v}$ label the diffusion coefficients of species $u$ and $v$, respectively. The parameter $\epsilon$ can be tuned at will, so as to modify the period of the Laplacian oscillations in time. We will in particular denote with $T$ the period obtained for $\epsilon=1$. As we shall prove, by forcing the network to oscillate at high frequencies, one can prompt a Turing-like instability. When the oscillations materialise as successive swaps between two (or more) network configurations, it is the frequency of the blinking that sets the patterns derive, also when Turing motifs cannot manifest for the model constrained on any of the considered (static) network snapshots. As a further condition, preliminary to the forthcoming development, we will assume that system (1) admits an homogeneous stable fixed point. In other words, we shall require that $(u_{i},v_{i})=(\bar{u},\bar{v})$ exists such that $\dot{u}_{i}=\dot{v}_{i}=0$, for all $i=1,\dots,N$.

Define the averaged Laplacian $\langle\mathbf{L}\rangle=\frac{1}{T}\int_{0}^{T}\mathbf{L}(t)dt$ and introduce the averaged reaction-diffusion system:

[TABLE]

A conventional linear stability analysis, calibrated so as to sense the spectral characteristics of the averaged network, can be performed and allows one to conclude on the possibility for patterns à la Turing in the framework of the time independent picture (2). Imagine that conditions are met so to enable a symmetry breaking instability of the homogeneous stable fixed point, as triggered by an external non homogeneous perturbation. Then, it can be rigorously proven that $\epsilon^{*}>0$ exists such that Eqs. (1) exhibit patterns of the Turing type, for $0<\epsilon<\epsilon^{*}$. Here, $1/\epsilon^{*}$ acts as an effective high frequency drive for self-organised spatial motifs to develop. The idea of the proof is sketched in the following, further details are relegated in the Supplementary Information (SI).

Label $\tau=t/\epsilon$ and introduce the compact notation $\vec{x}=(u_{1},\dots,u_{N},v_{1},\dots,v_{N})$. Then, by virtue of the theorem of averaging one can show that $\epsilon^{*}>0$ exists such that for all $0<\epsilon<\epsilon^{*}$ and $\tau=\mathcal{O}(1/\epsilon)$:

[TABLE]

where $\vec{y}(\tau)$ is the solution of the averaged system:

[TABLE]

subject to the initial conditions $\vec{x}~{}(0)~{}=~{}\vec{y}~{}(0)$. Since the reaction term, $F(\vec{x})=(f(u_{1},v_{1}),\dots,f(u_{N},v_{N}),g(u_{1},v_{1}),\dots,g(u_{N},v_{N}))$, is time independent, the average solely affects the diffusive part, organised in the block matrix $\mathcal{L}(t):=\left(\begin{smallmatrix}D_{u}\mathbf{L}(t)&0\\ 0&D_{v}\mathbf{L}(t)\end{smallmatrix}\right)$.

The solutions of, respectively, the time dependent system (with $\epsilon\in(0,\epsilon^{*}]$) and the averaged one, stay close for times $\mathcal{O}(1)$. Hence, if an injected perturbation destabilises the homogenous solution of (2), the same holds when the perturbation is applied to system (1). Remark however that the two systems behave similarly, only for a finite time span. We are in fact unable to prove an effective or exponential stability à la Nekhoroshev. Stated differently, we cannot guarantee that the smaller $\epsilon$ the longer the two systems agree. This in turn implies that the asymptotic patterns, which follows the initial instability and as obtained for respectively the averaged and the time dependent models, may, in principle, differ. We will return later on characterising the critical threshold $\epsilon^{*}$.

We shall here discuss a specific application (termed the twin network case) aimed at clarifying the conclusion reached above. To this end we will consider the Brusselator model, a widely studied model of chemical oscillators. This choice amounts to setting $f(u,v)=1-(b+1)u+cu^{2}v$ and $g(u,v)=bu-cu^{2}v$ where $b$, $c$ stand as free parameters. When the model is made spatially extended, the coaction of diffusion and reaction terms can disrupt the homogeneous fixed point, paving the way to the subsequent patterns derive. This occurs when the parameters are set so as to make the dispersion relation ($\max\Re\lambda$, the real part of the complex exponential growth rate $\lambda$) positive over a finite window of wave-numbers. When the diffusion takes place on a heterogeneous network, the dispersion relation is defined on a discrete support that coincides with the eigenvalues of the Laplacian operator. If the spectrum of the discrete Laplacian operator falls in a region where the continuous dispersion relation is negative, no instability can take place on the network, even if it can develop on a continuum support. This observation will be central for what follows.

Let us thus consider two simple networks, made of $N$ nodes arranged on a periodic ring, and label with $\mathbf{A}_{1}$ and $\mathbf{A}_{2}$ their respective adjacency matrices. The $N$ nodes are connected in couples (so we will assume for simplicity $N$ even), via symmetric links. In the first case (the network that is encoded in $\mathbf{A}_{1}$) nodes $2k-1$ and $2k$, for $k=1,\dots,N/2$ are tight together. In the other, the couples are formed by nodes $(2k,2k+1)$ for $k=1,\dots,N/2-1$, with the addition of pair $(N,1)$ (see Fig. 1 panel a). Both networks yield an identical Laplacian spectrum: two degenerate eigenvalues are found, $\Lambda^{1}=0$, with multiplicity $N/2$ (i.e. the number of connected components the network is made of) and $\Lambda^{N}=-2$, with multiplicity $N/2$. The parameters are set so that patterns cannot emerge for the Bussellator model defined on each of the networks introduced above. This is exemplified in Fig. 1 panel c.

We now introduce a periodically time varying network, specified by the following adjacency matrix $\mathbf{A}(t)$:

[TABLE]

where $\{r\}$ denotes the fractional part of the real number $r$ and $\gamma\in(0,1)$ is the fraction of occurrence of the first network (and thus $1-\gamma$ for the second one) during the period $T>0$. The Laplacian matrix $\mathbf{L}(t)$ can be computed accordingly. We get the averaged Laplacian matrix $\langle\mathbf{L}\rangle=\gamma\mathbf{L}_{1}+(1-\gamma)\mathbf{L}_{2}$ where $\mathbf{L}_{i}$ is the Laplacian matrix of the $i$–th network, $i=1,2$. As an example, assume $\gamma=0.3$, $T=1$ and compute the dispersion relation associated to the averaged network $\langle\mathbf{A}\rangle=\gamma\mathbf{A}_{1}+(1-\gamma)\mathbf{A}_{2}$. From inspection of Fig. 1 panel c), it is immediately clear that $\Lambda^{\alpha}$, the eigenvalues of the Laplacian $\langle\mathbf{L}\rangle$, fall in a region where $\Re\lambda_{\alpha}>0$. Hence, the Brussellator placed on the average network exhibits Turing patterns (see Fig. 1 panel d). Based on the above, we can find a positive $\epsilon^{*}$ such that for all $0<\epsilon<\epsilon^{*}$ the time dependent system (1) displays patterns. This is because this latter system is close enough to the averaged counterpart to be able to follow its orbits. On the other hand, no patterns can develop if $\epsilon$ is too large. These conclusions are clearly demonstrated in Fig. 2 panel a).

We shall now complement the analysis by elaborating on a recipe to estimate the critical threshold $\epsilon^{*}$. Assume $\epsilon=1$ and linearise system (1) close to the stable homogeneous equilibrium $(u_{i},v_{i})=(\bar{u},\bar{v})$, for all $i=1,\dots,N$:

[TABLE]

where $\delta X=(u_{1}-\bar{u},\dots,u_{N}-\bar{u},v_{1}-\bar{v},\dots,v_{N}-\bar{v})^{T}$ and $\mathbf{M}(t)=\partial_{X}F(\bar{u},\bar{v})+\mathcal{L}(t)$, where $\partial_{X}F(\bar{u},\bar{v})$ is the Jacobian of the reaction part evaluated at the homogeneous equilibrium. $\mathbf{M}(t)$ is $T$-periodic and piecewise constant, more precisely $\mathbf{M}(t)=\partial_{X}F(\bar{u},\bar{v})+\mathcal{L}_{1}:=\mathbf{M}_{1}$ if $\{t/T\}\in[0,\gamma)$ and $\mathbf{M}(t)=\partial_{X}F(\bar{u},\bar{v})+\mathcal{L}_{2}:=\mathbf{M}_{2}$ if $\{t/T\}\in[\gamma,1)$. We can thus solve exactly Eq. (5) over one period ${\delta X}(T)=e^{\mathbf{M}_{2}(1-\gamma)T}e^{\mathbf{M}_{1}\gamma T}\delta X(0)=:\mathbf{Q}\delta X(0)$ where the rightmost equality defines the monodromy matrix, $\mathbf{Q}$. Hence, for any $m\geq 1$, ${\delta X}(mT)=\mathbf{Q}^{m}\delta X(0)$ . The stability of the system is therefore determined by the spectral radius of the monodromy matrix, $\rho(\mathbf{Q})$. If the spectrum is entirely contained in the unit disk, $\delta X$ converges to zero and no patterns are allowed. Conversely, the applied perturbation can develop and eventually result in asymptotic stationary patterns. The same reasoning applies when $\epsilon$ is allowed to change freely. The monodromy matrix now reads $\mathbf{Q}_{\epsilon}=e^{\epsilon\mathbf{M}_{2}(1-\gamma)T}e^{\epsilon\mathbf{M}_{1}\gamma T}$ and one can conclude on the stability of the scrutinised system by computing its associated spectral radius $\rho(\mathbf{Q}_{\epsilon})$. Assume, in line with the above, that $\rho(\mathbf{Q})=\rho(\mathbf{Q}_{\epsilon})\rvert_{\epsilon=1}<1$. Then, the critical threshold $\epsilon^{*}$ can be quantified as:

[TABLE]

In Fig. 2 (b) the predictions of the theory are challenged versus numerical simulations, performed for the twin network setting, as introduced above, for different values of the system size $N$. Excellent agreement is found which testifies on the correctness of the proposed interpretative scenario. Interestingly, $\epsilon^{*}$ returns a different profile in $N$, depending on whether $N/2$ is even or odd. The explanation resides in the corresponding magnitude of the dispersion relations, as depicted in the inset of Fig. 2 (b). The instability is less pronounced when $N/2$ is odd, for small enough $N$. In this case, it is therefore necessary to push the modulation to the high frequencies even further to accompany the patterns crystallisation. A closed approximate formula for $\epsilon^{*}$ can be also derived following an approach that is discussed in the SI.

Summing up, patterns can be enforced in a reaction-diffusion system by switching periodically between distinct network configurations. Other examples aimed at illustrating our findings are discussed in the SI. In the remaining part of this Letter we generalise the results along a few possible avenues of investigations, delegating technical details and representative applications to the annexed SI.

We can straightforwardly adapt the above reasoning to the case of a periodic, continuously time varying network and provide a reliable estimate for the critical threshold $\epsilon^{*}$. In this case $\mathbf{M}(t/\epsilon)=\partial_{X}F(\bar{u},\bar{v})+\mathcal{L}(t/\epsilon)$ is a matrix whose elements are continuous functions of time. By using the Floquet-Magnus expansion Magnus1954 ; Casas2001 ; Blanes2009 one can write ${\delta X_{\epsilon}}(t)=e^{\mathbf{\Gamma}_{\epsilon}(t)}e^{t\mathbf{H}_{\epsilon}}\delta X_{\epsilon}(0)$, where $\mathbf{\Gamma}_{\epsilon}(t)$ is $\epsilon T$-periodic. Writing $\mathbf{H}_{\epsilon}=\sum_{n\geq 1}\mathbf{H}_{\epsilon,n}$, one can prove that $\mathbf{H}_{\epsilon,n}=\epsilon^{n-1}\mathbf{H}_{n}$. By evaluating the solution over one period $\epsilon T$ we obtain

[TABLE]

Since $\mathbf{\Gamma}_{\epsilon}(\epsilon T)=0$, the stability of the homogeneous solution is determined by the eigenvalues of $\mathbf{H}_{\epsilon}$, the so-called characteristic exponents of the system. For sufficiently small $\epsilon$, one can write $\mathbf{H}_{\epsilon}=\mathbf{H}_{1}+\mathcal{O}(\epsilon)$. The eigenvalues of $\mathbf{H}_{\epsilon}$ are hence close to those of $\mathbf{H}_{1}$. Using the Magnus series one can determine

[TABLE]

which is nothing but the modified Jacobian matrix that applies to the averaged system. If $\epsilon$ is sufficiently small, the periodic, continuously time varying system behaves as the average one. Computing the spectrum of (8) enables one to infer about the stability of the homogeneous solution. For a reaction-diffusion system hosted on network that is periodically modulated, with period $\epsilon T$ and $\epsilon<\epsilon^{*}$, patterns can arise, provided the average system can undergo a Turing like instability. Finally, we emphasise that the analysis readily extends to the case of non-periodically time varying networks. The only requirement is the existence of the generalised average $\langle\mathbf{L}\rangle=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\mathbf{L}(t)dt$. Relevant examples are addressed in the SI and include (i) quasi-periodic time dependent networks and (ii) networks constructed using aperiodic infinite words from a binary alphabet, as the Sturmian words.

In conclusion, we have here shown that, by properly tuning the network topology over time, one can drive the emergence of self-organised patterns, reminiscent of the Turing instability. The reaction-diffusion system endowed with a time varying support behaves as its average analogue, provided the network dynamics is made sufficiently fast. Hence, patterns are obtained for a model attached to a network that changes over time, if they are predicted to occur on the associated average system. This is a novel route to pattern formation that we conjecture relevant for all those applications where different species interact via a network that can adjust in time, as follows either an endogenous or exogenous drive.

Acknowledgments 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. The scientific responsibility rests with its author(s). D.F. acknowledges financial support from H2020-MSCA-ITN-2015 project COSMOS 642563.

Appendix A SUPPLEMENTARY INFORMATIONS

Appendix B Proof of the main result

The aim of this section is to provide a detailed derivation of our main result that we will cast in the form of a theorem. Let us first recall the definition of the model and then proceed to prove the theorem.

We assume, for the sake of simplicity, two different species that live and interact on a periodically time varying network. The number of nodes $N$ is kept fixed, while pairwise weighted links adjust in time. The network structure is embedded in a periodic time varying weighted adjacency matrix $A_{ij}(t)$, that results in time varying link strength, $s_{i}(t)=\sum_{j}A_{ij}(t)$. Further, we introduce the time varying Laplacian matrix defined $L_{ij}(t)=A_{ij}(t)-s_{i}(t)\delta_{ij}$. Observe that for all $t$ one has $\sum_{j}L_{ij}(t)=0$, namely $\mathbf{L}(t)$ is a row stochastic matrix.

Species diffuse across the network and interact on the same node as dictated by specific non-linear reaction terms. The model (equations (1) in the main body of the paper) reads:

[TABLE]

where $f$ and $g$ are non-linear fucntions and $D_{u}$ and $D_{v}$ denote the the diffusion coefficients associated, respectively, to species $u$ and $v$. The parameter $\epsilon$ sets the rate of the network evolution. For a subsequent use let us rewrite the above system in a more compact form by defining the $2N$–dimensional vector $\vec{x}=(u_{1},\dots,u_{N},v_{1},\dots,v_{N})$:

[TABLE]

where $F(\vec{x})=(f(u_{1},v_{1}),\dots,f(u_{N},v_{N}),g(u_{1},v_{1}),\dots,g(u_{N},v_{N}))^{T}$ and

[TABLE]

To study the emergence of patterns for model  (9), as follows a symmetry breaking instability of the Turing type, one cannot invoke the standard machineries, as eigenvectors and eigenvalues depend on time. To compute the dispersion relation following the canonical approach additional assumptions are to be enforced, as e.g. discussed in  Boccaletti2006 . To avoid this, we proceed with an alternative route which involves dealing with the so called theorem of averaging hereafter recalled for consistency (the interested reader can consult e.g. verhulst for further details on this topic):

Theorem 1** (Averaging in the periodic case)**

Let us consider the system

[TABLE]

where $f$ and $\partial_{x}f$ are defined, continuous and bounded in $[0,\infty)\times\mathbb{R}^{n}$ and assume $f(t,\cdot)$ to be $T$-periodic.

Let $\langle f\rangle(y)$ be the time average of $f(t,y)$, that is

[TABLE]

and let $y(t)$ be the solution of

[TABLE]

Then $x(t)-y(t)=\mathcal{O}(\epsilon)$ for $t=\mathcal{O}(1/\epsilon)$.

Based on the theorem above, and as discussed in the main body of the paper, we can relate the behaviour of system (9) to the behaviour of its averaged homologue, defined by replacing $\mathbf{L}(t)$ with the time averaged operator $\langle\mathbf{L}\rangle$ defined as:

[TABLE]

where $T>0$ is the period. As a side remark, we anticipate that our result holds true also if the time evolution of the network is not periodic, provided a generalised average exists, namely:

[TABLE]

More precisely we can state and proof the following theorem

Theorem 2

Let us define the averaged system 111$F(\vec{x})$ being time independent, the averaging affects only on the Laplacian matrix.:

[TABLE]

where $\vec{x}=(u_{1},\dots,u_{N},v_{1},\dots,v_{N})$, $\langle\mathcal{L}\rangle=\left(\begin{smallmatrix}D_{u}\langle\mathbf{L}\rangle&0\\ 0&D_{v}\langle\mathbf{L}\rangle\end{smallmatrix}\right)$ and $\langle\mathbf{L}\rangle$ is defined using Eq. (14). Assume moreover that the above system exhibits standard Turing patterns, namely that a stable homogeneous equilibrium $\bar{x}=(\bar{u},\dots,\bar{u},\bar{v},\dots,\bar{v})$ of system (15) can turn unstable, for an appropriate choice of the parameters involved, upon application of a non homogeneous perturbation.

Then there exists $\epsilon^{*}>0$ such that the fast varying version of Eq. (10)

[TABLE]

displays patterns of the Turing type for all $0<\epsilon<\epsilon^{*}$.

Proof Let us first introduce the rescaled time $\tau=t/\epsilon$ and the adapted variables $\vec{w}(\tau)=\vec{x}(t)\rvert_{t=\epsilon\tau}$, to rewrite Eq. (16) as (with the notation ${}^{\prime}=d/d\tau$):

[TABLE]

We can then apply theorem 1 and obtain

[TABLE]

for all $0<\epsilon<\epsilon^{*}$, where $\vec{y}(\tau)$ is the solution of the averaged system

[TABLE]

Back to the original variables, $\vec{x}$ and time $t$, the last statement can be rephrased as follows. Let $\vec{x}$ be the solution of Eq. (16) and $\vec{z}(t)=\vec{y}(t/\epsilon)$ the solution of Eq. (15) with the same initial conditions, then

[TABLE]

Namely the two above orbits stay close over a macroscopic time window. Hence, if a symmetry breaking instability occurs for system  (15), paving the way to the subsequent pattern derive, the same holds for the original system (16).

Appendix C On the estimate of the critical threshold

The aim of this section is to provide a detailed derivation of the formulae employed in the main text to estimate $\epsilon^{*}$. To this end we consider again the reference system written in a compact vectorial form:

[TABLE]

where $\vec{x}=(u_{1},\dots,u_{N},v_{1},\dots,v_{N})^{T}$, $F(\vec{x})=(f(u_{1},v_{1}),\dots,f(u_{N},v_{N}),g(u_{1},v_{1}),\dots,g(u_{N},v_{N}))^{T}$ and $\mathcal{L}(t)=\left(\begin{smallmatrix}D_{u}\mathbf{L}(t)&0\\ 0&D_{v}\mathbf{L}(t)\end{smallmatrix}\right)$.

Assume there exists $\bar{u},\bar{v}$ such that $\bar{x}=(\bar{u},\dots,\bar{u},\bar{v},\dots,\bar{v})^{T}$ is an homogeneous equilibrium of (19), and write $\delta\vec{x}=\vec{x}-\bar{x}$ to denote a tiny perturbation. Linearising system (19) close to the homogeneous solution yields:

[TABLE]

where $\partial_{x}F(\bar{x})$ is the Jacobian matrix of $F$ evaluated at $\bar{x}$.

In the rest of the section we shall assume that no Turing instability is allowed for $\epsilon=1$. In other words, the imposed perturbation fades away when $\epsilon$ is set equal to $1$, and system (19) converges to its homogeneous equilibrium. At variance, the time averaged system, which corresponds to the limiting case $\epsilon\approx 0$, can exhibit patterns of the Turing class. In the following, we shall write for short $\mathbf{M}_{\epsilon=1}=\mathbf{M}$. Our goal is to compute, in this setting, the critical threshold for the control parameter $\epsilon$ below which the fast varying version of the inspected system can also display Turing like patterns.

C.1 Temporal network of contact sequences

Let us consider a finite collection of networks indexed by $k\in\mathcal{K}$ and label $\left\{\mathbf{L}^{[k]}:k\in\mathcal{K}\right\}$ their associated Laplacian matrices. For fixed reaction terms, we then define, for every $k\in\mathcal{K}$, the matrix $\mathbf{M}^{[k]}=\partial_{x}F(\bar{x})+\mathcal{L}^{[k]}$, where $\mathcal{L}^{[k]}=\left(\begin{smallmatrix}D_{u}\mathbf{L}^{[k]}&0\\ 0&D_{v}\mathbf{L}^{[k]}\end{smallmatrix}\right)$. We further assume that for every linearised system involving a given $\mathbf{M}^{[k]}$, i.e. when the index $k$ is kept fixed, the zero solution is stable with respect to injection of a small non homogeneous perturbation $\delta\vec{x}$.

Let us now introduce a switching signal, $\sigma(t)$, that is, a piecewise constant function $\sigma:[0,\infty)\rightarrow\mathcal{K}$. This provides a direct map between time and indexes, thus enabling us to associate to each time its corresponding network

[TABLE]

Working along these lines, we generalise the construction of twin networks as discussed in the main body of the paper. For the twin networks, in fact, $\mathcal{K}=\{1,2\}$ and $\sigma(t)=1$ if $\mathrm{mod}(t,T)/T\in[0,\gamma)$ and $\sigma(t)=2$ otherwise. We are thus interested in the stability of the $1$-parameter family of linearised systems

[TABLE]

If the interval between any two consecutive discontinuities of $\sigma(t)$ is (on average) above a given threshold $\tau_{D}$, called the dwell time of the system, then the zero solution of Eq. (22) is stable (see for instance Morse_1996_dwell_time ).

Let $0=t_{0}<t_{1}<t_{2}<\dots$ be the times of discontinuities of $\sigma(t)$, and define $\tau_{k}$ by $t_{k}+\tau_{k+1}=t_{k+1}$ for $k\in\mathbb{N}$. The switching signal $\sigma(t)$ is constant over every interval $[t_{k},t_{k+1})$ and therefore the solution of Eq. (22) can be explicitly computed as

[TABLE]

We then get a discrete time system that is stable under injection of a small non homogeneous perturbation, if the joint spectral radius satisfies:

[TABLE]

Observe that even for a set composed by only two matrices, the joint spectral radius is cumberosme to compute Vankeerberghen2014_JSR . To allow a more straightforward computation of the critical threshold $\epsilon^{*}$, we proceed by making the simplifying assumption that the switching signal $\sigma(t)$ has a finite number of discontinuities $(n)$ and it is $T-$periodic, with $T=\sum_{i=1}^{n}\tau_{i}$.

As a result, we have that the solution of (23) at every integer multiple of periods $m\epsilon T$ is given by

[TABLE]

where the time-evolution operator over one period,

[TABLE]

where we have used that $\mathbf{M}_{\epsilon}^{\left[\sigma(t)\right]}=\mathbf{M}^{\left[\sigma(t/\epsilon)\right]}$ and thus $\mathbf{M}_{\epsilon}^{\left[\sigma(\epsilon t)\right]}=\mathbf{M}^{\left[\sigma(t)\right]}$ (see Eq. (22)).

The critical value of the acceleration parameter $\epsilon$ is thus given by

[TABLE]

where $\rho$ denotes the spectral radius.

C.2 Continuously varying networks

We now consider the more general case of a continuously time varying network, with Laplacian matrix $\mathbf{L}(t)$ whose entries are continuous integrable functions of time. The solution of the linearised system with $\epsilon=1$,

[TABLE]

can be written as

[TABLE]

if and only if the matrices $\mathbf{M}(t_{1})$ and $\mathbf{M}(t_{2})$ commute for any pair $t_{1},t_{2}$. Let us observe that, as $\mathbf{M}$ contains both a term related to the reactions and another related to the diffusion, the above assumption may not be satisfied, even if the Laplacian matrices associated to the networks do commute. The commutation between the Laplacian matrix and the reaction part is also needed. However one can always write:

[TABLE]

where $\mathbf{\Omega}(t)=\sum_{n=1}^{\infty}\mathbf{\Omega}_{n}(t)$, is known as the Magnus expansion Magnus_1954 ; Blanes2009_expansion-and-some-applications . From

[TABLE]

where the $B_{n}$’s are the first Bernoulli numbers ($B_{1}=-1/2$), $\mathrm{ad}_{X}\mathbf{A}=[\mathbf{X},\mathbf{A}]=\mathbf{X}\mathbf{A}-\mathbf{A}\mathbf{X}$ denote the matrix commutator and $\mathrm{ad}_{X}^{0}\mathbf{A}=\mathbf{A}$, $\mathrm{ad}_{X}^{k}\mathbf{A}=[\mathbf{X},\mathrm{ad}_{X}^{k-1}\mathbf{A}]$, one can iteratively compute all the terms of the series. For example,

[TABLE]

One can also prove Moan_2008-convergence-magnus-series the convergence of the above series provided $\int_{0}^{t}\left\|\mathbf{M}(\tau)\right\|_{2}\mathrm{d}\tau<\pi$.

As previously done when dealing with the switching network case, we assume for simplicity that the network evolution is periodic. Hence, the entries of $\mathbf{L}(t)$ are $T-$periodic continuous functions. The Floquet theorem ensures thus that the solution of (27) can be written as

[TABLE]

where $\mathbf{P}$ is $T$-periodic and bounded, and $\mathbf{F}$ does not depend on time. The asymptotic stability of the null solution is determined by the eigenvalues of $\mathbf{F}$, known as the characteristic exponents of the system. The Floquet-Magnus expansion Casas2001_floquet-exponential-perturbative allows one to write

[TABLE]

where the involved matrices are obtained as the sums of series 222The convergence for the series for $F$ follows from the above quoted condition $\int_{0}^{T}\left\|\mathbf{M}(\tau)\right\|_{2}\,d\tau<\pi$, whereas for the series related to $\mathbf{\Lambda}(t)$ we have the more restrictive sufficient condition $\int_{0}^{T}\left\|\mathbf{M}(\tau)\right\|_{2}\,d\tau<0.20925$.

[TABLE]

The matrix $\mathbf{\Lambda}(t)$ is $T$-periodic, with $\mathbf{\Lambda}(0)=\mathcal{O}$. Since $\delta\vec{x}(T)=\exp[\mathbf{\Omega}(T)]\delta\vec{x}(0),$ we have the identity:

[TABLE]

Let us now turn to considering the system (20) in the case $\epsilon<1$. We write $\mathbf{\Omega}_{\epsilon}(t)$ for the Magnus expansion corresponding to $\mathbf{M}_{\epsilon}(t)=\mathbf{M}(t/\epsilon)$, and apply the same notation to $\mathbf{\Lambda}(t)$ and $\mathbf{F}$. This leads to

[TABLE]

where $\mathbf{\Lambda}_{\epsilon}(t)$ is $\epsilon T$-periodic, and

[TABLE]

The stability is determined by the eigenvalues of the characteristic exponent $\epsilon T\mathbf{F}_{\epsilon}$, or equivalently, by looking at the sign of the real part of the eigenvalues of $\mathbf{\Omega}_{\epsilon}(\epsilon T)/(\epsilon T)$.

Let us observe that for every $n\geq 1$, we have $\mathbf{\Omega}_{\epsilon,n}(\epsilon t)=\epsilon^{n}\mathbf{\Omega}_{n}(t)$. In conclusion we obtain

[TABLE]

and, as a consequence, the stability of the zero solution of the fast-varying system ($\epsilon<1$) is determined by the sign of the real part of the eigenvalues of $\sum_{n=1}^{\infty}\epsilon^{n-1}\mathbf{\Omega}(T)/T$.

As remarked when discussing the implication of the averaging theorem, in the limiting case $\epsilon\rightarrow 0$, the averaged network can be used to determine the stability of the system. Indeed, if $\epsilon>0$ is small enough, the eigenvalues of $\sum_{n=1}^{\infty}\epsilon^{n-1}\mathbf{\Omega}_{n}(T)/T=\mathbf{\Omega}_{1}/T+o(\epsilon)$ lie on the same portion of the imaginary axis as the eigenvalues of

[TABLE]

The critical value of $\epsilon$ can be computed by imposing

[TABLE]

where $\sigma_{\gamma}$ is the spectrum of (40) for $\epsilon=\gamma$.

The computation of the first terms of the Magnus expansion, truncated as desired so as to reach the necessary precision on $\epsilon^{*}$, can be carried out by using the recursive formulas first obtained in Klarsfled_1989_recursive_terms :

[TABLE]

where the $\mathbf{S}_{n}^{(j)}$ follow from the recursion

[TABLE]

Working out these formulas explicitly,

[TABLE]

The above formulas are straightforward to implement, but prove numerically heavy to calculate since they involve multiple integrals of nested commutators. To overcome this limitation one can use the Time-stepping method which consists in performing a partition of the interval $[0,T]$ in $m$ subintervals of length $h$ and then compute the Taylor series for $\mathbf{M}(t)$ using the midpoint $t^{[k]}=(k-1)h+\frac{h}{2}=(k-\frac{1}{2})h$:

[TABLE]

Inserting this relation into the recursive formulas given above yields Blanes2000_high-order , up to order $5$:

[TABLE]

Here, we have used the simplified notation $[x_{1},x_{2},\ldots,x_{j}]=[x_{1},[x_{2},[\ldots,[x_{j-1},x_{j}]\ldots]]]$. Over each time interval $[t_{k-1},t_{k})$, we get as a viable approximation to the time evolution operator from $\delta\vec{x}(t_{k-1})$ to $\delta\vec{x}(t_{k})$, relative to the reference case $\epsilon=1$ :

[TABLE]

After appropriate truncation of the series, we have

[TABLE]

Note that here no integration of $\mathbf{M}(t)$ is required, so the computation is much faster to handle numerically. The analysis can be readily extended to the setting $\epsilon<1$. If one considers again $m$ subintervals 333Let us note that if $m$ is large enough, the convergence of the Magnus series on each subinterval is guaranteed. in $[0,\epsilon T]$, we then have

[TABLE]

As before, we only need to compute the terms of the Magnus series for the original system, that we then multiply by the appropriate integer power of $\epsilon$. In doing so, we can obtain an estimate for the critical value $\epsilon^{*}$ by evaluating the spectral radius of the matrix that appears in the right hand side of the latter equation.

C.3 A closed-form approximate expression for $\epsilon^{*}$

The aim of this section is to present an alternative method, directly inspired by the proof of the averaging theorem, to compute the critical threshold $\epsilon^{*}$. This procedure avoids dealing with the computation of the monodromy matrix and returns a closed formula for $\epsilon^{*}$, which proves adequate versus numerical simulations.

Digging into the proof of the averaging theorem one can appreciate that it relies on an invertible change of coordinates matching the two systems (Eq. (17) and Eq. (18)). Our idea is hence to quantify $\epsilon^{*}$ by determining the range of $\epsilon$ for which the sought invertibility can be achieved. The required change of variables is given by

[TABLE]

where $\vec{u}$ is explicitely given by

[TABLE]

where the terms involving the reaction part, $F$ is Eq. (19), cancel out because they do not depend explicitly on time.

To be able to invert relation (56) one has to require that its Jacobian

[TABLE]

is non singular, namely that $\det J\neq 0$. This is certainly true if $\epsilon=0$. We could hence determine an upper bound for $\epsilon$ as:

[TABLE]

By using the form for $\vec{u}$ given by Eq. (57) one can obtain

[TABLE]

and recalling the explicit form of $\mathcal{L}$ we get:

[TABLE]

In the following we apply this strategy to compute $\epsilon^{*}$ for the twin network setting introduced in the main body of the paper.

C.4 The critical threshold for the twin network

The analysis presented in the previous section C.3 can be pushed further in the case of switching twin–networks (see main body of the paper for technical details on the model formulation). In this case we can in fact exaclty compute $\langle\mathbf{L}\rangle=\gamma\mathbf{L}_{1}+(1-\gamma)\mathbf{L}_{2}$ and thus the right hand side of Eq. (59). In fact let

[TABLE]

then we have

[TABLE]

where $\{r\}$ denotes the fractional part of the real number $r$.

Let $\Lambda_{12}^{\alpha}$, $\alpha=1,\dots,N$, be the eigenvalues of $\mathbf{L}_{1}-\mathbf{L}_{2}$, then the roots of Eq. (59) depend on $t$ and are of the form (we assume to order the eigenvalues such that $\Lambda_{12}^{1}=0$)

[TABLE]

where

[TABLE]

which has a minimum at $t=\gamma T$. Hence Eq. (58) returns:

[TABLE]

where we set $\Lambda_{12}^{N}=\max_{\alpha}\lvert\Lambda_{12}^{\alpha}\rvert$, keeping in mind that the spectrum of $\mathbf{L}_{1}-\mathbf{L}_{2}$ contains both positive and negative eigenvalues.

Appendix D Another relevant application: the blinking network

Many interactions involving humans can be described by a binary contact lasting for a finite amount of time after which the two individuals get apart. During these events information can be mutually shared. If contacts are not mediated by electronic devices, viruses can pass from an individual to the other, hence triggering the epidemic spreading. As a first approximation, the dynamics of interaction within an large community can be broken down into a collection of successive pairwise exchanges, that extend over a finite window of time. Working in this setting, it is interesting to elaborate on the conditions (e.g. characteristic time of interactions) that discriminate between an homogenous or pattern like (in terms of virus load or information content) asymptotic equilibrium.

To this end we here consider a collection of networks, each made of $N$ isolated nodes with the exception of two nodes that are connected via an undirected link. Any two networks in this collection differ because a different link is activated. There are hence $N(N-1)/2$ different networks in total. All networks have the same Laplacian spectrum given by $\Lambda_{1}=0$ with multiplicity $N-1$ (i.e. the number of connected components the network is made of) and $\Lambda_{N}=-2$ with multiplicity $1$ (i.e. the minimal value the eigenvalues can assume in the case of complete bipartite network).

Let us now define a periodically time varying network $\mathbf{A}(t)$ as follows. Fix a period $T>0$ and divide the interval $[0,T]$ into $m=N(N-1)/2$ equal disjoint subintervals $\tau_{i}$, $\cup_{i}\tau_{i}=[0,T]$ and $\tau_{i}\cap\tau_{j}=\emptyset$ if $i\neq j$. Hence, $|\tau_{i}|=\tau=T/m$ (of course this assumption can be relaxed and subintervals with different length allowed for).

Assume to order the networks in the collection from $1$ to $N(N-1)/2$, then define $\mathbf{A}(t)$ by

[TABLE]

where $\mathbf{A}_{i}$ is the adjacency matrix associated to one of the networks of the collection. Then we periodically repeat the procedure, i.e. if $t>T$ then consider $(t\mod T)$ instead of $t$. Differently stated, during each time interval, a novel link is created between a newly selected pair of nodes and the connection active during the preceding time window deleted. Hence, at each time there is one and only one active link (see top left panel of Fig. 3 for an illustrative explanation). By making use of the aforementioned notation, we could also say $\mathbf{A}(t)=\mathbf{A}^{[\sigma(t)]}$, where $\sigma(t)=i$ if $t\in\tau_{i}$. Because of the assumption of equal length intervals $\tau_{i}$, the adjacency matrix of the averaged network is given by

[TABLE]

while the Laplacian matrix reads

[TABLE]

whose spectrum is $\Lambda_{1}=0$ with multiplicity $1$ and $\Lambda_{N}=-2/(N-1)$ with multiplicity $N-1$. One can thus choose the network size $N$ to ensure that the relation dispersion admits a positive real part, once evaluated in $\Lambda_{N}=-2/(N-1)$ (see bottom left panel of Fig. 3).

To carry out one test we consider again the Brusselator model running on top of a such time varying network. The reaction terms read therefore $f(u,v)=1-(b+1)u+cu^{2}v$ and $g(u,v)=bu-cu^{2}v$, where $b$, $c$ are parameters of the model fixed that we assign to be $b=8$ and $c=10$. The homogeneous equilibrium is thus $u_{*}=1$ and $v_{*}=0.8$. The diffusion coefficients are set to $D_{u}=3$ and $D_{v}=10$.

We also assign the remaining model parameters in such a way that each network of the collection ($\mathbf{A}_{i}$) cannot yield Turing patterns. This amounts to requiring $\lambda_{1}=\Re\lambda(\Lambda_{1})<0$ (the homogeneous equilibrium is stable) and $\lambda_{N}=\Re\lambda(\Lambda_{N})<0$ (see bottom left panel of Fig. 3).

Applying Theorem 2 one can find a positive $\epsilon^{*}$ such that for all $0<\epsilon<\epsilon^{*}$, the fast varying time dependent system (16) exhibits Turing patterns. This is due to the fact that, under this operating condition, system (16) is close enough to its averaged analogue to be able to follow its orbits. Conversely, system (16) does not exhibit Turing patterns if $\epsilon$ is too large (see bottom right panel of Fig. 3).

Appendix E Non-periodic networks

Our main result, Theorem 2, holds true also if the time varying network is not periodic: as mentioned earlier, we solely require the existence of the generalised average

[TABLE]

The aim of this section is thus to shortly discuss some examples of a non-periodic time-varying networks for which Eq. (64) is well defined.

The first example exploits the non commensurability of the periods of two periodic networks. More precisely, let $G_{1}(t)$, respectively $G_{2}(t)$, be a $T_{1}$, respectively $T_{2}$, periodic network. Then, assuming $T_{1}/T_{2}\not\in\mathbb{Q}$, the network $G(t)=G_{1}(t)+G_{2}(t)$ is quasi-periodic and one easily obtains that

[TABLE]

where $\langle\mathbf{L}_{i}\rangle$ is the average Laplacian matrix of the $i$–th network, computed over the period $T_{i}$.

Building on this setting one could generate interesting applications. Imagine that for all finite $t$, the network associated to $G(t)$ cannot develop patterns. Then we could speculate on the possibility that patterns arise in the averaged network in the limit $t\rightarrow\infty$. For instance, given two adjacency matrices $\mathbf{A}_{1}$ and $\mathbf{A}_{2}$ and two parameters $\gamma_{j}\in(0,1)$, $j=1,2$, we define

[TABLE]

then $G(t)=G_{1}(t)+G_{2}(t)$ is quasi-periodic. For all finite $t$ no patterns are allowed on $G(t)$ if they are not observed on neither $\mathbf{A}_{1}$, $\mathbf{A}_{2}$ nor $\mathbf{A}_{1}+\mathbf{A}_{2}$. However, the average limit network $\langle G\rangle=(\gamma_{1}+\gamma_{2})\mathbf{A}_{1}+(2-\gamma_{1}+\gamma_{2})\mathbf{A}_{2}$ could still yield patterns for a suitable choice of $\gamma_{j}$. We do not provide an explicit numerical evidence of this mechanism here, but instead turn to consider a different model of non periodic networks. For this further example, a numerical realization will be also given to support the conclusion of the analysis.

The final example that we are going to illustrate is based on aperiodic words (or binary strings) and follows the construction of a map associating to the letters (or entries) of the word a network. Let us consider for the sake of simplicity, two time periodic networks built on two different time windows of identical total duration:

[TABLE]

where $T>0$ is the period and $\gamma\in(0,1)$ a fixed parameter.

Given a string made of $w=(i_{n})_{n\in\mathbb{N}}$, $i_{n}\in\{0,1\}$, we can finally obtain a time varying network by selecting $P_{0}$ or $P_{1}$ according to the entries appearing in $w$. The Laplacian matrix $L(t)$ associated to this latter network will be aperiodic if the word $w$ also is. On the other hand, because of the definition of $P_{0}$ and $P_{1}$ (they are based on windows of the same length) one straightforwardly obtains: $\left\langle\mathbf{L}\right\rangle=\gamma\mathbf{L}_{1}+(1-\gamma)\mathbf{L}_{2}$.

Let us conclude by providing an explicit example (see also Fig. 5). Let us call the complexity of $w$ the application $p_{w}:\mathbb{N}\rightarrow\mathbb{N}$ which gives for every $n$ the number of subwords of length $n$ in $w$. Using the Morse-Hedlund theorem444A word $w$ is eventually periodic iff there exists $n\in\mathbb{N}$ such that $p_{w}(n)\leq n$. , an aperiodic word $w$ of minimal complexity is such that

[TABLE]

The words with minimal complexity are called Sturmian words 555Note that it is trivial to give examples of aperiodic words with higher complexity.. A well-known example is the Fibonacci word $w^{\textit{Fibonacci}}=w_{0}w_{1}w_{2}\dots$ where the $n$–th subword, $w_{n}$, is given by

[TABLE]

Observe that such word can be obtained using the following recursion

[TABLE]

where $s_{1}\cdot s_{2}$ is the concatenation of $s_{1}$ and $s_{2}$, from which the name Fibonacci follows. A numerical realization of the proposed scheme is reported in Figure 5: patterns appear in the time varying network also if they are formally impeded on each subnetworks components.