Critical synchronization dynamics of the Kuramoto model on connectome and small world graphs

TL;DR

This study investigates the synchronization behavior of the Kuramoto model on a large human connectome network, revealing critical-like dynamics and power-law synchronization durations that resemble brain activity patterns.

Contribution

It demonstrates the presence of critical synchronization dynamics on a large connectome network and explores the effects of network modifications, advancing understanding of brain-like criticality in complex networks.

Findings

Identified a crossover transition with power-law synchronization durations.

Observed exponents consistent with human brain experimental data.

Found similar critical behavior in networks with inhibitory interactions.

Abstract

The hypothesis, that cortical dynamics operates near criticality also suggests, that it exhibits universal critical exponents which marks the Kuramoto equation, a fundamental model for synchronization, as a prime candidate for an underlying universal model. Here, we determined the synchronization behavior of this model by solving it numerically on a large, weighted human connectome network, containing 804092 nodes, in an assumed homeostatic state. Since this graph has a topological dimension $d < 4$, 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 $\tau_t \simeq 1.2(1)$, away from experimental values for the brain. For comparison, on a large two-dimensional lattice, having additional random, long-range links, we obtain a mean-field value: $\tau_t \simeq 1.6(1)$. However, below the transition of the connectome we found global coupling control-parameter dependent exponents $1 < \tau_t \le 2$, overlapping with the range of human brain experiments. We also studied the effects of random flipping of a small portion of link weights, mimicking a network with inhibitory interactions, and found similar results. The control-parameter dependent exponent suggests extended dynamical criticality below the transition point.

Figures

Equations

$$⟨N_r⟩\sim r^d,$$

$$\dot{{\theta }_i}(t)={\omega }_{i,0}+K\sum _j{W}_{ij}\sin [{\theta }_j(t)-{\theta }_i(t)]$$

$$R(t)=\frac{1}{N}\left|\sum _{j=1}^N{e}^{i{\theta }_j(t)}\right|,$$

$$R(t,N)=N^{-1/2}{t}^{\eta }{f}_{\uparrow }(t/N^{\tilde{z}}),$$

$$R(t,N)=t^{-\delta }{f}_{\downarrow }(t/N^{\tilde{z}}),$$

$$t_k=1+1.08^k,$$

$${\eta }_{\mathrm{eff}}=\frac{\ln ⟨R(t_{k+3})⟩-\ln ⟨R(t_k)⟩}{\ln (t_{k+3})-\ln (t_k)}.$$

$${\tau }_t=1+\delta =1.5,$$

Full text

Introduction

Understanding the human brain, or in general neural systems is a great challenge of science, in particular the application of models and methods of statistical physics has been developing recently [1]. There are several types of whole brain models, ranging from continuous, integrate-and-fire models [2, 3] to discrete, activity spreading models [4, 5]. All of them are effective ones, trying to describe different features of neural functions measurable by neuroscience experiments. While different versions of integrating fire models are more detailed and use larger parameter space, simple activity spreading models try to capture basic features, like the emergence of power-laws (PL) of quantities via critical behavior [6]. Following experiments, these are the neuron activity avalanche size and duration distribution tails, before finite size cutoff. Criticality in these systems can be defined by the diverging correlation volume, as we tune a control parameter to a threshold value.

Criticality is an attractive hypotheses, because information processing and dynamic range is optimal [7, 8]. Neural activity avalanche measurements have found power-laws, which arise naturally close to a critical point of a phase transition [9, 10, 11, 12, 13]. The question of how a neural system would be tuned to this point has been debated. It was proposed to be by self-regulatory mechanisms [14] leading to self-organized criticality [15], or as the consequence of extended dynamical critical regions in spreading models [16, 17] in Griffiths Phases (GP) [18]. The measured scaling exponents have been found to be close to the mean-field transition values of discrete models [19].

It is also known, that individual neurons emit periodic signals [20], thus criticality may emerge by the collective behavior of oscillators at the phase synchronization transition point. However, not much is known about the dynamics of the synchronization or desynchronization process in these models [21, 22, 23]. Phase synchrony is essential for large-scale integration of information [24, 25], the role of the asynchronous state has remained more elusive [26]. Very recently theoretical analysis of the homogeneous 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 [27].

