The role of node dynamics in shaping emergent functional connectivity patterns in the brain

The contribution of structural connectivity to functional brain states remains poorly understood. We present a mathematical and computational study suited to assess the structure–function issue, treating a system of Jansen–Rit neural mass nodes with heterogeneous structural connections estimated from diffusion MRI data provided by the Human Connectome Project. Via direct simulations we determine the similarity of functional (inferred from correlated activity between nodes) and structural connectivity matrices under variation of the parameters controlling single-node dynamics, highlighting a nontrivial structure–function relationship in regimes that support limit cycle oscillations. To determine their relationship, we firstly calculate network instabilities giving rise to oscillations, and the so-called ‘false bifurcations’ (for which a significant qualitative change in the orbit is observed, without a change of stability) occurring beyond this onset. We highlight that functional connectivity (FC) is inherited robustly from structure when node dynamics are poised near a Hopf bifurcation, whilst near false bifurcations, and structure only weakly influences FC. Secondly, we develop a weakly coupled oscillator description to analyse oscillatory phase-locked states and, furthermore, show how the modular structure of FC matrices can be predicted via linear stability analysis. This study thereby emphasises the substantial role that local dynamics can have in shaping large-scale functional brain states.


ABSTRACT
The contribution of structural connectivity to dynamic and often highly variable brain states remains poorly understood.We present a mathematical and computational study suited to assess the structure-function issue.We treat a system of Jansen-Rit neural-mass nodes with heterogeneous structural connections estimated from diffusion MRI data provided by the Human Connectome Project.Direct simulations are used to determine the similarity of functional (inferred from correlated activity between brain areas) and structural connectivity matrices as a function of the parameters controlling dynamics of a single node, highlighting a non-trivial structure-function relationship in regimes that support limit cycle oscillations.To determine their relationship, we firstly determine the network instabilities that give rise to oscillations, and the set of false bifurcations that occur beyond this onset.In particular, we highlight that functional connectivity (FC) is inherited most robustly from structure when node dynamics are poised near a Hopf bifurcation, whilst near false bifurcations, structure only weakly influences function.Secondly, we develop a weakly coupled oscillator description to analyse oscillatory phase-locked states and, furthermore, show how the modular structure of the FC matrix can be predicted using a linear stability analysis.This study thereby emphasises that local dynamics can have a substantial role in shaping large-scale functional brain states.

AUTHOR SUMMARY
Patterns of oscillation across the brain arise because of structural connections between brain regions.However, the type of oscillation at a site may also play a contributory role.We focus on an idealised model of a neural mass network, coupled using estimates of structural connections obtained via tractography on Human Connectome Project MRI data.Using a mixture of computational and mathematical techniques we show that functional connectivity is inherited most strongly from structural connectivity when the network nodes are poised at a Hopf bifurcation.However, beyond the onset of this oscillatory instability a phase-locked network state can undergo a false bifurcation, and

