## Abstract

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.

## 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 structural connectivity only weakly influences functional connectivity. This highlights the important effect that local dynamics can have on large-scale brain states.

## 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 increasingly accepted that the understanding of brain function may be obtained from a network perspective, rather than by exclusive study of its individual subunits. 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 graph-theoretical properties of such large-scale networks have been well studied, highlighting key features including small-world architecture (Bassett & Bullmore, 2006; Liao, Vasilakos, & He, 2017), hub regions and cores (Oldham & Fornito, 2019; van den Heuvel & Sporns, 2013), rich club organisation (Betzel, Gu, Medaglia, Pasqualetti, & Bassett, 2016; Van Den Heuvel & Sporns, 2011), a hierarchical-like modular structure (Meunier, Lambiotte, & Bullmore, 2010; Sporns & Betzel, 2016), and economical wiring (Betzel et al., 2017; Bullmore & Sporns, 2012). The emergent brain activity that this structure supports can be evaluated by functional connectivity (FC) network analyses that describe patterns of temporal coherence in neural activity between brain regions. These highly dynamic patterns are widely believed to be significant in integrative processes underlying higher brain function (Van Den Heuvel & Pol, 2010; van Straaten & Stam, 2013), and disruptions in SC and FC networks are associated with many psychiatric and neurological diseases (Braun, Muldoon, & Bassett, 2015; Menon, 2011).

However, 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). Empirical studies suggest that while a 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 (Fox et al., 2005; Hutchison et al., 2013; Liegeois, Laumann, Snyder, Zhou, & Yeo, 2017; Preti, Bolton, & Van De Ville, 2017), motivating more recent consideration of *dynamic* (rather than time-averaged) FC networks, which have been proposed to more accurately represent brain function. An important example of SC-FC divergence is provided by resting-state networks, such as the ‘default mode network’ and the ‘core network’ (Thomas Yeo et al., 2011; Van Den Heuvel & Pol, 2010). These networks comprise brain areas that can be strongly functionally connected at rest (Van Den Heuvel & Pol, 2010), but can also temporally vary. Indeed, a neural ’switch’ has been proposed that facilitates transitions between resting-state networks (Goulden et al., 2014), and a theoretical study by Messé, Rudrauf, Benali, and Marrelec (2014) estimated that nonstationarity of FC contributes to over half of observed FC variance.

Theoretical studies deploying anatomically realistic structural networks obtained through tractography alongside neural mass models describing mean-field regional neural activity have 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, showing that FC reflects underlying SC on slow time scales, but significantly less so on faster time scales (C. Honey et al., 2009; C. J. Honey et al., 2007; Rubinov et al., 2009).

In the context of mean-field models, simulated (typically time-averaged) FC has been found most strongly to resemble SC 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). Similarly, 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 et al., 2019; Tewarie et al., 2018). These results provide a possible manifestation of the so-called critical brain dynamics hypothesis (Cocchi, Gollo, Zalesky, & Breakspear, 2017; Shew & Plenz, 2013). In Crofts, Forrester, and O’Dea (2016), both SC and FC 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. Deco, Kringelbach, Jirsa, and Ritter (2017) consider dynamic FC, with transient FC states described as metastable states, and in Deco et al. (2019), metastability of a computational model of large-scale brain network activity was used to predict which structures of the brain could be influenced to force a transition between states of wakefulness and sleep. Hansen, Battaglia, Spiegler, Deco, & Jirsa (2015) were also able to observe dynamic transitions between states resembling resting-state networks in a noise-driven, nonlinear mean-field model of neural activity.

In this paper, we adopt the mean-field neural mass approach and present a combined computational and mathematical study, which significantly extends the related works of Hlinka and Coombes (2012) and Crofts et al. (2016) to investigate how the detailed and rich dynamics of the intrinsic behaviour of neural populations, together with structural connectivity, combine to shape FC networks. Thereby, we provide a complementary investigation to many of the aforementioned studies which focus on the analysis of brain networks themselves, or those that employ statistical models, 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 FC and SC topology under systematic variation of model parameters associated with excitatory/inhibitory neural responses, and analysed by making use of techniques from bifurcation and weakly coupled oscillator theory.

## METHODS

### Neural mass model

