## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Relationships between cortical myeloarchitecture and electrophysiological networks

Edited by Marcus E. Raichle, Washington University in St. Louis, St. Louis, MO, and approved October 7, 2016 (received for review May 27, 2016)

## Significance

This paper identifies a significant relationship between cortical myeloarchitecture and functional connectivity in the human brain. This is a significant step toward understanding the role of myelin in shaping large-scale neural networks. Our results extend recent work showing that electrical activity promotes myelination and add significant weight to the argument that neural oscillations are a core mediator of brain connectivity.

## Abstract

The human brain relies upon the dynamic formation and dissolution of a hierarchy of functional networks to support ongoing cognition. However, how functional connectivities underlying such networks are supported by cortical microstructure remains poorly understood. Recent animal work has demonstrated that electrical activity promotes myelination. Inspired by this, we test a hypothesis that gray-matter myelin is related to electrophysiological connectivity. Using ultra-high field MRI and the principle of structural covariance, we derive a structural network showing how myelin density differs across cortical regions and how separate regions can exhibit similar myeloarchitecture. Building upon recent evidence that neural oscillations mediate connectivity, we use magnetoencephalography to elucidate networks that represent the major electrophysiological pathways of communication in the brain. Finally, we show that a significant relationship exists between our functional and structural networks; this relationship differs as a function of neural oscillatory frequency and becomes stronger when integrating oscillations over frequency bands. Our study sheds light on the way in which cortical microstructure supports functional networks. Further, it paves the way for future investigations of the gray-matter structure/function relationship and its breakdown in pathology.

The way in which integration of functionally specific brain regions supports ongoing cognition is one of the most important questions in neuroscience, and noninvasive in vivo imaging provides a tool to investigate this interregional connectivity in terms of both brain function and structure. Functional connectivity refers to statistical interdependencies between patterns of brain “activity” measured at separate cortical locations (1) and, even in the “resting state,” measured spontaneous brain activity defines nonrandom networks that are related to cognitive processes (2). The way in which these functional networks are supported by structural white-matter pathways is reasonably well understood (3). However, it is likely that the structure–function association extends to gray-matter morphology, for which fundamental understanding is lacking. Structural morphology of the cortex is known to vary significantly between individuals, and does so in an organized fashion. For example, individuals with a high cortical volume in Broca’s area typically exhibit high cortical volume in Wernicke’s area, reflecting a language network (4). Similar observations can be made between other associated cortical regions (5, 6). This is known as structural covariance (7⇓–9) and it allows the formation of matrices showing how structural properties of individual brain regions covary over subjects. In this paper, we assess structural covariance based upon cortical myeloarchitecture and probe its relationship to functional networks that are assessed based upon measured spontaneous brain activity.

Noninvasive mapping of cortical myeloarchitecture (10, 11) has grown in popularity in recent years, fueled by an increased interest in myelination-based parcellation (12). The general finding is that primary sensory cortices tend to be heavily myelinated whereas regions associated with multisensory integration are less myelinated (10, 13). More localized spatial changes in myelination also exist; for example, subtle subdivisions between regions are apparent in visual (14, 15), somatosensory (16), and auditory (17) cortices. There is evidence of cross-species changes in myelination (18), with the brains of nonhuman primates more heavily myelinated than those of humans. Although still unproven, findings suggest that cortical myelin may inhibit plasticity; for example, early sensory areas may require less plasticity, and therefore more myelin, whereas higher-order areas have less myelination, which might enable greater plasticity (10). Cortical myelination has been shown to change throughout development, and is not established fully until the third decade of life (19, 20). However, no studies have yet used myelin mapping to examine structural covariance or the link between a cortical myelin network and functional connectivity. Previous work does show a direct link between myeloarchitecture and function (21⇓–23). For example, a recent study (24) showed that optogenetic stimulation directed to increase neuronal firing in premotor cortex of mice promotes oligodendrogenesis and therefore myelination, thus providing a link between neuroelectrical activity and myeloarchitecture. Further work (25) suggests that the amount of cortical myelin in a region predicts the magnitude of electrophysiological responses. Taken together, the evidence converges to a hypothesis that, if neuronal firing acts to shape the spatial signature of myeloarchitecture, then networks reflecting the brain’s primary pathways of functional connectivity should be predictive of structural networks of intracortical myelin.

Magnetoencephalography (MEG) characterizes electrical activity in the brain via measurement of extracranial magnetic fields (26). The MEG signal from any one brain region is dominated by neural oscillations (rhythmic changes in electrical activity) that are observable in the 1- to 200-Hz frequency range. Evidence suggests that these oscillations represent an intrinsic process by which both short- and long-range functional connections in the brain are maintained. With this in mind, a growing body of work has begun to show that, via appropriate modeling of MEG data, networks of electrophysiological functional connectivity can be mapped (27⇓–29). The rich temporal complexity of MEG signals means that multiple ways to characterize functional connectivity exist (30). However, one of the most robust methods is amplitude envelope correlation (AEC), which probes temporal relationships between the envelope of oscillations, in a frequency band of interest, at spatially separate brain regions. Our previous work has shown that this intrinsic mechanism of (noninvasively measured) electrophysiological coupling underpins many of the observable resting state networks (RSNs) (27, 31). Given that these measurements represent the principal long-range functional connections in the brain, and given the evidence that electrical activity mediates myelination, we hypothesized that networks of oscillatory envelope correlation, measured between parcellated regions and in multiple frequency bands, would allow prediction of a network of structural covariance representing myeloarchitecture.

## Results

Fifty-eight volunteers (39 ± 12 y old, 27 male) took part in the study. Using ultra-high field MRI, we measured magnetization transfer (MT) (32) across the brain, which serves as a marker of myelin. These measurements were parcellated into 64 cortical regions according to the automated anatomical labeling (AAL) atlas (Table S1) and normalized by region volume to give myelin density estimates. Following this, assessment of structural covariance between all possible AAL region pairs allowed derivation of a network matrix showing the degree to which separate (AAL) regions exhibit similarities in their myeloarchitecture (see *SI Further Analyses* and Fig. S1 for an alternative methodology). The same individuals also underwent resting-state MEG acquisition. MEG data were parcellated according to the same atlas, and RSNs characterizing functional connectivity, between all possible AAL region pairs, were derived using AEC in five frequency bands. With the aid of a recently developed framework (31), we then characterized the relationship between the structural network, representing myelin density, and functional networks, representing the major pathways of electrophysiological communication.

