Structural and functional, empirical and modeled connectivity in the cerebral cortex of the rat.

Connectomics data from animal models provide an invaluable opportunity to reveal the complex interplay between structure and function in the mammalian brain. In this work, we investigate the relationship between structural and functional connectivity in the rat brain cortex using a directed anatomical network generated from a carefully curated meta-analysis of published tracing data, along with resting-state functional MRI data obtained from a group of 14 anesthetized Wistar rats. We found a high correspondence between the strength of functional connections, measured as blood oxygen level dependent (BOLD) signal correlations between cortical regions, and the weight of the corresponding anatomical links in the connectome graph (maximum Spearman rank-order correlation ρ=0.48). At the network-level, regions belonging to the same functionally defined community tend to form more mutual weighted connections between each other compared to regions located in different communities. We further found that functional communities in resting-state networks are enriched in densely connected anatomical motifs. Importantly, these higher-order structural subgraphs cannot be explained by lower-order topological properties, suggesting that dense structural patterns support functional associations in the resting brain. Simulations of brain-wide resting-state activity based on neural mass models implemented on the empirical rat anatomical connectome demonstrated high correlation between the simulated and the measured functional connectivity (maximum Pearson correlation ρ=0.53), further suggesting that the topology of structural connections plays an important role in shaping functional cortical networks.

1. Introduction

The study of brain connectivity from a network perspective (Newman, 2003; Strogatz, 2001) has become a promising framework to understand how action, perception, and cognition emerge from a dense ensemble of neural elements (Park and Friston, 2013). Leveraging advances in brain imaging and network science (Sporns, 2013; Sporns et al., 2005), recent approaches have focused on the topology and dynamics of large-scale projections linking anatomically distinct and functionally specialized brain regions (Bullmore and Bassett, 2011). The structure of these large-scale networks is thought to shape and constrain inter-regional interactions and computations.

Interactions between neuronal populations spanning brain-wide networks can be described from three different, but related, perspectives (Friston, 2011). Briefly, anatomical or structural connectivity (SC) refers to patterns of synaptic connections linking brain areas. Functional connectivity (FC) refers to statistical interdependence between activity time series recorded in different, often spatially remote, areas. Finally, effective connectivity (EC) refers to the influence that one neural system in one area exerts over another. While the interplay of these three modes of brain connectivity is not completely understood, some progress has been made combining anatomy with both resting-state and task-based functional connectivity (Hermundstad et al., 2013). Spontaneous or intrinsic neural activity (Cole et al., 2010; Fox and Raichle, 2007), as measured using fluctuations in the blood oxygen level dependent (BOLD) signals of resting-state functional MRI (rs-fMRI), has proven to be a useful technique for examining the extent to which structural patterns shape functional interactions between neural assemblies (Honey et al., 2010). Previous studies have shown that the presence of strong SC, as measured with diffusion-weighted MRI, between two areas increases the probability and strength of corresponding FC. Nevertheless, it has also been reported that strong FC may exist between areas with no (direct) anatomical connections, (Bowman et al., 2012; Damoiseaux and Greicius, 2009; Skudlarski et al., 2008), suggesting that indirect signaling and emergent dynamic processes make an additional strong contribution. For example, a study of SC and FC in macaque cortex supported the idea that functional interactions are strongly influenced by network-wide effects (Adachi et al., 2012).

Approaches simulating spontaneous cortical dynamics, in both humans and animal models, can yield invaluable insight into structure-function relationships by assessing the capacity of simulated dynamic patterns based on the SC scaffold to predict empirically measured BOLD signal correlations (Deco et al., 2011; Honey et al., 2007; Nakagawa et al., 2013). Connectomics data from animal models based on tract-tracing procedures allows in-depth characterization of SC networks. In contrast to MRI-based tractography, which provides coarse-grained undirected SC matrices, histological tracing technology yields highly resolved and directed connectivity information, hence providing important additional information for modelling cortical dynamics (van den Heuvel et al., 2016a). For instance, recent work relating the structural connectome of the mouse brain and the intrinsic BOLD signal dynamics within individual brain regions have shown the importance of considering both the weight and directionality of structural connections (Sethi et al., 2017). In addition, the mapping of functional connectivity networks in rodents provides an invaluable tool to understand neurological and psychiatric disorders from a more mechanistic way for translational research. (Gozzi and Schwarz, 2016; Pan et al., 2015). These experimental possibilities together with theoretical developments in network analysis are extending systems neuroscience from unimodal investigations of brain connectivity to a network-level understanding of structure-function interactions (Adachi et al., 2012; Diez et al., 2015; Goñi et al., 2014; Hsu et al., 2016; Mišić et al., 2016; Skudlarski et al., 2016; Stafford et al., 2014; Wang et al., 2015; Wirsich et al., 2016).

In this work, we examine the relationship between SC and FC in the rat cortical network. Using a detailed cortical SC matrix obtained from a carefully curated meta-analysis of published histological tracing data in rats (Bota et al., 2015), we first compare structural connections in the rat cortex with their corresponding spontaneous correlations extracted empirically from rs-fMRI data collected in a group of 14 Wistar rats. We then show the results of this comparison taking into account network-level effects by relating structural properties of brain connectivity to the functional modularity of rs-fMRI networks. Specifically, we study link reciprocity in both intra- and inter-modular connections as well as the structural motif frequency spectrum within functionally defined modules. Finally, we carry out computational simulations of neural mass models implemented on the empirical SC to generate a simulated FC matrix and compare it with the empirically measured FC. Overall, our results provide evidence on rs-fMRI BOLD signal correlations being constrained and shaped by the underlying structural connectivity patterns.

2. Materials and methods

2.1. Animals and MRI acquisition protocol

Experiments were carried out in a horizontal 7 T scanner with a 30 cm diameter bore (Biospec 70/30v, Bruker Medical, Ettlingen, Germany). The system had a 675 mT/m actively shielded gradient coil (Bruker, BGA 12-S) of 11.4 cm inner diameter. A 1H rat brain receive-only phase array coil with integrated combiner and preamplifier, no tune/no match, in combination with the actively detuned transmit-only resonator (BrukerBioSpin MRI GmbH, Germany) was employed. Data were acquired and processed with a Hewlett-Packard console running Paravision 5.1 software (Bruker Medical GmbH, Ettlingen, Germany) operating on a Linux platform.

For the rs-fMRI experiments, 14 Wistar rats were anesthetized with urethane (1.2 g/Kg). Anesthetized animals were placed in a custom-made animal holder with adjustable bite and ear bars, and positioned on the magnet bed. The animals were constantly supplied with 0.8 L/m O2 with a face mask and temperature was kept between 36.5 and 37.5 °C through a water heat-pad. The temperature, heart rate, SpO2, and breathing rate were monitored throughout the session (MouseOx, Starr Life Sciences, Oakmont, US).