One of the most fundamental models, showing phase synchronization is the Kuramoto model of interacting oscillators [28]. This is defined on full graphs, corresponding to the mean-field (MF) behavior [29], but as neural systems are not fully connected, we are interested in the phase synchronization transition in extended systems, where oscillators are located at graph points, possessing finite topological dimension $d$. This is defined by

[TABLE]

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).

Phase transition in the Kuramoto model can happen only above the lower critical dimension $d_{l}=4$ [30]. Below $d_{l}=4$ partial synchronization may emerge with a smooth crossover for strong coupling of oscillators, but a true, singular phase transition in the $N\to\infty$ limit is not possible. On higher dimensional full or random graphs the Kuramoto equation exhibits universal scaling dynamics of the phase order parameter [31, 32]. According to the theory of universality classes [19] this "simple" model can describe other, more complex models of the brain. Very recently it has been 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 [33]. As the consequence of quenched, purely topological heterogeneity an intermediate phase, located between the standard synchronous and asynchronous phases was found, showing "frustrated synchronization", metastability, and chimera-like states. This complex phase was investigated further in the presence of noise [34] and on simplicial complex model of manifolds with finite and tunable spectral dimension [35] as simple models for the brain.

Anatomical connections [36] and the synchronization networks of cortical neurons [37] indicate a small-world topology [38]. Here we will investigate the characteristic times, corresponding to synchronization or desynchroziation near the transition on a large human connectome graph (KKI-18) and compare it with results, obtained on 2d lattices with additional random, long range connections. The latter also exhibit small-world topology, because the 2d lattice has large clustering and the random, long range connections generate short path lengths among geometrically distant nodes. Our graphs are much larger than those considered before, allowing us to determine universal critical exponents that can be compared with experiments. Furthermore, we have heterogeneity in the intrinsic frequencies as well as in connection weights, which was found to be crucial in case of threshold model simulations [39]. Previously, extended discrete threshold model simulations of activity avalanches on KKI-18 did not support a critical phase transition [39]. It turned out the weight heterogenities were too strong to allow the occurrence of criticality. This means that only the strongly connected hubs played a role in the activation/deactivation processes and weak nodes just followed them. As this appears to be unrealistic and uneconomic in a brain of billions of neurons, an input sensitivity equilibration was assumed via variable, node dependent thresholds. This makes the system homeostatic and simulations proved the occurrence of criticality, as well as robust Griffiths effects [39, 40] in spreading models. Indeed, there is some evidence that neurons have a certain adaptation to their input excitation levels [41] and can be modeled by variable thresholds [42]. Very recently comparison of modeling and experiments arrived at a similar conclusion: equalized network sensitivity improves the predicting power of a model at criticality in agreement with the FMRI correlations [43].

Even more naturally, homeostasis can be achieved in real brains via inhibitory neurons [44, 45, 46, 47, 48], suppressing communications. This provides an alternative way for modifying the positive, undirected links of the KKI-18 graph to test the phase synchronization of the Kuramoto model in the presence of random, negative couplings.

Models and methods

We consider the Kuramoto model of interacting oscillators [28], with phases $\theta_{i}(t)$ located at $N$ nodes of networks, which evolve according to the dynamical equation

[TABLE]

Here, $\omega_{i,0}$ is the intrinsic frequency of the $i$-th oscillator, drawn from a Gaussian distribution with zero mean and unit variance and the summation is performed over other nodes, with connections described by the weighted adjacency matrix $W_{ij}$. The global coupling $K$ is the control parameter of this model, by which we can tune the system between asynchronous and synchronous states. We follow the properties of the phase transition through studying the Kuramoto order parameter defined by

[TABLE]

which is non-zero, above a critical coupling strength, $K>K_{c}$ tends to zero for $K<K_{c}$ as $R\propto\sqrt{1/N}$ or exhibits a growth at $K_{c}$ as

[TABLE]

in case of an incoherent initial state, with the dynamical exponents ${\tilde{z}}$ and $\eta$. In case of a coherent initial state it decays as:

[TABLE]

characterized by the dynamical exponent $\delta$. Here $f_{\uparrow}$ and $f_{\downarrow}$ denote different scaling functions.

