Directed functional and structural connectivity in a large-scale model for the mouse cortex

The communication between brain regions is often analyzed using structural and functional connectivity (Avena-Koenigsberger, Misic, & Sporns, 2018). The former refers to anatomical connections between brain regions generally quantified using tracer injections or diffusion magnetic resonance imaging (Ambrosen et al., 2020). The map of these connections is called “connectome” (Sporns, Tononi, & Kötter, 2005). Network measures are usually used to analyze the connectome, whereas nodes represent brain regions and edges refer to axonal projections (Bassett & Sporns, 2017; Bassett, Zurn, & Gold, 2018). Functional connectivity estimates brain communication from statistical relations between recorded brain signals (Avena-Koenigsberger et al., 2018; Reid et al., 2019). Particularly, directed functional connectivity methods use the concept of causality to infer both the intensity and the direction of the connections between brain regions (Bastos & Schoffelen, 2016). Even though there is some association between structural and functional connectivity, the relationship between them is not straightforward (Avena-Koenigsberger et al., 2018). While the former is practically static and composes the map of possible pathways for information flow between brain regions, the latter changes continuously and depends, for example, on the dynamical states of brain regions, noise, and strength of structural connections (Nunes, Reyes, & De Camargo, 2019).

During electrophysiological procedures, researchers typically record brain signals using electrodes positioned in different depths of brain regions. Even with the improvement in technologies for recording signals, it is usually possible to record signals only from a few areas compared with the number of sources of activity in the brain (Harris, Quiroga, Freeman, & Smith, 2016; Hong & Lieber, 2019; Schölvinck, Leopold, Brookes, & Khader, 2013). Thus, the functional connectivity analysis presents a problem because many unrecorded regions may indirectly influence other regions as common inputs (Bastos & Schoffelen, 2016; Reid et al., 2019; Sanchez-Romero & Cole, 2019). Therefore, the comparison between structural and functional connectivity becomes more complicated since spurious inferred causality relations can lead to misinterpretations of electrophysiological data.

Previous simulation studies evaluated the relation between directed functional connectivity and structural connections (Baccalá & Sameshima, 2001; Barnett & Seth, 2014; Mi, Cheng, & Zhang, 2014; Nunes et al., 2019). However, most of these studies used either autoregressive (Novelli, Wollstadt, Mediano, Wibral, & Lizier, 2019) or rate-based models (Mi et al., 2014) for the dynamics of each cortical area. These studies provided essential steps towards evaluating the reliability of causality measures. However, the time series obtained from autoregressive and rate models are distant from electrophysiological signals obtained in experimental laboratory conditions. Using spiking models, we can capture the dynamic of neuronal networks while generating simulated local field potential (LFP) signals from the synaptic currents. Also, most studies do not consider the impact of accessing only part of the activity in the brain.

In this work, we investigate the relationship between directed functional connectivity and structural connectivity in a large-scale network model of the cortex, derived from a cortical connectome of the mouse obtained using tracer injections (Gămănuţ et al., 2018). We used generalized partial directed coherence (GPDC), a frequency-domain method based on multivariate vector autoregressive (MVAR) models, which provides estimates of directed functional connectivity (Baccalá, Sameshima, & Takahashi, 2007; Sameshima & Baccalá, 2014). The mean correlation between the fraction of labeled neurons (FLN) and GPDC remained high (r > 0.6) even when only a few cortical areas were considered in the GPDC calculation, indicating that this causality measure provides reliable results in typical experimental conditions in which only recordings from a subset of areas are available.

The parameters used in this neuron model were previously reported and applied in some studies that modeled cortical neuronal populations (Barardi, Sancristóbal, & Garcia-Ojalvo, 2014; De Sancristóbal, Vicente, Sancho, & Garcia-Ojalvo, 2013; Sancristóbal et al., 2014).

