The effect of noise on the synchronization dynamics of the Kuramoto model on a large human connectome graph

We have extended the study of the Kuramoto model with additive Gaussian noise running on the KKI-18 large human connectome graph. We determined the dynamical behavior of this model by solving it numerically in an assumed homeostatic state, below the synchronization crossover point we determined previously. The de-synchronization duration distributions exhibit power-law tails, characterized by the exponent in the range $1.1<\tau_t<2$, overlapping the in vivo human brain activity experiments by Palva et al. We show that these scaling results remain valid, by a transformation of the ultra-slow eigen-frequencies to Gaussian with unit variance. We also compare the connectome results with those, obtained on a regular cube with $N=10^6$ nodes, related to the embedding space, and show that the quenched internal frequencies themselves can cause frustrated synchronization scaling in an extended coupling space.


Introduction
The organization of resting-state activity, i.e. the dynamics of the brain under the absence of external stimulation and no task condition, plays putatively a critical functional role given the fact that its maintenance requires a large part of the Email address: odor@mfa.kfki.hu (GézaÓdor) total brains energy budget [1,2]. There are nowadays empirical and computational evidences showing that the resting organization facilitates task-based information processing [3]. Resting brain networks as captured by Functional Connectivity (FC) maps very consistently show task-evoked activity such that individual differences in FC can predict individual differences in task-evoked regional activity [4,5,6]. From a mechanistic perspective, whole-brain models were able to demonstrate that resting-state organization conforms to a state of criticality that promotes responsiveness to external stimulation, i.e. resting-state organization facilitates task-based processing [7,8,9].
Neural activity avalanche measurements found size and duration distributions that can be fitted by power-laws before a size cutoff, which can arise naturally close to a critical point of a second order phase transition [10,11,12,13,14]. Criticality hypothesis has been advanced, because information processing, sensitivity, long-range and memory capacity is optimal in the neighborhood of criticality [15,16]. Criticality in models can be defined by the diverging correlation volume, as we tune a control parameter to a threshold value.
It has been debated how a neural system is tuned to criticality. At first selfregulatory mechanisms [17], leading to self-organized criticality [18] were proposed. Recently, it has been shown that as the consequence of heterogeneity extended dynamical critical regions emerge in spreading models [19,20] naturally. As real systems are mostly inhomogeneous and one must asses whether heterogeneity is weak enough to use homogeneous models for describing them. Heterogeneity can create rare-region effects of different relevancy [21]. They can generate so-called Griffiths Phases [22] in which scale-free dynamics appears over an extended region around the critical point with slowly decaying auto-correlations and burstyness [23]. This phenomenon was proposed to be the reason for the working memory in the brain [24]. Furthermore, in GP the susceptibility is infinite for an entire range of control parameters near the critical point, providing a high sensitivity to stimuli, beneficial for information processing.
As individual neurons emit periodic signals [25] it is natural to expect criticality in oscillator models at the synchronization transition point. Very recently analysis of Ginzburg-Landau type equations arrived at the conclusion that empirically reported scale-invariant avalanches can possibly arise if the cortex is operated at the edge of a synchronization phase transition, where neuronal avalanches and incipient oscillations coexist [26]. Several oscillator models have been used in biology, the simplest possible one is the Hopf model [27], which has been used frequently in neuroscience, as it can describe a critical point with scale-free avalanches, with sharpened frequency response and enhanced input sensitivity.
Indeed, Deco and colleagues developed a mesoscopic whole-brain model based on the Hopf model, which provides an excellent account of resting-state empirical FMRI [28] and MEG data [29]. Furthermore, the Hopf whole-brain model is able to be used in conjunction with higher-resolution functional parcellations that will increase model accuracy. The model consists of coupled dynamical units representing the cortical and sub-cortical brain areas from a given parcellation. The local dynamics of each brain area (node) is described by the normal form of a supercritical Hopf bifurcation, also called a Landau-Stuart Oscillator, which is the canonical model for studying the transition from noisy to oscillatory dynamics. The emerging global whole-brain dynamics results from the partial meta-stable entrainment of different clusters of brain areas that synchronize and build up different network micro-states.
Another complex model, describing more non-linearity 1 is the Kuramoto model [31,32], with phases θ i (t), located at N nodes of networks, according to the dynamical equationθ The global coupling K is the control parameter of this model, by which we can tune the system between asynchronous and synchronous states. The summation is performed over other nodes, with connections described by the weighted adjacency matrix W ij and ω i denotes the intrinsic frequency of the i-th oscillator. For simplicity we used for the g(ω i ) distributions Gaussian functions [33]. Earlier Eq. (1) was studied analytically and computationally on a human connectome graph network of 998 nodes and in hierarchical modular networks (HMN), in which moduli exist within moduli in a nested way at various scales [34]. As the consequence of quenched, purely topological heterogeneity an intermediate phase, located between the standard synchronous and asynchronous phases was found, showing "frustrated synchronization", meta-stability, and chimera-like states [35]. This complex phase was investigated further in the presence of noise [36] and on a simplicial complex model of manifolds with finite and tunable spectral dimension [37] as simple models for the brain.
We continued to investigate Eq. (1) on a large, weighted human connectome network, containing 804 092 nodes, in an assumed homeostatic state [33]. Homeostasis in real brains occurs via inhibitory neurons [38,39,8,40,41], here we modeled this by normalizing the incoming interaction strengths [42]. Recent experiments arrived at a similar conclusion: equalized network sensitivity allows critical behavior and produces model results, which reproduce measured FMRI correlations [43].
Since this graph has a topological dimension d < 4 [44], a real synchronization phase transition is not possible in the thermodynamic limit, still we could locate a transition between partially synchronized and desynchronized states. At this crossover point we observe power-law-tailed synchronization durations, with τ t ≃ 1.2(1), away from experimental values for the brain. Below the transition of the connectome we found global coupling control-parameter dependent exponents 1 < τ t ≤ 2, overlapping with the range of human brain experiments [14]. 2 Note however, that in case of the Kuramoto model the so-called spectral dimension is a more relevant factor influencing the scaling behavior [37].