We have also investigated the de-synchronization duration distributions by starting the system from fully synchronous or asynchronous states, near $K_{c}$ by measuring the time $t_{x}$ until $R(t_{x})$ first fell below the threshold value: $R_{T}=1/\sqrt{N}$, related to the synchronization noise in the incoherent phase (see Fig. 5). For this measurement we averaged over $\simeq 10^{4}$ runs, using independent random $\omega_{i,0}$ intrinsic frequencies and applied histogramming with increasing bin sizes: $\Delta t_{x}\propto t_{x}^{1.12}$ to estimate the probability distribution $p(t_{x})$.

The following graphs have been considered:

2d lattices with additional, random long-range connections such that $\langle k\rangle=5$ (2dll). We used periodic boundary conditions, simulating high dimensional graphs with supposedly mean-field behavior. 2. 2.

Weighted, symmetric large human connectome graph: KKI-18 [49] downloaded from the Openconnectome project [50]. 3. 3.

KKI-18, with 5% of the links turned to inhibitory: KKI-18-I

We applied the fourth order Runge-Kutta method from Numerical Recipes [51] and the boost library odeint [52] to solve Eq. (2) on various networks. Step sizes: $\Delta=0.1,0.01,0.001$ have been tested, but finally $\Delta=0.1$ precision found to be sufficient. Generally, the $\Delta<0.1$ precision did not improve the stability of the solutions, but caused smaller fluctuations due to the chaotic behavior of Eq. (2) which could be compensated by averages over many independent samples with different $\omega_{i,0}$. The criterion $\epsilon=10^{-12}$ was used in the RK4 algorithm and we parallelized the RK4 for NVIDIA graphic cards (GPU), by which we could achieve a $\sim\times 40$ increase in the throughput with respect to a single 12-core CPU. It is important to note, that first we verified the parallel GPU code, by comparing results on smaller sizes with those, obtained by the serial CPU program. Algorithmic and benchmark details will be discussed elsewhere [53].

We measured the Kuramoto order parameter with a fixed $K$, by increasing the sampling time steps exponentially

[TABLE]

which is a common method in case of PL asymptotic time dependences. In practice we estimate $t_{x}=(t_{k}+t_{k-1})/2$, where $t_{k}$ was the first measured crossing time. The initial conditions were generally $\theta_{i}(0)\in(0,2\pi]$ phases, with uniform distribution, describing fully disordered states. However, for comparison we also performed runs starting from the fully synchronized state: $\theta_{i}(0)=0$. Probability distribution tails were fitted using the least squares fit method above thresholds, fixed by visual inspection of the results. To see the corrections to scaling we determined the effective exponents of $R$ as the discretized, logarithmic derivative of Eq. (4) at these discrete timesteps $t_{k}$, near the transition point

[TABLE]

Here the brackets denote sample averaging over different initial conditions.

The KKI-18 graph has been downloaded from the Open Connectome project repository [50]. This network was generated from Diffusion Tensor Imaging (DTI) [54], approximating the structural connectivity of the white matter of a human brain. It comprises $N=\num{804092}$ nodes, connected via \num41523908 undirected edges, and several small sub-components, which were ignored here.  111Note, that keeping the sub-components, did not change the results within numerical accuracy.. Solving (2) on this graph allows running extensive dynamical studies on present day CPU/GPU clusters, large enough to draw conclusions on the scaling behavior without very strong finite size effects. These connectomes of the human brain possess \SI1mm^3 resolution, using a combination of diffusion weighted, functional and structural magnetic resonance imaging scans. They are symmetric, weighted networks, where the weights measure the number of fiber tracts between nodes. The large graph "KKI-18" used here is generated by the MIGRAINE method, described in [55]. They exhibit hierarchical levels by construction from the Desikan cerebral regions with (at least) two quite different scales. The graph structure can be seen on the Fig. 1, where the modules were identified by the Louvain algorithm [56], then the network of modules was generated and finally visualized using the Gephi tool [57]. This identified $144$ modules, with sizes varying between $8$ and $\num{35202}$ nodes.

In [49] it was found that, contrary to the small world network coefficients, these graphs exhibit topological dimension slightly above $D=3$ and a certain amount of universality, supporting the selection of KKI-18 as a representative of the large human connectomes available. This dimensionality suggests weak long-range connections, in addition to the $D=3$ dimensional embedding and warrants to see heterogeneity effects in dynamical models defined on them.