The characteristic decay time τ is 2 ms and 8 ms for excitatory and inhibitory synapses, respectively. When a presynaptic neuron j fires a spike at time t_j, g_syn,i,j is incremented by a synaptic weight w after the axonal delay d, which was set as 1 ms for all intra-areal connections (Sancristóbal et al., 2014). The value of w depends on the excitatory/inhibitory nature of the presynaptic and postsynaptic neurons. Furthermore, all neurons receive a background input given by a heterogeneous Poisson-process spiking activity with a rate of 7.3 kHz (Sancristóbal et al., 2014). The background input acts as an excitatory synaptic current. To add heterogeneity in our model, all synaptic weights w for recurrent connections and background input were taken from a Gaussian distribution (Table 1).

The mouse cortex’s large-scale network model is composed of 19 areas where a spiking neuronal population models each area with long-range and recurrent synapses. Parameters related to recurrent synapses were described in the previous session. Neurons from different areas are randomly connected with probability p_inter = 5%. The synaptic weights between cortical areas are based on the previously published anatomical connectivity dataset for the mouse cortex (Gămănuţ et al., 2018) obtained by retrograde tracer injections (Markov et al., 2014).

Long-range connections are excitatory, targeting either excitatory or inhibitory neurons with synaptic weight, $wlr,Ei$⁠, and $wlr,Ii$⁠, respectively. The index j represents the source area, i represents the target area, and N is the total number of simulated cortical areas. The axonal delay for long-range connections is given by the ratio between the inter-areal anatomical distance estimates between cortical areas and the conduction speed set as 3.5 m/s (Choi & Mihalas, 2019).

I_E,i accounts for both the local (within population) and global (inter-areal projections) excitatory synaptic currents, while I_I,i corresponds to the local inhibitory current. I_bkg,i is the synaptic current related to the background Poisson input. R represents the resistance of a typical electrode used for extracellular measurements, here chosen to be 1 MΩ (Sancristóbal et al., 2014). N_E is the number of excitatory neurons in each neuronal population.

The mean was subtracted from the simulated LFP signal. The resultant signal was filtered using a 1 kHz low-pass filter to avoid aliasing and downsampled to 1 kHz.

The MVAR model was estimated by the method of ordinary least squares (OLS; Hamilton, 1994). We used Akaike’s information criterion (AIC) to select model order (Supplementary Equation S1), choosing the order p ≤ 50 that had the minimum AIC (Supplementary Figure S6) value.

GPDC has values in the range from 0 to 1 and is invariant to scale, so the normalization of time series is unnecessary (Baccalá et al., 2007; Sameshima & Baccalá, 2014). Similar to other (directed) functional connectivity measures, unrecorded time series can lead to spurious estimates. Therefore, the reliability of estimates depends on the number of time series included in the estimates. For all analysis we used the peak GPDC value over all frequencies [0, $fs2$].

All simulations were performed using the simulator Brian2 (Stimberg, Brette, & Goodman, 2019) applying the exponential Euler method (Bower & Beeman, 2012) to integrate the differential equations with an integration step of 0.1 ms. Each simulation was 30 s long, generating sufficient data points to apply GPDC on the simulated LFP signals (Sommariva, Sorrentino, Piana, Pizzella, & Marzetti, 2019).

The large-scale network model of the mouse cortex contains 19 spiking neural populations with recurrent connections and excitatory long-range connections between populations, constrained by the directed and weighted structural connectome (Figure 1A and Figure 1B). The dynamical behavior of each simulated cortical area is predominantly asynchronous with transient spike synchronization (Palmigiano, Geisel, Wolf, & Battaglia, 2017; Uhlhaas et al., 2009) (Figure 1C), with the typical power spectral density (PSD) of LFP signals displaying a peak in the gamma band (Figures 1D and 1E) (Buzsáki & Wang, 2012). The firing rate of inhibitory neurons is 4.74 ± 0.11, higher than the excitatory neurons’ rate of 3.64 ± 0.42 (Figure 1F). Differences in population behavior are mostly due to inputs from other areas since we sample their parameters from the same distributions.