### Myelin Maps.

Fig. 1*A* shows the measured MT, averaged across subjects and plotted for all AAL regions. Red corresponds to high myelin density, blue indicates low myelin density, and gray shows regions where scan coverage was insufficient to gain an accurate estimate. High myelination was observed in primary sensory cortices, with statistical tests showing that this spatial variation was significant; for example, significantly (*P* < 0.001) higher myelination was found in primary motor cortex compared with the dorsolateral prefrontal cortex. A hemispheric division was also noted with significantly (*P* < 0.01) higher MT in left compared with right sensorimotor cortex. Fig. 1*B* maps the cortical variation of correlation between myelin density and handedness (the latter measured using the Edinburgh Handedness Inventory). A positive correlation denotes regions where right handers have more myelin than left handers; a negative correlation denotes regions where left handers have more myelin than right-handed individuals. Note a significant (*P* < 0.05) split in the polarity of the correlation between hemispheres.

### MT Structural Covariation.

The principle of structural covariance and the structural network are shown in Fig. 1 *C* and *D*, respectively. Fig. 1*C* shows correlation across subjects between myelin concentrations measured in left frontal inferior operculum (capturing Broca’s area) and left supramarginal gyrus (encompassing Wernicke’s area). Results show a significant (*P* < 0.001) positive correlation across subjects. Equivalent correlations can be derived between all possible AAL region pairs, and results are shown in Fig. 1*D*. Note that structural covariance between adjacent AAL regions is generally denoted by matrix elements close to the leading diagonal, whereas structural covariance between distal regions is represented far from the diagonal. The white boxes show structural covariance between homologous regions.

### MEG Networks.

MEG functional connectivity matrices are shown in Fig. 2. All networks are represented by a matrix similar to that in Fig. 1*D*. The 3D brains display all connections within 5% of the maximum value. It is clear that network structure differs between bands: Theta oscillations support connections in frontal occipital and parietal areas whereas alpha band-mediated connections are predominantly occipital. The beta band shows more widespread connectivity but is dominated by parietal and occipital connections. The low-gamma band was dominated by a sensorimotor network.

### The Relationship Between Structural Covariance and MEG.

Fig. 3*A* shows a “seed-based” structural covariance map (i.e., a single column in the matrix in Fig. 1*D*). Here, a seed region has been placed in right inferior parietal cortex; high values depict brain areas that show high structural covariance to the seed. A cross-hemispheric pattern is observable with high structural covariance between homologous regions. Fig. 3*B* shows the equivalent seed-based map calculated using beta band MEG data. Note the strong similarity between the functional (Fig. 3*B*) and the structural (Fig. 3*A*) networks. To generalize this relationship for all possible seed regions, we tested for correlation between the full MT matrix (representing all seeds; Fig. 1*D*) and the group-averaged functional connectivity matrices for all bands (Fig. 2). The resulting r^{2} values (bar chart in Fig. 3*C*) show that functional networks measured in the beta and low-gamma bands predict significantly the spatial pattern of structural covariance. Neither theta nor high-gamma bands showed a measurable relationship; the alpha band showed a trend. The inset images show seed regions for which structural covariance is best predicted by functional networks. Note that the structure/function relationship is strongest in parietal and occipital areas and weakest in the frontal lobes. It is noteworthy that the MEG-derived functional networks, particularly in the beta band, are driven in part by canonical RSNs (27), and by extension this suggests a significant relationship might also be found between the spatial signature of fMRI-derived RSNs and the structural network. This is indeed the case, and this significant relationship is shown in *SI Further Analyses*, *Functional Connectivity and Myeloarchitecture in RSNs*).

Although functional networks in individual frequency bands show significant correlation, it is likely that the structural network exists to support functional connectivity in all bands. For this reason, we sought to integrate the five MEG networks to test whether such combination could better predict structure than independent frequency bands. Two approaches were used. First, all five MEG matrices were combined in a linear weighted sum. Second, these same matrices were supplemented by nonlinear terms, formed based upon the square of each MEG matrix, and again a weighted sum derived. Importantly, the nonlinear terms have specific meaning: For any squared matrix, a single element, say [1, 2], represents the inner product of the connectivity profile of region 1 and region 2. Because this product is related to covariation, the matrix element will be high if the connectivity profile of region 1 to the rest of the brain overlaps with the equivalent connectivity profile of region 2. In this way the squared terms can be thought of as representing brain regions that share connections to similar areas. Fig. 4 *A*–*C* show connectivity matrices representing the structural network (Fig. 4*A*) and its prediction based upon linear (Fig. 4*B*) and nonlinear (Fig. 4*C*) combinations of MEG networks. These relationships, along with that for the best single frequency band, are further visualized in Fig. 4 *D* and *E*, which show “seed-based” structural covariance (top row) alongside equivalent maps made using the beta band (upper middle), the best linear combination (lower middle), and the best nonlinear combination (bottom). Seeds were placed in right lateral visual cortex and left superior frontal cortex in Fig. 4 *D* and *E*, respectively.

The relationship between structure and integrated functional connectivity is formalized in Fig. 4*F*, which shows r^{2} values representing correlation between the MT network and functional networks representing beta band only and linear and nonlinear predictions; the inset images show seed regions for which structural covariance is best predicted by the functional networks. Fig. 4 *G*–*I* show the r^{2} values plotted alongside their representative null distributions. In all cases, correlation between the structural network and MEG networks falls outside the null distribution. Note also that this relationship gets stronger when allowing integration over frequency bands.

## Discussion

Although recent years have seen significant progress in mapping the human connectome, the relationship between functional networks and cortical microstructure remains poorly understood. A recent study (24) in animals suggests that electrical activity promotes myelination. We therefore reasoned that, if functional networks represent major pathways of electrophysiological communication, then those pathways should shape myeloarchitecture and a significant relationship between functional connectivity and myelin should be observable. Our results support this, with significant correlation between the structural network and functional networks mediated by neural oscillations in the beta and low-gamma bands. This relationship became stronger when integrating MEG networks across frequency bands, suggesting that myeloarchitecture supports networks at all measurable electrophysiological timescales.

