Scale-freeness or partial synchronization in neural mass phase oscillator networks: Pick one of two?

Modeling and interpreting (partial) synchronous neural activity can be a challenge. We illustrate this by deriving the phase dynamics of two seminal neural mass models: the Wilson-Cowan firing rate model and the voltage-based Freeman model. We established that the phase dynamics of these models differed qualitatively due to an attractive coupling in the first and a repulsive coupling in the latter. Using empirical structural connectivity matrices, we determined that the two dynamics cover the functional connectivity observed in resting state activity. We further searched for two pivotal dynamical features that have been reported in many experimental studies: (1) a partial phase synchrony with a possibility of a transition towards either a desynchronized or a (fully) synchronized state; (2) long-term autocorrelations indicative of a scale-free temporal dynamics of phase synchronization. Only the Freeman phase model exhibited scale-free behavior. Its repulsive coupling, however, let the individual phases disperse and did not allow for a transition into a synchronized state. The Wilson-Cowan phase model, by contrast, could switch into a (partially) synchronized state, but it did not generate long-term correlations although being located close to the onset of synchronization, i.e. in its critical regime. That is, the phase-reduced models can display one of the two dynamical features, but not both.


Introduction
Characterizing the underlying dynamical structure of macroscopic brain activity is a challenge. Models capturing this large-scale activity can become very complex, incorporating multidimensional neural dynamics and complicated connectivity structures (Izhikevich and Edelman, 2008;Deco and Jirsa, 2012). Neural mass models, or networks thereof, that cover the dynamics of neural populations offer a lower-dimensional and therefore appealing alternative (Deco et al., 2008;Sotero et al., 2007;Ponten et al., 2010). To further enhance analytical tractability one may consider the (relative) phase dynamics between neural masses. We previously showed under which circumstances certain neural mass models can be reduced to mere phase oscillators (Daffertshofer and van Wijk, 2011;Ton et al., 2014)see also (Schuster and Wagner, 1990;Hlinka and Coombes, 2012;Sadilek and Thurner, 2015) and thus established a direct link between these two types of models. A minimal model describing phase dynamics is the Kuramoto network (Kuramoto, 1975), which in its original form consists of globally coupled phase oscillators. Generalizations of this model by adding delays and complex coupling structures result in a wide variety of complex dynamics (Acebr on et al., 2005;Breakspear et al., 2010). Even in its original form, however, the Kuramoto model is capable of showing non-trivial collective dynamics. A mere change of the (global) coupling strength can yield a spontaneous transition from a desynchronized to a synchronized state, i.e. the dynamics can pass through a critical regime.
Synchronization of neural activity plays a crucial role in neural functioning (Fries, 2009). In the human brain, synchronized activity can be found at different levels. At the microscopic level temporal alignment in neuronal firing is a prerequisite for measurable cortical oscillations (Pfurtscheller and Lopes Da Silva, 1999). However, it also manifests itself at the macroscopic level in the form of global resting state networks (Biswal et al., 1995;Deco et al., 2011). Synchronization properties are modulated under the influence of task conditions in, e.g., motor performance (van Wijk et al., 2012), visual perception (Singer and Gray, 1995), cognition (Fell and Axmacher, 2011) and information processing (Fries et al., 2007;Singh, 2012;Voytek and Knight, 2015). Epilepsy, schizophrenia, dementia and Parkinson's disease come with pathological synchronization structures (Schnitzler and Gross, 2005;Broyd et al., 2009). When aiming for a concise but encompassing description of brain dynamics, a macroscopic network model should capture this wide range of synchronization phenomena.
According to the so-called criticality hypothesis (Beggs, 2008), the human brain is a dynamical system in the vicinity of a critical regime. Its dynamics is located at the cusp of dynamic instability reminiscent of a non-equilibrium phase transition in thermodynamic systems (Chialvo, 2010;Beggs and Timme, 2012). The conceptual appeal of the critical brain lies in the fact that networks operating in this regime show optimal performance in several characteristics relevant to cortical functioning (Shew and Plenz, 2013). Critical dynamics often display power laws in multiple variables (Stanley, 1971) and have been observed in, e.g., size and duration distributions of neuronal avalanches (Beggs and Plenz, 2003) and EEG cascades (Fagerholm et al., 2015). Power-laws are also manifested as scale-free autocorrelation structures of the amplitude envelopes of encephalographic activity (Linkenkaer- Hansen et al., 2001;Palva et al., 2013). Very recently, long-range temporal correlations have been reported in fluctuations of the phase (synchronization) dynamics in neural activity (Botcharova et al., 2015b;Daffertshofer et al., 2018). The nature of these power-law forms in the correlation structure can be quantified by the Hurst exponent H (Hurst, 1951). Its value characterizes the correlations between successive increments of the signal, with H > 0:5 and H < 0:5, the first marking persistent behavior, i.e. positive, long-range correlations, in the time series, and the latter anti-persistent behavior, i.e. a negative auto-correlation.
In this study we considered the phase descriptions of two classical neural mass networks: the Wilson-Cowan firing rate and Freeman voltage model, both equipped with neurobiologically motivated coupling and delay structures. Coupling and delay structures were obtained from DTI data and the Euclidean distances between nodes, respectively. To anticipate, the two models lead to two qualitatively different phase synchronization dynamics.

Wilson-Cowan model
The first neural mass model we studied is the Wilson-Cowan model that describes the dynamics of firing rates of neuronal populations (Wilson and Cowan, 1972). We always considered properly balanced pairs of excitatory and inhibitory populations, E k ¼ E k ðtÞ and I k ¼ I k ðtÞ, respectively. We placed such pairs at every node of a network. Nodes were coupled to other nodes through the connections between excitatory populations given by a DTI-derived coupling matrix S kl forming a network of k ¼ 1; …; 90 nodes. We illustrate the basic structure of this network in Fig. 1 and the coupling matrix in Fig. 2, left panel. The dynamics per node k was cast in the form (1) where the coupling constants c EI , c IE , c EE , c II quantify the coupling strength within each (E/I) pair. The function Q ½x ¼ ð1 þ e Àx Þ À1 is a sigmoid function that introduces the thresholds θ E and θ I that need to be exceeded by the total input into neural mass k to elicit firing; the parameters a E and a I describe the slopes of the sigmoids. The delays τ kl were determined by conduction velocity and the Euclidean distance between nodes k; l. In the following, delay values are given in milliseconds.
Appropriate choices of the time constants μ k and external inputs q k guaranteed self-sustained oscillations in the alpha band (8-13 Hz). By randomizing the constant external inputs q k ; μ k across k we introduced heterogeneity in oscillation frequencies.