INTRODUCTION
Driven in part by advances in non-invasive neuroimaging methods that allow characterisation of the brain's structure and function, and developments in network science, it is becoming increasingly accepted that the understanding of brain function may be obtained from a network perspective, rather than by exclusive study of its individual sub-units.Anatomical studies using diffusion MRI allow estimation of structural connectivity (SC) of human brains, forming the so-called human connectome (Sporns, 2011;Van Essen et al., 2013) which reflects white matter tracts connecting large-scale brain regions.The graphtheoretical properties of such large-scale networks have been well studied, in both health and disease, highlighting key features including small-world architecture (Bassett & Bullmore, 2006;He, Chen, & Evans, 2007), hub regions and cores (van den Heuvel & Sporns, 2013), rich club organisation (Van Den Heuvel & Sporns, 2011), a hierarchical-like modular structure (Z.J. Chen, He, Rosa-Neto, Germann, & Evans, 2008;Meunier, Lambiotte, & Bullmore, 2010), and economical wiring (Betzel et al., 2017;Bullmore & Sporns, 2012).The emergent dynamic brain activity that this structure supports can be evaluated by functional connectivity (FC) network analyses, that describe the statistical patterns of temporal coherence in neural activity between spatially distinct brain regions.Time-series data from which such functional networks are constructed may be obtained by a range of methods, including electroencephalography (EEG), magnetoencephalography (MEG), and functional magnetic resonance imaging (fMRI).These highly dynamic FC patterns are widely believed to be significant in integrative processes underlying higher brain function (Van Den Heuvel & Pol, 2010;van Straaten & Stam, 2013) and, crucially, disruptions in both structural and functional brain networks have been shown to be associated with many psychiatric and neurological diseases such as epilepsy and schizophrenia (Braun, Muldoon, & Bassett, 2015;Menon, 2011).
However, though there is evidence to suggest certain network motifs organise functional connectivity (Adachi et al., 2011), the relationship between the brain's anatomical structure and the neural activity that it supports remains largely unknown (C.J. Honey, Thivierge, & Sporns, 2010;Park & Friston, 2013).In particular, the divergence between dynamic functional activity and the relatively static structural connections between populations is critical to the brain's dynamical repertoire and may hold the key to understanding brain activity in health and disease (Park & Friston, 2013), though current models have not yet been able to accurately simulate the transitive states underpinning cognition (Petersen & Sporns, 2015).Evidence from empirical studies suggests that while a direct structural connection between two brain areas is typically associated with a stronger functional interaction, strong interactions can nevertheless exist in their absence (Hermundstad et al., 2014;C. J. Honey et al., 2010); moreover, these functional networks are transient as the brain shifts between tasks (Fox et al., 2005).An important example of diverging SC-FC patterns is provided by resting-state networks, such as the so-called 'default mode network' and the 'core network' (Van Den Heuvel & Pol, 2010).These networks consist of distributed brain areas that are strongly functionally connected at rest (Van Den Heuvel & Pol, 2010), while it intuitively suggests the existence of direct structural connections between these areas; however, there is limited knowledge of how these interactions arise from the brain's neural structure (de Pasquale et al., 2012).Mišić et al. (2016) used multivariate statistical techniques to Network Neuroscience determine optimally associated structural and functional sub-networks.Importantly, the detected sub-networks typically consisted of functional relationships that arose from nonoverlapping sets of structural edges.
Theoretical studies deploying anatomically realistic structural networks obtained through tractography, alongside neural mass models describing mean-field regional neural activity have recently been used to further investigate the emergence of large-scale FC patterns (Breakspear, 2017;Deco et al., 2013;C. J. Honey, K ötter, Breakspear, & Sporns, 2007;Messé, H ütt, K önig, & Hilgetag, 2015;Ponce-Alvarez et al., 2015;Rubinov, Sporns, van Leeuwen, & Breakspear, 2009).These findings suggest that through indirect network-level interactions, a relatively static structural network can support a wide range of FC configurations.For example, a number of studies (C.Honey et al., 2009;C. J. Honey et al., 2007;Rubinov et al., 2009) have shown that functional networks reflect the underlying structural networks on slow time scales, but significantly less so on faster time scales.Moreover, simulated functional networks have been found most strongly to resemble the structural network when the dynamical system describing regional activity is close to a phase transition (Stam et al., 2016), and strong structure-function agreement is reported near Hopf bifurcations in Hlinka and Coombes (2012).These results providing a possible manifestation of the socalled critical brain dynamics hypothesis (Shew & Plenz, 2013).Furthermore, analysis of the dynamical systems underpinning neural simulations have shown to be a good fit to fMRI data when the system is near to bifurcation (Deco & Jirsa, 2012;Tewarie et al., 2018).In Crofts, Forrester, and O'Dea (2016), both structure and function are analysed together in a multiplex network, proposing a novel measure of multiplex structure-function clustering in order to investigate the emergence of functional connections that are distinct from the underlying structure.
In this paper, we present a combined computational and mathematical study, which significantly extends the related works of Hlinka and Coombes (2012) and Crofts, Forrester, and O'Dea (2016), to investigate the influence of structural connectivity, and the intrinsic behaviour of neural populations, on shaping FC networks.Thereby, we provide a complementary investigation to many of the aforementioned studies which focus on the analysis of brain networks themselves, by instead investigating the relationship between network structure and the emergent dynamics of these networks.Specifically, we consider synchrony between neural subunits whose dynamics are described by the neural mass model of Jansen and Rit (1995), and whose connectivity is defined by a tractography-derived structural network obtained from data in the Human Connectome Project (HCP) (Van Essen et al., 2013).Structure-function relations are interrogated by graph-theoretical comparison of functional and structural network topology under systematic variation of the model parameters associated with excitatory/inhibitory neural responses.Under these variations, we show that the degree to which structure drives function is related to the dynamical state of the nodes.Indeed, making use of techniques from bifurcation and weakly-coupled oscillator theory, we indicate how in some dynamical regimes the knowledge of nodal dynamics can explain the associated FC structure, without recourse to detailed structural information.In particular, we determine the way in which both true bifurcations and false bifurcations (for which one sees a qualitative change in the shape of an orbit, without a change of stability) contribute to the structure-function issue, and how they can drive this relationship in very different ways.When the node dynamics undergoes a Hopf bifurcation we see that function is largely inherited from structure, whilst for a false bifurcation the structural connectivity only weakly influences the functional connectivity.We show that this organisation of SC-FC correlation is reflected in the stability of phase-locked states in the weak-coupling Network Neuroscience limit and, moreover, that the eigenstructure of the linearised network model can be used to predict the modular structure of emergent FC networks.

