General scaling of maximum degree of synchronization in noisy complex networks

The effects of white noise and global coupling strength on the maximum degree of synchronization in complex networks are explored. We perform numerical simulations of generic oscillator models with both linear and non-linear coupling functions on a broad spectrum of network topologies. The oscillator models include the Fitzhugh–Nagumo model, the Izhikevich model and the Kuramoto phase oscillator model. The network topologies range from regular, random and highly modular networks to scale-free and small-world networks, with both directed and undirected edges. We then study the dependency of the maximum degree of synchronization on the global coupling strength and the noise intensity. We find a general scaling of the synchronizability, and quantify its validity by fitting a regression model to the numerical data.


Introduction
The emergence of collective and synchronous dynamics in large ensembles of coupled units is a ubiquitous phenomenon in nature and engineering. Its study has attracted much attention in a variety of fields, such as neuroscience, biology, physics, chemistry and social sciences [16,19]. For instance, there are proliferating indications that strong synchronization on large scales is related to pathological conditions of the human brain, e.g. epileptic seizures and Parkinsonʼs disease [24]. Subjects of current research include the comprehension of common properties of network synchronization in dependence on the individual node dynamics, the network topology, internode coupling types and the influence of various types of noise [3,4,17,18]. In the context of physiological networks and network medicine, for instance, a strong relationship between the topology of physiological networks and their physiological functions has recently been observed [5,6,12].
The focus of this study is the influence of uncorrelated white noise on the maximum degree of synchronization. While many studies investigated the phase transition associated with the onset of macroscopic synchronization [22,23], here we are particularly interested in the synchronizability in the aftermath of the phase transition. The effect of noise on the phasesynchronization of non-linear oscillators has for example been studied in [26], and it is known that white noise prohibits the capability of a system to achieve full synchronization by decreasing the maximum degree of synchronization [3]. It is not clear, however, how this decrease of synchronizability relates to noise intensity. Furthermore, it has been shown that the network topology has a great influence on the time-evolution of local patterns of synchronization on the path towards global coherence [11]. The question, though, whether the topology of the network has an influence on the maximum degree of synchronization, has not been answered yet.
In order to answer these questions, we develop a numerical simulation framework and study the dependency of the maximum degree of synchronization on the global coupling strength and the noise intensity. The framework incorporates three basic types of well known oscillators, namely the Fitzhugh-Nagumo model, the Izhikevich model and the Kuramoto phase oscillator. The oscillators are coupled by both linear and non-linear coupling functions. The coupling topologies include regular, random, small-world, scale-free and highly modular networks, with both directed and undirected edges.
We find a general scaling of the maximum degree of synchronization, and quantify its validity by fitting a regression model to the numerical data.
The paper is organized as follows. The models and coupling functions will be introduced in section 2, followed by a description of the network topologies in section 3. In section 4 the numerical simulation setup will be explained, and the measures for synchronization will be presented. The results are stated in section 5, and the concluding remarks are given in section 6.

The models
In the following, we will introduce three different oscillator models used in our simulations. The models are widely known and were chosen such that we would obtain a diverse set of oscillator models. The Kuramoto model is a ubiquitous phase oscillator model, simple enough to be mathematically tractable, yet sufficiently complex to display a large diversity of synchronization patterns. It is flexible enough to be adapted to many different contexts (e.g. biological models, associative memory models and laser arrays) [1]. The Fitzhugh-Nagumo model provides a simple yet basic representation of firing dynamics and has been broadly used as a model for cardiac cells and spiking neurons [10,15]. The Izhikevich model is a biologically plausible neuron model, capable of reproducing spiking and bursting behaviour of known types of neocortical and thalamic neurons [13]. The parameters of the Izhikevich model were chosen such that the model would reproduce a bursting behaviour, in order to further distinguish itself from the Fitzhugh-Nagumo model, thus expanding the diversity of the oscillator models included in this study.

Kuramoto model
The stochastic Kuramoto model [16] for N coupled and identical phase oscillators where θ i is the phase of the ith oscillator, ω π = 2 is its associated natural frequency, g the global coupling strength, 〈 〉 k the average graph connectivity (〈 〉 ≡ k L N 2 , with L denoting the total number of (weighted) links). ξ i stands for Gaussian uncorrelated white noise sources with expectation where D in will be referred to as the noise level, and M ij is the (weighted and/or directed) adjacency matrix of the simulated network. A weighted link is a link associated with a scalar value, quantifying properties as for instance the frequency of contact between actors in social networks, or the number of synapses connecting a pair of neurons in neural networks. Additionally, one distinguishes between directed and undirected graphs by whether the edges possess directional information or not, respectively.