Materials and methods
Now we extend this study by an additional, annealed zero centered Gaussian distributed noise term ξ(i) with unit variance and a coupling strength ṡ The effect of noise on synchronization of Kuramoto oscillators has been investigated by many works (see for example [45]). It was shown, that in high dimensional, meanfield models the noise can be neglected, it shifts the critical coupling rate K c only, but does not change the scaling behavior. Earlier we found the KKI-18 graph finite dimensional [44] and very heterogeneous, exhibiting a hierarchical modular topology. In low dimensional heterogeneous systems a large repertoire of attractors can exist, causing meta-stable states with different degrees of coherence and stability. Noise enables the system to jump into another close, more stable, attractor. As neurons work in a noisy background this study can be interesting for neuroscience. We follow the dynamical behavior of the system through studying the Kuramoto order parameter defined by 2 The topological (also called graph) dimension is defined by where N r is the number of node pairs that are at a topological (also called "chemical") distance r from each other (i.e. a signal must traverse at least r edges to travel from one node to the other).
which is finite, above a critical coupling strength K > K c , or tends to O(N −1/2 ) for K < K c . At K c , in case of an incoherent initial state it evolves as characterized by the dynamical exponents:z, η and f ↑ denotes a scaling function. Within the framework of the synchronization model we identified a spontaneous avalanche start times at t = 0 of the fully desynchronized initial condition and the end times t x , when R(t) returned back to 1/ (N), related to the synchronization value of independent oscillators. Note, that by changing the start times to the first up crossing value did not change the tail behavior we consider here. To obtain the avalanche duration probability distributions P (t x ) we performed ≃ 10 4 runs, with independent random ω i intrinsic frequencies and applied histogramming method with increasing bin sizes: ∆t x ∝ t 1.12 x . The following graphs have been considered: 1. 3D lattice with linear size L = 100 and periodic boundary conditions. 2. Weighted, symmetric large human connectome graph: KKI-18 [44] downloaded from the Open-connectome project [46].
To integrate the differential equation (1) we used a Graphics card (GPU) code, based on the fourth order Runge-Kutta method of Numerical Recipes [47] and the boost library odeint [48] on various networks.
Step sizes: ∆ = 0.1, 0.01, 0.001 have been tested and finally the ∆ = 0.01 precision found to be sufficient, for s ≤ 2. For larger noise amplitudes s this precision was not enough. In case of our initial attempt, with natural ω i < 0.02 we needed large s values to see a synchronization transition and even ∆ < 0.01 turned out to be insufficient, making the numerical analysis prohibitively time-consuming for available computer resources.
In general, the ∆ < 0.01 precision did not improve the stability of the solutions further, but caused smaller fluctuations due to the chaotic behavior of Eq. (1), which could be compensated by averages over many independent samples with different ω i . We used the ǫ = 10 −12 criterion in the Runge-Kutta algorithm and parallelized the solver for NVIDIA graphic cards, by which we could achieve a ∼ ×40 increase in the throughput on Tesla V100 GPU-s with respect to a single 12-core Intel Xeon Gold 6136 CPU. Algorithmic and benchmark details will be discussed elsewhere [49].
We determined the Kuramoto order parameter at a fixed K, by increasing the sampling time steps exponentially which is a common method at critical systems, where we expect power-law (PL) asymptotic time dependences. We estimated t x = (t k +t k−1 )/2, where t k was the first measured down crossing time. The initial conditions were generally θ i (0) ∈ (0, 2π] phases, with uniform distribution, describing fully disordered states. The probability distribution tails were fitted using the least squares fit method beyond a time, fixed by visual inspection of the results. To visualize corrections to PL scaling we determined the effective exponents of R as the discretized, logarithmic derivative of Eq. (5) were the brackets denote sample averaging over different initial conditions. We obtained the KKI-18 graph from the Open Connectome project repository [46]. This is based on the diffusion tensor image [50], approximating the structural connectivity of the white matter of a human brain. The graph version we downloaded in 2015 comprises a large component with N = 804 092 nodes, connected via 41 523 908 undirected edges and several small sub-components, which were ignored here 3 . This graph allowed us to run extensive dynamical studies on present day CPU/GPU clusters, large enough to draw conclusions on the scaling behavior without too strong finite size effects, hindering scaling regions. The large connectomes from [46] of the human brain possess 1 mm 3 resolution, obtained by a combination of diffusion weighted, functional and structural magnetic resonance imaging scans. These are symmetric, weighted networks, where the weights measure the number of fiber tracts between nodes. The KKI-18 graph is generated via the MIGRAINE method, described in [51]. It exhibits a hierarchical modular structure by the construction from the Desikan cerebral regions with (at least) two quite different scales. The graph topology is displayed on Fig. 1, in which modules were identified by the Leiden algorithm [52], and the network of modules generated and visualized using the pythonigraph library [53]. This identified 153 modules, with sizes varying between 7 and 35 332 nodes in the displayed case, however since this is a heuristic approach, these numbers vary by about 10 %. A recent experimental study has provided confirmation for the connectome generation used here [54]. This suggests that diffusion MRI tractography is a powerful tool for exploring the structural connection architecture of the brain.
In [44] it was found that, contrary to the small world network coefficients, these graphs exhibit topological dimension D = 3.05 [44], slightly above the embedding space and a certain amount of universality, supporting the selection of KKI-18 as a representative of the large human connectomes available.
To keep a local sustained activity requirement for the brain [55] and to provide a homeostatic state, we modified KKI-18 by normalizing the incoming weights of node i in [42]: W ′ i,j = W i,j / j∈neighb.of i W i,j at the beginning of the simulations. In reality such local homeostasis is the consequence of the competition of exhibitory and inhibitory neurons.