Neural mass model
We consider a network of interacting neural populations, representing a parcellation of the cerebral cortex, such that each area (node) corresponds to a functional unit that can be represented by a neural mass model, and with edges informed by structural connectivity.Neural mass activity is represented by the Jansen-Rit model (Jansen & Rit, 1995) of dimension m = 6, that describes the evolution of average post-synaptic potential (PSP) in three interacting neural populations: excitatory, pyramidal, and inhibitory.The Jansen-Rit model is described by three second order ordinary differential equations, which are commonly rewritten as six first order equations by adopting the notation (y 0 , . . ., y 5 ) for the dependent variables.The pairs (y 0 , y 3 ), (y 1 , y 4 ), and (y 2 , y 5 ) are associated with the activities of the pyramidal, excitatory and inhibitory populations respectively.Introducing an index i = 1, . . ., N to denote each node in a network of N interacting neural populations, we write the evolution of state variables as: Here f is a sigmoidal nonlinearity, representing the transduction of activity into a firing rate, and with the specific form The model is completed by the specifying the network interactions between excitatory cells in nodes i and j, encoded in a connectivity matrix with elements w ij , with an overall scale of interaction set by ε, and described in Structural and functional connectivity.The remaining model parameters, together with their physiological interpretations, are given in Table 1.
The Jansen-Rit model, defined by equation ( 1), can support oscillations that relate to important neural rhythms, such as the well known alpha, beta and gamma brain rhythms, and also irregular, epileptic-like activity (Ahmadizadeh et al., 2018).Moreover, the model is able to replicate visually-evoked potentials seen in EEG recordings (Jansen & Rit, 1995), from which FC may be empirically measured (Srinivasan, Winter, Ding, & Nunez, 2007).
In what follows, we consider the patterns of dynamic neural activity that arise under systematic variation of the model parameters A and B, these being chosen as the parameters of interest because they govern the interplay between inhibitory and excitatory activity, which would typically vary due to neuromodulators in the brain (Rich, Zochowski, & Booth, 2018).The choice of the remaining model parameter values are taken from Grimbert andFaugeras (2006), andTouboul, Wendling, Chauvel, andFaugeras (2011).It is known that a single Jansen-Rit node can support multi-stable behaviour which includes oscillations of different amplitude and frequency but, moreover, a network of these nodes can  1) and ( 2), along with physiological interpretations and values/ranges used in simulations.In particular, the values A and B, which modulate the strength of excitatory and inhibitory responses respectively, were chosen as the key control parameters for varying network activity.

Network Neuroscience
also exhibit various stable phase-locked states.A small amount of Gaussian noise is added to the global drive P on each node, in order to allow the system to explore a variety of these dynamical states.For direct simulations of the network we use an Euler-Murayama scheme, implemented in Matlab ® , with a fixed numerical time-step of 10 −3 and a noise standard deviation of 10 −1 .

