The Hopf whole-brain model and its linear approximation

Whole-brain models have proven to be useful to understand the emergence of collective activity among neural populations or brain regions. These models combine connectivity matrices, or connectomes, with local node dynamics, noise, and, eventually, transmission delays. Multiple choices for the local dynamics have been proposed. Among them, nonlinear oscillators corresponding to a supercritical Hopf bifurcation have been used to link brain connectivity and collective phase and amplitude dynamics in different brain states. Here, we studied the linear fluctuations of this model to estimate its stationary statistics, i.e., the instantaneous and lagged covariances and the power spectral densities. This linear approximation—that holds in the case of heterogeneous parameters and time-delays—allows analytical estimation of the statistics and it can be used for fast parameter explorations to study changes in brain state, changes in brain activity due to alterations in structural connectivity, and modulations of parameter due to non-equilibrium dynamics.

Whole-brain models are coupled stochastic dynamical systems in which nodes (i.e., brain regions) interact through couplings that represent anatomical connections estimated using diffusion imaging 1 , fiber tracing techniques 2 , or generative rules-such as the exponential distance rule 3 .Whole-brain models have proven to be useful to understand the emergence of correlations between neural populations or brain regions (or functional connectivity), as well as their spectral properties, in different brain states.In general, the ingredients of these models are a connectivity matrix between nodes, local node dynamics, noise, and, eventually, time-delays.Multiple choices for the local dynamics have been used depending on the studied behavior (e.g., network correlations, synchrony, metastability, etc.) and the data to be modelled (e.g., fMRI or M/EEG).Local node dynamics have been previously modelled using spiking networks 4 , conductance-based dynamics 1 , neural population dynamics 5 , neural mass models 6 , excitable systems 7 , phase oscillators [8][9][10] , and nonlinear oscillators 11 .In the present study, we examined the behavior of a network of nonlinear oscillators corresponding to a normal form of a supercritical Hopf bifurcation.This network model, first introduced by Matthews and Strogatz 12 to study collective behavior, is known as the Stuart-Landau model.It is a canonical model to study systems of coupled oscillators for which both the phase and the amplitude interact.The Stuart-Landau network has been used in diverse applications, from the study of coupled lasers 13 to neural networks 14 .In the context of neuroscience, this model is often referred as the Hopf model.In this model, as nonlinearities increase, isolated nodes transit through two qualitatively different dynamics: from damped oscillations to self-sustained oscillations.
The Hopf model has been used to study the link between brain structure and dynamics in resting-state conditions 15 and in different brain states, such as sleep 16 , low-level states of consciousness [17][18][19][20][21] , and psychedelic states 22 .Moreover, the Stuart-Landau model has been used to study the emergence of remote synchronization in human cerebral cortex 23 .Theoretical works have revealed sophisticate nonlinear emergent phenomena in the Stuart-Landau network such as oscillation and amplitude death 24,25 .Nevertheless, comparison of whole-brain models with resting-state neuroimaging data showed that the network operates in the simpler noisy-oscillation regime, suggesting that nonlinearities are small 15,26 .As we showed below, this case allows to strongly simplify the model to estimate the network statistics.This is important because, the Hopf model being a system of coupled stochastic differential equations, estimation of the network statistics (e.g., variances and covariances) requires extensive numerical simulations 27 , making often unpractical the exploration of a large part of the model's parameter space.
Here, we reviewed the Hopf model and derived network statistics using its linear approximation.The linearization allows analytical estimation of the statistics and can be used for fast parameter explorations without the need of extensive simulations.In order to facilitate future research, we have made the Matlab codes freely available online, allowing to perform the calculations for any connectome and for a large space of model parameters.
where η x and η y are uncorrelated white noises added to the real and imaginary parts, respectively.The variable z can also be written in polar coordinates, i.e., z = re iθ , where r = |z| = x 2 + y 2 1/2 is the module of z and θ = arctan y/x is its phase.Note that r ≥ 0 .In polar coordinates, we have r ṙ = x ẋ + y ẏ and r 2 θ = x ẏ − y ẋ .Thus, in absence of noise, Eq. (1) becomes: Equation (4) indicates that the phase evolves independently of r as θ(t) = ωt + ϕ , where ϕ is a constant phase.Clearly, a fixed point of Eq. ( 5) is r = 0 for which dr dt = 0 .The stability of the fixed point r = 0 depends on the parameter a , since deviations from r = 0 grow (i.e., dr dt > 0 ) if ar − r 3 > 0 and decrease (i.e., dr dt < 0 ) if ar − r 3 < 0 (Fig. 1A).For a < 0 , the solution r = 0 is stable as fluctuations around this point are attenuated.The eigenvalues of the system (2)-(3) are complex conjugates and equal to = a ± iω .For a < 0 , both eigenvalues have negative real part, indicating that the system relaxes to z = 0 with damped oscillations (see Fig. 1B), i.e., a spiral or focus solution.Note that, in this regime, addition of noise induces oscillations of the system.On the contrary, if a > 0 , r = 0 is unstable as fluctuations around it are amplified (Fig. 1A).In this latter case, a new fixed point appears given by r = a 1/2 , which is stable since fluctuations around it, r = a 1/2 + δr , are increased if δr < 0 , but decreased if δr > 0 .This solution is called a limit-cycle for which the system produces self-sustained oscillations with a constant amplitude and a constant angular frequency ω (see Fig. 1C).
In studies of whole-brain models, the brain signals (e.g.fMRI or MEG) are modelled by the real part of the state variables, i.e., x = Re(z).