The 3D lattice
In [33] we compared the connectome results with small-world graphs, generated from 2D lattices with additional random long-range links (2Dll). This was done as both the connectome and the 2Dll graphs exhibit small world properties and quenched topological disorder on top of the intrinsic heterogeneity. As the topological dimension of the connectome is slightly above the D = 3 embedding space [44], due to the long connections, now we have studied the growth of R(t) on the 3D lattice of linear size L = 100 by starting from states of oscillators with fully random phases and by averaging over 5000 − 10 000 internal frequency realizations up to t = 10 3 time steps. In this case we can test the effect of quenched heterogeneity of the ω i -s on the dynamical behavior, using a zero centered, unit variance g(ω) distribution. growth changes from convex to concave at K c ≃ 2.2, suggesting a crossover point there. This value is higher than what we obtained for the KKI-18 connectome [33]: K c = 1.7, as we have much lower connectivity now. The KKI-18 graph has an average node degree: k = 156, in contrast with the cubic lattice: k = 6, thus we need stronger global coupling to achieve synchronization. The initial behavior of growth at K c = 2.2 is like R(t) ∝ t 0.5 , with an exponent somewhat smaller, but close to the estimates obtained for the 2Dll and for the KKI-18 graphs: η = 0.6(1). This can be read-off from the local slopes before finite-size cutoff, shown in the inset of As for the connectome we determined the first crossing times t x , when R(t) of single realizations first fell below the threshold value 1/ (N) = 0.001. Following the usual histogramming procedure we obtained the de-synchronization distributions, which can be considered as avalanches, induced by a spontaneous synchronization and a relaxation process. The tail of the relaxation time distribution PDF at K c = 2.2 decays as ≃ t −1.33 (7) x (see Fig. 3), obtained by least squares power-law regression, applied for the tail region t x > 30. Below K c the curves shown on the figure provide good agreement with non-universal PL-s. However, for K < 1.7 the PL region seems to shrink. The ineffectiveness of the weak additive noise has also been demonstrated   in case K = 2.0 coupling, using a Gaussian annealed term of strength s = 1.
The emergence of frustrated synchronization by ω i heterogeneity is similar to the contact process with site disorder [21], where Griffiths Phase has been found in various spatial dimensions. To provide a counter-example we have also studied the 3D case with uniform oscillators and annealed noise. In the Appendix we show that in this case no extended scaling region, but a single critical point emerges.