Structural and functional connectivity
The structural connectivity was estimated using diffusion MRI data for 10 subjects obtained from the HCP (Van Essen et al., 2013).Briefly, we explain how this data is post-processed to derive connectomic data, though we direct the reader to Tewarie et al. (2019) and the references therein for a more detailed overview.60,000 vertices on the white/grey matter boundary surface for each subject (Glasser et al., 2013) were used as seeds for 10,000 tractography streamlines.Streamlines were propegated through voxels with up to three fibre orientations, estimated from distortion-corrected data with a deconvolution model (Jbabdi, Sotiropoulos, Savio, Gra ña, & Behrens, 2012;Sotiropoulos et al., 2016), using the FSL package.The number of streamlines intersecting each vertex on the boundary layer was measured and normalised by the total number of valid streamlines.This resulted in a 60,000 node structural matrix, which was further parcellated using the 78-node AAL atlas.This was used to describe connections between brain regions, providing an undirected (symmetric), weighted matrix whose elements w ij define the strengths of the excitatory connections in equations (1).To generate nodal inputs with commensurate magnitudes, the structural connectivity matrix was normalised by row so that afferent connection strengths for each node sum to unity.Noting that the matrices for SC and FC quantify entirely different relationships, the former being the density of tractography streamlines and the latter being the temporal correlation of neural activity between regions, we compare network topologies by thesholding and binarising such that only the top third of the weights (ordered Network Neuroscience by strength) are retained.This reduces the number of connections while maintaining the overall modular structure (Fig. 1).
In view of the non-linear oscillations supported by the network model given by (1), functional networks are obtained by computing a measure of pairwise mean phase agreement (MPA) between the simulated activities of nodes.Note that this differs from the commonly used metric of mean phase coherence (Mormann, Lehnertz, David, & Elger, 2000), which determines correlation strength in terms of the proclivity of two oscillators to phase-lock.This is unsuitable in our case since we consider situations where stable phase-locking naturally arises, and use relative phase differences to determine patterns of coherence.We choose y i = y 1 i − y 2 i as the variable of interest because its physical interpretation is the average potential of the pyramidal population, which is the main contributor to signals generated in EEG recordings (Teplan et al., 2002).Pairwise MPA measures the temporal average of the phase difference ∆φ jk (t) = φ j (t) − φ k (t), between two time-series indexed by j and k, where here the instantaneous phase φ j (t) is obtained as the angle of the complex output resulting from application of a Hilbert transform to the time-series, y j (t).The mean phase agreement of the time-series comprising M time-points t l (l = 1, . . ., M) is defined as: The matrix, R jk , that describes the strength of functional connection between brain regions, is binarised (keeping the top third of the strongest connections for consistency with SC) to obtain an unweighted functional connectivity network.Structure-function relations are assessed by computing the Jaccard similarity coefficient (Niwattanakul, Singthongchai, Naenudorn, & Wanapu, 2013) of the non-diagonal entries of the binarised SC and FC matrices.This describes the relative number of shared pairwise links between the two networks, providing a natural measure of structure-function similarity, ranging from zero for matrices with no common links to unity for identical matrices.

Bifurcation analysis
A natural way to analyse the neural mass dynamics is to use bifurcation theory.Here we construct diagrams for both the single-node and full network case.

Single node and network bifurcations
Bifurcations for a single node are readily computed using the software package XPPAUT (Ermentrout, 2002), using A and B as the parameters of interest.The result is a Hopf and saddle-node set in parameter space, which bounds a region of oscillatory solutions.This is shown in Fig. 2.
The corresponding diagram for the full network requires numerical analysis of a much higher dimensional system, described by N × m = 78 × 6 = 468 ODEs; this is computationally demanding, and so in Mathematical methods we develop a quasi-analytic approach by linearising the full network equations around a fixed point.The resulting equations can be diagonalised in the basis of eigenvectors of the structural connectivity, leading to a set of N equations, each of which prescribes the spectral problem for an m-dimensional system.Thus, each of these low dimensional systems can be easily treated without recourse to high performance computing.Moreover, this approach exposes the role that the eigenstructure of the structural connectivity matrix has in determining the stability of equilibria.We report the locus of Hopf and saddle-node sets in Fig. 4. Comparison of Figs 2 and 4 shows that the bifurcation structure of steady states for the full network is practically identical to that of the single node, highlighting the potential importance of single-node dynamics in driving SC-FC correlations.

Network Neuroscience
False bifurcations In Fig. 3 we consider in more detail the types of activity that the network model (1) supports.In particular, we observe that, under changes to parameter values within the oscillatory region (see highlighted parameter values in Fig. 2), the time-course of activity shifts from single-to double-peaked waves, which could have consequences for synchronisation of oscillations and, moreover, FC.The points of transition are known as false bifurcations since there is a significant dynamical change that occurs smoothly rather than critically.False bifurcations in a neural context have previously been seen as canards in single neuron models (Desroches, Krupa, & Rodrigues, 2013) as well as in EEG models of absence seizures (Marten, Rodrigues, Benjamin, Richardson, & Terry, 2009).In the latter case the false bifurcation corresponds to the formation of spikes associated with epileptic seizures (Moeller et al., 2008).
As illustrated in Fig. 3 the false-bifurcation transition is characterised by the change from a double-peaked profile (a) to a sinusoidal-like waveform (c) via the development of a point of inflection in the solution trajectory (b).Since this transition is not associated with a change in stability of the periodic orbit, these false bifurcations are determined by tracking parameter sets for which points of inflection occur.We refer the reader to Rodrigues et al. (2010) for details on methods for detecting and continuing false bifurcations in dynamical systems.The result of this computation is shown in Fig. 2. In the full network (not shown), this computation is more laborious, though gives very similar results to that for a single node.Subfigures in the upper row are plots of the timeseries solution, whereas the bottom row shows the trajectories of stable orbits in the (y, y ) plane.The chosen parameters lie at either side of the region where a smooth transition between activity types occurs, corresponding to a false bifurcation (see highlighted parameter values in Fig. 2).In (b), an inflection point occurs and is highlighted as a red star on the orbit.