Network model
The whole-brain dynamics are obtained by coupling the local dynamics of N Hopf nodes interconnected through a given coupling matrix C representing anatomical connections (Fig. 1C).In this study, to illustrate the method, we used a publicly available human diffusion MRI (dMRI) connectome from the Human Connectome Project (HCP) as the coupling matrix ( C ).The state variables of the network are given by the system of stochastic coupled nonlinear differential equations: where g (in s −1 ) represents a global scaling of the connectivity C and η j is uncorrelated white noise, i.e., �η j (t)� = 0 and �η j (t)η k (t ′ )� = σ 2 δ(t − t ′ )δ jk .Two versions of this model have been studied previously: the homogenous case for which the local bifurcation parameter is constant across nodes (i.e., a j = a) 15,16,28 and the heterogene- ous case for which nodes can have different local bifurcation parameters a j estimated from the data 15,20 .In both cases, ω j are estimated from the peak frequency of the data.
This model can be interpreted as an extension of the Kuramoto model to the case in which both the phase and the amplitude of the oscillators are allowed to vary and interact.In particular, the choice of the coupling function z k − z j promotes phase synchronization between coupled nodes.This can be seen by writing the deterministic system in polar coordinates: (1) represents a version of the Kuramoto model of phase oscillators for which couplings are modulated by the ratio of the amplitudes.The term sin θ k − θ j favors synchronization of nodes j and k , since an oscillator lagging behind another one ( θ k − θ j > 0 ) is sped up (a positive term sin θ k − θ j is added), whereas an oscillator leading another ( θ k − θ j < 0 ) is slowed down (a negative term sin θ k − θ j is added).In the case where the oscillations of the nodes are self-sustained (limit-cycles) and the couplings are weak, amplitude fluctuations are little compared to phase changes, and the system can be approximated by a Kuramoto model of phase oscillators interacting through couplings equal to C jk a k a j .In this study, however, we concentrated on the case of noisy oscillations (i.e., when nodes do not produce self-sustained oscillations).
Note also that the coupling function can have a stabilizing effect, since Eq. ( 6) without noise can be written as: żj = a j − gS j + iω j z j − z j 2 z j + g N k=1 C jk , where S j = N k=1 C jk is the strength of node j .In the case of S j > 0 , which is true in particular for positive connections C jk , the term −gS j < 0 contributes to the stability of the network.