Freeman model
The Freeman model (Freeman, 1975) describes the mean membrane potentials E k ; I k of neural populations. In analogy to (3) its dynamics per node can be given by Coupling structure of both neural mass networks. A proper balance between excitatory (E k ) and inhibitory (I k ) populations leads to self-sustained oscillations in the networksee also Appendix B.3. Self-coupling (c EE , c II ) is only present in the Wilson-Cowan network and, hence, plotted in gray. External inputs are indicated by q k . The coupling matrix S kl connecting the excitatory populations was based on structural DTI data.

Fig. 2.
Illustration of the DTI-derived structural connectivity matrix S kl and the matrix C F kl in (6). In the latter we incorporated the inter-pair coupling c EI and c IE together with the scaled structural connectivity K⋅S kl . The upper left block of C F kl has the same structure as S kl . The two diagonals represent the coupling strengths c EI and c IE .
where k ¼ 1; …; N with N being the number of excitatory populationsthe corresponding schematic is again given in Fig. 1. The sigmoid function Q ½x here covers the effects of pulse coupled neurons in the populations adjusted by the scaling parameter γ (Marreiros et al., 2008). The parameters α k ; β k represent mean rise and decay times of the neural responses in population k, which we here varied to introduce frequency heterogeneity in the system. Analogous to (1) appropriate parameter values guaranteed self-sustained oscillations in the alpha frequency band.
In (2) we separated the excitatory and inhibitory nodes to stress the similarity between both networks. For the sake of simplicity, however, we rather combine both equations in (2) with k; l ¼ 1; …; 2N. In (3), the delays τ kl between excitatory nodes are equal to delays τ kl in (2). Note that we assumed delays between the local pairs of excitatory and inhibitory nodes (½k l ¼ ½1; …; N; k þ N and ½k l ¼ ½l þ N; 1; …; N) to be negligible.

Phase description
The phase dynamics could be derived by transforming the systems to their corresponding polar coordinates around an unstable focus. We described its dynamics in terms of the periodic function A k cosðΩt þ ϕ k Þ, with A k denoting the amplitude, ϕ k the relative phase, and Ω the central frequency of the oscillation. We averaged the dynamics over one period 2π=Ω under the assumption that the characteristic time scale of the A k and ϕ k dynamics significantly exceeded this period, i.e. _ A k =A k ; _ ϕ k ≪jΩj. That is, the variables ϕ k ; A k evolved slowly enough to be considered constant within one period. Further, we assumed the time delays τ kl to be of the same order of magnitude as the period of oscillation, such that they could be captured by model-dependent phase shifts Δ kl in the phase dynamics. More details of the derivation are given in Appendix B, which builds on Daffertshofer andvan Wijk (2011), Ton et al. (2014). In consequence, both phase dynamics obeyed the form with ω k the natural frequency of the oscillator at node k, D kl the phase coupling matrix and Δ kl the time-delay induced phase shifts. For the Wilson-Cowan phase model (superscript WC) we obtained where Q ' denotes the derivative of Q resulting from a Taylor approximation around the points χ ð0Þ E;k , χ ð0Þ I;k ; their detailed expressions of the aforementioned unstable focus are given in Appendix B.1. By virtue of the definition of Q , one has Q ' ! 0. The definition of the parameters Λ k and ϖ k can also be found in B.1. In case of the Freeman phase model (superscript F) the expressions for ω k , D kl and Δ kl read here V ð0Þ l refers to the unstable focuswe again refer to B.1 for details. It is important to realize that due to the inhibitory coupling c EI , the coupling matrix C F kl and hence D F kl have both positive and negative entries. In contrast, for S kl and hence for D WC kl we have S kl ; D WC kl ! 0. This change in sign can have major consequences. For instance, the Wilson-Cowan phase model (5) resembles the Kuramoto-Sakaguchi model with phase lag Δ kl ; for τ kl ¼ τ ¼ 0 we have Δ kl 2 ð0;π=2Þ, and otherwise Δ kl 2 ðÀπ=2; π=2Þ since Ωτ kl 2 ð0; π=2Þ and because of our particular choice of parameters which yields Λ k ! 0. That is, a transition to full synchronization can occur if the coupling strength K exceeds a critical value K c . However, from (6) it follows that the left upper block of D F kl contains negative entries. Together with the π=2 phase shift, the Freeman phase dynamics is therefore more closely related to the repulsive cosine-variant of the Kuramoto network (Van Mieghem, 2009;Burylko et al., 2014). As will be shown below, this qualitative difference in dynamics led to profound contrasts in model behavior.

Simulations
Phase time series ϕ k ðtÞ were obtained by integrating the system (4) using either (5) or (6). We first determined the fixed points E k (Freeman) around which we observed stable oscillations. For the latter we also determined the characteristic frequency Ω and amplitudes A k . These initial estimates were achieved by a five second simulation of the systems (1)/(3) using an Euler scheme with time step Δt ¼ 1 ms. The choice for the Euler method was motivated by the implementation of delays in the coupling terms. Testing a more elaborated predictor/corrector algorithm (Shampine and Thompson, 2001) revealed little to no difference but required far more numerical resources.
Integration of the phase systems (4) was performed by means of an adaptive Runge-Kutta (4,5) algorithm with variable step size. Simulation time T ¼ 302 seconds matched the length of the available empirical time series, where we discarded the first two seconds of simulation to avoid transient effects due to the first random initial condition. To exclude effects of a specific natural frequency distribution of ω ð⋅Þ k , T was split into twenty bins with random duration T n > 14 s. The initial condition of bin n was matched to the last sample of bin n À 1. For each n a new set of parameter values q k ; μ k (Wilson-Cowan) or α k ; β k (Freeman) was chosen at random, under the constraint that the characteristic frequency Ω fell within the alpha band in all cases. The parameters μ k ; q k ; α k ; β k , and T n were drawn at random but the corresponding sets were kept equal across all simulation conditions, i.e. for all combinations of K and hτ kl i. We thus obtained twenty sets fω