Izhikevich model
The Izhikevich model [13] for N coupled oscillators ( = i N 1 ,..., ) is described by the ordinary differential equations with an after-spike resetting: 30, then is set to and is set to d According to [14], v i represents the membrane potential and u i a membrane recovery variable.

Fitzhugh-Nagumo model
The second neurological model we consider is the Fitzhugh-Nagumo model [9] for N coupled oscillators ( = i N 1 ,..., ). It is composed of the following differential equations: The variable v i represents the membrane potential and u i the recovery variable for the neuron membrane potential. Again, synaptic currents are transmitted via I i and ξ i in stands for Gaussian white noise sources as in equation (2) (2) and (3) with a noise intensity of = D 0.05 in is drawn in red. The two nullclines are drawn in green and blue. (b) The corresponding time-series of the membrane potential v(t) of the uncoupled Izhikevich neuron is drawn in green, the dashed black line represents the low-pass filtered time-series x(t) of the membrane potential (see equation (15)).

Coupling functions
For the Izhikevich and the Fitzhugh-Nagumo models, simulations are performed with both electrical and chemical coupling functions. Hence, I i in equation (4) and (7) takes the form: where ξ i ex stands again for Gaussian uncorrelated white noise sources with expectation i ex j ex ex ij but here it is accounting for extrinsic noise sources such as gap junctions and chemical synapses (e.g. synaptic release noise).
2.4.1. Electrical coupling. The electrical transmission, I el i , in equation (9), can be realized with a linear function of the form [20]: where g el represents the global electrical coupling strength, v i and v j stand for the membrane potentials of the post-synaptic and the pre-synaptic neurons respectively, and N is the number of neurons in the network. The local coupling strength between two connected neurons is obtained by the weighted adjacency matrix M of the simulated network.  (2) and (3) with a noise intensity of = D 0.05 in is drawn in red. The two nullclines are drawn in green and blue. (b) The corresponding time-series of the membrane potential v(t) of the uncoupled Fitzhugh-Nagumo neuron is drawn in green, the dashed black line represents the low-pass filtered time-series x(t) of the membrane potential (see equation (15)).

Chemical coupling.
There are several ways of modelling chemical synapses. In [20] for example, the authors use the approach of adding a first order dynamic for each synapse. A more effective computational implementation which still conserves the crucial properties is given by the following function [7] ∑ Γ Although taking into account several aspects of chemical synapses, some properties such as transmission delay are neglected in this model.

The networks
We have chosen six different networks from a broad spectrum of topologies, ranging from allto-all connectivity over regular, scale-free, small-world and modular networks to a random topology, including directed and undirected edges. We did so in order to thoroughly investigate the influence of the network topology on the maximal degree of synchronization and to test the generality of the scaling introduced below.
(a) An unweighted and undirected all-to-all network A, consisting of 256 nodes where every node is connected with any other node (global-coupling topology), as depicted in figure 3(a).
(c) figure 3(c) shows a computer-generated graph with two hierarchical levels of communities, as proposed in [2]. The undirected and unweighted network H consists of 256 nodes and is structured into two predefined hierarchical community levels. The inner communities consist of 16 nodes each and the outer communities consist of 64 nodes each. Each node has 13 links within its inner community, four links within its outer community and one more link with any other randomly chosen node in the network, adding up to a total of 1015 links in the entire network.
Furthermore, we have selected the real-world network of the somatic nervous system of the soil nematode C.elegans. The nervous system of C.elegans is the only one that has been almost completely mapped down to the synaptic level, and shares properties of small-world and scale-free networks [25]. The data is based on the most complete database to date, provided by [25]. It is composed of two adjacency matrices. The electrical synapse network G, undirected and weighted, connecting the 279 somatic neurons by a total of 887 (514, discarding weights) gap junctions and the chemical synapse network S, directed and weighted, connecting the neurons by a total of 6394 (2194, discarding weights) chemical synapses. The weight of a pair of connected neurons reflects the number of electrical (chemical) synapses connecting it.
(d) The simulations are performed on the combined network = + C G S, depicted in figure 3(d) (weights have been discarded in the Figure), which is simply the sum of the two adjacency matrices of the gap junction and the chemical synapse network. Gap junctions are thus treated as double-sided directed connections. This combined network consists of 279 nodes and 8168 (2990, discarding weights) directed connections.
(e) The unweighted realization of this network, as depicted in figure 3(e), will be referred to as U.
Furthermore, simulations are performed on a rewired surrogate network of the unweighted graph of C.elegans U, where only the degree-distribution is preserved. It is obtained by iteratively swapping randomly selected edges [21] of U. At each iteration, two links are chosen at random (( ↦ n n 1 2) and ( ↦ n n 3 4)) and rewired (( ↦ n n 1 4) and ( ↦ n n 3 2)), unless the respective new links do not already exist or introduce self-loops. Repeating this process often enough, the entire internal structure of the original network is destroyed, except the degree distribution.
(f) The random, directed and unweighted network, degree-matched to the unweighted combined network of C.elegans U, will be referred to as D. Its adjacency matrix is shown in figure 3(f).