We first compared the FLN values with the average GPDC over 10 simulations of the model. Most medium to strong connections from the structural connectome were also captured by the directed functional connectivity (Figure 2A and Figure 2B). We used the GPDC largest value (peak), but other approaches such as the average of GPDC over frequencies and area under the GPDC curve (Supplementary Figure S1) produced similar results.

Although the graph density of the structural connectome is 97% (Gămănuţ et al., 2018), most structural connections are weak, which leads to a prevalence of weak average GPDC values. Weak structural connections is a characteristic shared by connectomes from different mammals, with FLNs varying by several orders of magnitude, log-normally distributed (Buzsáki & Mizuseki, 2014; Gămănuţ et al., 2018; Markov et al., 2014; Theodoni et al., 2020). To evaluate the relation between structural and directed functional connectivity, we plotted GPDC values from 10 simulations against FLNs and fitted a linear model, obtaining the Pearson correlation r (Figure 2C). The scatterplot presents most points close to the origin due to the predominance of small values for the GPDC and FLN. The average Pearson correlation between FLN and GPDC is 0.74 (Figure 2D). We also verified that the average correlation between GPDC and FLN over bootstrap samples of 80 randomly chosen edges is 0.74 (Supplementary Figure S2). This correlation level is close to those obtained by other works that analyzed different structural connectomes using functional connectivity applied to empirical data (r ≈ 0.79; Hagmann et al., 2008) or firing rate models (r ≈ 0.73; Chaudhuri, Knoblauch, Gariel, Kennedy, & Wang, 2015).

The centrality of the cortical area seems to influence the variability of GPDC estimates over simulations. The variability of directed functional connectivity was measured by the coefficient of variation of GPDC (Figure 3A). The centrality, measured by the nodal in-strength (i.e., the sum of inward FLNs to a cortical area; Figure 3B), is positively correlated (r = 0.64) to the sum of the coefficients of variation (CVs) of the connections emerging from that area (source; Figure 3C). When the cortical area is considered the target of directed functional connectivity, the correlation with nodal in-strength is negative (r = −0.52) (Figure 3D). We performed the same analysis correlating the sum of coefficient of variation with eigenvector centrality (Supplementary Figure S3), and we obtained the same relationship, but with smaller Pearson correlation coefficients (r = 0.59 and r = −0.44). We should note that in both cases (source and target), the actual variability (standard deviation) increases with larger nodal in-strength values (Supplementary Figure S4).

We also investigated the relationship between the firing rate in a cortical area and the estimated activity that is arriving at this cortical area mediated by structural or directed functional connectivity pathways. The propagation of activity in the cortex is constrained by direct anatomical connections between areas and indirect paths (Vézquez-Rodríguez, Liu, Hagmann, & Misic, 2020), with the propagation of activity occurring mainly through the strongest long-range projections (Joglekar et al., 2018). The estimated activity mediated by FLNs is strongly correlated to the target areas’ firing rate (Figure 4A), while the correlation of estimated activity mediated by GPDCs and firing rates was 0.54 (Figure 4B). This indicates that GPDC estimates can be used to infer the propagation pathways, although less reliably than when using FLN values directly.

We analyzed the behavior of GPDC estimates when considering a reduced number of areas, reproducing typical experimental setups. We considered a visual and a frontoparietal cluster, each containing seven cortical areas (Gămănuţ et al., 2018) (Figure 5A). We evaluated the distribution of correlation between FLN and GPDC when GPDC estimates between all areas of each cluster are conditioned on the whole connectome, conditioned on the areas in each cluster, and using only pairwise (bivariate) estimates (Figures 5B and 5C). This analysis simulates the situations where an electrophysiologist has information only from a single cluster of cortical areas or a pair of areas. The highest correlations between the GPDC and FLN occurred when we conditioned GPDC to the whole connectome, followed by GPDC conditioned to the cluster area, and pairwise GPDC. Also, the correlation for the frontoparietal cluster was higher than for the visual cluster in all scenarios.