The Connectome graph Mapping ultra-slow self-frequency scales onto computationally feasible oscillations
FMRI measurements [56,28] found that in the human brain global phase synchrony of the BOLD signals evolves on a characteristic ultra-slow: < 0.01 Hz time scale. In modeling this, on smaller sized connectomes, using the noisy Kuramoto equation [56] and the Hopf model [28], the intrinsic frequencies were filtered in the 0.04 < ω < 0.07 Hz band. In both cases the temporal and spatial synchronization patterns were found to approximate well the empirical data. Here we extend this modeling for a large human connectome, which allows to test possible scaling behavior, with ω i -s of narrow spread (σ = 0.02) and typical mean value ω = ω i = 0.05.
First we show, that in the case of the Kuramoto model the narrow frequency band can be mapped onto the much wider σ = 1, ω = 0 Gaussian case, where we obtained recent results [33] on synchronization durations. The invariance over the mean frequency is obvious, as (1) is invariant to the global shift of a mean rotation frame ω → ω ′ . The oscillation-size dependence can also be gauged out by the following transformation: ω i → aω ′ i , t → (1/a)t ′ and K → aK ′ . Therefore, for small a values corresponding to empirical data, we can obtain the same results as for σ = 1 at rescaled, late times and using small global couplings. From a technical point of view this is very important, because we can recover the asymptotic results at ultra-slow frequencies using much shorter time scales in our simulations. Note, that in the case of the Hopf model the rotating frame transformation is also possible, but the width of the g(ω) distribution cannot be gauged out and acts as an independent parameter, causing a nontrivial phase structure of the oscillations [57]. Furthermore, a completely band filtered g(ω i ) spectrum with a flat top can be transformed onto the uniformly distributed natural frequencies. In this case it is known, that in the limit of infinite system size the Kuramoto model in the steady state undergoes a first-order phase transition [58]. At first-order transitions we don't expect critical scaling behavior and indeed our simulations for this case do not show PL de-synchronization duration tails (see Appendix).