Numerical simulation setup
In order to quantify the decrease in maximum synchronization of complex networks of oscillators due to Gaussian white noise and in dependence on the global coupling strength, we develop a simulation framework, as described in the following.
For the Izhikevich and the Fitzhugh-Nagumo model, the stochastic differential equations are solved by a standard first-step Euler method. Sufficient accuracy was achieved with a step size of Δ = t 0.1. Note that in all simulations, the intrinsic noise D in and the extrinsic noise D ex were set to equal values, therefore the noise level will simply be referred to as ≡ = D D D in ex . To avoid a priori synchronizations, initial conditions are drawn randomly from a uniform distribution and 50,000 iterations are calculated.
For the Izhikevich model, the generated time-series mainly contain two frequencies, a fast occurrence of spikes and a slow occurrence of bursts, as can be seen in figure 1(b). Since we are only interested in the synchronization of bursting activity, the following low-pass filter was applied to the signal = = = x v : (15) i t i t , 0 , 0 with a = 0.90. The time-series generated with the Fitzhugh-Nagumo model were smoothed by the low pass filter as well. To measure the synchronicity between two time-series, x i and x j , the Pearson correlation coefficient is calculated for the low pass filtered signal of each pair of neurons in the network. It is defined as where μ x ( ) is the mean value and σ x ( ) the standard deviation of the time-series. Because the initial conditions are chosen randomly, they are not necessarily close to the systemʼs attractor. The transient is therefore discarded, and the Pearson correlation coefficient is calculated from the th 20.000 iteration onwards. Furthermore, to eliminate random synchronizations, for each set of parameters (network, model, coupling method, global coupling strength and noise level) ten realizations are calculated and the arithmetic mean of all realizations = ∑ = R R ij r ij

