## Abstract

The complex network dynamics that arise from the interaction of the brain's structural and functional architectures give rise to mental function. Theoretical models demonstrate that the structure–function relation is maximal when the global network dynamics operate at a critical point of state transition. In the present work, we used a dynamic mean-field neural model to fit empirical structural connectivity (SC) and functional connectivity (FC) data acquired in humans and macaques and developed a new iterative-fitting algorithm to optimize the SC matrix based on the FC matrix. A dramatic improvement of the fitting of the matrices was obtained with the addition of a small number of anatomical links, particularly cross-hemispheric connections, and reweighting of existing connections. We suggest that the notion of a critical working point, where the structure–function interplay is maximal, may provide a new way to link behavior and cognition, and a new perspective to understand recovery of function in clinical conditions.

## Introduction

Significant efforts have been made to understand how the unique structural organization of the brain both constrains and enables its function and the resulting behavior. Using methods such as DTI/DSI (diffusion tensor imaging/diffusion spectrum imaging) to estimate the anatomical skeleton (Wedeen et al., 2005; Alexander and Lobaugh, 2007), and fMRI to measure the functional dynamics of whole-brain networks, it is possible to directly explore the relation between the two. Two regions that are anatomically linked tend to show a functional connection (Honey et al., 2009). The reverse is not true: a functional connection does not guarantee a direct structural link, primarily because of the indirect polysynaptic connections between regions. This relation is most evident in the study of resting-state networks (RSNs), where highly structured spatiotemporal patterns of intrinsic brain activity can be measured (Biswal et al., 1995; Damoiseaux et al., 2006; Fox and Raichle, 2007; Mantini et al., 2007; Brookes et al., 2011).

Several computational modeling studies have explained the emergence of RSNs by linking the fluctuating dynamics to the underlying anatomical structure (Honey et al., 2007; Ghosh et al., 2008; Cabral et al., 2011; Deco et al., 2011). One mechanism for the emergence of RSNs is the operation of the network at the edge of instability (Ghosh et al., 2008), when the parameters of the global system, such as interarea coupling, are such that the dynamical working point of the system operates at the edge of instability or criticality. At that particular critical point, the models maximally reproduce the empirical functional connectivity (FCs) (Deco and Jirsa, 2012). The global network dynamics at the critical instability amplifies the underlying structure of the anatomical connections and its interactions with the local dynamics. If we accept this point, then in principle, the empirical FC information, together with a computational model, could be used to estimate the structural connectivity (SC).

In this article, we define a self-consistent iterative procedure, with which we optimize the fit of the SC to the empirical FC. The basic idea was to redefine the strengths of the structural links successively according to the empirical FC. If, for a given SC framework, the model is not able to reproduce the FC between specific node pairs, then the strengths of those links can be redefined based on these FC values, and goodness-of-fit reassessed. The process is repeated iteratively, reducing the tolerance (*T*) of the difference between empirical and simulated FC links. The procedure continues until the fitting at the bifurcation is maximal (while ensuring that the added SC links are meaningful), giving the final enhanced SC matrix.

## Materials and Methods

#### Anatomical structure

##### Humans.

For the analysis of human data, we used a SC matrix derived from five healthy right-handed male human subjects using diffusion spectrum imaging (DSI) and tractography (Wedeen et al., 2005; Hagmann et al., 2008) based on a parcellation of 66 gray matter cortical areas. Each element of this SC matrix expresses the density with which two different brain areas are connected. That density is the normalized number of detected tracts linking those two brain areas. The neuroanatomical matrix was then averaged across the subjects.

##### Macaque.