*m*= 6, that describes the evolution of the average postsynaptic potential (PSP) in three interacting neural populations: pyramidal cells (

*y*

_{0}), and excitatory (

*y*

_{1}) and inhibitory (

*y*

_{2}) interneurons. These populations are connected with strengths

*C*

_{i}(

*i*= 1…4), representing the average number of synaptic connections between each population. The Jansen–Rit model is mathematically 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 therefore associated with the dynamics of the population average of PSPs and their temporal derivatives. The quantity of primary interest herein is

*y*=

*y*

_{1}−

*y*

_{2}, which is physiologically interpreted as the average potential of pyramidal populations and the main contributor to signals generated in EEG recordings (Teplan, 2002). 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:

*f*is a sigmoidal nonlinearity, representing the transduction of activity into a firing rate, and with the specific form

The model is identical to that presented in Jansen and Rit (1995) for a single cortical column, but is completed by specifying the network interactions as a function of average membrane potential of afferently connected pyramidal populations, encoded in a connectivity matrix with elements *w*_{ij} (described in ec2*Structural and Functional Connectivity*), with an overall scale of interaction set by *ε*. The remaining model parameters, together with their physiological interpretations and values (taken from Grimbert & Faugeras, 2006, and Touboul, Wendling, Chauvel, & Faugeras, 2011), are given in Table 1. A schematic ‘wiring diagram’ for the model indicating the interactions between different neural populations is shown in Figure 1.

Parameter . | Meaning . | Value . |
---|---|---|

C_{1}, C_{2}, C_{3}, C_{4} | Average number of synapses between populations | 135, 108, 33.75, 33.75 |

P_{i} | Basal extracortical input to main pyramidal excitatory populations | 120 Hz |

A, B | Amplitude of excitatory, inhibitory PSPs respectively | [2, 14] mV, [10, 30] mV |

a, b | Lumped time constants of excitatory, inhibitory PSPs | 100 s^{−1}, 50 s^{−1} |

ε | Global coupling strength | 0.1 |

w_{ij} | Coupling from node j to i | [0, 1] |

ν_{max} | Maximum population firing rate | 5 Hz |

v_{0} | Potential at which half-maximum firing rate is achieved | 6 mV |

r | Gradient of sigmoid at v_{0} | 0.56 mV^{−1} |

Parameter . | Meaning . | Value . |
---|---|---|

C_{1}, C_{2}, C_{3}, C_{4} | Average number of synapses between populations | 135, 108, 33.75, 33.75 |

P_{i} | Basal extracortical input to main pyramidal excitatory populations | 120 Hz |

A, B | Amplitude of excitatory, inhibitory PSPs respectively | [2, 14] mV, [10, 30] mV |

a, b | Lumped time constants of excitatory, inhibitory PSPs | 100 s^{−1}, 50 s^{−1} |

ε | Global coupling strength | 0.1 |

w_{ij} | Coupling from node j to i | [0, 1] |

ν_{max} | Maximum population firing rate | 5 Hz |

v_{0} | Potential at which half-maximum firing rate is achieved | 6 mV |

r | Gradient of sigmoid at v_{0} | 0.56 mV^{−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). It is known that a single Jansen–Rit node can support multistable behaviour, which includes oscillations of different amplitude and frequency, but, moreover, a network of these nodes can also exhibit various stable phase-locked states. A small amount of white noise is added to the extracortical input *P*_{i} on each node, in order to allow the system to explore a variety of these dynamical states: *P*_{i} + *dW*_{i}(*t*), where *dW*_{i}(*t*) is chosen at random from a Gaussian distribution with standard deviation 10^{−1} Hz and mean 0 Hz. For direct simulations of the network we use an Euler–Murayama scheme, implemented in Matlab, with a fixed numerical time step of 10^{−4}, which we have confirmed ensures adequate convergence of the method.

### Structural and functional connectivity