The fact that our functional networks are measured directly using electrophysiological imaging (as distinct from indirectly using hemodynamics) adds an extra (frequency) dimension to our study. In the past, neural oscillations were largely ignored in favor of measurements of evoked responses. However, the last decade has seen a surge of interest, showing these oscillations to be an integral feature of brain function. Recent work suggests that oscillations gate information flow in the cortex (33) and implies that oscillations are an intrinsic form of functional coupling (34). Measurable connectivity depends critically on the frequency band studied; indeed, this is shown in Fig. 2 with marked spatial differences between bands. The fact that the relationship between functional networks and myelin was strongest in the beta band is not surprising given that previous work (27, 31) has shown that beta oscillations mediate long-range connections in a large number of RSNs (see also *SI Further Analyses*, *Functional Connectivity and Myeloarchitecture in RSNs* and Fig. S2). The ubiquitous nature of beta-mediated connections is therefore the likely reason why the structural network correlates best with this band. However, it should be pointed out that an inherent problem with MEG is signal-to-noise ratio (SNR), which drops with increasing frequency; it is possible that, at high frequency, the drop in correlation between functional connectivity and myeloarchitecture reflects this SNR limitation. This said, the significant improvement in prediction when combining frequency bands does imply that oscillations across all timescales relate significantly to myelin structure. Interestingly, the addition of nonlinear terms also improved prediction; given that these additional terms represent brain regions that share common connections it is tempting to suggest that shared connectivity also affects myeloarchitecture.

Our estimation of myelin was based upon measurement of MT. Although serving as an efficient marker of myelination, it is important to note that there is no one-to-one relation between MT and myelin density. Further, if the present method was to be used in pathology care should be taken because the MT/myelin relationship may break down. Nevertheless, our results are in agreement with others, showing that the highest myelin concentration occurs in sensorimotor, auditory, and visual cortices. This finding supports an argument that myelin acts as a means to increase the speed of processing and inhibit plasticity in these primary sensory areas. On average, there was more myelin in left sensorimotor cortex than in right, and this likely reflects the fact that our cohort was biased toward right handers. In agreement with this, the hemispheric split shown in Fig. 1*B* suggests that this finding is reversed in left handers, implying that brain structure evolves to increase the speed of processing in the dominant hemisphere. Although such asymmetries in myelin have not been previously reported, associated asymmetries in function and structure have been shown previously both in humans (35⇓–37) and animals (38). A potential limitation of our methods is that myelin concentration was parcellated according to the AAL atlas. Although this parcellation allowed derivation of network graphs that could be compared with MEG, the regional parcellation afforded by the AAL atlas does not take into account known differences in myelin signature throughout the cortex. In future studies, a brain parcellation based upon myeloarchitecture (as distinct from cytoarchitecture) would be of high value (12). Although our study is unique in characterizing structural covariance based upon myelin density, previous work has shown that structural covariance networks can also be formed based on macroscopic properties such as cortical thickness (4), and further that such networks correlate with functional connectivity (39). Here, myelin content was normalized by region volume, meaning that our finding of correlation between a structural network and functional connectivity is not related trivially to the previous finding on cortical thickness. However, the question of how myelin density relates to cortical thickness remains, and this should be the topic of future work.

Given that the primary role of myelin is to increase the speed at which nerve impulses travel through neuronal pathways, it is intuitive that cortical myelination is shaped to support functional networks because this will act to maximize the efficiency of their formation. More speculative, however, is the process by which this structural support network evolves. Given that electrical activity on some particular pathway induces myelination, and given that RSNs emerge as early as the third trimester of gestation (40), it is tempting to argue that even before birth myeloarchitecture is being shaped by functional connectivity. Of course, the resulting structural changes would, in turn, refine functional connections and so the likelihood is that changes in myeloarchitecture and functional networks are linked intimately. This notion is supported by the fact that both myelination and functional connectivity change on a similar timescale throughout neurodevelopment. In the present study, our data preclude direct investigation of this interplay between structure and function. However, our study does pave the way for new investigations of this process via longitudinal studies of development. In addition, studies investigating how changing behavior alters both structure and function are becoming popular. This idea is not new (41), and recent evidence (42) suggests that both white-matter and gray-matter structure changes, even on a relatively fast timescale, when learning new skills. This research area would benefit from the use of tools presented here. Perhaps more importantly, our results have implications for future studies of disorders, in particular those involving demyelination or dysconnectivity. For example, gray-matter atrophy, lesions, and demyelination are a better correlate of physical disability and cognitive decline in multiple sclerosis than white-matter lesion load (43). In addition, inefficiency of functional networks has also been reported (39) in this disease. Our method might offer a means to link these findings. Similarly, severe psychosis has been linked with dysconnectivity (44), and recent work has begun to relate this to structural deficiencies (45). Once again, our work might offer a framework to link these abnormalities.

## Conclusion

We have probed the relationship between gray-matter myelination and electrophysiological networks, showing a significant correlation. This relationship is strongest for networks mediated by beta oscillations but becomes stronger when integrating across frequency bands, suggesting that myeloarchitecture supports connectivity across all bands. Our study sheds light on the way in which cortical microstructure supports functional networks, the latter being mediated by neural oscillations. Further, it paves the way for future investigations of the structure/function relationship and its breakdown in pathology.

## Methods

### Myelination and Structural Covariance.

Participants gave written informed consent and ethical approval was granted by the University of Nottingham Medical School Research Ethics Committee. MRI data were collected using a Phillips Achieva 7T system. A phase-sensitive inversion recovery T_{1} weighted image was acquired and used for MEG coregistration. MT data were obtained from z-spectra acquired using an MT-TFE sequence (46). MT imaging provides contrast based on the exchange of magnetization between free water and protons bound in macromolecules. Although it is likely that no one-to-one relationship exists, experimental and human studies have shown that MT is highly correlated with myelination, which is probably related to the high fraction of water in close proximity to myelin macromolecules. The procedure for extracting MT data from our imaging sequence has been described elsewhere (47). Briefly, z-spectra were corrected for B_{0} variation and fitted to a database of simulated spectra to extract myelination maps. To investigate structural covariation of myelin within the AAL regions, gray-matter-masked MT data were registered to the AAL atlas and a mean MT value was calculated for each region, for each participant. Due to confounding factors such as scanner drift, the modal value for each individual’s MT data was regressed from that individual’s regional values (7) (Fig. S3). Pearson correlation, measured across subjects, was used to quantify the relationship between MT values measured in AAL region pairs (Fig. 1*C*). These correlation values form elements of the structural matrix in Fig. 1*D*. An average MT map (Fig. 1*A*) was also generated by averaging regional MT values over participants. The relationship between myelination and handedness was probed via correlation between MT and handedness score.