For the analysis of macaque data, we used the CoCoMac database (http://cocomac.g-node.org), which is a compilation of published tracing studies (Stephan et al., 2001; Kötter, 2004). The brain areas of the monkeys were parcellated into 82 cortical regions (using the Regional Map of Kötter and Wanke, 2005). Each element of this SC matrix indexes the density of the connection between two brain areas. Those labels distinguish between “No,” “Weak,” “Moderate,” “Strong,” and “Unknown” connections. For the current application, Unknown connections were relabeled as zero, and unlabeled existing connections were assigned a rating of Moderate. Unlike the SC matrix from human DSI data, the SC matrix from CoCoMac is not symmetric.

#### Functional connectivity

##### Humans.

For the analysis of human subjects, we used a resting-state FC matrix obtained for the same five human subjects and using the same 66 cortical areas. The resting-state FC was calculated by measuring the corresponding fMRI BOLD signals during 20 min. An empirical FC matrix was obtained after regressing out the global signal and averaging across subjects (for details, see Honey et al., 2009).

##### Macaque.

For the analysis of the macaque data, we used a resting-state FC matrix obtained from six macaque monkeys using the 82 parcellated cortical areas adopted in the CoCoMac-based SC (Shen et al., 2012). Functional data were obtained from two 10 min “resting-state” scans, where the monkeys were lightly anesthetized (Hutchison et al., 2012). Head motion parameter estimates and BOLD estimates from white matter and cerebrospinal activity were used to regress out the physiological signals. Global regression was not applied. Matrices were averaged across the six monkeys.

#### Dynamic mean-field model

Here we used a dynamic mean-field (DMF) reduction of the detailed spiking/synaptic model (Deco and Jirsa, 2012), as implemented in The Virtual Brain platform (http://www.thevirtualbrain.org; Ritter et al., 2013; Sanz Leon et al., 2013). This DMF model is a mathematically consistent reduction of a spiking attractor network consisting of integrate-and-fire neurons with excitatory (NMDA) and inhibitory (GABA-A) synaptic receptor types (Brunel and Wang, 2001). These neurons are organized into an inhibitory population (20% of the neurons) and an excitatory population (80% of the neurons). The global brain dynamics of the network of interconnected local networks can be simply and consistently described by the following set of coupled differential equations:
where *H*(*x _{i}*) and

*S*denote the population rate and the average synaptic gating variable at the local cortical area

_{i}*i*,

*w*= 0.9 is the local excitatory recurrence,

*G*is a global scaling parameter which is optimized so that the system works at the edge of instability (see next section), and

*C*is the SC matrix expressing the neuroanatomical links between the areas

_{ij}*i*and

*j*. Parameter values for the input–output function are

*a*= 270,

*b*= 108 (Hz), and

*d*= 0.154 (s). The kinetic parameters are γ = 0.641/1000 (the factor 1000 converts the results to ms), τ

*= 100 (ms). The synaptic couplings are*

_{s}*J*= 0.2609 (nA) and the overall effective external input is

_{N}*I*

_{0}= 0.3 (nA). In Equation 1, υ

*is uncorrelated standard Gaussian noise, and the noise amplitude at each node is σ = 0.001 (nA).*

_{i}The BOLD fMRI signal was derived by first simulating the neuronal and synaptic activity for each region with the large-scale DMF model, and then using the Balloon-Windkessel model to generate the BOLD signal (Buxton and Frank, 1997; Mandeville et al., 1999).

#### Iterative self-consistency enhancement of SC based on FC empirical data

The enhancement of the original SC matrix is based on the addition of new links for pairs of nodes according to the corresponding FC between those nodes. The algorithm starts by first normalizing the maximal value of the empirical FC matrix (*FCe*) to 1. Then, the initial tolerance level is set to 1 and, starting with the original empirical SC matrix at the corresponding bifurcation (instability of the spontaneous state), the simulated FC matrix (*FCs*) is created. In doing this, the global scaling parameter *G* is adjusted such that the system works at the edge of the bifurcation where the spontaneous state loses stability (as shown by Deco et al., 2013, at that particular point the system optimally fits the FC and maximally reflects the underlying SC). Then, the tolerance is reduced (in our case by steps of 0.025, i.e., *T* = *T* − 0.025) and the simulated and empirical FC is checked for all connections. For all connections where |*FCe _{ij}* −

*FCs*| >

_{ij}*T*, the identified SC links are redefined to be equal to

*SC*= 0.15

_{ij}*FCe*if

_{ij}*FCe*> 0 and if

_{ij}*FCe*≤ 0,

_{ij}*SC*= 0, when it was zero in the original matrix; otherwise, if the link was not zero, it was scaled down to a minimum value of weights (0.0005). The starting weight of 0.15 was derived from the weight distribution in the DSI dataset, which was between 0.0005 and 0.18. Simulations were run with different starting values between 0.1 and 0.3 and yield the same final result. Thus, not only can new links be generated or reinforced, but existing links can also be weakened. If a strong link in the structural matrix generates a strong but inconsistent functional correlation, it could be reduced, but not completely eliminated, to better match a weak empirical functional correlation. For this new SC matrix, a new bifurcation is calculated and the corresponding BOLD signals and corresponding best-simulated FC are recalculated. The iterations continue for each new, lower tolerance level

_{ij}*T*=

*T*− 0.025 until

*T*= 0. By plotting the evolution of the correlation between those two matrices, we detect the level of tolerance at which the fit between the empirical and simulated FC is maximal. The argument is that new SC links that explain the FC are kept when the simulation is at the critical point. If the added SC links do not improve FC fitting, the new links are adding inconsistent information and therefore should not be considered.

## Results

Figure 1 schematizes the results for human FC data using the DMF approach. Figure 1*a* plots the bifurcation diagram describing the stationary states (attractors) for the DMF model as a function of the global scaling coupling factor *G*. For small values of *G*, only one spontaneous stable state exists, characterized by low firing rates in all cortical areas. At a critical value of *G* (*G* ≅ 1.5), a first bifurcation emerges and a new multistable state appears with multiple attractors and higher firing rates, while the spontaneous state is still stable. A second bifurcation appears at *G* = 2.5, characterized by a loss of stability in the spontaneous state. Just before this point, the simulations best match the empirical FC.

The DMF model demonstrates that as the global spontaneous state loses stability, the fluctuations of the ongoing activity of different cortical areas are structured and shaped by the underlying anatomical structure, yielding spatial correlations that best characterize the resting FC. At the critical instability, the spatial correlations between the different cortical areas simulated in the model (simulated FC) maximally reflect the underlying SC (Fig. 1*c*). In other words, the fact that the system operates at the brink of that instability is the main reason why noisy fluctuations around that point recapitulate the anatomical structure. The notion that the brain operates at the edge of criticality is not a new one (Beggs and Plenz, 2003; Haken, 2006; Kinouchi and Copelli, 2006; Aburn et al., 2012; Beggs and Timme, 2012), and what we reinforce here is that this critical point can be captured in a global attractor model of the whole brain (Deco and Jirsa, 2012).

With the DMF modeling framework in place, we then applied the iterative fitting to the human dataset. A problem with the SC matrix derived from DSI tractography is that many interhemispheric connections are absent. Figure 2*a* plots the fit between the empirical and simulated FC for different levels of tolerances, *T*. At the first iteration (*T* = 1), the results correspond to the best fit possible at the critical bifurcation with the original empirical SC matrix (compare Fig. 1*c*). For smaller values of *T*, links are added to the SC matrix and eventually the fitting is improved, meaning that the new links are self-consistently enhancing the SC and FC matrices. The maximal value observed in this subpanel (at approximately *T* = 0.2) corresponds to the optimally enhanced SC matrix. Figure 2*b* shows the percentage (with respect to the number of links in the empirical SC matrix) of new links added by the algorithm for successive tolerances *T*. Note that at the maximum of Figure 2*a*, the number of added missing intrahemispheric and interhemispheric connections is relatively very small (10% and 15%, respectively). Figure 2*c* shows the original and optimally enhanced SC matrices, whereas Figure 2*d* shows the original empirical FC and optimally simulated FC matrices. Most of the new links are interhemispheric between homologous areas (a comparable result has been obtained using a somewhat different formulation in the work of P. A. Robinson, S. Sarkar, G. M. Pandejee, and J. Henderson, unpublished observation). With the addition of these small numbers of extra links, we obtained a substantive improvement of the fitting of the FC matrix, from 0.4 (Fig. 1) to 0.75 (Fig. 2).

The fitting algorithm was next applied to an fMRI dataset taken from macaque monkeys (Hutchison et al., 2012; Shen et al., 2012). Because the link strengths were coarsely classified into just three categories, the quality of the fitting obtained with the original SC matrix was very poor (i.e., <0.1). Applying the iterative fitting algorithm, new links were added and the original link strengths were redefined based on the differences between the empirical and simulated FC matrices (Fig. 3). We observed that a relative low number of extra links improved the quality of the FC fitting notably (from 0.07 to 0.6), but the most significant contribution to fitting was the rescaling of link strengths from categorical to continuous values. After the enhancement, the weights are more exponentially distributed, as has been noted for node strength in human DSI data (Hagmann et al., 2008), and is consistent with scale-free architecture (Strogatz, 2001). We did an assessment of the relation between SC weights and distance and noted a weak but negative correlation (−0.15), which might be anticipated from recent macaque tracer studies (Markov et al., 2011). As more principles governing SC connection strengths become available, they could be incorporated into the optimization algorithm reported here.

Last, we used the DMF model to provide a stronger proof of concept that the procedure is indeed extracting a plausible reconstruction of the underlying SC. To do this, we started with an empirical, “ground truth” FC matrix (Fig. 2*d*, original FC) based on an optimal enhanced SC matrix (Fig. 2*c*, enhanced SC). To mimic the approximations to the SC matrix that may come from estimations using iterative DMF, we degraded the SC matrix by randomly erasing a percentage of links from 10% to 80% and modeled the FC using the iterative DMF procedure under conditions of criticality as before. Figure 4*a* shows, for different levels of degradation, the evolution of the fitting between the “empirical” and simulated FC (blue curve) as a function of tolerance. For each case we also plotted the fit between the corresponding recovered enhanced SC and the original SC (red curve). In all cases, we see that the maximum fit between the FCs and SCs is at a similar level of tolerance (dotted vertical line). Although the recovered connections may not be identical across degradations, the iterative DMF was able to recover a reasonable proportion of connections in the original SC matrix (*r* > 0.7). This is also shown in Figure 4, *b* and *c*, where the initial (triangle with dotted line of best fit) and best fits (circles with solid line) are plotted as a function of degradation level. For all degradation levels tested, the algorithm was able to achieve the same level of fit that was found using the original FC and original SC.

## Discussion

We have demonstrated in three different scenarios that it is possible to use a measure of brain network dynamics, FC, to infer the underlying SC. A central point here is that inference of SC from FC is possible only when the global parameters of the system are such that the system sits close and below the critical point of a bifurcation. Above or too far below this point, the relationship we present between structure and function breaks down. The most optimal working point is at the brink of the second bifurcation (Fig. 1), because it enhances its capacity to be variable by reducing the stability of existing states, thus offering an advantage through enhancing flexibility and sensitivity to incoming information. This finding reinforces the synergy of brain structure and function, and the proposition that FC, particularly as captured in RSNs, reflects the exploration of potential functional networks around the anatomical skeleton defined by SC (Deco et al., 2011).

The structure–function link we identified is valid only at the particular spatiotemporal scale explored here (i.e., low-frequency dynamics measured in BOLD fMRI). Indeed, recent empirical work on RSNs has noted that across long time windows, distinct functional networks form that relate to the structural connections. At shorter intervals, however, subnetworks form and dissolve as the full network evolves. The regions within each of the broad spatial patterns would move away from their core network to form other networks, and then transition to a new set of spatial patterns (Hutchison et al., 2013).

Apart from the exciting prospect of a new method to enhance connectomics, the notion of criticality as an operating point of brain function has significant implications. First, given the improvement in the SC matrix derived from human DSI data (Fig. 2), there is a great potential for using the approach presented here on individual tractography and FC data to suggest structural connections that are not well measured with DSI (for a complementary approach, see Hinne et al., 2012). This is not meant to supplant anatomical verification of the connection but, rather, to act as a pointer for further study. Second, there is a natural extension of our approach to clinical studies. Following a lesion, brain networks reorganize themselves structurally and functionally. In most cases, the lesioned brain network may not be at the critical operating point anymore. If our hypothesis is correct, then using the notion of the critical point may allow for the prediction of structural reorganization leading to functional and behavioral recovery.

## Footnotes

G.D. was supported by the European Research Council Advanced Grant DYSTRUCTURE (n. 295129), by the Spanish Research Project SAF2010-16085, and by the CONSOLIDER-INGENIO 2010 Programme CSD2007-00012. V.K.J. and G.D. are supported by FP7-ICT BrainScales. The research reported herein was supported by Collaborative Research Grant 220020255 from the James S. McDonnell Foundation. P.H. is supported by the Leenaards Foundation.

- Correspondence should be addressed to Anthony R. McIntosh, Rotman Research Institute of Baycrest Center, University of Toronto, Toronto, ON M6A 2E1, Canada. rmcintosh{at}research.baycrest.org