The structural connectivity was estimated using diffusion MRI data recorded with informed consent from 10 subjects, obtained from the HCP (Van Essen et al., 2013). Briefly, we explain how this data is postprocessed to derive connectomic data, though we direct the reader to Tewarie et al. (2019) and the references therein for a more detailed overview. Sixty thousand 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 propagated 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 enable a meaningful comparison between the network measures of SC and FC, the former reflecting the density of tractography streamlines and the latter that of correlated neural activity, we place them on a similar footing by thesholding and binarising, such that only the top 23% of the weights (ordered by strength) are retained; see Figure 2. Thresholding is a widespread technique for removing spurious connections that may not in fact be a realistic representation of brain connectivity. We note that our thresholding choice (that reduces the number of connections, while ensuring that the overall modular structure is unchanged) is commensurate with a recent study (Tsai, 2018), which employed DTI data averaged on the same brain atlas as used herein to consider thresholding approaches suitable to remove weak connections with high variability between (*n* = 30) different subjects. 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. This normalisation process permits some of the analysis that we undertake to help explain SC-FC relations (see ec3*Weakly Coupled Oscillator Theory*); however, we highlight that the results that we present herein are not crucially dependent on such a choice and so our conclusions generalise (see Supporting Information, section Mathematical Methods).

*y*

_{j}=

*y*

_{1j}−

*y*

_{2j}as the variable of interest because of its relation to the EEG signal, making it a good candidate to produce timeseries more readily comparable with empirical data. Pairwise MPC measures the average temporal variance 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 MPC of the time-series comprising

*M*time-points

*t*

_{l}(

*l*= 1, …,

*M*) is defined as:

Structure–function relations are assessed by computing the Jaccard similarity coefficient (Jaccard, 1912) of the nondiagonal 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.

Since the SC-FC correlation patterns of interest here arise naturally from global synchrony or patterns of phase locking of oscillatory node activity, the local stability of oscillatory node dynamics and of network (global or phase locking) synchrony is a natural candidate to explain the structures we observe. In the following subsections we consider bifurcation, false bifurcation, and weakly coupled oscillator theory approaches to address this.

### Bifurcation analysis

#### 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. We also observe a region of bistability bounded by fold bifurcations of limit cycles, in which the types of activity described in Figure 4A and 4C can both exist. This is shown in Figure 3. We refer the reader to Grimbert and Faugeras (2006), Touboul et al. (2011), and Spiegler, Kiebel, Atay, and Knösche (2010) for a comprehensive analysis of the bifurcation structure of the Jansen–Rit model.

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 the Supporting Information, section 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 eigenmodes of the structural connectivity matrix has in determining the stability of equilibria. We report the locus of Hopf and saddle-node sets for the network in Figure 5. Comparison of Figures 3 and 5 shows that the bifurcation structure of steady states for the full network is practically identical to that of the single node (even for moderate coupling strength—here, *ε* = 0.1), highlighting the potential importance of single-node dynamics in driving SC-FC correlations.

#### False bifurcations

In Figure 4 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 Figure 3), 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 Figure 4 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 Figure 3, where we observe the set of false bifurcations arising from the breakdown of two branches of fold bifurcations of limit cycles. In the full network (not shown), this computation is more laborious (and there is some delicacy in defining the bifurcation since the network coupling leads nodes to inflect at marginally different parameter values); however, we obtain very similar results to those obtained in Figure 3 for a single node (not shown).

### Weakly coupled oscillator theory

*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 on 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 the Supporting Information, section Mathematical Methods. Once this is known, the dynamics for the phases of each node in the network,

*θ*

_{i}∈ [0, 2

*π*), takes the simple form:

*Ω*= 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 Equation 1. For a given phase-locked state

*θ*

_{i}(

*t*) = Ω

*t*+

*ϕ*

_{i}(where

*ϕ*

_{i}is the constant phase of each node), local stability is determined in terms of the eigenvalues of the Jacobian of Equation 4, denoted by $H^$(

**Φ**) with

**Φ**= (

*ϕ*

_{1}, …,

*ϕ*

_{N})

^{⊺}, with components:

*ϕ*