Power-law behavior
We measured the amount of synchronization via the phase coherence, i.e. the modulus of the Kuramoto order parameter given as where ϕ k ðtÞ followed the dynamics (4). Next, we assessed the autocorrelation structure of RðtÞ by means of a detrended fluctuation analysis (DFA) (Peng et al., 1994). In DFA the cumulative sum of a time series yðkÞ, i.e. YðtÞ ¼ P t k¼1 yðkÞ, is divided into non-overlapping segments Y i ðtÞ of length T seg . Upon removing the linear trend Y trend i ðtÞ in segment i, the fluctuations F i ðT seg Þ corresponding to window length T seg are given by When these fluctuations scale as a power law, i.e. F i ðT seg Þ $ ðT seg Þ α , the fluctuations, and hence the associated autocorrelations, can be considered scale-free. The corresponding scaling exponent α resembles the Hurst exponent H (Hurst, 1951) and characterizes the correlation structure (the resemblance is proper if yðtÞ stems from a fractional Gaussian noise process). We assessed the presence of a power law in F i in a likelihood framework by testing this model against a set of alternatives (Ton and Daffertshofer, 2016). By applying the Bayesian information criterion (BIC) we could determine the model that constituted the optimal compromise between goodness-of-fit and parsimony (Burnham and Anderson, 2002). More details of the DFA and model comparison are given in A.1.
To determine the significance of the model results, we constructed surrogate data sets by generating 90 phase times series ϕ ðsurrÞ k ðtÞ, which equalled the number of excitatory nodes. The surrogate phase time-series consisted of random fluctuations around linear trends sampled from the ω F k distribution using the same T n partitions as in the model simulation conditions. We used a Wilcoxon rank-sum test to test the results from surrogate time series against simulated time series in a non-parametric way. For evaluation of the scaling exponents, we only incorporated those conditions that showed power-law scaling as assessed by the BIC.

Functional connectivity
We compared spatial correlation structures in terms of the functional connectivity matrices generated by both models with an empirically observed functional connectivity. For the latter we incorporated a previously published data set (Cabral et al., 2014;Daffertshofer et al., 2018). We refer to Cabral et al. (2014) for details concerning data acquisition and preprocessing of both the MEG and the DTI derived anatomical coupling matrix S that was used in the coupling matrices D ð⋅Þ kl given in (5&6). In brief, MEG of ten subjects was recorded in resting state conditions (eyes closed) for approximately five minutes. These MEG signals were beamformed onto a 90 node brain parcellation (Tzourio-Mazoyer et al., 2002), such that 90 time-series y k ðtÞ were obtained with a sampling frequency of 250 Hz. The signals y k ðtÞ were bandpass filtered in the frequency range 8-12 Hz and subjected to a Hilbert transform to obtain the analytical signal, from which the Hilbert phase could be extracted.
Using the phase time series from both MEG data, ϕ ðMEGÞ k ðtÞ, and phase time series generated by (4), ϕ ðsimÞ k ðtÞ, we calculated the pair-wise functional connectivity P ð⋅Þ via the pair-wise phase synchronization in the form of the phase locking value (PLV) (Lachaux et al., 1999). In its continuous form its matrix elements are defined as Note that for P ðMEGÞ we removed all individual samples that displayed relative phases in intervals I around 0 or AEπ, as for these samples true interaction and effects of volume conduction could not be disentangled (Nolte et al., 2004). The intervals I were defined as where Ω c is the center frequency (in this case 10 Hz) and F s the sampling frequency. In the simulations we calculated the PLV matrix according to (19) for each partition T n and afterwards averaged the thus obtained twenty PLV matrices.

Synchronization
Although both RðtÞ and P ð⋅Þ kl are synchronization measures, they measure two qualitatively different forms of synchronization, which is the reason why they offer resolution in either the temporal or the spatial domain, respectively. Functional connectivity P ð⋅Þ kl measures temporal alignment of two phase time series ϕ ð⋅Þ k ðtÞ, ϕ ð⋅Þ l ðtÞ by means of an averaging over time in (19), such that P ð⋅Þ kl provides resolution in the spatial domain, as indicated by the subscript kl. In contrast, from (17) it follows that calculating RðtÞ involves an average over k, i.e. over spatial coordinates, for each time instant t. This measure therefore provides resolution in time. That is, P ð⋅Þ kl measures temporal synchronization and offers spatial resolution, whereas for RðtÞ the opposite holds.
To gain more insight into the mechanisms responsible for the differential effects on synchronization behavior in the two models, we further considered the measures *

RðtÞ
that is, hRðtÞi is the temporal average of the order parameter and hP ð⋅Þ i corresponds to the average magnitude of pair-wise phase synchronization over the network.

Statistics
Model performance in terms of replicating spatial correlation structure was measured by calculating the Pearson correlation coefficient ρ between the lower triangular entries of P ðsimÞ and P ðMEGÞ , where we excluded the main diagonal to omit spurious correlations. Since the sampling distribution of the P ð⋅Þ entries are not normally distributed, we applied a Fisher z-transform before calculating the correlations.
Restricting ourselves to the lower-diagonal entries was sufficient due to the symmetry in the phase coherence measure. We also excluded the diagonal entries to avoid spurious correlations resulting from the trivial value P ð⋅Þ kk ¼ 1.

Power-law behavior
Only the Freeman phase dynamics generated power laws and thus long-range temporal correlations in the evolution of phase synchronization for a broad range of parameter values. In Fig. 3a we display the results as a function of coupling strength K and mean delay τ kl . Since for none of the parameter values the Wilson-Cowan phase model yielded power laws, we do not show the corresponding results for this model. The average value (AE SD) of the scaling exponents was α ¼ 0:56 AE 0:02 for the Freeman model, which is significantly different from the surrogate results α ¼ 0:501 AE 0:012 (p < 10 À4 ). In this average we only considered those realizations that were classified as power laws. This result qualitatively agreed with the observed value in MEG data (α ¼ 0:62, Daffertshofer et al., 2018), as both indicate persistent behavior and thus long-range temporal correlations. Fig. 3b provides examples of the log-log fluctuation plots for a single realization (K ¼ 0:7, hτi ¼ 9:4) for both the Freeman and the Wilson-Cowan based model (upper and lower panel, respectively). The latter clearly deviated from linearity indicating that it did not scale as a power law. There the model selection procedure assigned a piece-wise linear function (dashed gray in Fig. 3b) to the F i results confirming the deviation from linearity. Other parameter values yielded similar results for this model. The Freeman phase model yielded a power law with scaling exponent α ¼ 0:56; Fig. 3b, upper panel.