Weakly-coupled oscillator theory
Further insight into the phase relationship between nodes in a network can be obtained from the theory of weakly coupled oscillators (see, e.g., Hoppensteadt and Izhikevich (2012)).This technique reduces a network of limit cycle oscillators to a set of relative phases in a systematic way.The resulting set of network ODEs is (N − 1)-dimensional, as opposed to the (Nm)-dimensionality of the original system, and provides an accurate model as long as the overall coupling strength is weak (|ε| 1).This is because when all oscillators lie on the same limit cycle of a system, the interactions from pairwise-connected nodes can be considered as small perturbations to the oscillator dynamics.Moreover, the resulting set of network ODEs only depends upon phase differences and it is straightforward to construct relative equilibria (oscillatory network states) and determine their stability in terms of both local dynamics and structural connectivity.A method to construct the phase interaction function, H, for the network is provided in Mathematical Methods.Once this is known, the dynamics for the phases of each node in the network, θ i ∈ [0, 2π), takes the simple form: where Ω = 2π/T represents the natural frequency of an uncoupled oscillatory node with period T, and the second term determines phase changes arising from pairwise interactions between nodes.We emphasise that the T-periodic phase interaction function H(Ωt) = H(Ω(t + T)) is derived from the full system given by (1).For a given phase-locked state θ i (t) = Ωt + φ i (where φ i is the constant phase of each node), stability is determined in terms of the eigenvalues of the Jacobian of (4), denoted by H(Φ) with Φ = (φ 1 , . . ., φ N ) , with components: The globally synchronous steady-state, φ i = φ for all i, exists in a network with a phase interaction function that vanishes at the origin (i.e.H(0) = 0, which is not the case here), or for one with a row-sum constraint, ∑ j w ij = Γ = constant for all i, which is true for our specific structural matrix (for which Γ = 1).Note that the emergent frequency of the synchronous network state is given explicitly by Ω + εΓH(0).Using the Jacobian in (5), synchrony is found to be stable if εH (0) > 0 and all the eigenvalues of the graph Laplacian of the structural network, lie in the right hand complex plane.Since the eigenvalues of a graph Laplacian all have the same sign (apart from, in this case, a single zero value) then stability is entirely determined by the sign of εH (0).For example, for a globally coupled network with w ij = 1/N then the graph Laplacian has one zero eigenvalue, and (N − 1) other degenerate eigenvalues at −1, and so synchrony is stable if εH (0) > 0.
It is therefore useful to consider the condition εH (0) > 0 as a natural prerequisite for a structured network to support high levels of synchrony (without recourse to exploring the full Jacobian structure).A plot of εH (0) is shown in Fig. 4b.For completeness, however, the full Jacobian was also computed in order to account for the potential influence of detailed structure on the correspondence with the observed SC-FC agreement measured in simulations.To do this, the system given by (1) was integrated with ε = 0.001 to a (stable) Network Neuroscience phase-locked state, and relative phases computed.The eigenvalues of the Jacobian (eq.( 5)) were then computed, providing an indication of solution attractivity.The largest non-zero eigenvalue for each parameter choice is shown in Fig. 4(c).
It has been shown in Tewarie et al. (2018) that the eigenstructure of the structural connectivity matrix is predictive of the emergent functional networks arising after an instability from a steady state.The largest non-zero eigenvalue, which is related the most unstable eigenmode (or closest to instability), was found to be a good predictor of resultant FC by computing the tensor product of its corresponding eigenvector, v ⊗ v.Here we take this further by considering instabilities of the synchronous state.In this case the Jacobian (5) reduces to −εH (0)L ij and the phase-locked state that emerges beyond an instability of the synchronous state has a pattern determined by the a linear combination of eigenmodes of the graph Laplacian, since all eigenmodes destabilise simultaneously.It is known that the graph Laplacian can be used to predict phase-locked patterns (J.Chen, Lu, Zhan, & Chen, 2012) and has indeed been used to predict empirical FC from SC (Abdelnour, Dayan, Devinsky, Thesen, & Raj, 2018).Following from this, the eigenmodes of the Jacobian in ( 5) can be used as simple, easily computable proxy for the FC matrix.In Fig. 6 we compare the FC pattern from the (fully nonlinear) weakly coupled network with the linear prediction, to highlight its usefulness.FC in the weakly-coupled network is computed via MPA (equation ( 3)).For comparison, we use the tensor product sum, , which donates the k th eigenvector of the Jacobian for the synchronous state.These are weighted by their corresponding eigenvalues, λ k .