Noisy Kuramoto results
Having taken into account the transformation properties of the previous section we investigated the effect of annealed noise, by adding a time varying Gaussian distributed random numbers to the (1). Earlier in [33] we concluded that using step sizes ∆ < 0.1 in the Runge-Kutta-4 solver did not modify the results. Here we test the effect of ∆ again, as we add a stochastic noise on top of the chaoticity, inherent in the noiseless Kuramoto equation.
First we applied s = 1 noise to the KKI-18 homeostatic graph case [33]. As Fig. 4 shows slightly below the synchronization transition point: K = 1.4 < K c = 1.7 the Kuramoto order parameter growth curves fully agree with the noiseless case [33] one and the applied ∆ = 0.1, 0.01, 0.001 step sizes do not affect too much the growth regime. One can see a smaller saturation value for ∆ = 0.1 than for ∆ = 0.01, 0.001. Here each line corresponds to an average over ∼ 10 4 realizations, in which both the quenched and the annealed noise varies.
We have determined the first crossing times t x , when R fell below: 1/ (N) = 0.001094. Following the histogramming procedure we obtained the distributions p(t x ), which again show PL tails, characterized by exponents: τ t = 1.2(1) (see Fig. 5). The fitted τ t = 1.9(1) exponent is in the range of in vivo human neuro experiments: 1.5 < τ t < 2.4 [14]. Then we repeated the analysis for K = 1.3, s = 1 and found similar agreement with the noiseless results (see figure in the Appendix).
We have tested noise: s = 2 at K = 1.4 as shown on the inset of Fig. 4 and on 6, but again we just found insensitivity in the growth and scaling results. We have learned, that going to even stronger noise amplitudes requires smaller ∆-s to achieve numerical precision independence and it is questionable from neuroscience point of  view if it is worth to study such strong noises. We plan to study this later. In order to see the effect of the modules on the synchronization behavior in greater detail, we also computed the order parameter for each module separately. Fig. 7 shows the steady-state values for each module for different coupling parameters. While normal finite-size scaling can be observed below the transition point, larger modules start to synchronize above. We fitted PL-s for the sub-critical scaling, showing non-universal exponents in an extended K region.

Conclusion and discussion
Resting-state activity of the brain can be modeled by simple whole-brain models, exhibiting critical behavior at the edge of the synchronization transition. We have investigated the dynamical synchronization behavior of the Kuramoto model on a large, weighted human connectome network. The dynamical behavior of heterogeneous Kuramoto, especially for local interactions is a largely unexplored field, according to out knowledge. In case of identical oscillators heterogeneous phase lags or couplings have been shown to result in partial synchronization and stable chimera states [59,60,61,62] Realistic models of the brain, however, require oscillators [63] to be heterogeneous.
In particular, we focused on the effect of additive noise on the de-synchronization duration distributions. In Ref. [33] it was found that below the phase synchronization transition point these distributions exhibit non-universal power-law tails, with exponents overlapping with the empirical activity avalanche duration values for humans in vivo.
We found that weak Gaussian noises with amplitudes not larger than those of the quenched Gaussian self-frequencies do not affect the previous results within numerical precision. This means that time-dependent, thermal like noise does not destroy or alter the dynamical scaling behavior of this model.
Ref. [34] studied the synchronization behavior of the Kuramoto model on hierarchical modular networks with N = 4096 and human connectomes with N = 998 nodes. They found stretched exponential decay tails for the order parameter ρ(t) = 1 − R(t) in the absence of frequency heterogeneity, in agreement with analytic approximations using ω i = 0. They assumed that fixing all intrinsic frequencies to be identical does not decrease generality of the results. We have also tried to fit our p(t x ) results with stretched exponential functions, (see Appendix) but reasonable agreement could be found only for very low couplings, far below K c . So, we conclude that that the quenched disorder in the self-frequencies cause PL tails in the dynamical behavior of chimera-like states at the edge of criticality. These non-universal PL-s resemble to Griffiths Phase effects and the results obtained in case of the second order Kuramoto model for power-grids below the synchronization transition [64]. We have also detected module size scaling of R(t → ∞) below the transition point, with K-dependent exponents, which would be an interesting subject of further study.
We have shown that the empirical results with ultra-slow oscillations can be transformed onto zero mean Gaussian frequencies as the consequence of the Galilean symmetry of the Kuramoto equation. As we used the 4th order Runge-Kutta solver, we carefully confirmed that the step size ∆ = 0.01 is sufficient for the increased precision required by the additive noise. This sensitivity is related to the fast changes of the time dependent, additive noise.
Another point is the positiveness of the g(ω i ) distribution in the brain, as we don't expect neural oscillators 'rotating backwards'. This corresponds to the question of asymmetric distribution of natural frequencies, such that for g(ω i ) = 0 for ω i < 0. It has been shown that in case of uni-modal g(ω i )-s only the first derivative, the flatness of g(ω i ), matters here. Without a flat top, like an asymmetric triangle, one obtains the same universal critical behavior (β = 1/2) as for the original Kuramoto model with zero centered Gaussian [65]. Thus we expect the same dynamical behavior for an asymmetric, truncated Gaussian: g(ω i ) = 0 for ω i < 0 as we found by the symmetric Gaussian.
We have compared the results with those obtained on regular 3D lattices of similar size and found, that the heterogeneity of the self-frequencies already generate the non-universal scaling region, which can be called a frustrated synchronization phase, exhibiting chimeras [35]. At the synchronization transition point we found: τ t = 1.33 (7), slightly higher than in case of the connectome: τ t = 1.2(1) [33] and well below the mean-field value: τ t = 1.6(1) [33], so the topology plays a role mainly via the graph dimension. Determination of these exponent estimates depend strongly on the precise location of K c , which is harder to get in smaller systems, suffering finitesize cutoffs. Therefore the large systems considered here are justified. The local, 3D connected results can also be important for neuroscience, as local cortical connections appear before more distant ones at the creation of the circuit [66]. Thus this study can be interesting for understanding young brains, containing local connections mainly. Furthermore, 3D networks may occur in case of artificial intelligence neural systems.
The dynamical scaling behaviors have been found to be robust, supporting universality even if the Kuramoto model could be considered too simplistic to describe the brain. Note, however that in the weak-coupling limit equivalence of phase-oscillator and integrate-and-fire models was found [30], which may hold for the sub-critical region, where we observed the dynamical scaling. This provides a support for the edge-of-criticality hypothesis of oscillating systems near and below the synchronization transition point.
An interesting continuation of this work would be the study of the effect of phase shifts, caused by the finite signal propagation in the neural network or the introduction of a threshold, as in integrate-and-fire models, although by universality of critical systems we don't expect qualitative change in the scaling behavior.