### Resting-State MEG: Connectivity Analysis.

Three hundred seconds of eyes-open resting-state MEG data were acquired using a 275-channel CTF MEG system operating in third-order synthetic gradiometer configuration (sampling frequency of 1,200 Hz). Three head position indicator coils were placed at fiducial locations on the subject’s head and energized to facilitate continuous tracking of head location. Before acquisition, a 3D head digitization procedure was completed. Coregistration between MEG system geometry and individual brain anatomy was achieved by matching the digitized head surface to the equivalent surface extracted from an anatomical MRI. Functional connectivity was calculated between AAL regions. A scalar beamformer was used to obtain a single MEG signal representative of each region. These regional signals were frequency-filtered to the bands of interest and the confound of signal leakage mitigated using pairwise orthogonalization (48). The magnitude of the analytic signal, computed via a Hilbert transform, was used to generate the amplitude envelope of oscillations, for each regional time course in all frequency bands. Pearson correlation was then computed between envelopes for each region pair. In this way we generated a single adjacency matrix for each subject, and frequency band, representing whole-brain connectivity. These matrices were averaged over subjects.

### The Relationship Between MT and MEG.

We measured correlation (r^{2}) between the structural covariance matrix (Fig. 1*D*) and the group-averaged MEG networks (Fig. 2). To test statistically whether these correlation values were significant, we used a permutation test. A set of “pseudo-MEG matrices” were generated, each having spatial properties similar to the real matrices, but crucially they were not based on genuine data (*SI Methods* and Figs. S4 and S5). This approach accounts for the inherent spatial smoothness in the MEG-derived networks (*SI Methods*). Correlation between the structural matrix and 10,000 pseudomatrices yielded an empirical null distribution, and comparison of the genuine r^{2} value with the null distribution allowed computation of a *P* value.

Our supplementary information is split into three sections. In *SI Further Analyses* we describe additional analyses that, due to space constraints, were excluded from the main text. *SI Results* contains further visualization of primary results. *SI Methods* describes, in greater detail, our methodology.

## SI Further Analyses

### Cross-Subject Strength Covariance and Its Relationship with Myeloarchitecture.

#### Motivation.

In the main text we tested the hypothesis that a network of structural covariance, measured based upon myelination, could be predicted by functional connectivity measured using MEG. To this end, we used functional connectivity matrices averaged over subjects (MEG) to predict a structural covariance matrix formed via assessment of cross-subject correlation in myelin. Although highly instructive, this requires that a structural covariance matrix, based on intersubject variance, be compared with functional connectivity, based on intrasubject variance. A logical extension to this work is therefore to exploit the intersubject variance in MEG and perform a parallel analysis, to assess whether (*i*) cross-subject correlation in MEG offers a means to generate structured network graphs and (*ii*) whether such graphs also relate to myeloarchitecture.

#### Method.

With this in mind, we first reduced the dimensionality of the MEG connectivity matrices by measuring connectivity strength at each of the 64 AAL regions. Connectivity strength was defined as the summed connectivity between a single seed region, and all of the other 63 regions (i.e., connectivity strength is measured as the sum, in one dimension, of the MEG adjacency matrices shown in Fig. 2 – this is sometimes also termed “degree” and is shown schematically in Fig. S1*A*, *Left*). These measures generated 64 values of connectivity strength for each subject. Measurement of correlation in strength across subjects for all region pairs then generated a novel MEG metric, somewhat analogous to structural covariance, which we term “cross-subject strength covariance.” This analysis results in a single 64 × 64 adjacency matrix per frequency band, where a high value of strength covariance between regions A and B denotes that subjects with a high connectivity strength in region A also tend to have a high connectivity strength in region B. Note that although the metric is analogous to structural covariance in terms of the way in which it is computed it is based purely on measurement of brain function.

#### Results.

Fig. S1*A* shows connectivity strength, calculated for all brain regions and averaged over subjects, for all five frequency bands. Note the clear spatial patterns that differ between frequency bands, alpha and theta being predominantly occipital, beta more widespread, and low gamma dominated by the sensorimotor regions. Note also the significant variation in strength across bands, peaking in the beta band, as would be expected. Fig. S1*B* shows matrices representing cross-subject strength covariance. Interestingly, a high degree of structure is observed which, although obviously related to the “classical” functional connectivity matrices in Fig. 2, also differs markedly. Here, a clear and dominating feature, in particular in the beta band, is interhemispheric coupling (i.e., subjects with high connectivity strength in one region tend to exhibit similarly high connectivity strength in homologous regions of the opposite hemisphere). This is shown by the matrices and by the inset images. Images depict seed-based strength-covariance profiles (single columns in the matrices in the upper panel). Arrows mark the seed AAL region. Note that for frontal, parietal, temporal, and occipital seeds the highest strength covariance is in homologous areas of the opposite hemisphere. Fig. S1*C* shows the relationship between strength covariance and structural covariance based upon myelination. Note that this was measured in the same way as was described in *SI Methods* below. The bar chart shows correlation (r^{2}) between the structural network (Fig. 1*D*) and the strength covariance adjacency matrices (Fig. S1*B*). Statistical tests (based on pseudomatrices) show this correlation to be significant. The inset images show the seed regions that drive the relationship in the bar chart (i.e., red indicates a region whose structural connectivity profile and strength covariance profile is highly correlated).

#### Discussion.

The principal finding in the main text shows that brain areas that are highly functionally connected also exhibit cross-subject covariation in myeloarchitecture. Given this finding, it is perhaps unsurprising that networks representing cross-subject strength covariance also demonstrate a significant relationship with the myelin network (with similar r^{2} values). To explain this, consider the simplest possible case where a network graph is dominated by a single connection between regions A and B, and that the magnitude of this connectivity differs over individuals. In such a case, strength measures would be dominated by regions A and B, and, further, those strength measures would obviously covary over subjects because they are generated by the same single connection. It follows that strength covariance matrices should bear significant resemblance to standard functional connectivity matrices, and this is indeed the case (e.g., the four patches representing occipital connectivity that dominate the alpha band matrix in Fig. 2 also dominate the alpha band strength covariance matrix in Fig. S1*B*). This argument suggests that the significant relationship between strength covariance and structural covariance shown in Fig. S1*C* is secondary to the principal finding of our paper. What is compelling, however, is that as well as obvious similarities between functional connectivity and strength covariance, clear differences are also apparent (e.g., the interhemispheric connectivity in beta band observed in Fig. S1*B* is much less pronounced in classical functional connectivity measures). The novelty of the strength covariance measure, coupled with this clear demarcation, suggests that these measurements could provide new insights into brain architecture beyond what the widespread AEC (and related) functional connectivity measurements provide.