Functional connectivity & synchronization
As said, we quantified functional connectivity as pair-wise phase synchronization of the phase variables ϕ ð⋅Þ k followed either the dynamics (4) indicated by the superscript ðsimÞ or ðsurrÞ or the instantaneous Hilbert phase extracted from source-reconstructed MEG data, superscript ðMEGÞ. Model performance was quantified by means of the Pearson correlation coefficient ρ between the P ðsimÞ and P ðMEGÞ matrices. Maximal P ðMEGÞ -P ðsimÞ correlations were ρ ¼ 0:56 in both models for parameter values (K ¼ 0:8; hτ kl i ¼ 9:4) for the Freeman and (K ¼ 0:7; hτ kl i ¼ 9:4) for the Wilson-Cowan phase model (Fig. 4). This value is comparable to values reported in previous simulation studies (Cabral et al., 2014;Deco et al., 2013Deco et al., , 2014, but in contrast to the latter two studies no critical coupling strength at which model-data correlations collapse was found. The Freeman phase model appeared less sensitive to overall coupling strength than the Wilson-Cowan phase dynamics. We mentioned above that RðtÞ and P ð⋅Þ kl measure two qualitatively different forms of synchronization. That P ð⋅Þ kl and RðtÞ indeed constitute two different aspects of synchronization is reflected in the results. Whereas the qualitative difference in the phase coupling matrices D WC kl and D F kl did affect the autocorrelation structure in RðtÞ (Fig. 3b), i.e. the Wilson-Cowan model did not resemble power-law behavior while the The values on the x-axis are in milliseconds on a logarithmic scale (based on segment sizes of 10 4 to about 10 6:8 ms). On the vertical axis the expectation value of F i is shown that was determined via the corresponding probability densitiesp n ; see also A.1. The Wilson-Cowan phase model did not result in scale-free correlations for any of the parameter values, with typical results for the log-log fluctuation plots similar to the lower panel in Fig. 3b. The DFA result for the Freeman phase model (upper panel Fig. 3b) was classified as a power law with α ¼ 0:56this was significantly different from merely random noise when tested against surrogates.

Fig. 4. a:
Pearson correlation values between the Fisher-z transformed P S and P ðMEGÞ functional connectivity matrices for the Wilson-Cowan phase model as function of coupling K and mean delay hτ kl i (milliseconds). The colored shading codes the correlation values and correspond to the colorbar on the righthand side. The non-shaded area corresponds to the case in which the correlation was not significant. Fig. 4b: Similar to Fig. 4a but for the Freeman phase model. To respect weak coupling, the maximum coupling strength was set to K ¼ 0:7. Results were averaged over ten realizations for every parameter combination.
Freeman model did, it had only a minor influence on functional connectivity in that in both cases a similar maximum correlation with the empirical functional connectivity could be achieved (Fig. 4).
The averaged measures hRðtÞi and hP ð⋅Þ i served to quantify differential effects on synchronization behavior in the two models. In Fig. 5 we display the results as function of both delay and coupling together with the surrogate and MEG data values. As expected from the repulsive coupling in the Freeman phase model, phases were dispersed with hRðtÞi values significantly below surrogate values (p < 10 À4 ). In contrast, the Wilson-Cowan model resulted in a partially synchronized state, which better corresponded to empirical findings (gray solid lines). In accordance with the results on P ðsimÞ -P ðMEGÞ correlations, the qualitative difference between models in hRðtÞi did not transfer to pair-wise synchronization magnitude hP ð⋅Þ i. That is, both models yielded significantly larger hP ð⋅Þ i values than obtained for the surrogate data set (p < 10 À4 ). This was the case despite the spatial desynchronization of the network of Freeman models.
This striking result led us to further assess the dynamics of both phase models by supplementary simulations with considerably larger coupling strengths. These coupling values were beyond the weak coupling assumption, rendering the validity of the phase dynamics for these parameter values questionable; see also B.3. It did, however, provide additional insight into the dynamical properties of the phase model (4), especially with respect to the synchronizability of these networks. As shown in Fig. 6, the Freeman phase model did not allow for a (partially) synchronized state even for large coupling strength. In contrast, for sufficiently small delays the Wilson-Cowan phase model entered a fully synchronized state. This is consistent with the coupling structure of both models given by (6) and (5) respectively.
Although the degree of phase synchronization of the Freeman phase model was consistent with a repulsively coupled phase oscillator network, the inhibitory connections in C F kl made a direct comparison with the repulsive cosine variant of the Kuramoto network non-trivial. Nevertheless we expected these models to show similar behavior, since the inhibitory connections were rather sparse compared to excitatory ones (Fig. 2). To test this, we also considered an alternative: a Freeman model that only comprised the excitatory part D F kl , i.e. the left upper block of this matrix. Results are summarized in Appendix C. In a nutshell, these results indicate that the scale-free correlation structure displayed by the Freeman phase model is caused by (the nature of) the coupling between the excitatory units. Its dynamics can thus be understood by considering the phase dynamics (4)/(6) as a repulsively coupled Kuramoto network.

Discussion
We contrasted the phase dynamics derived from two seminal neural mass models, representing the integrated contribution of large numbers of neurons within populations. Neural mass models have been used for modeling a wide range of neural phenomena, ranging from the origin of alpha band oscillations and evoked potentials (Lopes Da Silva et al., 1974;Jansen and Rit, 1995) to the onset of pathological brain activity patterns such as epileptic seizures (Breakspear et al., 2006;Rodrigues et al., 2009). The phase reduction yielded phase oscillator networks that differed qualitatively in their coupling structure. Nevertheless, both models performed comparably well in the spatial domain as assessed by the P ðsimÞ -P ðMEGÞ correlations, i.e. they resulted in similar pair-wise  synchronization characteristics as featured by the experimentally observed data. A related finding has been reported by Mess e and coworkers (Mess e et al., 2014) who showed that model performance in terms of P correlations was relatively independent of nodal dynamics. Here it is important to realize that structural connectivity has a high predictive value for (empirical) functional connectivity (Bullmore and Sporns, 2009;Hlinka and Coombes, 2012;Ton et al., 2014), which here could also be confirmed by correlating P ðsimÞ kl and S kl (see Fig. D.1). That is, both models generated functional connectivity structures that were highly correlated with S kl . This same notion forms the basis for the general finding in RSN modeling studies that optimal model performance occurs near the critical point . The critical slowing down around the bifurcation point allows for a maximal reflection of S kl into functional connectivity (Deco et al., 2013). We showed that the reflection of structural into functional connectivity may occur in two models generating qualitatively different dynamics. This indicates that an inference about the dynamical regime, in particular regarding criticality, on basis of the P correlations alone is non-trivial; at least for the phase oscillator models considered here; see also Hansen et al. (2015) for are related conclusion. Whether this extends to more complex networks consisting of detailed neuronal models is beyond the scope of the current study.
In contrast to the pair-wise phase synchronization (PLV), the models differed qualitatively regarding the phase synchronization quantified by the phase coherence RðtÞ. In particular, only the Freeman phase model displayed scale-free autocorrelation structures observed in data, revealing complex characteristics in its dynamics. The values of the scaling exponents (α > 0:5) revealed the presence of long-range temporal correlations, which qualitatively agrees with the correlation structures reported in brain activity (Linkenkaer- Hansen et al., 2001;Palva et al., 2013;Botcharova et al., 2015b;Daffertshofer et al., 2018).
In the Kuramoto model (Kuramoto, 1975), critical coupling strength is the value of the coupling parameter K for which the desynchronized state loses stability and the system enters a (partially) synchronized regime (Strogatz, 2000). Here, synchronization is quantified by RðtÞ, and hence measures spatial synchronization in the network. Functional connectivity, however, is determined by the temporal alignment of, in the present study, phase signals ϕ l ðtÞ, ϕ k ðtÞ and thus reflects a fundamentally different form of synchronization. We showed that these forms of synchronization were affected differently by the generating dynamics: pair-wise synchronization largely agreed between models, whereas RðtÞ did not. The average value and the autocorrelation structure of RðtÞ were affected by the qualitative difference in coupling structure between models.
The Wilson-Cowan phase model displayed a transition into a fully synchronized state for sufficiently large coupling; see Fig. 6. Combined with the partial synchronization displayed in Fig. 5c and a, this indicates that the Wilson-Cowan model for K ¼ ½0:1; 0:7 is located at the onset of synchronization, i.e. in its critical regime. Although associated with critical dynamics (Stanley, 1971;Botcharova et al., 2015a), we did not observe power-law correlation structures in this model. Similar findings have been reported by Botcharova et al. (2015a) for phase difference time series Φ kl ðtÞ ¼ ϕ l ðtÞ À ϕ k ðtÞ, not only in case of the standard uniformly coupled Kuramoto network, but also for a more biologically plausible model incorporating a DTI derived coupling matrix and distance-related delays, see also Cabral et al. (2011). However, long-range temporal correlations were observed in Φ kl ðtÞ as well as RðtÞ in resting state brain activity in Botcharova et al. (2015b) and Daffertshofer et al. (2018), respectively. This suggests that the critical regime in Kuramoto-type networks has different properties compared to the dynamical regime of the resting brain, be the latter critical or not.
The desynchronized state for the repulsive coupling in the Freeman phase model (6) was consistent with various analytical results (Van Mieghem, 2009;Burylko et al., 2014;Pikovsky and Rosenblum, 2009); cf. Fig. 6. A desynchronized network, however, does not exclude complex dynamics as reflected in the presence of scale-free autocorrelations in the Freeman model. The topology of this model may be regarded as related to phase oscillator networks consisting of so-called conformists and contrarians studied by Hong and Strogatz (2011), Burylko et al. (2014). The latter showed that, even for small networks, a variety of complex dynamics including chaos may occur. A similar finding has recently been reported by Sadilek and Thurner (2015), who studied a two-layered Kuramoto network derived from the same Wilson-Cowan dynamics as considered here. They identified a chaotic region with the largest Lyapunov exponents arising at the boundary of synchronization, i.e. in the critical regime. By changing the value of the delay parameter range, this model could switch between a synchronized and desynchronized state through a bifurcation.
Despite the fact that the model in Sadilek and Thurner (2015) and the Wilson-Cowan phase model in the current study were derived from the same dynamics (1), both networks are quite different in their topology and in their delay structure. The model in Sadilek and Thurner (2015) contained an excitatory and inhibitory layer, whereas this was not the case in (4)/(5). The reason for this discrepancy is that Sadilek & Thurner described the oscillatory trajectory solely by the phase variable, whereas we also took amplitude into account; see also (Schuster and Wagner, 1990;Daffertshofer and van Wijk, 2011). As a consequence the reduction in dimensionality in the Wilson-Cowan phase description that we found when deriving (4)/(5) from (1), did not occur in Sadilek and Thurner (2015). As a consequence, the inhibitory connections in the neural mass dynamics were retained in the phase model in that study. A second distinction between both models is the delay structure. In both models (1) and (2) we regarded the delays between excitatory and inhibitory units to be negligible compared to those between excitatory units, as these connections represented long-range connections subject to finite conduction delays. In contrast, the delay parameter in Sadilek and Thurner (2015) quantified the delay between excitatory and inhibitory units and excitatory-excitatory delays were assumed to be zero.
With the two models considered here we could explain two profound phenomena observed in brain activity. The Freeman phase model generated the type of autocorrelation structures observed in brain activity, but its coupling structure resulted in a desynchronized network, i.e. low RðtÞ values, that did not agree with MEG recordings (see Fig. 5). Additionally it could not account for a transition into partially synchronized states, let alone the (pathological) fully synchronized one. In contrast, the Wilson-Cowan phase model could cover these synchronization phenomena, but it did not show the complex dynamics associated with (resting state) brain activity. The fundamental difference in coupling structure, combined with the analytical results discussed above, suggests that these dynamical properties are mutually exclusive for the models considered here.
We are left with the question, whether one of these models could be modified in such a way that it can exhibit both phenomena. First we have to admit that our DTI-based construction of anatomy and delays is a clear simplification of the 'real' structural connectivity. Adjusting this may have major consequences for the resulting phase dynamics. The aforementioned study by Sadilek and Thurner (2015) gives an indication for this, since they showed that a connectivity structure allowing for comparatively dense inhibitory connectivity yielded complex dynamics in the form of chaos. Interestingly, in other types of models inhibitory connections have been shown to be determinants in generating critical states (Mazzoni et al., 2007;Tetzlaff et al., 2010;Shew et al., 2011) and capacity of information transfer (Deco and Hugues, 2012;Shew and Plenz, 2013). However, these results reflected the dynamics within a neural population rather than the dynamics in the global cortical network considered here. Neurophysiological findings indicate that the long-range connections between areas are excitatory (Kandel et al., 2013;Douglas and Martin, 2004;Sotero et al., 2007) with inhibitory connections only providing local inhibitory feedback. From such a neurophysiological perspective we regard our coupling structure to be more realistic in the context of global cortical networks than the one in Sadilek and Thurner (2015). Thus, although incorporating inhibitory connections could potentially merge the dynamical properties of the Wilson-Cowan and Freeman phase descriptions, such a coupling structure would violate its neurophysiological plausibility and thus the appeal of deriving these networks from a neural mass dynamics. This is not to say that the network in Sadilek and Thurner (2015) is unrealistic from a neurophysiological point-of-view, but both the connectivity and the delay structure may be more representative of local cortical interactions than of global (cortical) networks considered here.
As an alternative one may extend the models to the stochastic regime, e.g., by adding noise to the firing rate or membrane dynamics. Dynamic noise is known for its capacity to alter the correlation structure of global outcome variables like the order parameter RðtÞ. Dynamic noise can also influence synchronization patterns and that not only by causing phase diffusion or shifting the critical point at which synchronization may emerge; in the case of common noise, it may even induce synchronization. A more detailed discussion of network dynamics under impact of random fluctuations, however, is far beyond the scope of the current study.
Delays in networks can lead to very complex dynamics. Since we considered the dynamics of the relative phases that were assumed to evolve slowly with respect to the oscillation frequency Ω, the delays between neural masses mapped to mere phase shifts in (4). Therefore a comparison of the networks in which delays explicitly influence the phase interactions, such as in Cabral et al. (2011) and the analytical results by Kim et al. (1997), Yeung and Strogatz (1999), Choi et al. (2000), cannot be readily made. In the case of delayed phase interactions, however, scale-free correlations could not be observed in a phase oscillator network incorporating a similar coupling scheme to the one employed here (Cabral et al., 2011;Botcharova et al., 2015a). Taken together, our findings suggest that phase oscillator networks without dense inhibitory coupling throughout the whole network, are not capable of showing the entire dynamic spectrum of resting state brain activity. Whether this limitation is posed by the phase oscillator network itself or the consequence of collapsing population dynamics onto a low-dimensional description in the form of a neural mass model remains to be seen.

Conclusion
We illustrated some challenges when deriving and interpreting the phase dynamics of neural mass models. As an example we employed networks of Wilson-Cowan firing rate models and networks of voltagebased Freeman models. The phase dynamics of these models differed qualitatively by means of an attractive coupling in the first and a repulsive coupling in the latter. While both phase dynamics did cover the functional connectivity observed in resting state activity, they failed to describe two pivotal dynamical features that have been reported in many experimental studies: (1) a partial phase synchrony with a possibility of a transition towards either a desynchronized or a (fully) synchronized state; (2) long-term autocorrelations indicative of a scale-free temporal dynamics of phase synchronization. The phase dynamics of the Freeman model exhibited scale-free behavior and the Wilson-Cowan phase model could switch into a (partially) synchronized state. However, none of the phase models allowed for describing both dynamical features in unison.
There is a range of possibilities to modify these models, e.g., by misbalancing excitatory and inhibitory units or by introducing delays that are biologically less plausible than the ones we chose. Alternatively, one may consider the phase dynamics further away from the onset of oscillations (Hopf-bifurcation) that limits analytic approaches. By either of these adjustments one may lose the direct link to the structural connectivity structure. In our example, neither of the phase dynamics can capture the full dynamical spectrum observed in cortical activity. We have to conclude that modeling phase synchronization and, in particular, inferring characteristics of its underlying neural mass dynamics require great care.

Acknowledgments
We would like to thank Mark Woolrich for his contribution to data acquisition and analysis and the fruitful discussion about interpretation of our findings.
We To assess the temporal character of ϕ ð⋅Þ k , we calculated the Kuramoto order parameter defined as where for all models k ¼ 1;…;N. That is, we only used the excitatory phases to calculate RðtÞ. In analogy with the procedure for empirical data discussed in Daffertshofer et al. (2018), we z-scored the RðtÞ time series, such that differences in scaling behavior could not be attributed to differences in the stationary statistics of the RðtÞ time series. We resampled RðtÞ to 250 Hz to match the sampling frequency of the data as well as to obtain an equally spaced time axis necessary for the detrended fluctuation analysis (DFA) (Peng et al., 1994) used to characterize the RðtÞ autocorrelation structure. To assess the presence of scale-free autocorrelations in RðtÞ, we used a modified version of the conventional DFA procedure. We shortly summarize this below; for a detailed explanation we refer to Ton and Daffertshofer (2016).
In line with the outline around Eq. (8), consider (the cumulative sum of) a time series YðtÞ, t ¼ 1; …; N that is divided into N=n non-overlapping segments Y i ðtÞ of length n with t ¼ 1; …; n. Upon removing the linear trend Y trend i ðtÞ in segment i, the fluctuations F i ðnÞ corresponding to window length n are given by A. Daffertshofer et al. NeuroImage xxx (2017) In the conventional DFA procedure one calculates the average fluctuation magnitudes We regarded the F i as a set of ⌊N=n⌋ realizations of the 'stochastic' variable F i and determined its probability density p n ðF i Þ. When these fluctuations scale as a power law, i.e. F i ðn⋅cÞ ¼ n α F i ðcÞ, we find that logðF i ðn⋅cÞÞ ¼ αlogðnÞ þ logðF i ðcÞÞ. Hence, under a transformation to logarithmic coordinatesñ ¼ logðnÞ, e F i ¼ logðF i Þ a power law appears as a linear relationship. To identify whether power-law scaling was present we fitted a set of candidate models f θ ðñÞ parametrized by the set θ. The linear model corresponding to power-law scaling was contained in this set, such that we could compare it against alternatives. For this comparison we defined the log-likelihood function as wherep n denotes the probability density p n transformed to the double logarithmic coordinate system. In (A.1) one evaluates for each n the probability densityp at model value f θ ðñÞ and defines the likelihood function as its product. The purpose of calculating L was to be able to use of the Bayesian Information criterion (BIC) defined as to compare different models f θ . In (A.2) M denotes the number of different interval sizes n, k the number of parameters in the model (the size of the set θ) and L max the maximized likelihood with respect to a particular model f ð⋅Þ θ . The model resulting in a minimal value of the BIC compared to alternative models was considered to be the optimal model; providing the optimal compromise between goodness-of-fit and parsimony (Burnham and Anderson, 2002). The set of candidate models was given by a combination of polynomial forms including the sought-for linear model. We further included an alternative exponential model fit as well as the form resulting from an Ornstein-Uhlenbeck process and, last but not least, a piece-wise linear model: The scaling exponent α was determined as the slope of the linear relationship f θ , i.e. α ¼ θ 2 in f 1 θ ðxÞ. When reporting mean α values, we only use those α values obtained in realizations for which the BIC indicated power-law scaling. We also calculated the finite-size corrected Akaike information criterion MÀkÀ1 which led to similar results (not shown). We determined F i for the range of interval sizes n ¼ [10, N=10], where N denotes the length of the time series, here amounting to 300⋅250 ¼ 7:5⋅10 4 samples.

Appendix B. Derivation of the phase dynamics
Appendix B.1. Wilson-Cowan phase dynamics The first model was derived from the Wilson-Cowan firing rate model in a regime where it displayed self-sustained oscillations. With E k and I k denoting firing rates of the excitatory and inhibitory populations, respectively, the dynamics read (Wilson and Cowan, 1972;Daffertshofer and van Wijk, 2011): Here μ k represents the neuronal membrane time constant, which slightly differed through k to introduce heterogeneity in the individual mass dynamics. The average value hμ k i ¼ 0:15 was chosen such that the characteristic oscillation frequency denoted by Ω fell in the alpha frequency band. The parameters c EI ; c IE denote coupling between the inhibitory and excitatory units within pair k with c EE ; c II indicating self-coupling within E k and I k respectively. The matrix S kl is a normalized DTI derived structural connectivity matrix describing the coupling between excitatory units in the network. It is scaled by overall coupling strength K. The parameter q k is an external input to node k. The sum of all input contributions is integrated in time when it exceeds the threshold θ E ;θ I . This integration is performed through the sigmoid function Q ½x. Note that both the external inputs q k and the contributions from the rest of the network mediated by S kl only act on the excitatory units, such that inhibitory feedback is only present locally (within pair k).
In deriving the phase dynamics of the Wilson-Cowan network we considered deviations from the fixed points E ð0Þ k , I ð0Þ k , that is E k ¼ E ð0Þ k þ δE k and I k ¼ I ð0Þ k þ δI k . Upon inserting this in (1) and combining this with the Taylor expansion for Q we obtain for the dynamics of the deviations, after restricting ourselves to n ¼ 1: Here the prime in Q ' denotes the first derivative of Q . By assuming oscillatory behavior close to the Hopf bifurcation point, which we guaranteed by an appropriate choice of parameter values, we could transform the system into Jordan real form and subsequently define polar coordinates: δE k ðtÞ ¼ A k cosðΩt þ ϕ k Þ and δI k ðtÞ ¼ A k sinðΩt þ ϕ k Þ. To derive the phase dynamics _ ϕ k we employed a combination of a slowly varying wave approximation and slowly varying amplitude approximation (Haken, 2004), or in brief averaging (Guckenheimer and Holmes, 2013). As discussed in Daffertshofer andvan Wijk (2011), Ton et al. (2014) this boils down to identifying two distinct time scales in the system: the oscillations with characteristic frequency Ω and the much slower dynamics of ϕ k . By this separation we could average over one period of the fast oscillations and only retain the slow dynamics _ ϕ k . For a detailed derivation we refer to Daffertshofer and van Wijk (2011), here we suffice with stating that the above approximations yielded equation (4) combined with (5). Note that from the corresponding Jacobian Due to the first order firing rate dynamics for each node, the phase dynamics (4)/(5) and in particular the matrix D WC kl is N Â N with N the number of excitatory populations. Hence, the structural connectivity matrix S kl and D WC kl have equal dimensions. Parameter values were chosen in such a way that the neural mass network (3) displayed alpha band oscillations and amounted to a E ¼ 1, a I ¼ 1, The other neural mass network that served as a basis for a phase reduction was a network of neural masses introduced by Freeman (1975), describing the dynamics of mean membrane potentials of neural populations as where k ¼ 1;…;2N. The parameters α k and β k represent mean rise and decay times of neural responses in population k. Before weighted by the coupling  Fig. 2, the delayed signals V l ðt À τ kl Þ are thresholded by the sigmoid function Q ½⋅. This function covers the lump sum effect of pulse coupled neurons l (Freeman, 1975;David et al., 2006), where the fraction in its argument may be interpreted as the cumulative distribution function of the normal distribution N ðV À θ; σ 2 Þ of the firing thresholds θ across the population. In (3) it is additionally scaled by a factor σ (Marreiros et al., 2008) to guarantee self-sustained oscillations.
To derive the expressions covering the phase dynamics of the Freeman neural mass network we used a very similar procedure as for the Wilson-Cowan model above. To this end we first re-cast the system (3) as a system of two first-order equations Analogous to the derivation of the Wilson-Cowan model we considered the dynamics of the deviations from the fixed points V k together with the Taylor expansion of Q ½⋅ given by (B.1) yielding These expressions can be considered equivalent to (B.2) and the subsequent procedure to obtain the phase dynamics also goes along the same lines. Thus by transforming ½δV k ; δU k →½A k cosðΩt þ ϕ k Þ À ΩA k sinðΩt þ ϕ k Þ and by applying the same approximations as in the derivation of the Wilson-Cowan phase dynamics, the phase dynamics for the Freeman model yielded the expressions given in (4) and (6). For a more detailed derivation we refer to Ton et al. (2014). For this model the parameter values amounted to c EI ¼ c IE ¼ 1, c EE ¼ c II ¼ 0, q k ¼ 20, θ ¼ 15, γ ¼ 250, α k 2 ½60; 80 and β k 2 ½165; 185. The parameters α k , β k were chosen randomly to introduce heterogeneity in the oscillation frequency in the network. Values were chosen such that self-sustained oscillations in the alpha band were guaranteed. Appendix B.3. A note on the phase reduction of neural mass models The reduction of high-dimensional and complex neuronal networks into phase oscillator networks is a field of ongoing theoretical research, see, e.g., (Ashwin et al., 2016;Hoppensteadt and Izhikevich, 1997;Ermentrout and Terman, 2010), but satisfactory answers that are both simple and at the same time mathematically sound are still sought after. The here presented method combining a slowly varying wave and slowly varying amplitude approximation (Haken, 2004), or in brief averaging (Guckenheimer and Holmes, 2013), clearly stands out for its direct transformation of the original parameters in the physiologically realistic neural mass models into the corresponding phase oscillator parameters. Using appropriate polar coordinates, the complex oscillatory dynamics of the underlying neural mass model is cast into a Kuramoto-Sakaguchi phase model. This reduction into the Kuramoto model (4) is inherent to our method, and does not need any heuristic argumentation why particular terms of the Fourier expansion can be discarded as in Schuster and Wagner (1990), Sadilek and Thurner (2015). The question that naturally arises is whether our phase reduction approach indeed captures a holistic picture of the true phase dynamics of the neural mass models.
A single node-specific Wilson-Cowan dynamics combined with the small coupling strengths K < 0:8 considered throughout the article calls for applying here the theory of weakly coupled neural networks (Hoppensteadt and Izhikevich, 1997). For small q k , two necessary conditions are fulfilled: First, each uncoupled Wilson-Cowan node is near a bifurcation through which oscillations are generated. Second, the coupling is considerably weak. To be more precise, for q k % À 0:35, each oscillator is near a supercritical Hopf bifurcation. With respect to the theory of weakly coupled neural networks, there is extensive literature on how to phase reduce coupled oscillatory systems, with each one close to a Hopf bifurcation (Schuster and Wagner, 1990;Hoppensteadt and Izhikevich, 1997;Sadilek and Thurner, 2015;Pietras and Daffertshofer). The reduction is usually achieved in two steps: a center manifold reduction followed by a phase reduction. Once the dynamics are cast onto the center manifold, leading to the Hopf normal form, phase reduction techniques aim at estimating the corresponding phase response curves prior to obtaining the actual phase interaction, or coupling function, which determines the phase dynamics. Interestingly, the method used in our article leads to exactly the same phase dynamics when starting from the Hopf normal form, see also (Pietras and Daffertshofer). However, properly deriving normal forms is usually a highly non-trivial task, and often circumvented with heuristic arguments. One method which can be applied straight-forwardly to stable-limit cycle oscillations away from generic bifurcations, is based on Malkin's theorem and implemented in Ermentrout's software XPP (Ermentrout, 2002;Pietras and Daffertshofer). A more thorough investigation contrasting different phase reduction techniques will be published elsewhere (Pietras and Daffertshofer). Here we note that while numerical phase reduction techniques can serve as a good reference to test for the validity of the parameterized phase models, their numerical strength also comes with the disadvantage that no proper mapping from the parameters of the original model into those of the phase reduced model is available. Indeed, for small q k near the Hopf bifurcation, we numerically found the coupling term to have a negligible cosine component in the non-delayed systemin good agreement with our reduction method, cf. (5).
All these considerations also apply to the Freeman neural mass model (3). Henceforth, we believe that up to now our method appears to be the best compromise that maps the physiologically realistic parameters in the Freeman model into an appropriate phase model.

Appendix C. A slightly alternative model
We expected the Wilson-Cowan model and the Freeman model to resemble similar synchronization behavior, because in our structural connectivity matrix the inhibitory connections were rather sparse compared to excitatory ones (Fig. 2). To test this, we considered an alternative: a Freeman model (FE) that only comprised the excitatory part D F kl , i.e. the left upper block of this matrix. The resemblance of the Freeman phase dynamics with a repulsively coupled network was confirmed by the FE model results. Fig. C.1 shows that the behavior of the FE model largely agreed with the Freeman phase model, despite the lack of inhibitory connections. Functional connectivity correlations with data were very similar in both models, as can be appreciated by comparing Fig. C.1a and Fig. 4, with maximal correlations of ρ ¼ 0:56 for K ¼ 0:8 and hτ kl i ¼ 9:4 in both models. The FE model also yielded scale-free autocorrelations in the phase synchronization dynamics (Fig. C.1b). The average scaling exponent was α ¼ 0:56 AE 0:02, which was almost equal to the Freeman phase model result. Again this was significantly different from surrogate values (p < 10 À4 ). The hP ðsimÞ i and hRðtÞi results in the FE model largely agreed with the Freeman model results as well. In particular, the FE model dynamics resulted in a spatially desynchronized network with, again, hRðtÞi significantly below surrogate values (p < 10 À4 ); see Fig. C.1c. Similar to the strong coupling case in the main text, the direct connection of the FE model with the underlying neural mass dynamics (3) was blurred by ignoring the inhibitory nodes. The results suggest, however, that the origin of the scale-free correlation structure displayed by the Freeman phase model is, in general, caused by nature of the coupling between the excitatory units. Its dynamics can thus be understood by considering the phase dynamics (4)/(6) as a repulsively coupled Kuramoto network.

Appendix D. Correlating functional and structural connectivity
We computed the correlation between P ðsimÞ kl and S kl to show that both models generated functional connectivity structures that were highly correlated with S kl ; see Fig Correlation values were averaged over ten realizations for each parameter combination. Note the similarity of this figure with Fig. 4 suggesting that the reflection of S kl is an important determinant in high functional connectivity correlations, cf. (Robinson, 2012).