To keep the local sustained activity requirement for the brain [4] and provide a homeostatic state, we modified KKI-18 by normalizing the incoming weights of node $i$ in [39]: $W^{\prime}_{i,j}=W_{i,j}/\sum_{j\in{\rm neighb.of}\ i}W_{i,j}$ at the beginning of the simulations.

In addition, it is well known that excitation is balanced globally by the inhibitory cells, which is assumed to be a $20\%$ fraction of all neurons. However, the vast majority of inhibitory connections are local as these cells have small dendritic and axonal trees [58]. Their range coincides roughly with the $1mm^{3}$ of voxels in the connectivity data of KKI-18. The tracts, obtained by DTI are in contrast mostly excitatory axonal fibers; middle and long range connections are made by axons of pyramidal cells, which are excitatory. But these axons still could target excitatory and inhibitory cells in a voxel’s area. To model this, we flipped the signs of weights of $5\%$ randomly selected links $i$ as $W^{\prime\prime}_{i,j}=-W^{\prime}_{i,j}$, creating thus a modified graph called: KKI-18-I. Such links are against local synchronization and can be considered as an inhibition mechanism of possible information retrieval mechanism via resonance [59].

Results

The 2dll graph

First we studied the growth of $R(t)$ on the 2dll model of linear size $L=6000$ by starting from states of oscillators with fully random phases and by averaging over $5000-10000$ $\omega_{i,0}$ realizations up to $t=10^{3}$. As Fig 2 shows power-law growth of synchronization emerge up to $t\simeq 100$ in the coupling region: $0.477\leq K\leq 0.478$. Following that the $R(t)$ curves veer up or down, depending on being super or sub-critical. Note, that for $t>\sim 800$ the curves begin to break down, owing to the finite size effect, when the growing correlation volume: $\xi\propto t^{\tilde{z}}$ exceeds the system size $N=L^{2}$. Looking at the effective exponents defined by (7) one can estimate the critical point: $K_{c}=0.4775(3)$, as we expect in the asymptotic $1/t\to\infty$ limit constant local slopes. Off-critical cases exhibit up or down veering curvatures. Here one can read-off the asymptotic value: $\eta=0.55(10)$, on the local slope inset of Fig 2, which is different from the MF value $\eta=0.75$, expected for the Kuramoto model [31] by scaling relations. The obvious discrepancy can be the consequence of very strong corrections-to scaling or some other quench disorder effect discussed further in [32]. Note, that in [31], in case of fully connected graphs, the $\eta=0.75$ exponent could hardly be seen, probably as the consequence of finite size and time limitations. While [31] achieved sizes $N\leq\num{819200}$, here we provide results for much larger systems, containing $N=\num{36000000}$ nodes, but without full topological order.

Instead of going into the details we just show that the relaxation time distributions, both for growth and for decay result in a robust fat tail behavior at the critical point, characterized by $\tau_{t}=1.6(1)$ (see Fig. 3). The evolution of single realizations, in case of growth runs are shown on Fig 5. The dashed line denotes the threshold at which the first passage time $t_{x}$ is measured. We have also measured $t_{x}$ in the case of fully coherent initial states in systems up to $t_{\mathrm{max}}=10^{4}$. Fig. 3 shows the tails of $p(t_{x})$ around the critical point for incoherent initial conditions for $L=6000$ and for coherent initial state with with $L=1000$. In the latter case the decay occurs following a long transient time, thus the time is divided by a factor $100$, but one can observe the same type of PL tails at $K\simeq 0.477$.

With this precision one cannot see a difference in the numerical scaling behavior at the critical point, but we note that the $\tau_{t}=1.6(1)$ estimate is slightly above that one could obtain using the scaling relation

[TABLE]

connecting dynamical exponents, see for example [19]. Note, that the $\delta=0.6(1)$ result is in agreement with those of the extensive decay results, presented in [32] and suggest that $\delta\simeq\eta$ on this scale, ruling out possible strong artifacts of the threshold value selection.

The Connectome graph

Normalized, positive weights