Linear approximation
Estimating the network statistics (e.g., the covariance matrix) of the system given by Eq. ( 6) requires long stochastic simulations, impeding the exploration of different model parameters.However, in the case of weak noise (7)   www.nature.com/scientificreports/and small non-linearities, one can estimate the statistics of the whole-brain network using a linear approximation that we describe in this section.
In the following, we use bold symbols to indicate column vectors and matrices.The dynamical system can be re-written in vector form as: ] T is the vector containing the strength of each node, i.e., S i = j C ij , and η = [η 1 , . . ., η N ] T represents a vector of uncorrelated noise.The symbol ⊙ is the Hadamard element-wise product, i.e., The superscript T denotes the transpose operator.
We studied the linear fluctuations δz around the fixed point z = 0 , which is the solution of dz dt = 0 (Fig. 2A).In the linearized system the higher-order terms (δz ⊙ δz)δz are discarded and only terms in the first-order in δz are kept.Using the real and imaginary parts of the state variables, the evolution of the linear fluctuations δu follows the stochastic linear equation: where the 2N-dimensional column vector δu = (δx, δy) = δx 1 , . . ., δx N , δy 1 , . . ., δy N T contains the fluctuations of real and imaginary parts.The 2N × 2N matrix A is the Jacobian matrix of the system evaluated at the fixed point: w h e r e F j = a j − x j 2 − y j 2 x j − ω j y j + g N k=1 C jk x k − x j f o r 1 ≤ j ≤ N ( r e a l p a r t s ) , a n d F j = a j − x j 2 − y j 2 y j + ω j x j + g N k=1 C jk y k − y j for N + 1 ≤ j ≤ 2N (imaginary parts).By evaluating the partial derivatives at the fixed point, the Jacobian matrix can be written as a block matrix: We considered the heterogenous model for which the parameters a and ω were drawn from normal distributions N (a 0 , �a) and N (ω 0 , �ω) , respectively, with means a 0 and ω 0 , and standard deviations a and �ω .The connectivity matrix C was given by the HCP structural connectivity in a parcellation with N = 1000 nodes (Schaefer parcellation).We numerically calculated the eigenvalues of the Jacobian matrix for different values of a 0 and the global coupling g (normalized by the 2-norm of the connectivity matrix C ) and we evaluated the stability of the origin.The origin is stable if Re( max ) < 0 , where max is the eigenvalue with largest real part.Note logarithmic scale in the x-axis.Grey: the origin is unstable, i.e., Re( max ) > 0 .Blue: the origin is stable, Re( max ) < 0 , and a j < 0 for all nodes.Light blue: the origin is stable, Re( max ) < 0 , and a j > 0 for at least one node.Parameters: a = 0.2 ; �ω = 0.1 × 2π .(B) Proportion of positive bifurcation parameters ( a j > 0 ), for g/�C� = 0.7.
, where diag(v) is the diagonal matrix whose diagonal is the vector v .As shown below, the Jacobian matrix determines the statistics of the linear system.Note that the Jacobian depends on all the parameters of the model.Given an initial condition δu(0) at t = 0 , the general solution of a stochastic linear system such as Eq. ( 10) is given by 29 : where W is an 2N-dimensional Wiener process, σ is the noise amplitude, and e tA is the exponential matrix defined as: where I is the identity matrix.The right-hand side of Eq. ( 13) is the sum of the deterministic behavior plus a stochastic integral representing the diffusion due to noise.
The linearization is only valid if the origin z = 0 is a stable solution of the system, i.e., if all eigen- values of A have negative real part.Note that, in complex representation, the Jacobian writes A = diag a + iω − gS + gC = diag(a + iω) − gL , where L = S − C is the Laplacian matrix of the network.It is known that the Laplacian matrix is positive semidefinite: the eigenvalues nonnegative and µ 1 = 0 30 .Let j be the eigenvalues of A , the origin is asymptotically stable if Re( max ) < 0 , where max is the eigenvalue with largest real part.In the case of homogeneous local bifurcation and intrinsic frequency parameters, i.e., diag(a + iω) = (a + iω)I , the eigenvalues of A relate to those of −gL and we have Re( max ) = a − gµ 1 = a .Thus, in this case, the origin is stable if a < 0 .For the heterogenous case, however, there is not a direct expression for Re( max ) which depends on the contribution of the matrices diag(a + iω) and −gL , and stability needs to be evaluated numerically.For the HCP coupling matrix and the heterogeneous case, we found that the stability of the origin fixed point increases as a function the global coupling g and that, for sufficiently large g , the origin is stable even if a j > 0 for some nodes (Fig. 2A).Indeed, for strong coupling and close to instability, the majority of nodes can have a j > 0 while the origin remains stable (Fig. 2B).In other words, the focus solutions of single nodes can be unstable by themselves, but are stabilized by network interactions-as observed in simpler oscillator networks 31 .