_{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 Equation 5, synchrony is found to be stable if

*ε*

*H*′(0) > 0 and all the eigenvalues of the graph Laplacian of the structural network,

*ε*

*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 Figure 5B. 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 Equation 1 was integrated with *ε* = 0.001 to a (stable) phase-locked state, and relative phases computed. The eigenvalues of the Jacobian (Equation 5) were then computed, providing an indication of solution attractivity. The largest nonzero eigenvalue for each parameter choice is shown in Figure 5C.

*v*⊗

*v*. Here we take this further by considering instabilities of the

*synchronous*state. In this case the Jacobian Equation 5 reduces to −

*ε*

*H*′(0)𝓛

_{ij}and the phase-locked state that emerges beyond 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 (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 Equation 5 can be used as simple, easily computable proxy for the FC matrix when the system is poised at a local instability. In Figure 7 we compare the FC pattern from the (fully nonlinear) weakly coupled network with a linear prediction to highlight its usefulness. blackIn this case, MPC (Equation 3) is not ideally suited for our study because it struggles to discern between phase locking and complete synchrony, yet we consider situations where stable phase locking naturally arises. Therefore, FC in the weakly coupled network is computed via the new metric of mean phase agreement (MPA), whereby patterns of coherence are determined by a temporal average of relative phase differences:

*v*

_{k}= ($vk1$, …, $vkN$), which denotes the

*k*

^{th}eigenvector of the Jacobian for the synchronous state. These are weighted by their corresponding eigenvalues,

*λ*

_{k}, and we include the

*N** unstable eigenmodes.

## RESULTS

Figure 5 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 Figure 5A the Jaccard similarity between SC and FC is computed from direct numerical simulations of the Jansen–Rit network model (Equation 1). Beyond the onset of 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 fold bifurcations of limit cycles and false bifurcations of a single node (in black), reproduced from Figure 3. Indeed, it would appear that these mathematical constructs are natural for organising the behaviour seen in our in silico experiments. We reiterate that we have confirmed that the organising SC-FC features that we here identify are not crucially dependent on the binarisation, thresholding and normalisation procedure, described in ec2*Structural and functional connectivity* and are qualitatively similar under variation of coupling strength (see Supporting Information, section Mathematical Methods); moreover, results obtained via MPC and of MPA are indistinguishable (data not shown). In Figure 5B we show a plot of *H*′(0). Recall from ec3*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 Figure 5A, highlights that when synchrony is unstable (*ε**H*′(0) < 0) SC only weakly drives FC. Moreover, this instability region coincides with the region of bistability and the false bifurcation, stressing the important role of these bifurcations for understanding SC-FC correlation.

Of course, there is a much finer structure in Figure 5A 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 Figure 5C. The most stable (nonsynchronous) phase-locked states occur in the neighbourhood of the false bifurcations, as well as in the region of bistability and 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 (Figure 5B).

To test the predictive power of the weakly coupled theory, in Figure 6 we compare the emergent FC structure obtained from direct simulations of the Jansen–Rit network model (Equation 1) against direct simulations of the weakly-coupled oscillator network (Equation 4). For the former, the phases required to compute the mean phase agreement (Equation 7) are determined from each time series by a Hilbert transform; in the latter case, the phase variables from Equation 4 are employed directly. Since the weakly coupled reduction of the Jansen–Rit model is deterministic, these computations were ran in the absence of noise (*dW*_{i} = 0 for all nodes). As expected, we find excellent agreement between the modular FC structure in the case for very weak coupling, with this agreement reducing with increasing *ε*, as quantified by a reduction in Jaccard similarity (from 0.98 in panel A to 0.65 in C). This is a manifestation of the network moving from a dynamical regime that can be well described by the weakly coupled reduction (Equation 4) to one where stronger network interactions dominate. Since an analogous theory does not exist for stronger coupling, we do not consider here how SC-FC relations arise from network dynamics within a strongly coupled framework. Moreover, through the instability theory of the synchronous state we can construct a proxy for the FC as described in ec3*Weakly coupled oscillator theory*. In Figure 7 we compare simulated FC with that predicted by $R^$ (Equation 8; i.e., using the unstable eigenmodes of the Jacobian at synchrony), for parameter values that lie just beyond the onset of instability of the globally synchronous state and near the false bifurcation set (see Figure 5A and 5B). We observe that the key features of the FC are captured by the eigenmode prediction; indeed the (weighted) Jaccard similarity coefficient between predicted and simulated FC (both scaled to [0, 1]) is calculated to be 0.82. 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 SC and FC, and the utility of bifurcation theory and phase oscillator reduction techniques (that are naturally positioned to explain the generation of patterns of synchronous node and network activity) 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 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 SC strongly drives FC 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. Messé et al. (2014) considered statistical models to predict FC from SC (in particular, a spatial simultaneous autoregressive model (sSAR), whose parameters can be estimated in a Bayesian framework) and found, interestingly, that simpler linear models were able to fare at least as well. More recently, Saggio, Ritter, and Jirsa (2016) were also able to make predictions of FC from empirical SC data (and vice versa) using a simple linear model. Since the only free parameter of their model for SC is the global coupling strength, results from this method are efficient and computationally inexpensive. We have not attempted to reproduce empirical data here, but we have shown that similar predictions can be made using bifurcation theory and network reduction techniques; such an approach allows us to consider in more detail, and explain, the influence of the rich neural dynamics supported by the Jansen–Rit model on SC-FC relationships. Nevertheless, it is important to note that the FC structures we are concerned with are averaged over long time scales and therefore represent a static FC state, as opposed to dynamic FC (as discussed in Introduction). Use of such static FC networks as a clinical biomarker is widespread; however, subject variability in FC means that their predictive power is restricted to group analyses (Mueller et al., 2013). To capture the rich dynamic FC repertoire exhibited in empirical resting-state data, for example, the distinct hierarchical organisation in switching between FC states (Vidaurre, Smith, & Woolrich, 2017), will require alternative approaches. One such approach is dynamic causal modelling, as employed in Goulden et al. (2014) and Van de Steen, Almgren, Razi, Friston, and Marinazzo (2019) for empirical data.

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, fold of limit cycle, 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 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.

A natural extension to the work presented here would be the inclusion of conduction delays, characterised by Euclidean or path-length distances between brain regions, which are certainly important in modulating the spatiotemperal coherence in the brain (Deco, Jirsa, McIntosh, Sporns, & Kötter, 2009). These would manifest as constant phase shifts in the weakly coupled reduction of the model (Ton, Deco, & Daffertshofer, 2014). For strongly coupled systems the mathematical treatment of networks with delayed interactions remains an open challenge. Recent work in this vein by Tewarie et al. (2019) focusses on the role of delays in destabilising network steady states, and techniques extending the Master Stability Function to delayed systems (Otto, Radons, Bachrathy, & Orosz, 2018) may be appropriate for treating phase-locked network states.

In summary, the findings reported here suggest that there are multiple factors which give rise to emergent FC. While structure clearly facilitates FC, the degree to which it influences emergent FC states is determined by the dynamics of its neural subunits. 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.

## AUTHOR CONTRIBUTIONS

Michael Forrester: Investigation; Writing - Original Draft. Stephen Coombes: Supervision; Writing - Review & Editing. Jonathan Crofts: Supervision; Writing - Review & Editing. Stamatios Sotiropoulos: Data curation; Writing - Review & Editing. Reuben O’Dea: Supervision; Writing - Review & Editing.

## FUNDING INFORMATION

Michael Forrester, Engineering and Physical Sciences Research Council (http://dx.doi.org/10.13039/501100000266), Award ID: EP/N50970X/1.

## TECHNICAL TERMS

- Diffusion tensor imaging:
A magnetic resonance imaging technique that measures the diffusion of water in tissue in order to produce axonal fibre tract images.

- Structural connectivity:
A pattern of anatomical links between between distinct brain regions. This is often measured using diffusion tensor imaging.

- Functional connectivity:
A pattern of statistical dependencies between between distinct brain regions. This is often measured using coherence or correlation measures between time-series.

- Neural mass model:
A phenomenological model for the activity of a neuronal population cast as a system of ordinary differential equations.

- Hopf bifurcation:
The appearance or disappearance of a periodic orbit through a local change in the stability properties of an equilibrium point under parameter variation.

- Human Connectome Project:
A consortium research project to build a human connectome for structural and functional connectivity, and facilitate research into brain disorders.

- Weakly coupled oscillator theory:
A reduction of a network of weakly interacting limit cycle oscillators to a system of fewer dynamical equations that describes the evolution of relative phases between nodes.

- Jaccard similarity:
An index for comparing members for two sets of binary data to see which are shared and which are distinct. The higher the index, the more similar the two sets.

- False bifurcation:
A qualitative change in the shape of a periodic orbit without a change in stability, say from a small to a large amplitude oscillation, that occurs under a very small parameter variation.

- Saddle-node bifurcation:
A local bifurcation in which two equilibria of a dynamical system collide and annihilate.

## REFERENCES

## Competing Interests

Competing Interests: The authors have declared that no competing interests exist.

## Author notes

Handling Editor: Petra Ritter