We extended the analysis to evaluate the effect of cluster size on GPDC correlation to FLN. We used cluster sizes ranging from 3 to 15 areas. We created 150 random clusters sampled from all areas in the connectome for each cluster size and computed the Pearson correlation for the GPDC (a) conditioned on the whole connectome, (b) conditioned on the cluster areas, and (c) evaluated using pairwise data. For cases (a) and (b), the Pearson correlation increases, and the standard variation decreases as we increase the cluster size (Figure 6), showing that it is advantageous to include more areas in the GPDC calculation. Surprisingly, the correlation between structural and directed functional connectivity when using simulated signals from a few cortical areas (blue dots) is similar to using signals from the whole cortex (black dots), with most points showing statistically different results. The bivariate GPDC (Figure 6) had a statistically significant lower average Pearson correlation for all cluster sizes with four or more areas, indicating that these measures are affected by interference from ignored signals.

Our results shed light on the relationship between structural and directed functional connectivity in circumstances similar to those faced by electrophysiologists. They indicate that the reliability of directed functional connectivity estimates and their relationship with structural connectivity depends on the number of areas considered. Nevertheless, the GPDC conditioned on few cortical areas had similar results to the GPDC conditioned on all areas, providing evidence that it is possible to obtain a reasonable relationship between structural and directed functional connectivity in electrophysiological experiments even with signals from few areas.

Previous studies evaluated the relationship between structural and functional network connectivity strength on electrophysiological data (Straathof, Sinke, Dijkhuizen, & Otte, 2019), with some using undirected functional connectivity measures (Grandjean, Zerbi, Balsters, Wenderoth, & Rudin, 2017; Stafford et al., 2014). But in electrophysiology studies, researchers do not have access to signals from unrecorded areas and have only estimates of structural strengths from tracers. Using large-scale network models solves this problem, as the researcher has access to all variables in the system, allowing a better understanding of the obtained functional connectivity results.

The relationships between structural and functional connectivity have been largely unexplored through large-scale network models (Bansal, Nakuci, & Muldoon, 2018), and the existing models use neural mass descriptions (rate models) to describe each area’s activity (Honey et al., 2009; Novelli & Lizier, 2021). However, information propagated between brain regions can be characterized not only by the rate code but also by the temporal code (Bieler et al., 2017; Hahn, Ponce-Alvarez, Deco, Aertsen, & Kumar, 2019; Kumar, Rotter, & Aertsen, 2010; Luczak, McNaughton, & Harris, 2015; Seth, 2015), and hypotheses are pointing to spike-timing and spike coherence as essential components of cortical communication (Hahn, Bujan, Frégnac, Aertsen, & Kumar, 2014; Palmigiano et al., 2017; Tiesinga & Sejnowski, 2010). Spiking neuronal populations have richer dynamical behaviors than rate models and better resemble cortical activity; through spiking neuronal networks it is possible to investigate the consequences of spike synchronization (Palmigiano et al., 2017), model different approaches for the propagation of information (Hahn et al., 2019; Joglekar et al., 2018), and generate simulated LFP signals from the synaptic currents, which better resemble biological LFP signals (Mazzoni et al., 2015). Moreover, our obtained correlations are in the same range as the studies using more complex electrophysiological data (Straathof et al., 2019).

The centrality of a cortical area affects the variability of GPDC estimates in different ways when such area is examined as the target or source of functional connections. Strong functional connectivity generally occurs between areas with direct structural connections (Honey et al., 2009), and network measures applied to structural connections can help predict the resting-state functional connectivity (Goñi et al., 2014). However, as far as we know, no previous work has indicated that the variability of directed functional connections could be partially explained by centrality measures applied to structural connectivity. We also noticed that synchronization is strongly correlated to the centrality of the node (Supplementary Figure S5). So it is likely that stronger long-range connections targeting an area increase the synchronization of spikes in this area, and the increased synchronization changes the variability in directed functional connectivity. Indeed, it was observed in previous work that synchronization has an important role in directed functional connectivity (Palmigiano et al., 2017).