Network statistics: covariances
In the following, we derive the network statistics of the linear system.The network mean activity (first order statistic) is trivial since fluctuations around the origin z = 0 have null mean.A first interesting statistic is the covariance of the fluctuations around the origin, i.e., C v = �δuδu T � , where the superscript T denotes the trans- pose operator.For a stochastic linear system such as Eq. ( 10), the motion equation of the covariance matrix C v is given as: where Q n = �ηη T � is the covariance matrix of the noise.For uncorrelated noise, Q n is diagonal, i.e., Q n = σ 2 I .The derivation of Eq. ( 15) is based on Eq. (10) which can be formally written as: dδu = Aδudt + dW , where W is an 2N-dimensional Wiener process with covariance �dWdW T � = Q n dt .Using Itô's stochastic calculus, we get to calculate the evolution of the covariance d δuδu T .Since �δudW T � = 0 , taking the expectations and keeping terms in first order of the differential dt (since dt 2 can be made arbitrarily small), we obtain: The stationary covariance matrix can be obtained by solving dC v dt = 0 , which leads to the following algebraic equation: Equation ( 16) is an algebraic Lyapunov equation that has a unique solution provided that A is asymptoti- cally stable.The Lyapunov equation can be solved using the eigen-decomposition of the Jacobian matrix.Let A = V DV −1 , where D is a diagonal matrix containing the eigenvalues of A , denoted i , and the columns of matrix V are the eigenvectors of A .Multiplying Eq. ( 16) by V −1 from the left and by the conjugate transpose of V −1 , noted V − † , from the right we get: where the matrix M is given as: A fast, stable numerical solution of Eq. ( 16) can be obtained using the MatLab function lyap.m that uses the Bartels-Stewart method 32 based on the Schur decomposition of the matrix A.  www.nature.com/scientificreports/Moreover, knowledge of the Jacobian matrix and the stationary covariance gives the stationary lagged covariances of the state variables, defined as C v (τ ) = �δu(t + τ )δu(t) T � .Using the general solution of the system given by Eq. ( 13), we get: where C v (0) = C v is the covariance matrix (i.e., zero-lag).The lagged covariance has been used to described the temporal structure of whole-brain activity 33 .

Network statistics: power spectral densities
In the frequency domain, the power spectral density (PSD) of fluctuations around the fixed point is also determined by the Jacobian matrix.Taking the Fourier transform F of Eq. ( 10), we get: where δ u(ν) and η(ν) are the Fourier transforms of δu(t) and η(t) at frequency ν , respectively.Using the relation δ u = −(A + i2πνI) −1 η , we get the cross-spectrum of the linear fluctuations: The real part of the cross-spectrum (also called co-spectrum) represents the simultaneous covariance at frequency ν .Its imaginary part (called quadrature spectrum) is the covariance of time-series lagged by a phase π/2 at frequency ν .At each frequency ν , the PSDs of the nodes, φ j (ν) , are given by the diagonal terms of ψ(ν) and the coherence between nodes, γ jk (ν) , is given by the normalized cross-spectrum, i.e., γ jk (ν) = ψ jk (ν)/ φ j (ν)φ k (ν) 34 .For uncorrelated noise, the PSD is given as: The Fourier transform is also a useful tool to study the system in the case of time-delays.Consider the nonlinear Hopf network with delayed interactions: where τ jk represents the time-delay of the interaction between nodes j and k .For simplicity, one can assume that τ jk is given by the Euclidean distance between nodes j and k divided by a constant transmission velocity v .Delayed-interactions can be treated in the Fourier space, since the change of variable t ′ = t − τ jk leads to Using the linear approximation and the Fourier transform, we get: where Ŵ is the matrix containing the delays, i.e., Ŵ jk = τ jk , the elements of C ⊙ e i2πŴ are C jk e i2πτ jk , and B is the block matrix given by: .From the cross-spectrum ψ we can obtain the PSD of each node (i.e., diagonal terms), the lagged-covariances (i.e., the inverse Fourier transform of the cross-spectrum), and the covariance matrix C v by integrating the real part of ψ over frequencies: In summary, in the linear approximation, the stationary instantaneous and lagged covariance matrices, the cross-spectrum, and the PSDs of the model can be obtained through algebraic operations including the Jacobian matrix.This can be done both in the homogeneous and the heterogeneous cases, and also in the presence of time delays.