### Functional Connectivity and Myeloarchitecture in RSNs.

#### Motivation.

The main text shows that a significant relationship exists between MEG-derived functional connectivity and a structural network based upon myeloarchitecture. However, our previous work (27) has shown that similarly derived MEG functional connectivity metrics, in part, represent the electrophysiological basis of RSNs commonly identified by fMRI. It thus follows that these canonical RSNs (e.g., sensorimotor, frontoparietal, default mode, etc.) should predict our network of structural covariance.

#### Method.

Ten commonly observed canonical RSNs were selected from a previously published fMRI study (49); these included the sensorimotor network, two visual networks, three frontoparietal networks, the default mode network, a bilateral frontal network, a bilateral parietal network, and an occipitoparietal network. To translate these volumetric weighted RSN maps into AAL space, we first multiplied each network map by a set of binary masks depicting each individual AAL region. We then summed (over voxels) the weightings for each RSN within each AAL region and normalized by the number of voxels in each region. This generated a single 1 × 64 vector, per RSN, showing the extent to which each AAL region is involved in that network. These vectors are visualized in Fig. S2*A*, *Upper*. Following this, to generate a matrix description of each RSN, the outer product of each 1 × 64 vector with itself was calculated; this generated 10 weighted adjacency matrices (one for each network), which are shown in Fig. S2*A*, *Lower*. Finally, to test our hypothesis that canonical networks would predict structural covariance, we constructed a general linear model:** S** represents the vectorized structural covariance network (Fig. 1

*D*),

**comprised 10 columns representing the vectorized RSN matrices, parameters in**

*G***, and**

*S***represents error in the fit. This approach allows a means to explain the structural network as a linear weighted sum of fMRI-derived RSNs. To quantify the success of the fit, Pearson correlation between the model,**

*e***, was calculated.**

*S*To test this GLM fit statistically, a permutation approach was used based on pseudomatrices. Ten thousand iterations of the algorithm were run and, on each iteration, the design matrix, ** G**, was used to predict a new pseudomatrix. Each pseudomatrix was constructed based upon the structural covariance matrix. Again, Pearson correlation between the model and the pseudomatrix was calculated, and these values were used to generate an empirical null distribution. Comparison of the real correlation with the null distribution allowed derivation of a

*P*value to determine the significance of the relationship between fMRI RSNs and structural covariance. In addition to Pearson correlation, for each iteration of the null distribution we measured the magnitude of the fitted parameters,

Finally, for completeness, this same procedure was used to model each of the five MEG adjacency matrices using the RSNs; this provided a means to replicate previous findings in interrogating which fMRI-derived RSNs are represented by MEG functional connectivity measurements.

#### Results and discussion.

Fig. S2*B* shows the results of predicting structural covariance based on RSNs. The left-hand side of the plot shows the structural covariance matrix (left-hand matrix) and the fit based upon a linear combination of RSNs (right-hand matrix). Note the visual similarity. The center panel shows the value of Pearson correlation between the model and the data (in yellow) and the associated empirical null distribution (in blue). Note that the real value of correlation falls outside the null distribution, indicating the significance of the fit and confirming our hypothesis above. The bar chart in the lower right-hand panel shows the parameters, *A* to the structural network; note that all RSNs contribute, with the largest contribution coming from the sensorimotor and frontoparietal networks. Our statistical test also conformed this; the network maps, inset in Fig. S2*B*, show the three networks that contributed most. ** indicates a significant contribution following correction for multiple comparisons; * indicates a trend (i.e., a *P* value less than 0.05 but that did not survive FDR correction). The fact that the largest contributions come from networks predominantly centered around the parietal and occipital regions (i.e., visual, sensorimotor, and frontoparietal networks) is in good agreement with results in the main text, which show the strongest relationship between functional connectivity and myeloarchitecture exists in the occipitoparietal regions.

Finally, Fig. S2*C* shows the equivalent analyses applied to each of the five MEG adjacency matrices (i.e., we test whether the MEG connectivity data can be predicted using RSN network maps). Results show a significant relationship for all five frequency bands. The strongest correlations were observed in the alpha, beta, and low-gamma bands (bar chart in Fig. S2*C*), and this is in good agreement with previous work (27). The inset images show networks that contribute maximally to each band. As would be expected, most networks are represented in the beta band, including visual, sensorimotor, and frontoparietal.

## SI Results

### The Effect of Common Mode Rejection.

To protect against confounding factors, such as scanner drift, the common mode was regressed from MT data (7). However, Fig. S3 shows the MT maps from the main text, made using uncorrected data. Fig. S3*A* shows correlation between MT and handedness. Positive correlations show regions where myelin is higher in right handers. Negative correlations show regions where myelin is higher in left handers. Fig. S3 *B*–*D* show seed-based structural covariance maps with seed locations in parietal, occipital, and frontal cortices, respectively.

### Linear and Nonlinear Prediction of MT.

Fig. S5 shows validation of our linear and nonlinear fits to the structural covariance matrix. Fig. S5*A* shows the five MEG matrices (from Fig. 2, Fig. S5*A*, *Left*) and their squares (Fig. S5*A*, *Right*). Note that for our linear fit only the five original matrices (Fig. S5*A*, *Left*) are combined, whereas the nonlinear combination incorporates a sum of all 10 matrices.

