{'fraction_non_alphanumeric': 0.033677811550151975, 'fraction_numerical': 0.0052684903748733535, 'mean_word_length': 4.577757685352622, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stable resonances and signal propagation in a chaotic network of coupled units

We apply the linear response theory developed in \cite{Ruelle} to analyze how a periodic signal of weak amplitude, superimposed upon a chaotic background, is transmitted in a network of non linearly interacting units. We numerically compute the complex susceptibility and show the existence of specific poles (stable resonances) corresponding to the response to perturbations transverse to the attractor. Contrary to the poles of correlation functions they depend on the pair emitting/receiving units. This dynamic differentiation, induced by non linearities, exhibits the different ability that units have to transmit a signal in this network.


I. INTRODUCTION
Currently, there is considerable research activity in network dynamics. This is clearly motivated by the wide expansion of communication media (mobile phones, Internet, multimedia, etc.), but also by the growing interest in network modeling of biological processes (neural networks, genetic networks, ecological networks ...). A large part of these studies focuses on topological properties of the underlying graph. However, in many cases, the nodes of the networks are units behaving in a non linear way. For example, in a communication network a relay regenerates (amplifies) weak signals, but it has a finite capacity and saturates if too many signals arrive simultaneously. A neuron has a non linear response to an input current, a gene expression is determined by a non linear function of the regulatory proteins concentration, etc.. These non-linearities might modify the network abilities in a drastic way. For example, a relay may have a high graph connectivity ("hub"), but the dynamics drives it close to its saturation point, so that it has a weak reactivity to the changes in the inputs coming from the other units and a poor capacity to transmit information. Consequently, the information is transmitted via other units, possibly weaker links, and, in this regime, these units become temporary "hubs" though they may have a low graph connectivity, while the main hub is decongested. In biological networks similar effects may arise. For example, the capacity of a neuron to transmit a specific excitation strongly depends on its state, determined itself by the overall currents coming from afferent neurons. This suggests us that the mere study of the graph topological structure of a network with non linear units is not sufficient to capture the full dynamical behavior. However, there are relatively few studies which analyze the joint effect of topology of the network and non-linearity. Nevertheless, these networks are dynamical systems with a large number of degrees of freedom, and so dynamical systems theory and statistical mechanics provide a powerful framework to state problems in a well-defined way and to propose solutions.
In this paper, we analyze the following situation. We consider a network composed by a set of N units receiving and transmitting signals. At each time step t the unit i receives a bench of signals coming from each unit connected to it, and it emits, at time t + 1, a signal which is a sigmoid function of the global input [see eq. (2)]. In the model studied below, the global dynamics has generically a chaotic attractor, provided that the non linearity of the transfer function is sufficiently large [see Section II]. In spite of the presence of chaos it is possible to analyze how a periodic signal of weak amplitude, superimposed upon a chaotic background, is transmitted in the network. However, as discussed above, this analysis requires the consideration of the network structure as well as nonlinear effects.
The main tool we use for this investigation is the linear response theory developed by D. Ruelle [1] for hyperbolic dynamical systems (e.g. dissipative systems with a chaotic attractor) in a non equilibrium steady state. This theory allows us to compute explicitly the variation of the average value of a generic observable, induced by a time dependent signal of weak amplitude. Indeed, provided that the amplitude of the signal is sufficiently small (but finite), this variation is a linear function of the signal and a linear response operator is explicitly given in terms of the dynamic evolution. In our case, this operator has a simple expression (see eq. (6)). The effects of a periodic signal emitted by a unit on a receiving unit is characterized by the Fourier transform of the linear response, called susceptibility in the sequel (see section IV). This gives us a frequency response curve (see Indeed, in this theory, the linear response operator is the sum of two contributions. There is a regular term, corresponding to the response to perturbations "parallel" to the attractor (more precisely locally projected along the unstable manifold). This term is actually a correlation function [2] and, consequently, it obeys classical relations such as the Fluctuationdissipation theorem. The poles of its Fourier transform are called Ruelle-Pollicott resonances [3] or "unstable" poles. They give the rate of mixing of the chaotic system or equivalently, the relaxation rate to equilibrium for a perturbation "on" the attractor. These poles are independent of the observable. Therefore, in our case, they are independent of the pair emitting/receiving unit (see Fig 2). When focusing on the response to real frequency one observes therefore resonance peaks common to all pair of units, and these peaks are also present in the Fourier spectrum of the corresponding correlation function.
The second term corresponds to the response to perturbations locally projected along stable manifolds, namely transverse to the attractor. Therefore, it exists only in the dissipative case. It does not obey fluctuation-dissipation theorem and has drastically different properties than the first term. In particular its poles ("stable" poles) are expected to be distinct from the unstable poles. In this paper, we indeed exhibit such stable poles.
To the best of our knowledge, this is the first example where these poles are explicitly exhibited, though their existence was theoretically proved. Moreover, we show numerically that the stable poles depend on the pair emitting/receiving unit (see Fig. 3). When focusing on the response to real frequency one observes therefore specific resonances peaks (see Fig. 1). This shows that a unit receiving a periodic signal emitted from another unit may respond in a specific way to this signal, the amplitude depending both on the signal frequency and on the emitting unit. Note that according to the discussion above this effect cannot be observed by studying correlation functions.
The paper is organized as follows. In section II we introduce the model and discuss its properties. The section III recall briefly the main results of Ruelle linear response theory used in this paper. An explicit computation of the linear response is performed. It shows the explicit contributions of the network topology and of the non linearity in a signal propagation. In section IV we compute numerically the frequency response curve and discuss the different resonance peaks. The poles of the complex susceptibility for a few pairs of units are computed and compared in the section V. Our main conclusions are then drawn.

II. MODEL
Consider the following dynamical system, originally proposed in the context of Neural Networks [see [4,5,6] and references therein]. The output signal is a function of the weighted sum of the signals arriving at i at time t and is given by: The weights J ij 's may be positive (excitatory), negative (inhibitory) or zero (no direct link between j and i). They are in general non symmetric (J ij = J ji ). Thus, the matrix of weights, J, defines an oriented graph such that there is a link from j to i if and only if J ij = 0. The global dynamics can also be written as: where where the parameter g controls the non linearity. In terms of input/output ratio, a unit amplifies weak signals (if g > 1), but with a limited capacity: f "saturates" if the local field is too strong, and the variations of the output signal are all the weaker as the local field is big. Thus, the capacity of i to retransmit a signal emitted from k does not only depend on the weight J ik but also on the state of saturation of i when it receives the signal coming from k. Note also that the Jacobian matrix Therefore, in this model, the dynamical contraction is closely related to the saturation of the sigmoid transfer function.
In order to emphasize the effects of the nonlinearity and minimize the effect of the network topology, one may assume that the network is fully connected and that the J ij 's are drawn randomly with respect to some smooth distribution (uniform, Gaussian, etc . . . ).
As an example, one may fix the average value [J ij ] = 0 and the variance [J 2 ij ] = 1 N (to ensure the correct normalization of the local field with the size N). This example is interesting because the system (2) exhibits a wide variety of dynamical regimes (static, periodic, quasi periodic, chaotic). More precisely, it has been shown in [5] that it generically exhibits a transition to chaos by quasi periodicity when g increases. Note that the same transition occurs if the network is sparse [4] with K > 2 neighbors (K can be random) chosen at random, provided the variance of the J ij 's scales like 1 K . However, we do not address this case in this paper since we want to minimize the effect of the network structure. Note also that this type of transfer function allows dynamical regimes where several attractors coexist. It has been indeed shown in [5,6] that, adding a threshold θ to the local field, there exists a region in the parameter space g, θ where two attractors coexist. This region can be analytically computed. However, in the present paper, the parameters are located in a region where there is only one attractor and all initial conditions converge to this attractor.
Let us now assume that the non linearity is large enough so that the global dynamics has a chaotic attractor (with all Lyapunov exponents bounded away from zero and at least one positive Lyapunov exponent).
We now add a signal of small amplitude ξ(t) to the output of some units. Then the evolution of the perturbed system, denoted byũ, is given by : Note that the formalism introduced below accommodates the generalization where ξ(t) depends also on u(t), but we do not consider this case here.
We want to investigate the capacity of the network to transmit signal ξ(t) superimposed upon the chaotic background. This is a complex problem since after a few time steps the total signal arriving at time t at k includes the sum of contributions corresponding to different paths followed by ξ, with different time delays. Moreover, along a path the signal can be . Finally, the dynamics being chaotic, after a sufficiently long time the signal is distorted by the nonlinearities and scrambled by mixing.
To tackle this problem we analyze how the differenceũ(t) − u(t) between the perturbed and unperturbed dynamics behaves on average as a function of ξ(t). When ξ(t) is small enough, and in spite of the initial condition sensitivity intrinsic to chaotic systems, it can be shown that this difference is a linear functional of ξ(t). This is the content of the linear response theory developed by D. Ruelle [1] for chaotic and dissipative [10] system. Some aspects of this theory are briefly recalled in the next Section.

III. LINEAR RESPONSE THEORY.
The unperturbed dynamical system (2) has a strange (chaotic) attractor for sufficiently large g. Usually, strange attractors carry a natural probability measure called the Sinai-Ruelle-Bowen (SRB) measure [7]. If one prepares the system (2) with an initial macrostate µ having a uniform density (i.e. µ(du) = du), corresponding to selecting typical initial conditions, then, provided that the limit exists, the SRB measure is the asymptotic macrostate ρ = lim t→+∞ G t µ where G t µ is the image of µ under the t-th iterate of G. The SRB measure has several remarkable features which make it "natural" [11]. One of its most important property for practical purpose is the following: If A is some observable (a smooth function of u), its average with respect to ρ, is equal to the time average along typical trajectories. This means that "ensemble average" and time average are equivalently for typical trajectories. This is especially useful for numerical computations (see next Section).
Applying a time dependent perturbation ξ(t) to the system induces time dependent changes in the statistical averages. More precisely, the natural extension of the SRB measure defined above is a time dependent SRB measure ρ t . It is given by the (weak) limit lim s→+∞G t . . .G t−s µ. The corresponding average will be denoted by < > t .
It has been established in [1] that a linear response theory exists for uniformly hyperbolic diffeomorphism [12]. In our framework, this means that, provided that ξ(t) is sufficiently small, and for any smooth observable A, the variation < A > t − < A > is proportional to ξ(t) up to small non linear corrections. In other words, ρ t is differentiable with respect to the perturbation. The derivative is called the linear response.
The theory developed by Ruelle allows one to compute the linear response, for general perturbations depending both on time t and state u, and for a general observable A. In our context, however, where the considered observables are simply the variables of systems (2)- (3), the linear response has a simple form, which can be written as: where χ(τ ) represents the averaged Jacobian matrix: for τ ≥ 0. Otherwise χ(τ ) = 0 (which is consistent with the requirement of causality).
A remarkable consequence of Ruelle theory is that χ(τ ) is a bounded function for all τ ≥ 0. In particular, it does not diverge exponentially fast, despite the presence of a positive Lyapunov exponent. As discussed below, this is essentially a consequence of exponential mixing.
In what concerns network dynamics, equation (5) is interpreted as giving the average response of unit i of the system when the network is submitted to weak signal ξ(t). In particular it is seen that if only one unit j is perturbed at time t = −1 by a kick of amplitude ǫ [that is ξ(t) = ǫe j δ(t + 1) with the Kroenecker symbol δ and the j-th unit vector e j ], then ǫ χ ij (t) gives precisely the mean response of unit i at time t. This suggests to define the susceptibility of the network as the Fourier transform of χ ij (t), namely: This matrix function will be numerically computed and studied in the next Section. We conclude the present Section by analyzing further the structure of χ ij (τ ) in the case of dynamical system (1). Here one can decompose χ ij (τ ) as : The sum holds on each possible paths γ ij (τ ), of length τ , connecting the unit k 0 = j to the unit k τ = i, in τ steps. One remarks that each path is weighted by the product of a topological contribution depending only on the weight J ij and a dynamical contribution. Since, in the kind of systems we consider, functions f are sigmoids, the weight of a path γ ij (τ ) depends crucially on the state of saturation of the units k 0 , . . . , k τ −1 at times 0, . . . , τ − 1. Especially, though a signal has many possibilities for going from j to i in τ time steps, some paths may be "better" than some others, in the sense that their contribution to χ ij (τ ) is higher. One can decompose the response function (6) into two terms. The first one is obtained by locally projecting the Jacobian matrix on the unstable directions of the tangent space. This term will be named the"unstable" response function. It corresponds to linear response of the system to perturbations locally parallel to the local unstable manifold (roughly speaking perturbations "on" the attractor). One can show that the linear response associated with this type of perturbation is in fact a correlation function, as found in standard fluctuationdissipation theorems [1]. Hence, as usual for correlation functions of a chaotic system, it decays exponentially (because of mixing) and the decay rates are associated with the poles of its Fourier transform. More precisely, these exponential decay rates correspond to the imaginary part of the complex poles of the unstable part of the susceptibility (8). Thus they will be called "unstable" poles. More generally, it can be shown that these poles are also the eigenvalues of the operator governing the time-evolution of the probability densities (which we denoted above as G t µ), the so-called Perron-Frobenius operator [3]. Therefore, these poles, whose signatures are visible in the peaks of the modulus of the correlation functions, do not depend on the observable, though some residues may accidentally vanish for a given observable.
The second term [13] is obtained by locally projecting the Jacobian matrix on the stable directions of the tangent space. It corresponds to the response to perturbations locally parallel to the local stable manifold (namely transverse to the attractor). Therefore, it is Hence, both stable and unstable terms are exponentially damped, ensuring the convergence of the series (5), but for completely different reasons. Moreover, the stable and unstable part of the linear response have drastically different properties. As a matter of fact, the stable part is not a correlation function and it does not obey the fluctuation-dissipation theorem. In particular, the unstable poles and stable poles are expected to be distinct.
In this paper, we give for the first time an evidence of this theoretically predicted effect.
Moreover, we show that the stable poles indeed allow to distinguish the units in their capacity to transmit a signal.
For this we first numerically compute the susceptibility (7) for real values of ω. The computation is based on the following remark. Let us consider perturbations ξ (1) (t) = ǫe j cos(ωt) and ξ (2) (t) = −ǫe j sin(ωt) and let u (1) , u (2) denote the variables of the corresponding per-turbed systems: Then it follows from (5) that : Note that the time dependent average response to periodic perturbation is therefore periodic.
The linear response at time t is an infinite sum corresponding to contributions of time delayed signals following different paths. Since the signal is sinusoidal the terms in this sum may interfere in a constructive way [but exponential damping prevent the series to diverge, ensuring the existence of a linear response].
Sinceχ ij (ω) is independent of t, then it is equal (for ω = 0) to the time averagê The time-dependent averages < u (k) i > t involve an average over initial conditions in the distant past. One can reasonably assume that the above average over t makes the average over the initial conditions unnecessary. Then one may write : where the u (k) i (t) (k = 1, 2) are obtained by iterating maps (9). This provides a straightforward way to compute the susceptibility, where most of the computing time goes into computing the orbits u (k) (t).
As an example, we performed a numerical computation of the dynamical system (2) where we take a fixed realization of J ij 's, with N = 8 units. There is a quasi-periodic transition to chaos as g increases. The system is studied for g = 3.5 corresponding to a positive Lyapunov exponent λ 1 = 0.158, while the second one is λ 2 = −0.183. The system is therefore weakly hyperbolic (all Lyapunov exponents bounded away from 0).
The functionχ(ω), the Fourier transform of the matrix (8), has been computed with a resolution δω = π 2048 ≈ 1.53 × 10 −3 . The average is done with 26214400 samples. We did several runs where we varied the length T of the time average in (12). We checked that the global structure is the same. In particular the amplitude of the susceptibility |χ(ω)| does not depend on T (see note [12]). Also the fluctuations decrease like 1 √ T according to the central limit theorem.
In Fig. 1 we have plotted the modulus of the susceptibilitiesχ 33 ,χ 45 ,χ 71 . Comparing these functions, one remarks that there are thin peaks essentially located at the same frequencies, with different heights. Moreover, these frequencies are harmonics of a fundamental frequency (ω 0 ∼ 0.166). This is expected from the frequency locking in the quasi-periodic transition preceding chaos. Some of these frequencies are also present in the Fourier spectrum of the correlation functions but with a smaller amplitude and some peaks are indistinguishable from the background. Instead, all harmonic peaks are revealed in the susceptibility spectrum.
But we also note that for many peaks, the width varies strongly from a pair ij to another.
This means that the resonance strength depends on which unit is excited and which unit responds. In particular, some peaks are very thin, corresponding to an accurate resonance while some others are broad. In terms of poles, this means that the imaginary part are distinct and consequently the corresponding poles are different [see next section]. Finally there are additional peaks strongly dependent on the pair ij.
Thus, a simple glance to Fig. 1 tells us that the frequency response of a unit i to the excitation emitted by a unit j strongly depends on the pair i, j. As discussed above, and numerically shown below, this difference comes from the stable part of the linear response.
Consequently, the specificity of the response is revealed only if one consider perturbations transverse to the attractor. [Note that, generically, the signal is a perturbation having local projections both on local stable and unstable spaces.]

V. UNSTABLE AND STABLE POLES.
Resonances correspond to poles in the complex plane. As a matter of fact, the position of the maximum of the peak corresponds to the real part of the pole, its width is related to its imaginary part, and the value of the maximum is related to the residue. From this observation, we developed an algorithm to estimate the residue width and locations of the poles. Let ω 0 = ω r + iω i be a simple pole ofχ and A its residue. If one multipliesχ by a phase factor e iψ then the real and imaginary part rotate continuously, without changing   the modulus. If the pole is close enough to the real axis then there exists a phase ψ such that, on the real axis, the real part has a characteristic Lorentzian shape symmetric with respect to ω r while the imaginary part is antisymmetric. Then a nonlinear curve fitting allows us to determine A, ω r , ω i . Once a local analysis has roughly determined the poles, a global nonlinear fit (Levenberg-Marquardt [9]) allows us to localize the poles with a better accuracy.
In fig. 2 we have plotted the real and imaginary part of the poles of several correlation functions. One notices that all pair of units have poles at the same value of ω, within the error bars. We have also plotted in Fig. 3 the modulus of the susceptibilitiesχ 33 ,χ 36 ,χ 63 (left column) and the corresponding poles (right column) with the poles of the correlation functions. As expected from Fig. 1 we observe common poles (unstable poles) but also distinct poles (stable poles) that, moreover, strongly depend on the pair receiving/emitting unit.
Finally, note that some poles are very close to the real axis. Since their imaginary part is related to the coherence time of the response to a kick, this tells us that the response to a pulse may subsist on quite a bit long times, though the underlying dynamics is chaotic.
[Recall however that the linear response measures variations of average value of observables].

VI. CONCLUSION.
This paper gives an example of network dynamics where the nonlinearity induces particularly prominent effects that cannot be anticipated by the mere analysis of the graph topology. In particular we exhibit a dynamic differentiation in the capacity that a unit has to transmit information. We also argue on theoretical grounds, and numerically show (see Fig. 2) that the dynamics differentiation is not revealed by correlation functions. It is purely an effect of the dynamics transverse to the chaotic attractor that must be handled with the proper tool. We show that the linear response gives quite a bit more information than the correlation function, provided that its computation takes into account the singularity of the SRB measure transversally to the attractor. This is the case with Ruelle linear response theory and this opens the perspective for developing an extension of statistical mechanics for the analysis of networks dynamics with nonlinear units.