Comparison with stochastic simulations
We compared the predictions of the linear approximation against the statistics obtained using stochastic simulations of the nonlinear model.The coupling matrix was given by the human dMRI connectome from HCP, with N = 1000 nodes.The model parameters a and ω were drawn from normal distributions N (a 0 , �a) and N (ω 0 , �ω) , respectively, with means a 0 and ω 0 , and standard deviations a and �ω .We simulated the system for T = 3 min after letting it reach the stationary regime and we used n = 100 realizations of the system with different random initial conditions.
We used the linear approximation to study the fluctuations around the origin.We first examined the predictions of the linear approximation when the stability of the origin is strong ( Re( max ) < −1 ).In this case, the approximation accurately estimates the covariances (Fig. 3A), the auto-and cross-covariances (Fig. 3B,C), and the PSDs (Fig. 3D,E).To study the accuracy of the prediction as a function of the origin's stability, we varied the local bifurcation parameter a 0 in the homogeneous case (i.e., a = 0 ).This analysis, that requires to simulate the system for different parameters a 0 , was done using a subsampled of the network, with N = 250 nodes (see Methods).As the origin loses stability, nonlinear terms become non-negligible, it is thus expected that the linear approximation fails close to Re( max ) → 0 .We quantified the goodness of the prediction through two measures: (i) the R-squared value ( R 2 ) of the correlation between covariances obtained from numerical simulations C sim v and those obtained with the linear approximation C lin v , and ii) the relative error ( E ) between the matrices using www.nature.com/scientificreports/ the Frobenius norm: E = �C sim v ��C lin v �/�C sim v � (Fig. 3F).We found that the linear approximation accurately estimates the covariances ( R 2 > .99 and E < 0.1 ) for Re( max ) < −0.15.
We also evaluated the predictions of the linear approximation in the case of time-delays.The delay-coupled Hopf model has been recently studied using numerical simulations 35 .In this case, the interaction delays between nodes can be approximated using the Euclidean distance between brain regions divided by a transmission velocity v .Here, we used the distances from the HCP data, which yield an average distance between nodes equal to 79 mm.The intrinsic frequencies were chosen from a normal distribution centered on ω 0 2π = ν 0 = 1 Hz and with standard deviation equal to �ω 2π = 0.2 Hz.For this example, we chose a transmission velocity v such that the aver- age transmission delay D is of the same order of the average intrinsic period of the network, i.e., D ∼ ν 0 −1 .The parameters a were drawn from the normal distribution N (a 0 , �a) , with a 0 = −1 and a = 0.3.As previously, we simulated the system for T = 3 min after letting it reach the stationary regime and we used n = 100 realizations of the system with different random initial conditions.The linear approximation accurately approximates the PSD of the nodes (Fig. 4A).Moreover, integration of the cross-spectrum, obtained using the linear approximation, gives an accurate prediction of covariances (Fig. 4B).

