Collective fluctuations in networks of noisy components

Collective dynamics result from interactions among noisy dynamical components. Examples include heartbeats, circadian rhythms, and various pattern formations. Because of noise in each component, collective dynamics inevitably involve fluctuations, which may crucially affect functioning of the system. However, the relation between the fluctuations in isolated individual components and those in collective dynamics is unclear. Here we study a linear dynamical system of networked components subjected to independent Gaussian noise and analytically show that the connectivity of networks determines the intensity of fluctuations in the collective dynamics. Remarkably, in general directed networks including scale-free networks, the fluctuations decrease more slowly with the system size than the standard law stated by the central limit theorem. They even remain finite for a large system size when global directionality of the network exists. Moreover, such nontrivial behavior appears even in undirected networks when nonlinear dynamical systems are considered. We demonstrate it with a coupled oscillator system.


Introduction
Understanding fluctuations in dynamically ordered states and physical objects, which consist of networks of interacting components, is an important issue in many disciplines ranging from biology to engineering. When each constituent component of a system is noisy due to, e.g., thermal fluctuations, it generally occurs that the entire system collectively fluctuates in time.
Such collective fluctuations may be advantageous or disadvantageous in functioning of the systems depending on situations. For example, reduction in noise is likely to improve information processing in retinal neural networks [1][2][3][4]. Precision of biological circadian clocks [5][6][7][8] may be improved by reduction in collective fluctuations (i.e., fluctuations in collective activities). On the other hand, maintaining a certain amount of fluctuations in an ordered state is advantageous for stochastic resonance [9] and Brownian motors [10].
Despite the relevance of collective fluctuations in a variety of systems, theoretical frameworks that formulate collective fluctuations are missing. The central limit theorem states that, if the dynamical order is simply the averaged activity of noisy components, the standard deviation of the collective fluctuation would decrease with the number N of noisy components as N −1/2 . However, scaling is unclear in systems of interacting components. Clarifying the property of collective fluctuations in such systems will give us insights into the mechanisms and design principles underlying the regulation of noise in, for example, living organisms and chemical reactions, and also into possible controls of fluctuations in collective dynamics.
In this study, we analyze an ensemble of components subjected to independent Gaussian noise that interact on general networks, including complex networks and regular lattices. We first consider a linear dynamical system, which can be regarded as linearization of various systems, such as networks of periodic or chaotic oscillators [11,12], the overdamped limit of elastic networks [13], a consensus problem treated in control theory [14]. We show that collective fluctuations are determined by the connectivity of networks. It turns out that the scaling N −1/2 is the tight lower bound, which is obtained for undirected networks. General directed networks yield a slower or nonvanishing decay of collective fluctuations with an increase in N. We then argue such nontrivial behavior appears even in undirected networks when nonlinear systems are considered. In particular, we show that linearization of coupled nonlinear oscillator systems on undirected media yields linear dynamics on asymmetric networks, such that the slow decay of the collective fluctuation is relevant.

Model and analysis
Consider a network of N components obeyinġ where x i is the state (or the position) of the ith component, √ D i is the intensity of noise, ξ i is the independent Gaussian (generally colored) noise, and w ij is the intensity of coupling and can be also regarded as originating from the Jacobian matrix of underlying nonlinear dynamical systems such as coupled oscillator systems that we consider later. We allow negative weights and asymmetric coupling; w ij can be negative or different from w ji . Equation (1) is a multivariate Ornstein-Uhlenbeck process [15,16].
For convenience, we represent Eq. (1) aṡ is the asymmetric Laplacian defined by [17,18]. L always has a zero eigenvalue with the right eigenvector u ≡ (1 . . . 1) ⊤ , i.e., Lu = 0. This eigenvector is associated with a global translational shift in state x and corresponds to the fact that such a shift keeps Eq. (1) invariant. We assume the stability of the ordered state represented by x 1 = . . . = x N in the absence of the noise (i.e., D i = 0 for all i); the system relaxes to the ordered state from any initial condition. This is equivalent to assuming that the real parts of all the eigenvalues of L are positive except for one zero eigenvalue, i.e., 0 ≡ λ 1 < Reλ 2 ≤ . . . ≤ Reλ N . This is a nontrivial condition for general networks with negative weights. However, for networks with only non-negative weights, i.e., w ij ≥ 0 (1 ≤ i, j ≤ N), this property holds true when the network is strongly connected or all the nodes are reachable by a directed path from a single node [17,19,20].
We are concerned with collective fluctuations in dynamics given by Eq. (1). To quantify their intensity, we decompose x as We call y(t)u the collective mode. In the absence of noise, the dynamical equation for y(t) is given bẏ Therefore, y(t) is a conserved quantity of the dynamics. The remainder mode ρ(t) is associated with relative motions among the components. Because of the stability assumption, ρ(t) asymptotically vanishes with characteristic time (Re λ 2 ) −1 . Therefore, all the values of eventually go to the same value y that is determined by the initial condition, i.e., y = vx(0).
In the presence of noise, we obtaiṅ Because ξ i is the independent Gaussian noise, this equation reduces tȯ where ξ(t) is the Gaussian noise having the same statistical property as that of each ξ i (t).
Thus, y(t) performs the Brownian motion with effective noise strength σ and is unbounded.
The remainder mode ρ(t) fluctuates around zero because of its decaying nature. Therefore, the long-time behavior of x i (t) = y(t) + ρ i (t) is approximately described by a single variable y(t) for any i. We denote as σ, which depends on the structure of the network, the intensity of collective fluctuations. σ can be calculated for given network.
In practice, the average activity of the population,x ≡ N i=1 x i /N, but not the activity at individual nodes, may be observed. Becausex = y + N i=1 ρ i /N and N i=1 ρ i /N can be neglected in a long run, σ also characterizes the fluctuations ofx.

