Topologically determined optimal stochastic resonance responses of spatially embedded networks

We have analyzed the stochastic resonance phenomenon on spatial networks of bistable and excitable oscillators, which are connected according to their location and the amplitude of external forcing. By smoothly altering the network topology from a scale-free (SF) network with dominating long-range connections to a network where principally only adjacent oscillators are connected, we reveal that besides an optimal noise intensity, there is also a most favorable interaction topology at which the best correlation between the response of the network and the imposed weak external forcing is achieved. For various distributions of the amplitudes of external forcing, the optimal topology is always found in the intermediate regime between the highly heterogeneous SF network and the strong geometric regime. Our findings thus indicate that a suitable number of hubs and with that an optimal ratio between short- and long-range connections is necessary in order to obtain the best global response of a spatial network. Furthermore, we link the existence of the optimal interaction topology to a critical point indicating the transition from a long-range interactions-dominated network to a more lattice-like network structure.


Introduction
In the last decade, the theory of complex networks has become a hot topic in social, biological and information sciences [1,2]. The so-called small-world [1] and scale-free (SF) [3] networks have been found within such diverse structures as food webs [4], computer networks [5], social networks [6], cellular and metabolic networks [7] and in the nervous system [8,9]. Such networks are characterized by some exceptional properties. For small-world networks, it is known that on the one hand, the average path between vertices is relatively short, like in a random network, but on the other hand, the clustering coefficient indicating the cliquishness of a typical neighborhood is large, like in a regular graph. Furthermore, in many complex networks, the connectivity of the vertices is distributed according to power law, thus indicating the SF character of the network. The most prominent feature of the SF networks is the presence of a few highly connected vertices, i.e. hubs, whereas there are many nodes with a few connections.
In order to understand the emergence of the SF networks, the mechanism of preferential attachment was introduced [3], in which new vertices connect preferentially to existing wellconnected ones. As an alternative, the varying fitness model was introduced [10]- [12], where fitness values are assigned to individual vertices and links are established between them with a probability depending on their fitness values. Thus, the vertex fitness mechanism leads to networks whose topological properties are determined by the characteristics of the vertices. In particular, it was shown that the SF distributions and even non-SF distributions of fitness values lead to networks with power-law degree distributions [12]. Nevertheless, many real networks are embedded in metric space and the establishment of links is also a function of the Euclidean distance between the vertices [13]- [15]. The most prominent examples of such genuine networks are the Internet [13], mobile communication [16], traffic networks, power grids and the nervous system [9].
Stochastic resonance (SR) is a phenomenon in which the response of a nonlinear system to an external sub-threshold signal is optimized with the addition of noise (for a review see [17]). Even after almost three decades after its discovery, this contradictory observation is still a vibrant avenue of research for several scientific disciplines ranging from physics [18] to biology [19] and chemistry [20]. In parallel with studies of isolated systems with relatively small numbers of degrees of freedom, the research focus has been shifting to coupled and spatially 3 extended systems with many degrees of freedom. Those studies revealed numerous examples of nontrivial collective dynamic behavior in various setups (for a review see [21]), ranging from noise-induced propagation [22], spatiotemporal SR [23] and spatial coherence resonance [24] to array-enhanced SR [25], system size SR [26], diversity-induced resonance [27], etc.
With the advances in network science, the frontier of interest has shifted from regular networks to ensembles with complex interaction topologies. Studies on the SR phenomenon on networks with a small-world topology in general emphasize the advantages of added shortcut links reflected by an enhanced collective response and a greater synchronization [28]. Similar investigations incorporating the SF network models also emphasize the advantages of complex interaction topologies [29]. Nevertheless, in contrast to previous investigations, where mostly topological properties of networks in connection with the SR phenomenon have been studied, we also take into account spatial aspects of the networks. For this purpose, we make use of the spatially embedded vertex fitness network model [13,15], where a fitness value is assigned to each element of the network and vertices are connected with the probability depending on the vertex fitness and their spatial localization. With a control parameter we can then adjust the relative influence of those two factors on the connectivity of individual vertices. In this manner, the network topology can smoothly be altered between an SF network with dominating longrange connections originating mostly from a few hubs and a network where, in general, only adjacent oscillators are connected. Our results indicate that besides an optimal noise intensity, there is also a best possible network configuration at which the greatest signal amplification via SR is achieved. Most remarkably, for different distributions of excitability, the optimal network topology is always between a highly heterogeneous network with mainly long-range connections and a strong geometric regime. We treat two types of oscillators constituting the network, one bistable and the other excitable, thus indicating that our findings are independent of the inherent dynamics driving the nodes. Moreover, we can also identify the value of the network control parameter giving the optimal network structure.