RESULTS
Fig. 4 shows plots in the (A, B) parameter space highlighting our studies on the combined influence of SC and node dynamics on FC.The region bounded by the bifurcation curves, obtained via a linear instability analysis of the network steady state, is where the network model supports oscillations as well as phase-locked states.In Fig. 4(a) the Jaccard similarity between structure and function is computed from direct numerical simulations of the Jansen-Rit network model (1).Beyond the onset of an oscillatory instability (supercritical Hopf bifurcation) the emergent phase-locked network states show a nontrivial correlation with the SC.This varies in a rich way as one traverses the (A, B) parameter space, showing that precise form of the node dynamics can have a substantial influence on the network state.The highest correlation between SC and FC coincides with a Hopf bifurcation of a network equilibrium (shown as a solid white line), whilst a band of much lower correlation coincides with the false bifurcation of a single node (in black), reproduced from Fig. 2. Indeed, it would appear that these mathematical constructs are natural for organising the behaviour seen in our in silico experiments.
In Fig. 4(b) we show a plot of H (0). Recall from Weakly-coupled oscillator theory that a globally synchronous state (which is guaranteed to exist from the row-sum constraint) is stable if εH (0) > 0. Comparison with Fig. 4(a), highlights that when synchrony is unstable (εH (0) < 0) SC only weakly drives FC.Moreover, this instability region coincides with the false bifurcation, stressing the important role of false bifurcations for understanding SC-FC correlation.Of course, there is a much finer structure in Fig. 4(a) that is not predicted by considering either the bifurcation from steady state, or the weakly-coupled analysis of synchronous states, and so it is illuminating to pursue the full weakly coupled oscillator analysis for structured networks.The eigenvalues of the Jacobian, corresponding to more general stable phase-locked states, can be used to give a measure of solution attractivity.The largest eigenvalue is plotted in Fig. 4(c).The most stable phase-locked states occur in the neighbourhood of false bifurcations, as well as along the existence border for oscillations, defined by a saddle node bifurcation.Furthermore, apart from near false bifurcations, stronger stability of the general phase-locked states corresponds with stronger stability of global synchrony (Fig. 4(b)).

Network Neuroscience
To test the predictive power of the weakly-coupled theory, in Fig. 5 we compare the the emergent FC structure obtained from direct simulations of the Jansen-Rit network model (1) against direct simulations of the weakly-coupled oscillator network (4).For the former, the phases required to compute the mean phase agreement (equation ( 3)) are determined from each timeseries by a Hilbert transform; in the latter case, the phase variables from equation ( 4) are employed directly in (7).As expected, we find excellent agreement between the modular FC structure in each case for very weak coupling, with this agreement reducing as ε is increased.Moreover, through the instability theory of the synchronous state we can construct a proxy for the FC as described in Weakly-coupled oscillator theory.In Fig. 6, we compare simulated FC with the eigenmode predictor (eq.( 7), where v k is the eigenvector of associated with the k th largest non-zero eigenvalue of −εH (0)L), choosing parameter Network Neuroscience values that lie just beyond the onset of instability of the globally synchronous state (in this case, near the false bifurcation set shown in Fig. 4a,b).This is a much more efficient way of simulating an emergent FC pattern, since it does not require brute-force forward integrations of the model, which may take a long time to converge.All of these results highlight the strong impact that nodal dynamics can have on the correlation between structure and function, and the usefulness of bifurcation theory and network reduction techniques to provide insight into how SC-FC correlations are organised across parameter space.