The firing rate of cortical areas is explained by the estimated activity flow, as proposed by Cole et al. (2016). When using GPDC as an estimate of structural connections, the correlation between actual and estimated activity in the target area decreases to 0.54. This indicates that directed functional connectivity can be used to estimate the activity flow. Although it is less reliable than when using the actual structural connection strengths, researchers may only have access to directed functional measures.

The relationship between structural and directed functional connectivity is the largest when GPDC is conditioned to all areas in the connectome and decreases as we reduce the number of areas. Gămănuţ et al. (2018) identified six clusters in the mouse connectome (prefrontal, frontal, parietal, cingulate, temporal, and visual) based on the same approach used to investigate the macaque cortex (Ercsey-Ravasz et al., 2013). We evaluated the relationship between GPDC and FLNs in the visual cluster and in a combination of the prefrontal, frontal, and parietal clusters, which we called frontoparietal. We did not use the other clusters, which had a small number of regions. The average correlation was in the range of correlation obtained for random clusters, with r = 0.76 for the frontoparietal and r = 0.50 for the visual cortex. This indicates that within anatomical clusters the relationship between GPDC and FLNs does not change in relation to randomly selected areas. These results also show that GPDC estimates provide statistical information on structural connections even when considering only a few areas. However, when considering individual connections, there can be large differences between GPDC estimates and actual structural connection strengths.

Our large-scale network model has some limitations. First, modeled neuronal population parameters are drawn from the same distributions with activity in the gamma band range (Figure 1). The activity of cortical areas in mice occurs in multiple frequency ranges (Sherman et al., 2016; Tsurugizawa, Djemai, & Zalesky, 2019) and the relationship between structural and functional connectivity depends on the frequency (Vezoli et al., 2020). A second limitation is that we do not model changes in network states, which are known to influence functional connectivity (Stitt et al., 2017). Some studies in computational neuroscience have already explored multistability and temporal patterns of functional connectivity (Deco & Jirsa, 2012; Golos, Jirsa, & Daucé, 2015; Orio et al., 2018). Finally, we considered only cortical areas in our large-scale network model, excluding subcortical areas, which have a more complex dynamic (Erö, Gewaltig, Keller, & Markram, 2018; Knox et al., 2018). Future studies can overcome these limitations by creating richer spiking network models, with different operating frequencies and evolving neuronal dynamics. These models are difficult to create but would allow one to compare functional connectivity values with structural connection strength in more dynamic settings.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00206.

Ronaldo V. Nunes: Conceptualization; Investigation; Methodology; Project administration; Software; Validation; Visualization; Writing – original draft. Marcelo B. Reyes: Conceptualization; Methodology; Supervision; Writing – review & editing. Jorge F. Mejias: Conceptualization; Methodology; Supervision; Writing – review & editing. Raphael Y. de Camargo: Conceptualization; Methodology; Supervision; Writing – review & editing.

Ronaldo V. Nunes, Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (https://dx.doi.org/10.13039/501100002322), Award ID: Finance Code 001.

Spiking models:

A network model composed of spiking neurons (nodes) connected by synapses (edges). A spiking neuron is a simplified neuron model, which generates discrete spike events. It is also referred to as a spiking neuronal population model.

Local field potential (LFP):

A transient electrical activity in the extracellular medium resulting from ionic flows in multiple neurons.

Multivariate vector autoregressive (MVAR) model:

Autoregressive models permit predicting future values in time series from past values. Multivariate vector models extend them to work with multiple interdependent time series.

Power spectral density (PSD):

Description of the distribution of the power of a signal in terms of its frequencies.

Coefficient of variation (CV):

The ratio between the standard deviation and the mean. It measures the relative dispersion around the mean.