Parameter exploration and data fitting
Finally, we studied how well the linear approximation predicts the correlations of resting-state (rs-) fMRI signals.
For this, we analyzed rs-fMRI signals from the HCP, from 1003 participants.First, we calculated the correlation matrix (or functional connectivity, FC) averaged over participants, in the parcellation with N = 1000 nodes.Sec- ond, we computed the FC for the heterogenous linearized Hopf model constraint by the HCP dMRI connectivity matrix.Finally, we compute the correlation between FC matrices obtain from the data and from the linearized Hopf model.The model parameters a and ω were drawn from normal distributions N (a 0 , �a) and N (ω 0 , �ω) , respectively, with a = 0.2 , ω 0 = 2π , and �ω = 0.1 .Note that, here, the local parameters a j and ω j were taken from normal distributions and were not fitted/optimized using the data as in previous work 20 .We evaluated the fitting of the empirical FC in the parameter space a 0 , g , for varying mean local bifurcation parameter and global coupling (Fig. 5).We found that, for this particular example, the best fit of the FC was obtained when the coupling was high enough with respect to the norm of the connectivity matrix (i.e., g �C� ∼ 1-100).The fitting values are similar to what was found with previous numerical simulations with the same parcellation 36 .In that previous work, however, long-range connections were added to the connectivity, which improve the fit.Also, in previous studies 15,16 , a narrow band-pass filter was applied to the fMRI signals, thus making the signals strongly oscillatory, which might explain the fit increase close to the onset of self-sustained oscillations.

Discussion
Using a linear approximation, we have derived network statistics of the Hopf whole-brain model.The linearization allows analytical estimation of the stationary instantaneous and lagged covariance matrices, the crossspectrum, and the PSDs of the model.This can be done in the most general form of the model, namely in the delay-coupled heterogeneous case.The linearization provides good estimates of these quantities as soon as nonlinear terms do not dominate (as it is the case sufficiently close or beyond the bifurcation).This occurs when the origin is stable.Exploration of the parameter space, for which the origin destabilizes and dynamics are strongly nonlinear, could be treated using approximations more sophisticated than the linear approximation, for example, using higher-order phase reduction 37 .
Synchronization among brain regions has been studied in multiple previous studies using different neuroimaging techniques 10,18,[38][39][40][41][42] .The present model is a canonical model to describe, at a phenomenological level, the synchronization of oscillators with phase and amplitude interactions, previously used to study large-scale brain dynamics 15-18, 20, 22 .However, the neuronal/synaptic mechanisms underlying the brain's large-scale synchronization are not fully understood.Noisy oscillations around a fixed point can be understood using more realistic, yet still simple, models composed of interconnected excitatory and inhibitory neural populations such as the Wilson-Cowan model 43 or the stabilized supralinear network 44 , for which linear fluctuations can be studied in light of the biological interpretation of the different parameters.It is worth noting that the linear fluctuations around the fixed point are rich in structure, as shown here by their structured covariance and cross-spectrum which are determined by local dynamics, network interactions, network stability, time-delays, and noise propagation.Even richer dynamics could emerge in the case of strongly nonlinear dynamics, which might be the subject of future research.
There are several applications of the present framework.The estimated network statistics can be used to track changes in the brain state, e.g., in the case of low-level states of consciousness [17][18][19][20][21] , anesthesia 20,45 , sleep 16 , etc., or to evaluate the effect of lesions in the connectome [46][47][48][49] .
We here tested the model predictions using rs-fMRI data, but the model can be used to approximate MEG data in different frequency bands 28,35 .The slow time scale of fMRI signals allows to neglect the effect of conduction delays between the different brain regions, which are orders of magnitude faster-tens of milliseconds 50 -than the periods of the model oscillators, and treat the interactions as instantaneous.In the case of MEG data, however, delayed interactions can have an important effect for sufficiently fast frequency bands.Thus, the linear approximation of delay-coupled Hopf whole-brain model derived here can represent a valuable tool to study the PSDs and cross-spectrum of MEG, which are well-established methods for FC analysis in the frequency domain [51][52][53][54] .
Furthermore, recent studies suggest that dynamics out of equilibrium are relevant to describe the wholebrain [55][56][57] .The present model can be used to track non-stationarities by assuming that changes in parameters are sufficiently slow relative to the time it takes for the system to reach equilibrium 58 .In this way, using the linear model to fit the stationary statistics of the system measured in short time windows, it is possible to infer the change in network parameters over time.
Finally, since the goal of the present study was to derive the linear statistics of the model, rather than fitting functional data, we used models for which the local parameters were not estimated from the data, opposite to previous studies 15,20 .Future research could combine the present linear approximation with algorithms to optimize the parameters of the model, such as the N local bifurcation parameters.This can be achieved using genetic algorithms applied to infer optimal local parameters 59 , allowing to compare the learned parameters in different brain states, disorders, or across aging, for example.The use of the linear approximation would allow to estimate the parameters for large N.
In all the above applications, one would need to systematically verify that the origin of the model is a stable fixed point and that the real part of the leading eigenvalue does not approach zero.