DISCUSSION
In this paper, we investigate the degree to which the dynamical state of neural populations, as well as their structural connectivity, facilitates the emergence of functional connections Network Neuroscience in a neural-mass network model of the human brain.We have addressed this by using a mixture of computational and mathematical techniques to assess the correlation between structural and functional connectivity as one traverses the parameter space controlling the inhibitory and excitatory dynamics and bifurcations of an isolated Jansen-Rit neural mass model.Importantly, SC has been estimated from HCP diffusion MRI datasets.We find that structure strongly drives function when the system is close to a Hopf bifurcation, whereas in the neighbourhood of a false bifurcation, this drive is diminished.These results emphasise the vital role that local dynamics has to play in determining FC in a network with a static SC.In addition, we show that a weakly-coupled analysis provides insight into the organisation of SC-FC correlation features across parameter space, and can be exploited to predict emergent FC structure.
The modelling work presented here is relevant in a wider neuroimaging context-for example, epilepsy is often considered to be caused by irregularities in synchronisation (Lehnertz et al., 2009;Mormann et al., 2003;Netoff & Schiff, 2002).It is noteworthy that the changes in synchrony patterns that we observe arise from local dynamical considerations as opposed to large scale structural ones.In the Jansen-Rit model, the bifurcations organising emergent FC take the form of Hopf, saddle and false bifurcations.False bifurcations have received relatively little attention in the dynamical systems community (a notable exception being the work of Marten et al. (2009)), although our results indicate that they may Network Neuroscience be significant for understanding how 'synchronisability' of brain networks is reduced during seizures.This phenomena was reported in Schindler, Bialonski, Horstmann, Elger, and Lehnertz (2008), which also found that synchronisability increases as the patient recovers from seizure state.
In summary, the findings reported here suggest that there are multiple factors which give rise to emergent FC.While structure clearly facilitates functional connectivity, the degree to which it influences emergent functional states is determined by the dynamics of its neural sub-units.Importantly, we have shown that local dynamics has a clear influence on SC-FC correlation, as does network topology and coupling strength.Our combined mathematical and computational study has demonstrated that a full description of the mechanisms that dictate the formation of FC from anatomy requires knowledge of how both neuronal activity and connectivity are modulated and, moreover, exposes the utility of bifurcation theory and network reduction techniques.This work can be extended to more complex neural mass models such as that derived in Coombes and Byrne (2019), to further explore the relationship between dynamics and structure-function relations in systems with more sophisticated models for node dynamics.