Appendix
In this Appendix we show results for the synchronization duration distribution of Kuramoto on the KKI-18 model at K = 1.3. As Fig. A.1 shows dependence on the numerical integration step size ∆ is negligible and the results agree with the former noiseless case. A least-squares fit for the tails results in τ t = 2.1(1), which is in the experimentally measured range for activity avalanche durations.
We have also tried to test in the noiseless case whether the P (t x ) distributions would follow stretched-exponential decay instead of PL-s. As Fig. A.2 shows the − ln(P (t x ))) curves do not exhibit straight lines on the log-log plot, except perhaps for K = 1.2 only, which is far from K c . The method has been tested by considering 3D lattices of size L = 100. of uniform oscillators and annealed noise with s = 5 amplitude. In this case we don't have any heterogeneity and an order-disorder phase transition emerges by increasing the coupling, which should belong to the XY criticality in 3D [67]. As one can see on the attached graph (Fig. A.3) we obtain a PL at the phase transition point K = 0.05(1) and fast decays below it. Here the fitted slope is τ t = 1.1 (3). This agrees with the expectation, coming from scaling relations and known data for the critical dynamics of the XY with model-A dynamics in 3D: β = 0.347(1), ν ⊥ = 0.670(1) [68], Z = 2 [69], τ t = 1 + β/(ν ⊥ Z) = 1.25 (2). Note, that for K > 0.05 the distributions become singular, characterized by 1/t x decay, before a finite size cutoff. Here we also show duration distribution results for the Kuramoto model with uniform g(ω i ), which undergoes a first-order phase transition. We do not expect critical scaling behavior and indeed our simulations for this case do not show PL de-synchronization duration tails. The p(t x ) curves saturate, followed by a sharp decay for different K-s as shown on Fig. A.4.

Competing interest
The authors declare that they have no competing interests.

Consent for publication
Not applicable.

Ethics approval and consent to participate
Not applicable.

Funding
The work has been performed under the Project HPC-EUROPA3 (INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action under the H2020 Programme; in particular, G.Ó. gratefully acknowledges the support of Center for Brain and Cognition Theoretical and Computational Group Universitat Pompeu Fabra / ICREA Barcelona and the computer resources and technical support provided by BSC Barcelona. Support from the Hungarian National Research, Development and Innovation Office NKFIH (K128989), the Initiative and Networking Fund of the Helmholtz Association via the W2/W3 Programme (W2/W3-026) and the Helmholtz Excellence Network DCM-MatDNA (ExNet-0028) is acknowledged.

Availability of data and materials
The codes and the graphs used here are available on request from the corresponding author.

Authors' contributions
G.Ó designed, coordinated this research, drafted the manuscript and performed computations and data analysis. J. K. developed GPU codes and participated in data analysis. G. D. contributed to the manuscript, conceivied of the study and participated in research coodination. The authors read and approved the final manuscript.