The core models of coupled oscillators
For the sake of generality, we consider two prototypical nonlinear systems: one bistable and the other excitable. First we consider an ensemble of coupled bistable overdamped oscillators, whose dynamics is driven by Langevin equations of the form where dynamical variables x i describe the state of the ith oscillator located at the ith site, D defines the strength of Gaussian noise with zero mean and autocorrelation ξ i (t)ξ j (t ) = δ i j δ(t − t ), the sum signifies the coupling and it runs over all units, whereby ε i j = 1 if the unit i is coupled to unit j, while otherwise ε i j = 0. g is the coupling constant and A i is the amplitude of the external periodic forcing with frequency ω. As the second system, we examine coupled excitable elements. In particular, we make use of a two-dimensional (2D) iterated map, the so-called 'Rulkov map', describing neuronal dynamics as proposed by Rulkov [30], where u i (n) and v i (n) are the slow and fast dynamical variables, respectively, and are considered as dimensionless variables, n is the discrete time index and α, χ and σ are the system's parameters. By choosing χ = σ = 0.001 1, we ensure that the dynamics of u i (n) is much faster than that of v i (n). Furthermore, by setting α = 1.95 < 2, we ensure that the system exhibits a single excitable steady state. In the case of minor deviations from the excitable steady state, the system remains quiescent. Only substantial perturbations push the system into the firing state, which is characterized by large-amplitude spikes. Although the dynamics of excitable systems is qualitatively different from the one characterizing the bistable oscillator, both setups have proven to be suitable candidates for the realization of the SR.

External forcing
The external periodic forcing is of the so-called sub-threshold type, which means that in the absence of noise the external driving can neither provoke hopping from one stable state to the other in the case of bistable oscillators, nor trigger large-amplitude excitations in the case of excitable oscillators. Moreover, the sub-threshold amplitude of the external forcing A i is assumed to follow a power-law distribution Here the values of A i are assigned deterministically as shown below [15], where i = 1, 2, . . . , N and N signifies the number of units. In order to ensure that all oscillators are subjected to the sub-threshold driving, we bound values A i by setting the maximum allowed value A i (max), so that in both cases none of the units can provoke a response without the support of noisy perturbations, whereby the characteristics of the distribution remain unchanged.

Constructing the network of coupled oscillators
The generation of an interaction network of coupled oscillators is based on the vertex fitness network model [11] and the spatially embedded vertex fitness network model [15]. The amplitudes of external forcing A i are used for characterizing the fitness of each vertex and the network connectivity is determined as shown below, where l i, j signifies the Euclidean distance between the ith and jth oscillators in a d-dimensional space and δ is the control parameter that characterizes the topology of the network. In particular, for δ 1, the connectivity is predominantly dependent on the fitness of the vertices and consequently the long-range interactions dominate in the network structure. On the other hand, for δ 1, the distances between oscillators become the key constraint that impacts the topology of the network, so that mainly short-range interactions exist between the oscillators. The ith and jth cells are connected if the connectivity exceeds the threshold value i→ j > . Given a fixed set of N units, we can set the number of links (or mean vertex degree k ) in the network by choosing the appropriate . Note that the model given by equation (5) is similar to the probability function employed by Morita [15], where for simplicity reasons the distance l i, j was defined by the L max norm (supremum norm), while in our calculations, the usual Euclidean L 2 norm is used. In figure 1  In order to quantify the network properties, we calculate the corresponding degree distributions. As shown in the lower row of figure 1, the degree distribution for δ = 0.4 indicates an SF structure of the network, and for δ = 5, the degree distribution obeys a Poisson distribution, as found in random geometric graphs. For intermediate values of δ, the degree distribution is smoothly altered between these two extremes. In particular, for δ = 2, the fraction of long-range connections is, in comparison with the case where δ = 0.4, obviously decreased. However, there still exist some well-connected vertices, i.e. hubs, but their extent is declining due to the increased number of connections between the neighboring units. As δ is further increased (δ = 3.8), short-range interactions begin to dominate in the network and the degree distribution becomes more and more similar to the Poisson distribution, as typical for δ 1. Note that for other distributions of the fitnesses of vertices, qualitatively similar networks appear as a function of δ, e.g. for δ close to 0, we obtain a network with a few hubs, which are connected with the majority of the nodes via long-range links. However, in this case, the degree distribution does not exactly follow the power law. Namely, according to Servedio et al [12], an SF network can be produced by arbitrary distributions of vertex fitnesses, but in this case one has to find a suitable linking probability function, which is in general different from those used in the present paper, even at δ = 0. For that reason and because of the investigation of the role of the scaling exponent β in the continuation of the paper, we use a power-law distribution of forcing amplitudes.