General properties
We assume for simplicity that i . It is straightforward to extend the following results to the case of heterogeneous D i . The vector v is uniform, i.e., w ji are indegree and outdegree, respectively [18]. Undirected networks satisfy this condition. In this case, we obtain σ = N −1/2 , which agrees with the central limit theorem.
The normalization condition N i=1 v i = 1 guarantees that σ ≥ N −1/2 for any v. Therefore, undirected networks are the best for reducing collective fluctuations. In the case of directed or asymmetrically weighted networks, v i is generally heterogeneous, and σ > N −1/2 . We will show later that this is also the case for nonlinear systems on undirected networks. When the weight w ij is nonnegative for any i and j, the Perron-Frobenius theorem guarantees that v i is nonnegative for all i [21]. In this case, we obtain The case σ = 1 is realized by a feedforward network, in which a certain component i 0 has no inward connection (i.e., k in i 0 = 0). Then, v i 0 = 1 and v i = 0 for i = i 0 , which yields σ = 1 irrespective of N; the collective fluctuations are not reduced at all with an increase in N. When negative weights are allowed, some elements of v may assume negative values. Then, σ may be larger than 1, in which case collective fluctuations are larger than individual noise.
We note that σ 2 is the so-called inverse participation ratio [22]. σ −2 can be interpreted as the effective number of components that participate in collective activities; the remaining components are slaved.

Directed scale-free networks
We demonstrate our theory by using some example networks. First, we consider directed scalefree networks, schematically shown in Fig. 1 wherev Therefore, This approximation is sufficiently accurate for uncorrelated networks [18]. and where · is the ensemble average. When γ out < 2, a winner-take-all network is generated [23,24], and there exists a node i such that k out . When γ out ≥ 2, the extremal criterion results in the maximum degree increasing with N as N 1/(γout−1) (γ out ≥ 2) in many networks [24,25]. Then, we obtain [26] and Therefore, we obtain The fairly heterogeneous case γ out < 2, in which the average outdegree diverges as N → ∞, effectively yields a feedforward network. The case γ out ≥ 3, where the second moment of the outdegree converges for N → ∞, reproduces the central limit theorem. The latter result is shared by the directed version of the conventional random graph. The case 2 ≤ γ out < 3 yields a nontrivial dependence of σ on N. In Fig. 2(a), we compare the scaling exponent β, where σ ∝ N −β obtained from the theory (solid line; Eq. (16)) and numerical simulations of the configuration model [23,27] with the power-law degree distribution with minimum degree 3 (open circles). The fitting procedure is explained in Fig. 2

Directed lattices
The second example is the directed one-dimensional chain of N nodes depicted in Fig. 1(b).
, and w j,i = 0 (j = i − 1, i + 1). For this network, by solving vL = 0, we analytically obtain and Values of σ for various ǫ and N are plotted by solid lines in Fig. 3(a). Interestingly, for ǫ = 1, The results for the two-dimensional lattice depicted in Fig. 1(c) are plotted by solid lines in Fig. 3(b). To confirm our theory, we also carried out direct numerical simulations of Eq. (1) with Gaussian white noise for these directed lattices. The results indicated by circles in Fig. 3 indicate an excellent agreement with our theory.
A similar result is obtained for the Cayley tree (see Appendix C).