T2-weighted anatomical images were collected using a rapid acquisition relaxation enhanced sequence (RARE), applying the following parameters: field of view (FOV) = 40 × 40 mm; 15 slice thickness = 1 mm; matrix size = 128 × 128; effective echo time (TEeff) = 56 ms, repetition time (TR) = 2 s, and a RARE factor of 8. The B0 field distribution in a large voxel (40 × 40 × 40 mm3) containing the whole head was acquired (FieldMap). Briefly, the brain was localized with T2-weighted RARE sequence, and first- and second-order shims adjusted with MAPSHIM application in a sufficiently large voxel containing the whole brain. Functional MRI acquisition was performed using a GE-EPI sequence in 30 coronal slices applying the following parameters: FOV = 25 × 25 mm; slice thickness = 0.5 mm; matrix size = 50 × 50; segments = 1; FA = 60°; echo time (TE) = 15 ms; repetition time (TR) = 2000 ms (300 samples per run, 10 min), rendering an isotropic voxel of 0.5 × 0.5 × 0.5 mm3. Between one and three runs were acquired from each animal. T2-weighted anatomical images with exactly the same geometry were collected using a RARE sequence using the following parameters: FOV = 25 × 25 mm; 30 slices; slice thickness = 0.5 mm; matrix size = 200 × 200; TEeff = 56 ms; TR = 2 s; RARE factor = 8.

2.2. Preprocessing of MRI data