Network bifurcations of equilibria
Each node of the network is described by the m-dimensional Jansen-Rit model, with m = 6.There are N = 78 nodes in the network.Analysing bifurcations of network equilibria requires finding a set of m × N eigenvalues from the linearised system.Defining y i ∈ R m as (y 1 i , . . . ,y m i ) allows us to write the system of first-order ODEs given by (1) for the network as: where and B = (0, 0, 0, 0, AaP, 0) , K(y) = (0, 0, 0, 0, f (y 1 − y 2 ), 0) , with Here we have introduced the 3 × 3 identity matrix I 3 , and the 3 × 3 zero matrix 0 3 .The network steady state y i = y i , for i = 1, . . ., N, is defined by setting the left hand side of (8) to zero.We now linearise (8) by setting y i (t) = y i + u i (t), where u i (t) is a small perturbation.This gives, Network Neuroscience where DK(y) ∈ R m×m is the Jacobian of K(y).The only two non-zero entries of this matrix are given by [K(y)] 5,2 = f (y 1 − y 2 ) = −[K(y)] 5,3 .It is now useful to define DF i ≡ [M + L f (τy i ) τ] and DG j ≡ εaADK(y j ), so that DF i is the Jacobian which describes the the intra-mass dynamics of node i and DG j is the Jabobian for the effect of the inter-mass interactions with node j.Then we may write ( 12 where U = (u 1 , . . ., u N ) , and ⊗ denotes the tensor product.This system can be simplified by considering the eigenvalues of the connectivity matrix w ∈ R N×N (with components w ij ).We introduce a matrix of normalised eigenvectors, E, and a corresponding diagonal matrix of eigenvalues, Λ = diag(µ 1 . . .(15) The system (15) is in a block diagonal form and so it is equivalent to the set of decoupled equations given by Solving E = 0 for λ produces a set of eigenvalues which can be tracked to determine bifurcations.Since stability requires the real part of all eigenvalues to be negative, if one of these eigenvalues crosses the imaginary axis the solution can undergo either a saddle-node bifurcation (Re λ = 0 = Im (λ)) or a Hopf bifurcation (Re λ = 0, Im (λ) = 0).

Network Neuroscience
To investigate the nature of phase-locked oscillatory states in the Jansen-Rit network, it is appropriate to use weakly-coupled oscillator theory.For a recent review see (Ashwin, Coombes, & Nicks, 2016).This gives rise to the set of equations ( 4) where the phase interaction function H is determined in terms of two quantities.The first is the so-called phase response or adjoint Q ∈ R m , that describes the response of an attracting limit cycle to a small perturbation.This can be computed by solving the adjoint equation.Here y(t) is a T-periodic of the Jansen-Rit node model and , denotes a Euclidean inner product between vectors.The second ingredient comes from writing the physical interactions in terms of phases rather than the original state variables.This is easily done by writing y i (t) = y(θ i /Ω).The phase interaction function is then obtained as The adjoint equation is readily solved numerically by backward integration in time (Williams & Bowtell, 1997), whilst the integral in ( 19) can be evaluated using numerical quadrature.

Figure 1 .
Figure 1.The original structural matrix (a) is derived from DTI data taken from the Human Connectome Project database and parcellated on to a 78-region brain atlas.This is thresholded and binarised to keep the top third strongest connections (b) and normalised by row so that ∑ N j=1 w ij = 1 for all regions i) in (c).

Figure 2 .
Figure 2. Two-parameter bifurcation diagram in the (A, B) plane in the single-node case of the Jansen-Rit system of equations (1).Red lines are Hopf bifurcations, black are false bifurcations and blue lines represent saddle points.Parameters are taken from Table 1.The three coloured points at B = 22, A = 7.0, 7.7, 9.0 indicate parameter values at which we observe distinctly different dynamics as shown in Fig. 3.

Figure 3 .
Figure3.Activity profiles of y = y 1 − y 2 , the potential of the main population of pyramidal neurons for a node in the Jansen-Rit network, with B fixed at 22 and (a) A = 9.0;(b) A = 7.7; (c) A = 7.0.Subfigures in the upper row are plots of the timeseries solution, whereas the bottom row shows the trajectories of stable orbits in the (y, y ) plane.The chosen parameters lie at either side of the region where a smooth transition between activity types occurs, corresponding to a false bifurcation (see highlighted parameter values in Fig.2).In (b), an inflection point occurs and is highlighted as a red star on the orbit.

Figure 4 .
Figure 4. (a) Jaccard similarity coefficient between SC and FC when the network supports an oscillatory solution, averaged over 30 realisations of initial conditions chosen at random.Warmer colours indicate greater SC/FC correlation.Here we have superimposed the bifurcation diagram for the network steady state, which shows the oscillatory region being bounded by Hopf/saddle-node sets in solid/dashed white lines respectively.Boxes are Bogdanov-Takens points and false bifurcations are indicated by a black line; (b) The value of H (0) in the A, B-plane.When this value is positive/negative, the globally synchronised solution is stable/unstable (if it exists); (c) The largest non-zero eigenvalue of the Jacobian for the full weaklycoupled oscillator network (equation (5)), calculated at a stable phase-locked state.More negative values indicate a stronger stability.

Figure 5 .
Figure 5.Comparison of FC patterns from averages of realisations of the weakly-coupled oscillator model with corresponding Jansen-Rit simulations at A = 5, B = 19, computing averages over 600 realisations with initial conditions chosen at random.We do this for (a) ε=0.01;(b) ε=0.1;(c) ε=1 to show how the weakly-coupled theory becomes less predictive for stronger coupling strengths.

Figure 6 .
Figure 6.Comparison of (a) a FC prediction given by the a linear combination of eigenmodes of the weakly-coupled oscillator system, given by tensor products of eigenvectors of the SC network's graph Laplacian (b) a direct simulation of the Jansen-Rit network model and using (3), calculated for A = 6, B = 18, which lies near the existence border for stable synchronous solutions (see Fig. 4(b)).
µ N ), such that wE = EΛ.Imposing the change of variables V = (E ⊗ I m ) −1 U transforms (13) to d dt V = (E ⊗ I m ) to establish that for any block diagonal matrix A, formed from N matrices of size m × m, that (E ⊗ I m ) −1 A(E ⊗ I m ) = A.Moreover, using standard properties of the tensor operator, (E ⊗ I m ) −1 (w ⊗ I m ) = (E −1 w) ⊗ I m = (ΛE −1 ) ⊗ I m = (Λ ⊗ I m )(E −1 ⊗ I m ).
It is convenient to write the dynamics for a single uncoupled Jansen-Rit node in the form ẏ = F(y), with F, y ∈ R m .Using the notation above we have explicitly that F(y) = My + B + L f (τy).The adjoint is given by the T-periodic solution

Table 1 .
Parameters in the Jansen-Rit model, given by equations (