Fig. S5*B* shows our rate-of-improvement analysis. In Fig. S5*B*, *Right*, the yellow dashed line represents the improvement in r^{2} afforded by moving from the beta band only to first linear and then nonlinear fits, in real data. The blue lines show equivalent gradients for the null distribution, which was based upon pseudomatrices. In Fig. S5*B*, *Left*, these data are further visualized with the gradient for real data (yellow) plotted against the gradients for pseudodata (blue). Note that the null distribution shows positive gradients, meaning that the introduction of more terms in the model allows a better fit to the data, even in the absence of any genuine effect. However, the gradient for real data falls outside this empirical null, showing clearly that the improvements observed are in excess of what would be expected statistically. Fig. S5*C* shows the results of our cross-validation procedure for the linear combination (Fig. S5*C*, *Upper*) and nonlinear combination (Fig. S5*C*, *Lower*). Subjects were split into two groups, A and B; in both cases, the blue curve shows the r^{2} values for the estimation of MT data from group A, using linear and nonlinear combinations of MEG matrices, with parameters also derived from group A. Conversely, the red curve plots the distribution of the cross-validation r^{2} values; in other words, it shows the estimation of MT data from group B, using linear and nonlinear combinations of MEG matrices from group B, but with parameters derived independently from group A. The yellow curves show the null distributions (sham r^{2} values). Note that, as would be expected, a slight drop in r^{2} is noted when using independent parameters from a different subject group. However, importantly, this drop is slight and the median correlation falls outside the null distribution in both cases. This analysis shows clearly the robustness of the fits. Finally, Fig. S5*D* shows the relative contributions of the original and the squared terms to the nonlinear fit. The matrix on the left shows the five linear terms combined and the matrix on the right shows the five nonlinear terms combined. The sum of these two matrices generates the result in Fig. 4*C*. The bar chart shows the matrix norms corresponding to the linear and squared terms. Note that the nonlinear contribution is smaller.

### AAL Regions.

Table S1 lists the 78 AAL regions used for analysis in the main text. Following MT preprocessing, 14 regions (marked with an asterisk) were removed from further analysis.

## SI Methods

### Participants.

Seventy-seven healthy participants were recruited to the study. Nineteen participants were excluded because of either unsuitable MEG or MT data. Following preprocessing, a final cohort of 58 participants were included [mean age of 39 ± 12 y (mean ± SD), 27 male]. Participants completed an online screening form to assess health and lifestyle; this included the Edinburgh Handedness Inventory (50). MRI and MEG data were acquired in all subjects, in two scanning sessions. To avoid the effects of tissue magnetization, MEG data were acquired first.

### Resting-State MEG Data Acquisition.

MEG data were acquired using a 275-channel CTF MEG system (MISL) operating in third-order synthetic gradiometer configuration, at a sample frequency of 1,200 Hz. The resting-state paradigm comprised a 5-min recording during which participants fixated on a small centrally positioned red circle presented on a gray background that was back-projected onto a screen placed ∼40 cm from the participant’s face. Subjects were seated in the MEG system. Three head position indicator coils were placed at fiducial locations (nasion, left, and right preauricular points) on the subject’s head. These coils were energized periodically to facilitate continuous tracking of the head location throughout data acquisition. Before MEG acquisition, a 3D head digitization procedure was completed for each individual (Polhemus Inc.), characterizing head shape relative to the fiducial markers. Coregistration between MEG system geometry and individual brain anatomy images was achieved by matching the digitized head surface to the equivalent surface extracted from the participant’s individual anatomical MRI.

### Resting-State MEG Connectivity Analysis.

MEG data were down-sampled to 600 Hz and divided into 30 10-s epochs. All epochs were inspected visually and those containing excessive interference or head movement (>5 mm from starting position) were discarded.

Functional connectivity was calculated between 78 discrete cortical regions, defined based on the AAL atlas (51). To assess the time evolution of electrical activity within each region a scalar beamformer was used (52). Each AAL region was split into 4-mm cubic voxels and the beamformer-estimated time course of electrical activity was derived for each voxel. For beamforming, we calculated the data covariance matrix within a 1- to 150-Hz frequency window and a time window spanning the whole experiment (53). Regularization was applied using the Tikhonov method with a regularization parameter equal to 5% of the maximum eigenvalue of the unregularized matrix. The forward model was based upon a dipole approximation (54) and a multiple local sphere head model (55). Dipole orientation was determined using a nonlinear search for optimum SNR. Beamformer time courses for all voxels in any one region were sign flipped where necessary to account for the arbitrary polarity introduced by source orientation estimation. Following this, time courses were combined in a weighted average, where the weighting was defined by a 3D Gaussian function (full width at half maximum of 17 mm) centered on the region’s center of mass (i.e., those voxels further from the center of the AAL region were down-weighted) (56). This procedure allowed a single beamformer-derived electrophysiological time course for each AAL region.

Regional signals were filtered into the theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), low-gamma (30–70 Hz), and high-gamma (70–100 Hz) frequency bands. Within each band, the confound of signal leakage was mitigated using a pairwise orthogonalization scheme described previously (48, 57). Following leakage correction, a Hilbert transform was used to derive the analytic signal, and the absolute value of the analytic signal was computed to give the amplitude envelope of oscillations, for each time course, in all five frequency bands. Pearson correlation was then computed between the envelopes for each region pair; this was done separately for each 10-s epoch and results averaged over epochs within each participant. In this way we generated a single 78 × 78 adjacency matrix for each subject, and each frequency band, representative of connectivity between all pairs of AAL regions. These connectivity matrices were averaged over subjects within each band to generate five group level connectivity matrices, representing frequency specific amplitude envelope connectivity across the cortex. It should be pointed out that this AEC method has been used extensively in previous work (for reviews see refs. 29 and 58) and represents an intrinsic mode of electrophysiological coupling (34).

### MRI Data Acquisition.