Evaluating the response of the network to the external forcing
In order to quantify the global response of the system to the external forcing, we calculate the Fourier coefficient Q of the mean field, where the coefficients are defined as Furthermore, we also define the local Fourier coefficients Q i as  (7) and (8), n defines the number of oscillation periods used. Note that s i in equations (6) and (8) denotes x i for the case of bistable oscillators, whereas it signifies u i when excitable oscillators are analyzed. Since the Fourier coefficients characterize the amplitude of the input frequency ω in the output signal, they represent a suitable and commonly used measure for the quantification of the SR in the system. Furthermore, in order to examine the correlation between the individual oscillators in the system, we calculate the correlation matrix R, whose ijth element is defined as where s i has the same meaning as in equations (8), ands i and σ s i are the mean value and the standard deviation of the time series s i (t) of the ith oscillator, respectively. The correlation coefficient R i j thus indicates the relationship between the dynamics of the ith and the jth oscillator.
The employed system parameters used throughout this study were g = 0.1, ω = 0.01, A i (max) = 0.08 for the network system of bistable oscillators and g = 0.005, ω = 0.004, A i (max) = 0.012 for the network of excitable oscillators (Rulkov map). The number of nodes was N = 100, unless otherwise stated, and the average degree of the network was set to k = 5. All calculations were averaged at least over 50 independent realizations, in order to diminish statistical fluctuations originating from the stochastic dynamics and from randomness incorporated into network generations.

Results
We analyze the collective response of coupled bistable and excitable oscillators as a function of the network structure. For this purpose, we calculate the Fourier coefficient Q for different noise intensities D and different network configurations characterized by different δ. The results are presented in figure 2: the upper row for the ensemble of bistable oscillators and the lower row for the excitable system; in both cases three different values of the scaling exponent β were considered. It can be observed that besides an optimal noise intensity, there is also a best possible network configuration for which the system's collective response is at best correlated with weak external driving. The comparison between the upper and lower rows of figure 2 indicates that this observation does not depend on the inherent dynamics of elements constituting the network. Furthermore, for different values of β, qualitatively similar results are obtained, since in all cases an intermediate network configuration exists at which the greatest synchrony between the forcing and the response of the system is achieved. However, for higher values of β, the optimal network structure is attained at lower values of δ.
Next we investigate the role of the system size expressed by different values of N . Similar to figure 2 (N = 100), in figure 3 Q as a function of D and δ is shown for N = 200 (left column) and N = 300 (right column) and for both the cases of the inherent dynamics of units constituting the network. Evidently, in all cases the optimal responses coincide with those inferred by N = 100 (see the left column of figure 2), which confirms that the reported phenomenon is qualitatively independent of the system size. Importantly, however, the overall values of Q decrease as the system size is increased, which is a consequence of the lower average forcing amplitudes to which all oscillators are subjected. This reflects the fact that on the one hand the values A i are arranged in descending order of i, as given by equation (4), and on the other hand, we are obliged to fix the maximal amplitude A i (max) in order to ensure that all oscillators are subjected to sub-threshold forcing. Therefore, lower amplitudes of periodic excitations result in a decreased resonant response of individual oscillators and hence a decreased collective response Q is attained as N is increased.
To gain more insights into the observed phenomenon, we show in figure 4 the responses of individual units Q i as a function of δ at optimal noise intensity. In the left panel, the results for the case of bistable oscillators are shown, whereas in the right panel the same analysis for coupled excitable elements is performed. In both cases, we use β = 2 as a representative example, but analogous interpretations are attained irrespective of the scaling exponent characterizing the distribution of driving amplitudes. At this point it should be emphasized that the amplitudes of periodic forcing are decreasing with increasing i, as given by equation (4), so that accordingly the responses of individual units are expected to decrease with increasing i as well. However, the decrease in the resonant response reflected by values of Q i obviously depends on the network structure. There exists a most favorable δ, at which the resonant responses of individual units decrease the most slowly with increasing i. These results corroborate the findings already described in figure 2 that indeed an optimal global response is achieved in the intermediate regime between an SF network and a network with only short-range interactions.
In order to explore how the existence of an optimal topologically determined SR response is related to dynamical correlation among oscillators, we calculate the correlation matrices R for different network configurations. To provide insights into the relation between the topological and dynamical features of the network, we show in figure 5 characteristic connectivity matrices ε along with the correlation matrices R for both cases of the inherent dynamics of units constituting the network. The optimal noise intensity and three different network configurations are considered: a network with dominating long-range connections and a power-law degree distribution (δ = 0.4); a random geometric network (δ = 5); and a network whose topological properties are between those two extremes (δ = 2) and that was found to ensure the best synchrony between the external driving and the mean field. Evidently, for intermediate values of δ, the greatest overall correlation between oscillators is achieved, which once again confirms that the optimal response is attained if both long-and short-range interactions constitute the structure of the network.
For the quantification of this visual assessment we further calculate the average correlation coefficient The results shown in figure 6 indicate that indeed the best average correlation is achieved for intermediate values of δ and, in addition, the maximum roughly coincidences with the optimal δ at which also the greatest signal-to-noise ratio was attained. Apparently, the optimal SR response can be associated with the network topology, which ensures the best transmission of weak forcing across the network.
A more advanced analysis of this optimal global response of the system being realized at a given threshold, at a critical (optimal) δ, in a phase-transition manner, is presented in the appendix.