Results
We are interested in the influence of the noise level D and the global coupling strength g on the mean correlation 〈 〉 R for the Izhikevich and the Fitzhugh-Nagumo model, and the order parameter O for the Kuramoto model, respectively. For the sake of convenience, the order parameter for a simulation with the Kuramoto model, O , will be referred to as 〈 〉 R as well, The following analysis of the simulations is based on the interpretation of the mean correlation as a function of the noise level and the global coupling strength: 〈 〉 = 〈 〉 g D R R ( , ). In total, 3064 simulations were performed, each corresponding to a specific model, coupling method, network, global coupling strength and noise level. They are divided into 17 subsets of data as shown in table 1.
In figures 4(a) and 4(b), the dependency of the average synchronization 〈 〉 R on the coupling strength g and the noise level D is shown for the all-to-all network A of electrically coupled Fitzhugh-Nagumo neurons and for the degree-matched random network of C.elegans, D, of Kuramoto phase oscillators, respectively. The result is a two-dimensional surface in threedimensional euclidean space, where every grid point on the surface represents the mean of ten realizations of a simulation and is associated with the corresponding mean correlation 〈 〉 R . Lines of constant noise level D are projected on the ( 〈 〉 g R )-plane (green lines), lines of constant coupling strength g are projected on the ( 〈 〉 D R )-plane (red lines) and lines of constant mean correlation 〈 〉 R are projected on the (gD)-plane (blue lines). Comparing the shape of the surfaces described by the correlation function 〈 〉 = 〈 〉 g D R R ( , ), it becomes apparent that they closely resemble each other (note that the input data, the global coupling strength g and the noise level D were rescaled to the range [0,1]). For constant noise levels D the mean correlation 〈 〉 R resembles a sigmoid function, where the steepness of these curves declines with an increasing noise level, and the inflection point moves towards larger values of g. The lines of constant coupling strength (red lines) resemble sigmoid curves as well, and for increasing values of the coupling strength, the steepness of the sigmoids declines, and the inflection point moves towards higher values of D. The blue lines, representing intersection lines of different (gD)-planes with the surface described by the correlation function, show how the coupling strength scales with the noise for fixed values of the mean correlation. The scaling seems to be of the form ∼ β g D , with β ≈ 1 or slightly larger. Remarkably, the same behaviour is observed independently of the models, coupling methods and networks that the numerical simulations were conducted with (not all figures shown). In order to quantify the deviation of the numerical data from the observed functional dependency, we now introduce two regression models, which will be fit to the numerical data by a least squares method, and the normalized root-mean-square deviation (NRMSD) will serve as a measure for the difference between the observed values of 〈 〉 R and the values implied by the regression models. The NRMSD is defined as the square root of the mean square error (MSE) normalized by the range of observed values where the sum in (18) is taken over all data points of a simulation setup, and ⋆ R g D w ( , , ) n n ( ) ( ) is the model output for a given g n ( ) and D n ( ) . The optimal model parameters are given by the vector w, which is retrieved by a gradient descent method.
The first regression model is of rather low complexity, based only on the observation that the dependency of the coupling strength on the noise level for fixed values of 〈 〉 R is close to Table 1. Division of all numerical simulations. The first column states the model of the simulations, followed by the underlying network, the coupling method and the number of simulations. The fifth column quotes the normalized root-mean-square deviation (NRMSD) of a fit to the linear model ⋆ R g D lin w ( , , ) , l as described in equation (20), and the last column quotes the NRMSD of a fit to the non-linear model ⋆ R g D nonlin w ( , , ) , nl as described in equation (21). In the bottom row, the total number of simulations is stated, followed by the mean NRMSD over all simulation setups (±one standard deviation) for the linear model and the non-linear model. are the model parameters, which again were fit separately for each simulation setup.
In table 1, the NRMSDs for both regression models and all simulation setups are summarized. The average NRMSD of all setups for the non-linear regression model has, as expected, a very low value of = ± NRMSD 4.8% 1.3% nonlin and is therefore capable of fitting the data very well. With a value of = ± NRMSD 7.6% 1.8% lin , the average NRMSD for the linear model is about 1.6 times larger, but considering the simplicity of the model, its capability of fitting the data is surprisingly good.

Discussion
In summary, we performed numerical simulations on a variety of oscillator models in a broad spectrum of network topologies, coupled by both linear and non-linear coupling functions. The models range from the Fitzhugh-Nagumo model in a continuously spiking state, over the Izhikevich model in a periodically bursting state to the Kuramoto model of identical phase oscillators. The network topologies include scale-free, small-world, regular, random and highly modular networks, with both directed and undirected edges.
We were interested in the maximum degree of synchronization, in dependence on the global coupling strength g and the intensity of the white noise sources D. We found common characteristics independent of the oscillator model, network and coupling type. For fixed values of the noise intensity, we found a sigmoidal dependency of the synchronizability on the global coupling strength, and for fixed values of the global coupling strength, the dependency of the synchronizability on the noise intensity seems to be sigmoidal too. Furthermore, the scaling of the noise intensity with the global coupling strength for fixed values of the average synchronization of a network seems to be close to linear.
Introducing allowed us to quantify the deviation between the numerical results and the proposed sigmoidal dependency for each simulation setup in terms of the normalized root-mean-square deviation (NRMSD), as stated in equation (19). Given the simplicity of the regression model, and the diversity of oscillator models, networks and coupling methods, the consistently small NRMSDs ( = ± NRMSD 7.6% 1.8% lin ; see , capable of reproducing a non-linear relationship between noise intensity and coupling strength for fixed values of the mean correlation as well as a rising slope of that relation for increasing values of the mean correlation, only diminishes the average NRMSD to = ± NRMSD 4.8% 1.3% nonlin . We have thus shown that the maximum degree of synchronization can be approximated quite adequately by a simple 2-dimensional sigmoidal function with a linear relation between noise intensity and global coupling strength for a given synchronizability, reminiscent of the known linear relation between the noise intensity and the critical coupling strength ϵ λ = + D 2( ) c (as in our study we consider identical oscillators, we have λ = 0) for the mean-field Kuramoto model in the thermodynamic limit [3].
It is left to future work to assess whether this approximation can be derived analytically, at least for the mean-field Kuramoto model, and if it holds true in the case of non-identical oscillator models. Other issues to explore would be the influence of the oscillator model, the coupling method and the network topology on the parameters of the regression model.