Neuroimaging ethics
The Washington University-University of Minnesota (WU-Minn HCP) Consortium obtained full written informed consent from all participants to study procedures and data sharing outlined by HCP, and research procedures and ethical guidelines were followed in accordance with Washington University institutional review board approval.

Functional MRI data
In this study we analyzed publicly available rs-fMRI data from the Human Connectome Project (HCP), from 1003 participants.The participants were scanned on a 3 T connectome-Skyra scanner (Siemens).The rs-fMRI data was acquired for approximately 15 min, with eyes open and relaxed fixation on a projected bright cross-hair on a dark background.The HCP website (https:// www.human conne ctome.org/) provides the details of participants, the acquisition protocol and preprocessing of the functional data.

Parcellation
Schaefer and colleagues created a publicly available population atlas of cerebral cortical parcellation based on estimation from a large data set (n = 1489) 60 .They provide parcellations of regions of interest (ROIs) available in surface spaces, as well as MNI152 volumetric space.We used the Schaefer parcellation with 1000 areas and estimated the Euclidean distances from the MNI152 volumetric space 60 and extracted the timeseries from HCP using the surface space version.Finally, for the analysis presented in Fig. 3F, we subsampled the connectivity by choosing only 250 ROIs.This allowed us to simulate the stochastic nonlinear dynamical system for a large amount of repetitions, initial conditions, and varying parameters.

Structural connectivity using dMRI
Structural connectivity was estimated from diffusion spectrum and T2-weighted imaging data from 32 participants from the HCP database, scanned over 89 min.Acquisition parameters are described in detail in the HCP website 61 .The freely available Lead-DBS software package (http:// www.lead-dbs.org/) provided the preprocessing which is described in detail in Horn and colleagues 62 .Standardized methods in Lead-DBS were used to produce the structural connectomes for the Schaefer parcellation scheme 60 .The connectivity weight C ij = C ji was given by the number of fibers connecting two brain regions.To have values between 0 and 1, we normalized the weights by dividing them by the largest value, i.e., max(C).

Statistics and reproducibility
The goodness of the linear prediction of fMRI FC was given by the Pearson correlation between the vectorized FC averaged over all subjects and the model FC for all combinations of parameters a 0 , g (Fig. 5).Stochastic numerical simulations were performed using Euler's method, with a simulation step size equal to 0.001 s and 0.005 s in the absence and presence of delays, respectively.The system was simulated for T = 3 min after letting it reach the stationary regime after 20 s; the stochastic simulations were repeated n times with differ- ent random initial conditions ( n = 100).For the subsampled system of 250 nodes (Fig. 3F) we used: T = 10 min and n = 200.The PSDs of simulated time-series were estimated using the fast Fourier transform.MATLAB (R2021a) software was used to perform all analyses and to simulated the model.Numerical simulations were performed in a 50-nodes computer cluster (Intel® Xeon® E5-2684 at 2.1 Ghz, 256 GB RAM, 1 TB disk).

Figure 1 .
Figure 1.Hopf model: single-node and network dynamics.(A)The fixed points of a Hopf node have modules which are the roots of ṙ = ar − r 3 .For a < 0 , the solution r = 0 is stable since deviations from r = 0 are attenuated (i.e., ṙ < 0 ).On the contrary, if a > 0 , r = 0 is unstable as fluctuations around it are amplified (i.e., ṙ > 0 ).In this latter case a new fixed point appears given by r = a 1/2 , which is stable since fluctuations around it, r = a 1/2 + δr , are increased if δr < 0 , but decreased if δr > 0 .The arrows indicate the direction of flow and are given by the sign of ṙ .(B) Single-node dynamics for a < 0 .The system relaxes with damped oscillations from the initial condition (white circle) to the origin of the complex plane.Insets: top: in the absence of noise ( η = 0 ) the oscillations die out; bottom: in the presence of noise ( η = 0 ) the oscillations are noise-driven. (C) Single-node dynamics for a > 0 .The system produces self-sustained oscillations.Insets: top, deterministic system; bottom, stochastic system.(D) Network model.The whole-brain network is composed of N Hopf nodes interconnected through anatomical connections.Here, we used dMRI connectivity from the Human Connectome Project (HCP), in a parcellation with N = 1000 nodes. (E) Example dynamics for five nodes of the network.Parameters: a j = −0.5 (homogeneous); g = 1 ; ω j = 10 rad.s -1 ; σ = 0.3.

Figure 2 .
Figure 2. Linear stability of the origin.(A)We considered the heterogenous model for which the parameters a and ω were drawn from normal distributions N (a 0 , �a) and N (ω 0 , �ω) , respectively, with means a 0 and ω 0 , and standard deviations a and �ω .The connectivity matrix C was given by the HCP structural connectivity in a parcellation with N = 1000 nodes (Schaefer parcellation).We numerically calculated the eigenvalues of the Jacobian matrix for different values of a 0 and the global coupling g (normalized by the 2-norm of the connectivity matrix C ) and we evaluated the stability of the origin.The origin is stable if Re( max ) < 0 , where max is the eigenvalue with largest real part.Note logarithmic scale in the x-axis.Grey: the origin is unstable, i.e., Re( max ) > 0 .Blue: the origin is stable, Re( max ) < 0 , and a j < 0 for all nodes.Light blue: the origin is stable, Re( max ) < 0 , and a j > 0 for at least one node.Parameters: a = 0.2 ; �ω = 0.1 × 2π .(B) Proportion of positive bifurcation parameters ( a j > 0 ), for g/�C� = 0.7.

Figure 3 .
Figure 3.Comparison with numerical simulations.(A) Comparison between variances and covariances obtained using numerical simulations and the linear approximation.The black line indicates the identity line.(B,C) Autocovariances (B) and lagged covariances (C) for numerical simulations (black trace) and the linear approximation (red dotted trace) for three example nodes (B) and pairs of nodes (C).(D) PSD for six example nodes and their linear predictions (solid lines).The frequency was normalized by the average intrinsic frequency ν 0 = ω 0 /(2π) . (E) Comparison between the peak frequencies (normalized by ν 0 ) obtained using numerical simulations and the linear approximation.The black line indicates the identity line.Model parameters for panels (A-E): a 0 = −1 ; a = 0.3 ; g = 3 ; ω 0 = 2π ; �ω = 0.2 × 2π ; σ = 0.01 .(F) Accuracy of the prediction for different values of Re( max ) .The origin is stable for Re( max ) < 0 .We quantified the goodness of the prediction through the R-squared value ( R 2 ) of the correlation between covariances obtained from numerical simulations and those obtained with the linear approximation.In the analysis presented in panel (F) we used a subsample of the network, i.e., N = 250 nodes.Model parameters: a = 0.3 ; g = 3 ; ω 0 = 2π ; �ω = 0.2 × 2π ; σ = 0.001.

Figure 5 .
Figure 5. FC prediction in parameter space.Correlation between FC matrices obtain from the data and the linearized Hopf model, for varying mean local bifurcation parameter and global coupling.Grey: the origin is unstable, i.e., Re( max ) > 0 .Between the horizontal line and the grey zone, the nodes can have a j > 0 while the origin remains stable.Note the logarithmic scale of the x-axis. https://doi.org/10.1038/s41598-024-53105-0