In case of the KKI-18 graph first we determined the crossover point via the inflexion condition, which separates up (convex) and down (concave) veering curves of the growth runs (see Fig 4). As we can see, the transition is much smoother than what we obtained in the 2dll graph. The lower inset of Fig 4) shows the steady state values $R(t\to\infty)$ as the function of $K$, with a very low level of synchronization above the transition. This smooth crossover behavior is not surprising, as the topological dimension of this graph is: $d=3.05<d_{l}=4$ [49]. This behavior is in agreement with PET and fMRI studies, which suggest that the magnitude of activity change from rest to task is rather small. Looking at the shapes of the $R(t)$ curves and the saturation level of the corresponding local slopes we can estimate this crossover at: $K_{c}=1.65(5)$, with an effective scaling exponent $\eta_{\mathrm{eff}}\simeq 0.6(1)$. This exponent value is smaller than the $\eta=0.75$ MF value for the Kuramoto model, but close to our 2dll graph results.

Having determined the transition point we run the numerical solver at control parameter values below $K_{c}$, by starting with thousands of random initial states and measuring the first crossing times $t_{x}$, when $R$ fell below: $1/\sqrt{(}N)=0.001094$ (see Fig. 5).

Following a histogramming procedure, with PL growing bin sizes in $t_{x}$, we obtained the distributions $p(t_{x})$, which exhibit PL tails, characterized by the exponents $1<\tau_{t}<2$ (see Fig. 6). Here the almost $\sim 1/t$ decay at $K_{c}\geq 1.7$ marks synchronized phase, with singular behavior; at $K_{c}$ one obtains $\tau_{t}\simeq 1.2$.

The $\tau_{t}=1.2(1)$ is out of the range of neuro experiments: $1.5<\tau_{t}<2.4$ [13], but a good agreement/overlap can be found in the sub-threshold region. To investigate the effects of inhibition as in [39] here we also studied the modified KKI-18 graph, in which we flipped a small fraction of weight links randomly.

Inhibitory weights

We repeated the analysis of the previous section for the KKI-18-I graph, possessing $5\%$ inhibitory link fraction symmetrically : $W^{\prime\prime}_{ij}=W^{\prime\prime}_{ji}=-W^{\prime}_{ij}$. First we located the transition point as shown on Fig. 7.

The crossover to synchronization occurs at $K_{c}=1.9(1)$, slightly higher than in case of the KKI-18 network. The tails of the $p(t_{x})$ probability distributions exhibit PL-s with $1<\tau_{t}\leq 2$ in the $1.4<K<1.8$ region. These exponent values overlap the range of experiments (see Fig. 8).

The large variation of $\tau_{t}$ below the transition, suggests the strong disorder may cause some GP effects, but in the lack of true phases we can only claim resemblance with recent results on power-grid networks [60], where we pointed out a relation to a phenomena called frustrated synchronization [33, 35].

We have redone this analysis for another random $5\%$ inhibitory link graph sample and arrived to the same results. In [39] the robustness of GP in case of the threshold model dynamical behavior has been tested by the random neglection of 20% of links in one direction. Here we considered the situation of the neglection of all links in one direction: $W^{\prime\prime}_{ij}=-W^{\prime}_{ij}$, $W^{\prime\prime}_{ji}=0$. Even in this extreme anisotropic case one case find a similar extended scaling region below a very smooth crossover as shown in the Supplementary material.

Finally, we considered graphs with $5,10,20\%$ inhibitory node assumptions, by flipping signs of (out or in) link weights of these randomly selected sites. Note, that nodes here represent big bunches of neurons. We show the duration distribution results here for one of the $5\%$ inhibitory node case, when signs of out links are reversed. Below the synchronization transition point, which is at $K_{c}=1.7(1)$ we can find again a region : $1.35<K<1.7$, where PL tailed de-synchronization duration durations emerge as before (see Fig. 9), characterized by the exponents $1<\tau_{t}<2$.

We have also studied the synchronization behavior on graphs with fully inhibited node cases of $5,10,20\%$, by flipping the weight signs of in-links. Even without weight normalization this produces crossovers at $K_{c}\simeq 0.18$ with $\eta_{\mathrm{eff}}\simeq 0.23(5)$, different from the $\eta\simeq 0.6$ value we found up to now. The tails of the $p(t_{x})$ probability distributions exhibit PL-s with well inside the range of experiments (see Supplementary material).