MRI data were collected using a Phillips Achieva 7T system, equipped with a whole head volume transmit coil and a 32-channel receive head coil. A high-resolution anatomical phase-sensitive inversion recovery [PSIR; field of view (FOV) = 240 × 216 × 160 mm^{3}, 0.8-mm isotropic resolution] image was acquired and used for MEG coregistration (discussed above). MT data were obtained from z-spectra acquired using an MT-TFE sequence (46). Each point on the z-spectrum was acquired using a saturation-prepared 3D TFE sequence. The saturation consisted of a train of *n* = 20 Gaussian windowed sinc rf pulses of bandwidth BW = 200 Hz, 30 ms long, repeated every T = 60 ms (50% duty cycle), with a phase increment between each pulse, and a spoiler gradient applied at the end of the train to remove any residual transverse magnetization. Z-spectra were acquired by varying the off-resonance saturation frequency from −5 kHz to 5 kHz, including a scan at 50 kHz off-resonance for normalization. This was repeated for the three nominal B_{1sat} values of 0.33, 0.65, and 1.09 μT (in B_{1rms}). The imaging readout was a volume acquisition with a readout train of 410 gradient echoes, TE/TR/FA = 2.7 ms/5.8 ms/8°, FOV = 192 × 192 × 60 mm^{3}, 1.5-mm isotropic image resolution, low-high k-space acquisition, and a SENSE factor (RL) of 2. The 3D volume acquisition required five repetitions of this cycle. Using this 3D nonsteady-state approach, a 15-point z-spectrum (plus an additional point for normalization) was acquired in 8 min (24 min total for the three powers). The amplitude of readout pulses was modulated to avoid large variations in signal at the start of the TFE train; there were 2 ramped rf pulses before acquisition, followed by 4 ramped rf pulses at the start of the acquisition, with the remaining pulses at constant flip angle alpha = 8° (B_{1read}). Additionally, a B_{0} map (double echo method, FOV = 252 × 255 × 100 mm^{3} at 3 × 3 × 2 mm^{3} voxels) was acquired for B_{0} correction (1 min), along with a whole-head B_{1} map (dual TR method, at 20 ms and 120 ms TR, FOV 205 × 180 × 132 mm^{3} at 3.2 × 4 × 4 mm^{3} voxels, 3 min) and a T_{1} map [dual readout PSIR data, FOV = 240 × 216 × 160 mm^{3} at 0.8-mm isotropic resolution, SENSE factor (RL) of 2.2 and (FH) of 2, 10 min], giving a total scan time of 38 min.

### MRI Data Analysis and MT Mapping.

The procedure for extracting high-fidelity MT data from the CEST imaging sequence has been described in detail in a recent paper (47). All data were coregistered onto the MT space using a high contrast-to-noise ratio (hCNR) image created by averaging across z-spectra dynamics acquired at the highest B_{1} amplitude. Each z-spectra dataset was motion corrected [using MCFLIRT in FSL (59)] and registered onto the hCNR image. Similarly, B_{0}, B_{1}, and PSIR scans were registered onto the hCNR image. All z-spectra were then corrected for B_{0} fluctuations (commonly observed at ultra-high field strength). The B_{0} corrected z-spectra were then fitted to a database of simulated spectra to extract myelination (MT) maps using a four-pool model with 15 off-resonance values sampled (47). This method includes B_{1} and T_{1} correction, because it is important to obtain unbiased MT maps as described in ref. 47.

To investigate covariation of myelin within the same 78 cortical AAL regions used for MEG, MT data were first segmented to separate gray and white matter. Gray-matter masks were generated from a high-resolution anatomical image (PSIR) using SPM. These were then thresholded at a high probability value. This, in effect, meant a highly conservative gray-matter mask that minimizes any partial volume effects due to the relatively low resolution of the MT data. Gray-matter-masked MT data were then registered to the AAL atlas. A mean MT measure was calculated for each region, for each participant, creating 78 values of MT per person. Due to limited FOV and intersubject variability, 14 AAL regions were either missed or poorly characterized by the MT acquisition and so were removed from any subsequent analyses. The remaining 64 regions were averaged across subjects to give a single map, showing average MT in each region. Due to potentially confounding factors in the MRI such as scanner drift, the common mode for each individual’s MT data was calculated and regressed from that individual’s regional MT values (7). Following this, regional MT values were correlated across subjects to generate a structural matrix. In everything that follows we assume that these MT maps, and the associated structural network, are representative of cortical myeloarchitecture.

To confirm that our myelin maps were related to human physiology we assessed the relationship of myelin with handedness to test a hypothesis that right handers would exhibit greater myelin content in the left hemisphere, whereas left handers would exhibit greater myelin content in their right hemisphere. Given the right-handed bias in our cohort of 58 participants (24 reported being strongly right handed), we first tested statistically for a difference in myelination, across the subject group, between the left and right motor cortices, using a paired *t* test. Following this, we assessed correlation across subjects between MT and handedness, the latter being assessed by the Edinburgh handedness inventory. The questionnaire was scored such that strongly right-handed individuals have a score of +12, strongly left-handed individuals have a score of −12, and ambidextrous individuals fall between these limits. Thus, for any one brain region, a positive gradient between MT and handedness would indicate more myelin for right handers. Likewise, a negative gradient would indicate more myelin for left handers. We hypothesized that a negative gradient would dominate the right hemisphere whereas a positive gradient would dominate the left hemisphere. To test this statistically, first the difference in mean (across regions) gradient for each hemisphere was measured. This value was then tested against an empirical null distribution. To generate the null, the gradient difference was computed in the same way; however, the gradients for each region were defined following a randomization of the order of the handedness scores across subjects. In other words, we reasoned that, if there was no relationship between handedness and myelination, the handedness scores could be switched randomly between subjects with no effect on the slope difference detected. Hence, the real gradient difference would fall within the null distribution.

### The Relationship Between MT and MEG in Individual Frequency Bands.

To assess a relationship between the myelin structural network and the MEG functional networks we probed correlation between network matrices. All matrices were first vectorized and redundancy was removed, so that the leading diagonal (which represents only within-region correlation) did not contribute to the measurement. Following this, linear correlation was measured between the structural network matrix and all five MEG functional connectivity matrices (representing each of the five frequency bands). This yielded five values of correlation (represented by r^{2}). In testing whether these correlation values are statistically significant, a significant problem arises. Specifically, given the logical ordering of brain regions in the AAL atlas and the finite spatial resolution of MEG, functional connectivity matrices exhibit inherent smoothness, meaning that separate connections between pairs of brain regions are not independent. This means that simply assuming that the number of degrees of freedom in the matrix is equivalent to the number of measured connections would vastly inflate the prescribed significance of the correlation. To eliminate this problem we used a permutation test. A set of “pseudo-MEG matrices” was generated, each having spatial properties similar to those of the real matrices but, crucially, not based on genuine data. This procedure has been described recently by Tewarie et al. (31); briefly, to obtain the pseudomatrices we first performed an eigenvalue decomposition of the real MEG-derived matrices. Each eigenvector was then randomized using a phase-based technique (60). Reconstruction of the matrix postrandomization yielded a pseudomatrix similar in mathematical structure to the genuine adjacency matrices but not reflecting genuine MEG-derived functional connectivity (Fig. S4). Correlation between the structural matrix and 10,000 iterations of the pseudomatrix (with a different phase distribution on each iteration) was measured and an empirical null distribution derived. Comparison of the r^{2} value from real MEG connectivity matrices with the empirical null from the pseudomatrices allowed computation of a *P* value. Significant correlation was assigned if the *P* value was less than 0.01 (i.e., a threshold of 0.05 corrected for multiple comparisons across bands).