Summary
In this paper, we have studied the SR phenomenon in spatial networks, where the connectivity is a function of the location of oscillators and the amplitude of external forcing. By varying a control parameter, the network topology was smoothly altered between an SF network with a high inherent degree of inhomogeneity of individual oscillators and a random geometric network where no substantial discrepancies in the degrees of individual nodes are found. Our results show that besides an optimal noise intensity, there is also a most favorable interaction topology at which the best correlation between the response of the network and the imposed rhythm is achieved. For various distributions of amplitudes of external forcing, an optimal topology was always found in the intermediate regime, where long-and short-range connections exist between the oscillators. Here, units that are subjected to the highest forcing amplitudes act as local hubs that govern the dynamics of interconnected neighborhood, and long-range links primarily represent connections between those hubs (see the second panel from left in figure 1). For comparison, in the SF network (left panel in figure 1), the majority of the units are connected exclusively to a few hubs, which operate as global pacemakers with a very large node degree, and thereby influence very different parts of the network. It seems that the entraining of the network is less efficient if there are large discrepancies in the node degrees. On the other hand, in a network with mainly short-range interactions (see the right column of figure 1), there are no significant differences in node degrees, and the leading oscillators remain unexpressed and hence fail to efficiently dictate the rhythm of the network. Apparently, a suitable number of hubs and an optimal ratio between short-and long-range connections are necessary for ensuring the greatest global response of the whole network. The results presented here are also in agreement with our previous study, where the influence of interaction topology on the regularity of noiseinduced oscillations in an ensemble of excitable oscillators has been analyzed [34]. These results have revealed that besides an intermediate noise intensity, there is also a best possible network configuration at which the noise-induced oscillations are most regular. The reported optimal network architecture approximately coincides with those in the present paper, thus indicating that a most favorable topology exists for various dynamical processes in spatially embedded complex networks. Furthermore, from the results presented in the appendix, we observed that the critical point of the oscillators dynamics giving the optimal global response coincides with δ c = d/(β − 1), the critical point of the network structure, signaling the change in network organization from a long-range dominated (δ < δ c ) to a more compact one (δ > δ c ). We hope that our results presented here will help in further clarifying the impact of the spatial organization of coupled oscillators on their global dynamic response.  between nodes, we find [33] From this, we can compute the mean of C, where (d) is the d-dimensional solid angle for example 2π in 2D space. It is interesting to note that when the parameter δ approaches the critical value δ c = d/(β − 1), the cost of links C drops to zero. Expanding C in series to the first order as δ → δ c gives C ∝ (δ − δ c ) 1/(β−1) .
The critical point δ c = d/(β − 1) denotes the transition from the long-range dominated network to a more lattice-like network structure. By plotting the links length fluctuation against δ(β − 1), we observe that the curves obtained with different exponents β collapse to a single one ( figure A.1, right panel). We further explore the dependence of the resonant responses on δ(β − 1) for different values of β and for both types of oscillators. For a better comparison, we normalize the Fourier coefficients Q, so that they are bounded in the unit interval and denote it as Q norm . The results showing the normalized resonant responses at optimal noise intensity are presented in figure A.2. It can be observed that for all values of β, the curves reach their maximums at roughly the same value of (β − 1) δ = d = 2. This finding indicates that indeed the optimal global response is realized at the critical point δ c .