Data preprocessing within runs was carried out using FSLv5.0 (FMRIB Software Library, http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) (Jenkinson et al., 2012) and MATLAB 2014a (The MathWorks, Inc., Natick, MA, United States, https://www.mathworks.com/). Once images were converted to NIfTI (Neuroimaging Informatics Technology Initiative, http://nifti.nimh.nih.gov/) data format, the original voxel size, (x,y,z), was scaled up by a factor of 10. This step is very common when analyzing rodent data to accurately apply the same algorithms (largely those involving spatial transformations) as in human analyses (Kalthoff et al., 2013; Pan et al., 2015).

The very first volume of fMRI data was used as reference across runs of the same subject for head motion correction, brain segmentation and co-registration. As suggested by (Kalthoff et al., 2011), head motion correction was applied to each individual slice and restricted to (coronal) in-plane translations (x, y) and rotation (z) to reduce signal fluctuations related to respiration in anesthetized rats. After applying motion correction and brain segmentation (Smith, 2002), global intensity normalization was set to 1000 and spike detection was performed through using DVARS measure (Power et al., 2012). Note that DVARS is highly dependent on the particular dataset (Power et al., 2014); we therefore did not select an absolute threshold and instead considered as outliers those temporal points exceeding the 75th percentile + 1.5*IQR (interquartile range). None of the runs were discarded since the number of spikes was below 30 (out of 300 samples) in all cases (15 ± 4.6 spikes per run, mean ± SD), hence ensuring a minimum length of 9 min to estimate functional interactions. By using a nuisance regression model, each voxel was corrected for: (1) the three rigid body parameters (translation in x and y, and rotation in z) previously computed and their derivatives (backward difference); (2) a single regressor per spike with a “b0, f0” window (Satterthwaite et al., 2013) (i.e. neither preceding nor following samples were used); (3) the global signal and its derivative (backward difference); and (4) two regressors modelling the mean and a linear trend. A Band-pass filtering (nonlinear high-pass filter of = 50 , and Gaussian linear low-pass filter of = 2 ) was applied to retain those frequencies ranging from 0.01 and at around 0.1 Hz (the resulting spectrum of this filter is shown in Supplementary Fig. 1). Spatial smoothing was not applied to avoid introducing spurious high correlations between a node and its neighbors (Fornito et al., 2010). Finally, filtered rs-fMRI data were co-registered to the brain-extracted T2-weigthed image using a rigid body transformation and then normalized to a rat template described elsewhere (Schwarz et al., 2006). Normalization was carried out through using an affine deformation (Jenkinson et al., 2002; Jenkinson and Smith, 2001) and the template was resampled to match the original resolution of functional images (i.e., 0.5 × 0.5 × 0.5 mm3).

2.3. Rat connectome and definition of brain areas

We made use of a directed anatomical network coming from a systematic curation of the primary neuroanatomical literature for the rat (for specific details as to annotation and collation methodology, see ref. (Bota et al., 2015)). In particular, the dataset was originally composed of 73 cortical gray-matter regions that were defined according to the Swanson-04 hierarchical nomenclature for the rat central nervous system (Swanson, 2004). Connections in the rat cortical association macro-connections (RCAMs) matrix are not yet defined bilaterally and are encoded by means of eight different ordinally arranged weight categories: not present, very weak, weak, weak/moderate, moderate, moderate/strong, strong, and very strong.

A critical point in neuroimaging research is the use of a common framework for localization of brain structures, which allows comparison among results coming from different studies. Nomenclature of the rat connectome project is available in the Swanson space (SwS) (Swanson, 2004) and anatomical connectivity between regions is directly given in matrix form, that is, no Analyze or NIfTI images are provided. On the other hand, to the best of our knowledge, the most complete atlas available in Analyze format (easily converted to NIfTI format) is that developed by Schwarz et al. (2006), where brain areas are given in the Paxinos&Watson space (PWS) (Paxinos and Watson, 1998). Therefore, the first task before carrying out network analysis was to establish a correspondence of brain regions between both spaces with the aim of aligning functional correlation matrices with the anatomical network. To this end, we carefully inspected both atlases as follows: (1) a particular reference brain region from SwS was localized in its corresponding coronal plane; (2) we sought the coronal slice/s in PWS containing the previous region, and (3) by taking into account the spatial distribution of surrounding areas as well as the relative anterior and posterior planes, we identified which region or group of regions from PWS matched with the reference region.

Since not all of the original structures contained in the RCAMs matrix were available in the NIfTI atlas as a single mask, some of them were grouped to cover the whole brain cortex (Supplementary Fig. 2 schematizes the process carried out for computing these new connections). Thus, we created a new pattern of cortical connections for the secondary visual area (VISs), which was composed of the original anterior latero-lateral area (VISlla), anterolateral area (VISal), rostrolateral area (VISrl), intermediolateral area (VISli), laterolateral area (VISll), mediolateral area (VISlm), posterolateral area (VISpl), and anteromedial area (VISam). The tenia tecta cortical region (TT) was composed of the dorsal and ventral parts (TTd and TTv, respectively). Interactions of the retrosplenial area ventral part (RSPv) with other areas were obtained from the original ventral part (RSPv), ventral part zone a (RSPv-a), and ventral part zone b/c (RSPv). The hippocampal region was partitioned into dorsal and ventral parts as follows. The dorsal part of the hippocampal region (HIPd) was composed of CA1 dorsal part (CA1d), CA3 dorsal part (CA3d), CA2 dorsal part (CA2d), dentate gyrus dorsal part (DGd), and induseum griseum (IG). Whereas the ventral part of the hippocampal region (HIPv) was composed of CA1 ventral part (CA1v), CA3 ventral part (CA3v), CA2 ventral part (CA2v), and dentate gyrus ventral part (DGv). Note that the current RCAMs matrix does not distinguish between dorsal and ventral parts for CA2, CA3 and DG. To overcome this limitation, and with the aim of including the whole hippocampal region in our analysis, structural patterns contained in the rat connectome for CA2, CA3 and DG were assumed to be equal for both ventral and dorsal parts. Categorical weight for the same region in different parts (e.g., between CA2d and CA2v) were assumed to be very strong. A new pattern of structural connections was also created for the cortical amygdalar nucleus (COA), which was formed from the anterior part (COAa), posterolateral part (COApl), and posteromedial part (COApm). The secondary auditory areas (AUDs) was composed of the ventral and dorsal auditory areas (AUDd and AUDv, respectively). The medio-ventral part of the orbital area (ORBmv) included the ventral and medial parts (ORBv and ORBm, respectively). The presubiculum (PRE) and parasubiculum (PAR) were grouped to form the PreParaSubiculum region (PREPAR). The basolateral amygdalar nucleus (BLA) was composed of the anterior and posterior parts (BLAa and BLAp, respectively). Finally, structural interactions between the basomedial amygdalar nucleus (BMA) and the rest of nodes were obtained from the anterior and posterior parts (BMAa and BMAp, respectively).

Definition of nodes is a very important step for brain network analysis. In this work, two different brain parcellations were used, which were based on the information available in the NIfTI atlas. Table 1 presents the 50 brain regions from SwS and their counterpart regions from PWS used for simulations (see subsection 2.8). Robust BOLD signals were not consistently available for the above 50 ROIs in our dataset due to variations in the field of view across rats. Thus, a second parcellation was used in which a reduced set of 32 ROIs (Table 1) were selected for the comparison of empirical and simulated FC, as well as to investigate intra- and between-module structural patterns.

Table 1

Cortical brain regions used as network nodes in this work. All regions listed were used for simulations, whereas those marked with an asterisk were removed for ensuring consistent extraction of BOLD signals across rats.

Swanson structure	Abbreviation	Paxinos&Watson structure	Abbreviation
Primary somatomotor area	MOp	Primary motor cortex	M1
Secondary somatomotor areas	MOs	Secondary motor cortex	M2
Primary somatosensory area	SSp	- Primary somatosensory cortex	S1
		- Primary somatosensory cortex, barrel field	S1BF
		- Primary somatosensory cortex, dysgranural region	S1DZ
		- Primary somatosensory cortex, forelimb region	S1FL
		- Primary somatosensory cortex, hindlimb region	S1HL
		- Primary somatosensory cortex, jaw region	S1J
		- Primary somatosensory cortex, jaw region, oral surface	S1JO
		- Primary somatosensory cortex, trunk region	S1Tr
		- Primary somatosensory cortex, upper lip region	S1ULp
Supplemental somatonsensory area	SSs	Secondary somatosensory cortex	S2
Visceral area	VISC	Granular insular cortex	GI
Infralimbic area	ILA	Infralimbic cortex	IL
Gustatory area	GU	Dysgranular insular cortex	DI
Main olfactory bulb*	MOB	- Glomerular layer of the olfactory bulb	Gl
		- Granular cell layer of the olfactory bulb	GrO
		- Internal plexiform layer of the olfactory bulb	IPl
		- Mitral cell layer of the olfactory bulb	Mi
		- External plexiform layer of the olfactory bulb	EPl
Accessory olfactory bulb*	AOB	- Accessory olfactory bulb	AOB
		- External plexiform layer of the accessory olfactory bulb	EPlA
		- Glomerular layer of the accessory olfactory bulb	GlA
		- Granule cell layer of the accessory olfactory bulb	GrA
		- Mitral cell layer of the accessory olfactory bulb	MiA
Anterior olfactory nucleus	AON	- Anterior olfactory nucleus, dorsal part	AOD
		- Anterior olfactory nucleus, external part	AOE
		- Anterior olfactory nucleus, lateral part	AOL
		- Anterior olfactory nucleus, medial part	AOM
		- Anterior olfactory nucleus, posterior part	AOP
Tenita tecta	TT	- Tenia tecta, layer 1	TT1
		- Tenia tecta, layer 2	TT2
		- Tenia tecta, layer 3	TT3
Piriform area	PIR	- Piriform layer	Pir
		- Region external to piriform layer	Pir/ext
		- Region internal to piriform layer	Pir/int
		- Piriform cortex	PirCtx
		- Cortex amygdala transition zone	CxA
Postpiriform transition area*	TR	Amygdalopiriform transition area	APir
Nucleus of lateral olfactory tract*	NLOT	- Nucleus of the lateral olfactory tract	LOT
		- Nucleus of the lateral olfactory tract, layer 1	LOT1
		- Nucleus of the lateral olfactory tract, layer 2	LOT2
Cortical amygdalar nucleus*	COA	- Anterior cortical amygdaloid nucleus	ACo
		- Posterolateral cortical amygdaloid nucleus	PLCo
		- Posteromedial cortical amygdaloid nucleus	PMCo
Primary auditory area	AUDp	Primary auditory cortex	Au1
Secondary auditory areas	AUDs	- Secondary auditory cortex, dorsal area	AuD
		- Secondary auditory cortex, ventral area	AuV
Secondary visual area	VISs	- Secondary visual cortex, lateral area	V2L
		- Secondary visual cortex, mediomedial area	V2MM
		- Secondary visual cortex, mediolateral area	V2ML
Primary visual area	VISp	- Primary visual cortex, binocular area	V1B
		- Primary visual cortex, monocular area	V1M
Anterior cingulate area, dorsal part	ACAd	Cingulate cortex, area 1	Cg1
Anterior cingulate area, ventral part	ACAv	Cingulate cortex, area 2	Cg2
Prelimbic area	PL	Prelimbic cortex	PrL
Orbital area, lateral part*	ORBl	Dorsolateral orbital cortex	DLO
Orbital area, medio-ventral part	ORBmv	- Medial orbital cortex	MO
		- Ventral orbital cortex	VO
Orbital area, ventrolateral part	ORBvl	Lateral orbital cortex	LO
Agranular insular area, dorsal part	AId	Agranular insular cortex, dorsal part	AID
Agranular insular area, ventral part	AIv	Agranular insular cortex, ventral part	AIV
Agranular insular area, posterior part	AIp	Agranular insular cortex, posterior part	AIP
Retrosplenial area, dorsal part	RSPd	Retrosplenial granular B cortex	RSGb
Retrosplenial area, lateral agranular part	RSPagl	Retrosplenial agranular cortex	RSA
Retrosplenial area, ventral part	RSPv	Retrosplenial granular A cortex	RSGa
Posterior parietal association areas	PTLp	Parietal association cortex	PtA
Temporal association areas*	TEa	Temporal association cortex	TeA
Ectorhinal area*	ECT	Ectorhinal cortex	Ect
Perirhinal area*	PERI	Perirhinal cortex	PRh
Entorhinal area, lateral part*	ENTl	Lateral entorhinal cortex	LEnt
Entorhinal area, medial part, dorsal zone*	ENTm	Medial entorhinal cortex	MEnt
PreParaSubiculum*	PREPAR	- Presubiculum	PrS
		- Parasubiculum	PaS
Postsubiculum	POST	Postsubiculum	Post
Subiculum, dorsal part	SUBd	Subiculum, dorsal part	DS
Subiculum, ventral part*	SUBv	Subiculum, ventral part	VS
Hippocampal region, dorsal part	HIPd	- Hippocampus posterior, dorsal part	HCpd
		- Hippocampus fronto dorsal	HCfd
		- Dentate gyrus, dorsal part	DGd
		- Field CA3 of hippocampus, dorsal part	CA3d
		- Indusium griseum	IG
Hippocampal region, ventral part	HIPv	- Hippocampus posterior, ventral part	HCpv
		- Field CA3 of hippocampus, ventral part	CA3v
		- Dentate gyrus, ventral part	DGv
Claustrum	CLA	Claustrum	Cl
Endopiriform nucleus, dorsal part	EPd	Dorsal endopiriform nucleus	DEn
	EPv	Ventral endopiriform nucleus	VEn
Endopiriform nucleus, ventral part*
Lateral amygdalar nucleus*	LA	- Lateral amygdaloid nucleus, dorsolateral part	LaDL
		- Lateral amygdaloid nucleus, ventrolateral part	LaVL
		- Lateral amygdaloid nucleus, ventromedial part	LaVM
Basolateral amygdalar nucleus*	BLA	- Basolateral amygdaloid nucleus, anterior part	BLA
		- Basolateral amygdaloid nucleus, posterior part	BLP
		- Basolateral amygdaloid nucleus, ventral part	BLV
Basomedial amygdalar nucleus*	BMA	- Basomedial amygdaloid nucleus, anterior part	BMA
		- Basomedial amygdaloid nucleus, posterior part	BMP
Posterior amygdalar nucleus*	PA	- Amygdalohippocampal area, anterolateral part	AHiAL
		- Amygdalohippocampal area, posteromedial part	AHiPM

2.4. Construction of functional networks

The mean time courses within each of the 32 regions of interest (ROIs) were extracted and converted to z-scores (i.e. fMRI BOLD time series were centered and scaled to have zero mean and unit variance). Since the collation process of cortical associations in the current RCAMs ((Bota et al., 2015)) do not differentiate between hemispheres, voxels within ROIs were also combined across hemispheres (Supplementary Figs. 3, 4 and 5 present seed-based correlations analysis in three different brain regions uncovering bilateral networks in the dataset). That is, the time course from a given region, say the primary somatomotor area, was obtained as the average between signals from right and left primary somatomotor areas. Importantly, to create a homogeneous parcellation by ensuring that all regions had the same size, we adopted the procedure carried out in (Alexander-Bloch et al., 2012), where ROIs were eroded until all of them were exactly the same volume, 1 mm3 here (0.5 mm3 in each hemisphere). Before computing functional interactions, time samples marked during preprocessing were “scrubbed”. Functional networks were then estimated for each single run by computing the Pearson correlation coefficient. Next, raw correlations were converted to Fisher’s z-values and connectivity matrices of each run were averaged within subject. Finally, FC matrices across the 14 subjects were averaged and a group network was obtained and then used for comparative analysis with the RCAM matrix and simulated FC matrices.

2.5. Modularity

Community structure detection involves the partition of a network into “modules” or “clusters” wherein nodes are highly connected to each other and only sparsely connected with nodes of different modules (Newman and Girvan, 2004). In neuroscience applications, community detection allows the grouping of neural elements (e.g. brain regions) in both anatomical and functional networks into distinct modules (Sporns and Betzel, 2016). Despite the fact that community detection is conceptually straightforward, its application is methodologically challenging as indicated by the large number of algorithms and approaches dealing with it (Fortunato, 2010; Fortunato and Hric, 2016). Whereas it is also possible to perform graph partition with overlapping modules, we only consider in this work the detection of non-overlapping modules.

One popular approach for community detection attempts to maximize a quality function commonly known as the modularity function. Let W be the functional matrix, then the problem of community detection can be formalized as:

In the previous expression, represents the actual weight of the connection between node i and j, whereas refers to a specified null network model so that nodes within communities are internally more connected than expected by chance. The term ∈ [1, …, ] stands for the assignment of the node i to the k module, and the Kronecker function (, ) is equal to unity when node i and j belong to the same community, = , and zero otherwise. Thus, only nodes belonging to the same community contribute to maximize . The modularity function is typically modified by including a constant factor before summation with the aim of setting the maximum of the quality function to one. This last step has no influence on the final partition.

The precise value of the null model relies strongly on the nature of the network being analyzed. Traditionally, the most popular null model is the so-called configuration model, represented by the following expression: with and . However, it has recently been shown that approaches using Eq. (1) can lead to biased results when correlation matrices are considered (Bazzi et al., 2016; MacMahon and Garlaschelli, 2015). We therefore made use of a different null model known as the constant Potts model (Traag et al., 2011). This approach is named constant because the weight matrix in expression (1) is compared to a tunable parameter gamma :

By varying gamma, one is able to reveal community structure at different scales (given by the number of communities and their size), hence mitigating the problem known as the resolution limit (Fortunato and Barthélemy, 2007).

In this work, we used 109 different gamma values ranging linearly from −0.02 to 0.25 in increments of 0.0025 and, for each of them, we ran the Louvain algorithm 10,000 times (Blondel et al., 2008). After optimizing , we computed the mutual similarity (as quantified by using the z-score of the rand index, or z-Rand) over all pairs of partitions within a given setting of (Traud et al., 2011). We then obtained a consensus partition for each (Bassett et al., 2013; Lancichinetti and Fortunato, 2012).

2.6. Reciprocity

Link reciprocity is a network measure that allow assessing the tendency of node pairs to form mutual connections between each other, and whose definition has been extended for weighted networks (Squartini et al., 2013). Briefly, a real-world network can be decomposed into a symmetric (reciprocated) part and an asymmetric (non reciprocated) part. The weighted reciprocity, r, is then computed as the ratio between the total reciprocity weight and the total weight of the network. Afterwards, this quantity is scaled relative to the average weighted reciprocity derived from a random null model, (), as follows:

Thus, ρ indicates the tendency of the network to reciprocate (ρ > 0) or to avoid reciprocation (ρ < 0). In this work, ρ was computed by generating 5000 random networks preserving the in- and out-degree sequences as well as the total strength of the real network. As our main interest was to compare weighted intra-modular (a)symmetries compared to those taking place between modules, the measure ρ was computed separately for each type of links. Put simply, by imposing the partition of functional modules on the weighted anatomical matrix, we could study whether modularity obtained from rs-fMRI data showed differences in link reciprocity between intra- and inter-modular projections.

2.7. Network motifs

The dynamic behavior of a complex system relies critically on the underlying interconnections and how nodes are linked to form specific subgraphs or network motifs (Milo et al., 2002). The existence of these “building blocks” can be recognized in many directed real-world networks such as biological and technological networks, as compared with randomized networks. Motif analysis has been applied in a variety of brain networks to identify motifs that are significantly increased in frequency over various null models (Sporns and Kötter, 2004; van den Heuvel et al., 2012; Varshney et al., 2011). In contrast to link reciprocity, which is a measured quantifying second-order topological properties, networks motifs examine higher-order connectivity patterns.

The total number of different subgraphs or network motifs depends on the number of nodes being considered. For instance, for structural motifs of size = 3 (third-order topological properties) there are 13 different possible alternatives or classes through which three nodes can be linked (Table 2). Previous studies have addressed the problem of modularity by clustering subgraphs or higher-order connectivity patterns (Arenas et al., 2008; Benson et al., 2016), rather than considering dyadic interactions. Interestingly, this framework, wherein community structure detection and network motif approaches are related, can reveal new insights into the topological organization of complex networks. In this work, we were interested in the higher-order structural organization within functionally defined modules. Based on the so-called motif modularity (Arenas et al., 2008), we calculated the ratio of the number of occurrences (summarized in the motif frequency spectrum) for a given motif class restricted to within communities, , to the motif frequency spectrum computed over the whole network for the same class, :

Table 2

The 13 different structural motif classes which can be obtained for motifs of size = 3, along with the number of functional instances that each can generate.

Structural motif class	Number of functional motifs	Structural motif class	Number of functional motifs
	1		10
	1		9
	1		10
	3		10
	4		24
	3		54
	4

We calculated the normalized motif frequency spectrum for structural motifs of size = 3. In this way, we could study whether modules obtained from functional data are associated with a specific pattern of intra-modular anatomical connections. Binary motif analysis was performed for three connectivity levels of SC: (1) keeping all connections, (2) removing very weak links, and (3) removing both very weak and weak/moderate links.

The observed value was statistically compared against an ensemble of random networks (5000 samples) preserving the number of incoming edges, outgoing edges and mutual edges of each node of the actual network (Milo et al., 2002). Generating null models preserving mutual edges allows us to ascertain that higher-order patterns does not simply emerge as a result of lower-order properties (e.g., link reciprocity), but the abundance of networks motifs inside communities is a feature of real-world networks. Thus, we calculated for each reference network, and we eventually obtained a null distribution that allowed us to statistically evaluate the fraction of every motif class within communities using the − :

2.8. Simulations with The Virtual Brain

Finally, to investigate further the degree to which rat structural connectivity can predict functional connectivity, simulations of rat spontaneous neural dynamics were performed using the open-source neuroinformatics platform The Virtual Brain (TVB). TVB, amongst other features, allows for the modelling of the dynamics of large-scale brain networks by giving users access to a library of Python-based models, which range from complex models based on nonlinear differential equations with dozens of parameters to simpler linear or oscillatory models. TVB simulations take structural data about a brain - structural connectivity, spatial coordinates of ROIs, and tract length distances between ROIs - and output simulated neural activity time series for each ROI, using a given model type and set of simulation parameter values (Sanz-Leon et al., 2015). Raw simulated activity time series were passed through a Balloon-Windkessel model to generate synthetic BOLD responses; the period of the BOLD monitor was set to 2000 ms, reflecting the TR of fMRI data (Friston et al., 2000; Sanz Leon et al., 2013). The BOLD hemodynamic response function was created using a Volterra kernel with the default parameters supplied by TVB. Numerical integration during the simulations was performed with a stochastic variant of Heun’s method with a step size of 0.1 ms, with a noise term drawn randomly from a Gaussian distribution added. The standard deviation of the Gaussian was determined by the noise parameter, and was never more than 1% of the state variables to which it was added in each integration step. Raw BOLD time series were subjected to global signal regression to match the processing steps carried out on the empirical time series.

In order to take advantage of as much available structural information as possible, while keeping the simulation results in the same parcellation scheme as the empirical FC results, the full structural network based on the 50-region SwS atlas was used as input for the simulations. Once a simulation was complete, a 50-region simulated resting state FC matrix was calculated based on the Pearson correlation between the simulated BOLD time-series for each ROI pair. Comparison with the empirical FC matrix (using raw Pearson correlations) was carried out only considering the 32 regions that were common in the empirical and simulated data by computing the pairwise Pearson correlation between the upper triangular portions of the two matrices.

In order to test the stability of simulation results, multiple runs of simulations with the same parameters and structural data were performed while varying the initial random seed. FC matrices derived from 40 simulations of length 10 s with identical parameters but varying random seeds never had pairwise correlations lower than 0.96, and variability declined with length of the simulations - 10 simulations of length 180 s all had pairwise correlations above 0.99.

2.8.1. Parameter space exploration

In order to maximize the similarity between the functional connectivity matrix derived from simulated time series and the empirical FC matrix, two model parameters were optimized: one scaled global coupling, and another represented the axonal conduction velocity, which controls the time delay of signals between regions. To automate the exploration of the parameter space of multiple TVB models, we adapted the open-source Python utility “Distributed Evolutionary Algorithms in Python” (DEAP) (Fortin et al., 2012), which allows for the implementation of evolutionary optimization algorithms. We used the package to optimize the values of simulation parameters within given ranges. This approach was inspired by BluePyOpt, an open-source software created by members of the BlueBrain project team which implements DEAP to optimize model parameters; however, BluePyOpt is most useful in, and primarily intended for, the optimization of neuron-level models, so using DEAP itself as opposed to BluePyOpt proved simpler for TVB simulations (Van Geit et al., 2016). DEAP implementation allowed us to stochastically search a specified parameter space for any given model, with the aim of finding parameter values that allow TVB models to generate realistic simulated FC as measured by pairwise correlation.

DEAP optimization was performed on all available TVB models, with the exception of two models designed to simulate epileptic activity. Initial parameter ranges were derived from the TVB source code. The evolutionary algorithm implemented using DEAP operated as follows. For each model type, a collection of 600 different combinations of values for each of the global (model and coupling) parameters being optimized were initialized by randomly selecting values within a given range for each parameter - these combinations are analogous to a population of individuals. Individuals with high fitness, i.e. whose parameter value combinations produced simulated FC matrixes with high similarity to the empirical FC matrix, were copied, parameters were randomly swapped between individuals, and some individuals had their parameter values reassigned to random values within that parameter’s range – processes analogous to natural selection, mating, and mutation respectively. One iteration of selection, mating and mutation of the various individuals collectively forms a “generation”. Optimization was performed over the course of 36 successive generations, and the simulated time period of every simulation corresponded to 15 min of real time.

3. Results

3.1. Characterization of SC and rs-fMRI-derived FC

Categorical weights of the RCAMs matrix (Bota et al., 2015) were encoded between 0 (not present) and 7 (very strong) (Fig. 1a, left panel). Out of the 32 selected regions defined for the cortex, 49.8% of the possible pair-wise structural connections were not present, 7.2% were very weak, 8.6% weak, 6.7% weak/moderate, 10.2% moderate, 5.9% moderate/strong, 8.8% strong, and only 2.8% very strong (Fig. 1a, middle panel). These values confirmed the overall sparsity of the anatomical network. We next derived the FC within the same set of anatomical regions, computed as the correlation between the average BOLD signals in each of the 32 cortical regions (Fig. 1b, left panel). As expected, pair-wise functional interactions between the 32 regions included negative correlations, with FC values following a slightly right skewed and unimodal distribution with a pronounced peak around zero (Fig. 1b, middle panel).

Fig. 1

Structural and empirical functional networks of the rat cortex. (a) RCAMs matrix built from (Bota et al., 2015) (left panel), histogram of the weight distribution (middle panel), and structural connections as a function of the Euclidean distance between pairs of ROIs (right panel). Color scale represents the categorical weights of structural links (0, not present; 1, very weak; 2, weak; 3 weak/moderate; 4, moderate; 5, moderate/strong; 6, strong; 7, very strong). (b) Functional connectivity between pairs of ROI time courses (left panel), probability density estimate of the Fisher’s z-values (middle panel), and pair-wise functional connections as a function of the Euclidean distance between ROIs (right panel). Color scale represents the strength of functional interactions between pairs of nodes. The same color scale for structural and functional networks is used throughout the work. Each dot represents a particular functional connection between pairs of nodes. Euclidean distances are measured between centers of mass.

Although different studies have revealed that wiring minimization itself cannot account for all topological features of anatomical and functional networks (Betzel et al., 2016; Véertes et al., 2012), the essential role of physical distances between pairs of nodes is widely recognized. We plotted the pair-wise SC and FC values as a function of distance (approximated as Euclidean distance) and found that both SC and FC values tend to be higher between brain regions separated by short distances (Fig. 1, right panels).

With the aim of investigating the relationship between the SC and the rs-fMRI-measured FC in the rat cortex, we first calculated the association between anatomical and functional matrices (Fig. 2). The Spearman correlation quantified to what extent the strength of a structural connection predicted the corresponding functional connection between two regions. Using all possible value pairs in the SC connectome (i.e., present and not present connections) and the FC matrix, we obtained a Spearman (rank-order) correlation of = 0.48 ( < 2 × 10−29) between the structural and functional networks. The Spearman correlation did not vary after removing all node pairs for which an anatomical connection was found to be not present, = 0.48 ( < 2 × 10−19). This result reveals a strong positive linear association between SC and FC matrices. In addition, we also calculated the association between both networks while controlling the effect of the Euclidean distance. In this case, the Spearman correlation dropped to = 0.35 ( < 3 × 10−15) and = 0.33 ( < 2 × 10−9), relative to the total or zero weight-corrected SC matrix, respectively, supporting the role of geometric distances between region pairs in brain network organization. Functional connectivity patterns without global signal regression and the relationship with the SC matrix are shown in Supplementary Fig. 6.

Fig. 2

Pair-wise functional interactions as a function of the underlying structural connections. Box and whiskers plots representing the median (horizontal red line), first, 1, and third, 3, quartiles (blue box), as well as the lower and upper (black) whiskers of the functional connections. Those values greater than 3 + 1.5*(3 − 1) or lower than 1 − 1.5*(3 − 1) were considered as outliers (red crosses).

3.2. Relationship between functional modules and structural network-level patterns

Partitioning a network or graph into components offers the ability to uncover modules at different spatial scales. Indeed, each module or component can be interpreted as a subnetwork itself. Fig. 3a shows the number of communities extracted from the functional data along with the value of the modularity function, which represents the quality of the obtained consensus partition, and the consistency of the partitions as measured by the similarity of the modules generated in 10,000 iterations of the partition algorithm (see subsection 2.5). Given that different community structures might be equally feasible and plausible, we investigated the distribution of motif classes within functionally define modules at different scales of resolution, rather than targeting a unique partition. We examined structural patterns at = −0.0125 (Fig. 3b), = 0.005 (Fig. 3c), and = 0.02 (Fig. 3d); these parameter settings were selected for several reasons. First, these scales were located within intervals of values where the number of communities remained relatively constant, hence yielding partitions with two ( = −0.0125), three ( = 0.005) or five ( = 0.02) communities. Second, given a number of communities, the selected scales were the most stable partitions, as indicated by a high mean and a low variance of the z-Rand index. Finally, we did not explore functional partitions beyond = 0.02 because the resulting partitions yielded numerous modules of very small size, which makes motif analysis impractical. For example, at = 0.035, the cortex was partitioned into eight modules, with three modules of three nodes, one module of two nodes, and another singleton module.

Fig. 3

Community structure. (a) Different properties were computed at various levels (Louvain algorithm run 10,000 times): number of communities, modularity function (Q), measuring the quality of the obtained consensus partition, and mean and variance of the z-score of the rand indices (z-Rand), measuring the similarity over all pairs of partition within scales. Partitions below −0.015 and above 0.24 were unstable. For clarity, only scales in the range [−0.015, 0.15] are shown, and Q was rescaled between 0 and 10. (b) Functional partition obtained at = −0.0125 (two communities). (c) Functional partition obtained at = 0.005 (three communities). (d) Functional partition obtained at = 0.02 (five communities).

At the most coarse scale ( = −0.0125), the cortex was partitioned into two modules (Fig. 3b). The wider module (M1) contained most regions of the sensory-motor cortex including the primary somatomotor (MOp), somatosensory (SSp, SSs), visceral sensory-motor (VISC, ILA), gustatory (GU), olfactory (AON, TT, PIR), auditory (AUDp, AUDs), and visual areas (VISs, VISp). Association areas such as orbital (ORBmv, ORBvl), agranular insular (AId, AIv, AIp) and posterior parietal (PTLp), together with the cortical subplate, namely, claustrum (CLA) and the dorsal part of the endopiriform nucleus (EPd) were grouped in the same module. On the other hand, the second module (M2) brought together all parts of the hippocampal formation included in this analysis (HIPv, HIPd, SUBd, POST) as well as retrosplenial (RSPv, RSPagl, RSPd), prelimbic (PL), and anterior cingulate (ACAv, ACAd) areas. Along this association structures, only one sensory-motor region was included, namely, the secondary somatomotor region (MOs). As the resolution parameter was increased ( = 0.005), the cortex was split into three functional modules (Fig. 3c). The previous module M1 was divided into two different modules, whereas M2 remained largely intact, with only the infralimbic area relocated from M1 into M2. Finally, at the finer scale ( = 0.02), a total of five modules were obtained (Fig. 3d). In this case, the original module M1 was partitioned into three communities and M2 into two. Again, the infralimbic area was relocated with module M2.

The next part of our analysis focused on investigating structural connections and their relationship with functional modularity. Regardless of the resolution scale, both intra- and inter-modular links showed positive reciprocity, indicating that brain regions tend to be mutually connected within and between modules. However, when comparing reciprocity between both types of links, connections of the RCAMs matrix within communities showed a stronger reciprocity compared to projections connecting nodes from different modules: 0.506 vs 0.275, 0.534 vs 0.304, and 0.575 vs 0.326 when the brain cortex was partitioned into two, three, and five modules, respectively.

We next investigated the density of each of the 13 potential motif classes associated with structural motifs of size 3 within functionally identified modules and for the selected scales (Fig. 4). We obtained a very high – for motif class 13 for all partitions and for all connectivity levels investigated (when keeping all connections, when removing very weak links, and when removing both very weak and weak/moderate links), indicating that this maximally densely connected structural motif is significantly enriched within functional communities. Comparable results were obtained for motif class 12 for scales 0.005 (Fig. 4b) and 0.02 (Fig. 4c), and for −0.0125 (Fig. 4a) after deletion of very weak links and of both very weak and weak/moderate links. A few other motif classes were significant at > 3:7( < 0.0001), for example motif class 10 and 11 for scales 0.005 and 0.02 and only when removing both very weak and weak/moderate links. Regardless of the considered scale and connectivity level, the significance level of all other motif classes was ≤ 3:7.

Fig. 4

Network motif analysis. Functionally obtained partitions are imposed over the RCAMs matrix (upper row). Statistical evaluation of each of the 13 structural motifs classes (lower row) is shown for different binarized versions of the RCAMs matrix, namely, when keeping all connections (threshold = 1), removing very weak links (threshold = 2), and removing both very weak and weak/moderate links (threshold = 3). (a) Motif analysis at = −0.0125. (b) Motif analysis at = 0.005. (c) Motif analysis at = 0.02.

3.3. Dynamic simulations based on the rat connectome

Simulations of spontaneous rat cortical activity produced patterns of simulated FC that were significantly correlated with empirical FC, reaffirming that the topology of the rat cortex SC predicts and shapes functional relationships. To achieve realistic model dynamics we identified optimal model parameters, in this case the conduction velocity and a global coupling parameter that scales the strength of structural connections. Optimal values were identified by conducting stochastic optimizations through an evolutionary algorithm powered by the DEAP python package (see subsection 2.8.1).

We subjected a broad range of computational models available in TVB to joint optimization of their conduction velocity and global coupling parameters through DEAP. DEAP performed automatic parameter optimizations, with only parameter ranges as input. The highest performing model available in TVB, the Generic 2D Oscillator model, produced a simulated FC matrix with a Pearson correlation of = 0.53 with the empirical FC matrix ( < 10−70) when the conduction velocity parameter was set to 3.06 m/s and the global coupling parameter had a value of 0.3734 (Fig. 5). These simulation results implementing a generative model for FC based on the network of structural connections and simple biophysical node dynamics demonstrate are consistent with the notion that a large proportion of the observed patterning of FC is shaped by the underlying SC matrix.

Fig. 5

Comparison between experimental and simulated functional connectivity (FC) matrices. (a) Simulated FC network. (b) Empirical FC network. (c) Probability density estimate of the simulated FC matrix. (d) Scatterplot showing the linear dependency between experimental and simulated FC matrices.

4. Discussion

The study of structure and function at the large-scale in rodent animals offer the possibility of understanding the mechanisms underlying psychiatric and neurological disorders for subsequent translation research (Gozzi and Schwarz, 2016; Jonckers et al., 2015). Nevertheless, knowing and describing the healthy rodent brain connectivity is the first step for that endeavor. The organizational principles of the rat and mouse structural networks have been previously described (Bota et al., 2015; Oh et al., 2014; Swanson et al., 2016; van den Heuvel et al., 2016b). Comparisons between the mouse connectome and rs-fMRI connectivity have revealed that both networks are intimately correlated (Sethi et al., 2017; Stafford et al., 2014). In this work, we have carried out comparisons between the cortical rat connectome (Bota et al., 2015) and the corresponding functional network obtained from rs-fMRI data. The main findings of this study reveal that: (1) in agreement with previous studies in humans and various animal models, the level of functional interaction between two anatomically connected brain areas in resting-state is significantly predicted by the strength of the underlying structural connection; (2) fMRI BOLD signal fluctuations in rat brain cortex can be robustly partitioned into functional modules or clusters; (3) reciprocity of links connecting nodes from the same functionally defined modules is stronger than those linking regions from different modules; (4) within functional communities under a control condition and in resting-state, densely interconnected structural motifs (i.e. class 12 and 13 for motifs of size = 3) are significantly enriched; and (5) simulations of the intrinsic neural activity based on the rat connectome can reproduce much of the empirical functional patterns obtained empirically.

4.1. Role of higher-order patterns in shaping cortical functional modular organization

Comparing SC and FC without taking into account network-level effects only supplies a partial view of the inherent complexity of the brain. Therefore, network approaches are increasingly gaining prominence to understand the SC-FC coupling (Wang et al., 2015). Here, we have first shown that the 32 cortical brain regions considered in this work disclose a distinct community structure, partitioned into two, three or five communities, at = −0.0125, = 0.005 and = 0.02, respectively. Previous studies have investigated the modular structure of functional connections over the whole rat brain under different states of consciousness (D’Souza et al., 2014; Liang et al., 2012). Other authors have focused on specific resting-state networks (e.g. the default mode network), revealing functional subcomponents and how interactions between and within these submodules can be modulated in age-related neurocognitive disorders (Hsu et al., 2016). Here, we have restricted our analysis to the brain cortex. Future work may address additional brain regions, for example including the association connectome of the rat cerebral nuclei (Swanson et al., 2016).

After evaluating the modular organization of the cortex, we first observed a stronger weighted link reciprocity within functionally defined modules compare to inter-modular links. We then observed a high ratio of structural motif classes 12 and 13 within functional modules. These findings are consistent across resolution scales. In particular, whereas motif class 12 is overrepresented at scale −0.0125 after removing very weak links and both very weak and weak/moderate links, motif class 12 at scales 0.005 and 0.02, and motif class 13 at all resolution scales were strongly significant for all thresholds applied. This indicates that the functional coherence of FC communities may in part reflect an underlying aggregation of densely connected anatomical motifs. Note that motif analysis was carried out on three binarized versions of rat connectome since we were interested in the count of motif classes rather than in the strength of the connections. Additionally, the abundance of these anatomical subgraphs can not be trivially explained by lower-order features because the generated network null models also preserved the number of mutual edges of individual nodes (Milo et al., 2002). The important role of motifs in shaping the SC-FC coupling has been previously observed in the macaque cortex (Adachi et al., 2012). Nevertheless, whereas these authors computed motifs over the whole network, motif analysis in this work was restricted to within modules and to small motifs of size 3. For graphs with a greater size, a richer repertoire of larger structural motifs could be investigated (e.g., network motifs of size 4 and higher (Sporns and Kötter, 2004)). There has been an increasing interest in detecting coherent groups of nodes forming specific network motifs (Arenas et al., 2008; Benson et al., 2016). In this framework, one can specify in advance a particular motif class of interest to guide the clustering process. When domain-specific knowledge is not available, it is also possible to analyze which type of motif organizes and shape the modular structure of complex networks.

The term “motif” is widely applied for describing and understanding recurring circuits of interactions that take place on real-world networks (Alon, 2007). For instance, in the context of large-scale brain networks, as in the present report, structural motifs consist of a set of nodes and potential pathways supporting communication, whereas functional motifs refer to specific combinations of connections that can be activated within structural motifs (Sporns and Kötter, 2004). That is, both structural and functional motifs are defined for directed anatomical networks (Table 2). In other settings, for example when approaching the functioning of the brain from a computational perspective (Turkheimer et al., 2015) or in the context of activation spreads recorded with different modalities (Frostig et al., 2008; Mohajerani et al., 2013), the concept of motif is related but not completely the same as that employed in this work. The common message behind this general concept is that the brain is a highly hierarchical and complex system wherein organizational principles are repeated across space and time (Turkheimer et al., 2015).

4.2. Modelling how structure shapes function

The high levels of similarity between empirical FC and FC derived from simulations using multiple computational models lend additional support to the finding that structural connectivity both constrains and enables functional connectivity in mammalian brains (Honey et al., 2009). Structural connectivity data derived from tract tracing experiments is generally considered to more closely approximate “ground truth” than connectivity derived from DTI streamlines, with lower levels of both type I and type II errors (Chen et al., 2015; Dauguet et al., 2007). The high quality of the structural data, combined with the optimization of model parameters through DEAP, likely contributes to the high levels of FC prediction seen across different nonlinear models (max = 0.53, 33.5% of variance explained). This level of correlation is somewhat higher than that typically achieved with nonlinear models using SC data derived from DTI images (Falcon et al., 2016; Messé et al., 2015). Importantly, the structural network derived from tract tracing experiments is directed, in contrast with SC obtained from DTI, which may allow a neural mass model to produce a richer set of dynamics through asymmetric connections between regions and more accurately defined pathways of communication between regions (Knock et al., 2009). Our findings from the simulation model support the notion that a (simple) generative model based on structural connectivity and neural mass dynamics can account for a significant portion of the observed patterns in empirical FC.

There are several limitations to the simulation portion of this work. First, while the structural data used here is presumably of higher quality than DTI, it is incomplete in that many parts of the rat brain are not included. For example, the omission of sub-cortical structures such as the cerebral nuclei or the thalamus may limit the potential range of dynamics captured in the simulation. Furthermore, interhemispheric connections are not included in the SC matrix, which may be important for simulating more realistic neural dynamics. Second, our estimate of the anesthetized rat hemodynamic response function may not accurately represent rodent BOLD responses. The BOLD forward model available in TVB is intended for use with primate data, and does not incorporate possible regional variations in the hemodynamic response function (Sanz-Leon et al., 2015). Third, the stochastic nature of the DEAP based parameter optimization may result in the algorithm converging on local optima. These limitations can be addressed in future work as additional parts of the rat connectome become available (e.g. (Swanson et al., 2016)), with more detailed physiological recordings of rodent fMRI responses, and further improvements in simulation and optimization algorithms.

4.3. Construction of brain graphs

A very critical step for construction and analysis of brain network is that of dividing the brain into discrete regions. The absence of a well-established and accessible Swanson atlas in NIfTI format for parcellation of fMRI data required us to perform a manual matching procedure between cortical regions of Swanson space (SwS) (Swanson, 2004) and Paxinos&Watson (PWS) (Paxinos and Watson, 1998). This entails several caveats. First, although Swanson regions can also be expressed in stereotaxic coordinates to correspond to Paxinos & Watson regions, coronal planes from both spaces do not always exactly match with respect to bregma reference. Therefore, taking into account the spatial distribution of adjacent planes is fundamental to achieve an accurate alignment of brain regions between different spaces. Second, a common problem when using different spaces or atlases (not only restricted to SwS and PWS) concerns nomenclature of brain regions, as different spaces may use different nomenclature to label the same brain regions. In this case, we have used the Swanson nomenclature after carefully inspecting both SwS and PWS (see subsection 2.3 and Table 1). Finally, some cortical regions were not available in the NIfTI atlas as a single mask (e.g., CA1 and CA2). To overcome these limitations, several areas from the anatomical matrix were aggregated together by averaging the original categorical weight. In addition, some cortical regions were discarded to ensure a representative time course in each node from functional networks. All 32 ROI signals included in our analysis were averaged across hemispheres, had exactly the same volume (i.e., 1 mm3 in total, 0.5 mm3 in each hemisphere), and were extracted from exactly the same voxels across rats. The consequences of these methodological steps might inevitable mask relevant feature of individual nodes. For instance, parahippocampal regions, such as perirhinal and entorhinal cortices, were not included in our analysis and, however, these areas have shown to be critical for cortico-hippocampal integration in mice (Bergmann et al., 2016). Despite these considerations, the main research goal of the present report was to assess the empirical SC-FC coupling in the rat brain cortex with a special emphasis in comparisons at regional scale (functional modularity and structural network motifs), rather than the analysis at more local scale (Ferezou et al., 2007; Kaiser, 2011).

4.4. Relevance of data preprocessing: global signal regression

An important issue in fMRI investigations is that of controlling for head motion and physiological artifacts (mainly those related with respiratory and cardiac cycles) when estimating functional interactions between brain areas, especially in task-free paradigms (Van Dijk et al., 2012; Yan et al., 2013). Even though anesthesia can reduce rough head movements, variance explained by motion parameters can substantially improve the specificity of FC in anesthetized rodents (Kalthoff et al., 2011). Numerous techniques have arisen with the aim of correcting or minimizing variations in the BOLD signal as a consequence of these artifacts: scrubbing (Power et al., 2012), spike regression (Satterthwaite et al., 2013), motion parameter regression (Friston et al., 1996), Independent Component Analysis (ICA) (Pruim et al., 2015), as well as its combined use with global signal regression (GSR). The inclusion of the global signal as a nuisance regressor is a widely used but controversial preprocessing step. How negative correlations should be interpreted after GSR is very challenging given the nature of the BOLD signals (Fox et al., 2009; Murphy et al., 2009). In our dataset, Spearman correlation demonstrated a significant linear dependence between SC and empirical FC after controlling for the global signal. Whether or not other pre-processing steps would further increase the coupling between BOLD signals and the underlying SC is an open question that would require a separate study comparing the effects of different preprocessing pipelines (Ciric et al., 2017). The removal of global signal variance has also applied in a previous study relating structure and function in the mouse brain (Stafford et al., 2014). It is important to acknowledge that the application of GSR has shown to be effective to reduce the relationship between connectivity and motion, but accentuates distance-dependent effects. On the other hand, temporal censoring techniques (e.g., scrubbing or spike regression) do appear to be suitable to reduce distance-dependent effects, but using additional degrees of freedom in return (Ciric et al., 2017). In this work, we preprocessed rs-fMRI data combining these techniques and ensuring a minimum temporal window of 9 min to estimate functional interactions between node pairs.

4.5. Effect of anesthesia on functional connectivity

Another very important consideration when examining fMRI data in rodents, and other animal models, has to do with the effect of the anesthetic agent on brain physiology (Pan et al., 2015). Due to the multi-compartmental origin of fMRI signals (Moreno et al., 2013), interference of anesthetics on the neural computations, as well as on the neurovascular coupling itself are expected. Advantages and disadvantages of a number of anesthetic compounds, as well as mixtures of them, have been reported in the literature (Paasonen et al., 2016; Williams et al., 2016). It is well established that the degree of functional inter-hemispheric coupling is dose-dependent, and previous studies have uncovered bilateral functional networks in the rodent brain in a variety of anesthetic agents at lower doses (Gozzi and Schwarz, 2016). For instance, whereas the use of urethane anesthesia at higher levels (1.5 g/Kg) has compromised functional inter-hemispheric coupling in mice, cortico-cortical correlations improves when lower doses (1.2 g/Kg) are used (Grandjean et al., 2014). In this work, where subjects were urethane-anesthetized at 1.2 g/Kg, BOLD signals from homotopic ROIs were averaged across hemispheres because bilateral networks were uncovered using different seeds (Supplementary Figs. 3, 4 and 5). We note, though, that for group studies where subtle inter-hemispheric differences could be of relevance, extraction of ROIs should be performed separately. Finally, it is important to highlight that the significant match found in our study between FC and SC could be favored by a decrease in the repertoire of functional configurations that the system can visit under anesthesia, as previously reported in monkeys (Barttfeld et al., 2015). Therefore, the remaining functional interactions under anesthesia would most likely reflect the strongest, hard-wired, connections in the brain, this is, the structural connectivity.

5. Conclusions

We have proved in this work the structure-function relationship at large-scale of the rat by directly comparing the structural and functional connections spanning the brain cortex (maximum Spearman rank-order correlation between was type of networks was 0.48), and relating second-order properties (reciprocity) and higher-order structural connectivity pattern (network motifs) with functional communities. We highlight here the importance of densely connected structural motifs in shaping the community structure of functional networks in resting-state. The accuracy of computational models of the spontaneous neural activity based on the rat structural connectome (maximum Pearson correlation between the empirical and functional FC was 0.53) further supports the idea that structural connectivity is coupled to and shapes functional connectivity in cortical networks.