### The Relationship Between MT and MEG in Linear and Nonlinear Combinations.

It is likely that, if our hypothesis that functional connectivity predicts myeloarchitecture is true, then the relationship would not be confined to any single frequency band. Rather, it is likely that the structural network develops based on a combination of electrophysiological activity across all bands. For this reason we used a recently developed mathematical framework (31) to assess whether linear or nonlinear combinations of frequency bands could predict the structural network better than any single MEG band. A three step approach was adopted:

• First, a single-frequency model was used to assess the predictive value of each individual MEG frequency band (as described above).

• Second, all five MEG frequency bands in linear combination were used to predict the structural network.

• Finally, a multivariate truncated Taylor series was adopted allowing for both linear terms (i.e., the regional MEG matrices in Fig. 2) and squared terms (i.e., the square of each of the five MEG matrices in Fig. 2), 10 terms in all.

For all three steps above, a goodness-of-fit measure (r^{2}) was used to assess the efficacy of the model. Importantly, the squared terms have specific meaning because they represent covariance between the neuronal connectivity profiles of separate regions. In other words, for any nonlinear term, a matrix element, say [1, 2], will be high if the connectivity profile of region 1 to the rest of the brain overlaps with the connectivity profile of region 2 to the rest of the brain. This physical interpretation helps explain why their inclusion might improve the prediction of structural matrices; if two regions share electrophysiological interactions with similar areas, then it is likely that their myeloarchitecture may evolve together.

Our hypothesis was that as the model got more complex (i.e., moving from single frequency band to first linear and then nonlinear integration) we would see a better fit to the structural covariance matrix, and that this would be due to the fact that myeloarchitecture supports networks in all frequency bands. However, a confound here is that, when introducing more complexity into a model, the goodness of fit for that model is likely to increase, simply because one tries to fit data using more terms. We therefore tested whether the improvements observed were merely the result of increased model complexity. To do this, the rate of improvement of the fit with increasing model order was calculated and tested statistically. The gradient of the slope representing improvement was first measured using real data. Second, a null distribution of model improvement was constructed. Here, five pseudomatrices generated based upon the five MEG matrices in Fig. 2 were used. These five sham matrices were combined (both linearly and nonlinearly) and a fit to the MT data derived, alongside an associated r^{2}. A gradient, representing improvement in fit to the MT data with increasing model complexity, was then derived. This procedure was repeated 10,000 times (different realizations of the pseudomatrices on each iteration). The 10,000 values of gradient then represented an empirical null distribution; this was compared with the value from real data.

Finally, we performed a cross-validation to test whether the prediction of the structural network based on multiple frequency bands was robust. If the estimated parameters contain meaningful information, then applying the same set of parameters on an independent dataset would result in a significantly better prediction of the structural network than an equivalent prediction based on pseudo MEG networks. To this end, the cohort of 58 subjects was split into two equal groups, A and B (29 subjects in each). For each group a structural network and average MEG networks were estimated. First, the structural network in group A was predicted based on linear and nonlinear combinations of the MEG networks from group A, resulting in an r^{2} value. Second, parameters from this prediction were used to map the MEG networks to the structural network in group B to obtain a “cross-validation” r^{2}. Third, the same parameters from group A were used to map pseudo MEG networks obtained from group B to the structural network of group B*,* which yielded in a “sham” r^{2} value. This procedure was repeated 10,000 times, with the 58 subjects split differently on each iteration. This resulted in three r^{2} distributions (original, cross-validation, and sham). The parameters of the initial prediction were said to contain meaningful information if the median of the distribution of cross-validation r^{2} values could be found outside the 95% area under the curve of the sham r^{2} distribution.

## Acknowledgments

This work was funded by Medical Research Council (MRC) New Investigator Research Grant MR/M006301/1, MRC Partnership Grant MR/K005464/1, and MRC Doctoral Training Grant MR/K501086/1.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: matthew.brookes{at}nottingham.ac.uk.

Author contributions: D.K.J., K.D.S., P.G.M., P.A.G., and M.J.B. designed research; B.A.E.H., P.K.T., O.E.M., and N.G. performed research; P.K.T., O.E.M., N.G., P.A.G., and M.J.B. contributed new reagents/analytic tools; B.A.E.H., P.K.T., O.E.M., N.G., and M.J.B. analyzed data; and B.A.E.H. and M.J.B. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: MEG and MT adjacency matrices have been deposited at www.figshare.com and can be found with DOI: 10.6084/m9.figshare.4056441.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1608587113/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵.
- Beckmann CF,
- DeLuca M,
- Devlin JT,
- Smith SM

- ↵.
- Meier J, et al.

- ↵
- ↵.
- Andrews TJ,
- Halpern SD,
- Purves D

- ↵.
- White LE, et al.

- ↵.
- He Y,
- Chen ZJ,
- Evans AC

- ↵
- ↵.
- Mechelli A,
- Friston KJ,
- Frackowiak RS,
- Price CJ

- ↵
- ↵
- ↵
- ↵.
- Glasser MF,
- Van Essen DC

- ↵
- ↵.
- Sereno MI,
- Lutti A,
- Weiskopf N,
- Dick F

- ↵
- ↵.
- Dick F, et al.

- ↵.
- Glasser MF, et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Gibson EM, et al.

- ↵
- ↵.
- Cohen D

- ↵.
- Brookes MJ, et al.

- ↵
- ↵.
- O’Neill GC,
- Barratt EL,
- Hunt BAE,
- Tewarie PK,
- Brookes MJ

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Driver ID, et al.

- ↵
- ↵
- ↵
- ↵.
- Doria V, et al.

- ↵.
- Hebb DO

- ↵
- ↵
- ↵
- ↵.
- Iwabuchi SJ,
- Liddle PF,
- Palaniyappan L

- ↵.
- Mougin O,
- Clemence M,
- Peters A,
- Pitiot A,
- Gowland P

- ↵.
- Geades N, et al.

*Magn Reson Med*10.1002/mrm.26459. - ↵
- ↵.
- Filippini N, et al.

- ↵
- ↵
- ↵.
- Yoshimoto T,
- Kotani M,
- Kuriki S,
- Karibe H,
- Nakasato N

- Robinson S,
- Vrba J

- ↵
- ↵
- ↵
- ↵.
- Brookes MJ, et al.

- ↵
- ↵
- ↵
- ↵