Oscillator dynamics
As an application of our theory to nonlinear systems, we examine noisy and rhythmic components. As a general, tractable, yet realistic model, we consider a network of phase oscillators [11,28,29], whose dynamical equation is given bẏ where φ i ∈ [0, 2π) and ω i are the phase and the intrinsic frequency of the ith oscillator, respectively, A ij is the intensity of coupling, and f (·) is a 2π-periodic function. We assume that, in the for small noise intensity is tested by carrying out direct numerical simulations of Eq. (19) with i is satisfied in the directed one-and two-dimensional lattices, as shown in Fig. 4(a) and (b), respectively.
When there is some dispersion in ψ i in a phase-locked state, the relation σ ≈ N −1/2 may be violated even in undirected networks. This is because the effective weight is generally asymmetric (i.e., w ij = w ji ) unless f (·) is an exact odd function. In reality, f (·) is usually not an odd function [11,[29][30][31]. As an example, we consider target patterns (i.e., concentric traveling waves), which naturally appear in spatially extended oscillator systems [11,32]. We carry out direct numerical simulations of Eq. (19) on the two-dimensional undirected lattice with linear length √ N = 50, f (φ) = sin(φ − α) + sin α, and α = π/4. Such a function may be analytically derived from a general class of coupled oscillators [11], and it approximates a variety of real systems [29,31,33]. We set ω i = ω 0 + ∆ω (∆ω ≥ 0) for 4 × 4 pacemaker oscillators in the center and ω i = ω 0 for the other oscillators, where ω 0 is arbitrary and set to 1. A target pattern is formed when there is sufficient heterogeneity in the intrinsic frequency [11]. A region with high intrinsic frequency acts as a pacemaker. A snapshot for ∆ω = 0.3 is shown in Fig. 5(a).
As observed, the radial phase gradient is approximately constant, which makes the effective network similar to the directed two-dimensional lattice depicted in Fig. 1(c). Therefore, as shown in Fig. 5(b), v i calculated numerically decreases almost exponentially with the distance from the center. We find that the dependence of σ on N, shown in Fig. 5(c), is similar to that for directed lattices.
We emphasize that the network is undirected (i.e., A ij = A ji ). We have also theoretically confirmed that our results are valid for the continuous oscillatory media under spatial block noise, which models chemical reaction-diffusion systems (see Appendix D).

Conclusions
In summary, we have obtained the analytical relationship between collective fluctuations and the structure of networks. In undirected networks, the fluctuations decrease with the system size N as N −1/2 ; this result agrees with the central limit theorem. In general directed networks, the collective fluctuations decay more slowly. For example, in directed scale-free networks, we obtain N −β with 0 < β < 1/2. In networks with global directionality, the fluctuations do not vanish for a large system size. We have also demonstrated that such nontrivial dependence appears even in undirected networks when nonlinear systems are considered. We have focused on systems of nonleaky components. Results for coupled leaky components will be reported elsewhere.
Our results are distinct from earlier results demonstrating the breach of the central limit theorem due to heavy-tailed noise [34] or the correlation between the noise in different elements [35,36].
Finally, because our theory is based on a general linear model, it can be tested in a variety of experimental systems. An ideal experimental protocol is provided by photo-sensitive Belousov-Zhabotinsky reaction systems, in which the heterogeneity, noise intensity, and system size can be precisely controlled by light stimuli [32]. Experiments with coupled oscillatory cells, such as cardiac cells and neurons under an appropriate condition, would be also interesting.

Acknowledgments
We thank Istvan Z. Kiss

Appendix A: Derivation of the collective mode
To derive the collective mode y(t)u, we note that there exists a nonsingular matrix P such thatL ≡ P −1 LP is its Jordan canonical form [17,21]. We assume thatL 11 = λ 1 = 0 and Under the variable change (y y r ) ≡ P −1 x ∈ R N , where y ∈ R and y r ∈ R N −1 , the coupling term is transformed into −L(y y r ). Then, in the absence of the dynamical noise, y = N i=1 v i x i is a conserved quantity, which is the collective mode. ρ(t) in Eq. (3) is given by P r y r , where P r is the N by N − 1 matrix satisfying P = (u P r ).

Appendix B: Collective fluctuations in regular lattices with arbitrary dimensions
Consider a directed two-dimensional square lattice with a root node. As depicted in Fig. 1(c), the edges descending from the root node and those approaching the root node in terms of the graph-theoretic distance are given weight 1 and ǫ (0 ≤ ǫ ≤ 1), respectively. We define layers such that the layer ℓ (≤ ℓ max ) is occupied by the nodes whose distance from the root node is equal to ℓ. Layer 0 contains the root node only, and layer ℓ (≥ 1) contains 4ℓ nodes. We consider the lattice within a finite range specified by ℓ ≤ ℓ max . Note the difference from the case of the one-dimensional chain examined in the main text ( Fig. 1(b)), where the root node is located at the periphery of the chain. However, the scaling of σ is not essentially affected by this difference.
The symmetry guarantees that the four nodes in layer 1 have the same value of v i . Consider a node in Fig. 1(c) that is labeled 2 and adjacent to two nodes labeled 1. There are four such nodes. The equation in vL = 0 corresponding to this node is given by (2 + 2ǫ)v 2 = 2ǫv 1 + 2v 3 .
The other four nodes labeled 2 in Fig. 1(c) yield a different equation (1 + 3ǫ) Similarly, we obtain (2 + 2ǫ)v ℓ = 2ǫv ℓ−1 + 2v ℓ+1 for all but four nodes in layer ℓ. The other four nodes satisfy (1 + 3ǫ)v ℓ = ǫv ℓ−1 + 3v ℓ+1 . Despite this inhomogeneity, v ℓ ∝ ǫ ℓ satisfies all these equations. By counting the number of nodes in each layer, the properly normalized solution is The difference between the one-and two-dimensional cases lies in the number of nodes in each layer, which affects the normalization of v ℓ and hence the value of σ. In the limit of a purely feedforward network, σ is independent of the system size, i.e., lim ǫ→0 σ = 1. In the case of undirected networks, the central limit theorem is recovered, i.e., lim ǫ→1 σ = N −1/2 . In the limit of infinite space, we obtain For a general dimension d, layer 0 has a single root node, and layer ℓ (1 ≤ ℓ ≤ ℓ max ) has From Eq. (S.24), we obtain . (S.26) Note that lim ǫ→0 σ = 1 and lim ǫ→1 σ = N −1/2 . In the limit ℓ max → ∞, Eq. (S.25) becomes

Appendix C: Collective fluctuations in the Cayley tree
Consider a Cayley tree with degree k and a specific root node. We assume that the maximum distance from the root node is equal to ℓ max . The edges descending from the root node and those approaching the root node are assigned weight 1 and ǫ, respectively. The exact value of v i in layer ℓ, denoted by v ℓ without confusion, is obtained via By solving Eq. (S.29), we obtain exists only when ǫk < 1, and it is equal to Appendix D: Target patterns in continuous media under spatial block noise We show that our results for the coupled oscillator system in the d-dimensional lattice are also valid for that in the continuous Euclidean space. We assume that Gaussian spatial block noise is applied. This type of noise has been used in experiments [32].
We consider the d-dimensional nonlinear phase diffusion equation given by where r ∈ R d is the spatial coordinate, ω > 0 is the intrinsic frequency, ν > 0 is the diffusion constant, and µ > 0 is the coefficient of the nonlinear term [11]. The term s(r) represents the localized heterogeneity, which is positive near the origin and vanishing otherwise.
The synchronous solution corresponding to the target pattern is written as φ(r, t) = Ωt + ψ(r), where Ω and ψ(r) satisfy Ω = ω + ν∇ 2 ψ + µ (∇ψ) 2 + s(r). (S.34) Let x(r, t) be a small deviation from the target pattern defined by x ≡ φ−(Ωt+ψ). Linearizing Eq. (S.33) using x(r, t), we obtain ∂ t x(r, t) = Lx, where the linear operator L is given by We define the inner product as Now, we introduce the perturbation to Eq. (S.33) as follows: where ξ(r, t) represents a weak perturbation to the target pattern. Similarly to Eq. Let us assume that ξ(r, t) is the Gaussian spatial block noise characterized by where ℓ is the vector index for the block R d (ℓ). Using Eq. (S.41), Eq. (S.40) is transformed dr v(r). (S.44) From Eq. (S.43), we find that the intensity of the collective fluctuation is given by Note that v ℓ satisfies the normalization condition as follows: dr v(r) = dr v(r) = 1.