Conclusion and discussion

Brain experiments support evidence for power-law distributed activity avalanches. These have been explained mainly by discrete, threshold type of models, showing exponents close to the MF values, in agreement with neuro measurements. Oscillatory activity is widespread in dynamic neuronal networks [20]. The Blue brain project [61] suggests that the cortical dynamics operates at the edge of a phase transition between an asynchronous phase and a synchronous one with emerging oscillations [62]. A recent MF theory showed the emergence of scale-free avalanches at the edge of synchronization [27].

An extension, taking into account network heterogenities of a large human connectome is provided here, within the framework of the Kuramoto model. We determined the phase synchronization transition points and provided characteristic time exponents, which describe synchronization or de-synchronization events. We found good agreement with the neuro-experimental values in an extended range below the crossover point. We investigated the case, when $5,10,20\%$ fraction of the randomly selected links or nodes are set inhibitory. The obtained $\tau_{t}$ duration exponents have been found to be invariant for these fractions and are in the range of experiments. This conclusion has already been derived within the framework of discrete threshold models [39, 40].

At the transition point the Kuramoto model on these homeostatic KKI-18 network exhibits: $\tau_{t}=1.2(1)$, well below the results for the 2dll graph, $\tau_{t}=1.6(1)$, which is expected to be a system of MF interactions. However, even for this MF like model the dynamical exponents were found to be slightly away from the Kuramoto MF values [31]. This can be the result of enormous corrections to scaling, or due to quenched heterogeneity effects of the 2dll graph. The details of this problem is discussed in [32]. The effective growth exponent on the investigated connectome networks is $\eta\simeq 0.6$, near to the 2dll graph results, except for the node inhibited case, where $\eta_{\mathrm{eff}}\simeq 0.25(5)$ has been found.

Although in the KKI-18 graph the topological dimension is below $d_{l}=4$, a crossover behavior can clearly be identified. Around this smeared transition we found scale-free de-synchronization ’avalanche’ tails, like in case of the dynamical criticality of GP, pointing out relation to possible frustrated synchronization effects [33, 35]. Modules of the connectome graph enhance rare-region effects or frustrated synchroniztion domains.

As it was discussed in [39] such coarse grained connectomes suffer possible sources of errors, like unknown noise in the data generation; underestimation of long connections; radial accuracy, influencing endpoints of the tracts and hierarchical levels of the cortical organization; or transverse accuracy, determining which cortical area is connected to another. Still important modifications, such as inhibitory links, directedness, or random loss of connections up to $20\%$ confirmed the robustness of dynamical scaling, suggesting that fine network details may not play an important role. It is also important, that the PL tail in the weight distribution is similar to what was obtained by a synaptic learning algorithm in an artificial neural network [63]. A very recent experimental study has provided confirmation for the connectome generation used here [64]. Those results suggest that diffusion MRI tractography is a powerful tool for exploring the structural connectional architecture of the brain.

From a neuroscience point of view one may find the Kuramoto model too simplistic to describe the brain. However, at least in the weak coupling limit equivalence of phase-oscillator and integrate-and-fire models has been found [65], that may hold for the sub-critical region, where we observed the dynamical scaling. At criticality universal scaling is expected to hold and thus the exponents obtained for this model can well describe the synchronization transition of other, more complex models. We believe that investigating a basic model is an important and necessary first step for neuroscience as this can provide a representative for a whole class of more real ones.

An interesting continuation of this work would be the study of the frequency entrainment of oscillators, which can exhibit real phase transition at $d=3.05$ of the KKI-18 graph, or consideration of more complex models or graphs than the ones investigated here. Another open point could be the determination of avalanche sizes in these synchronization processes. The codes and the graphs used here are available on request from the corresponding author.

Acknowledgments

We thank R. Juhász for the useful comments and Wesley Cota for generating and providing us Fig. 1. We gratefully acknowledge computational resources provided by NIIF Hungary, the HZDR computing center, the Group of M. Bussmann and the Center for Information Services and High Performance Computing (ZIH) at TU Dresden via the GPU Center of Excellence Dresden. We thank S. Gemming for support. Support from the Hungarian research fund OTKA (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.