Extracting Dynamical Understanding From Neural-Mass Models of Mouse Cortex

1. Introduction

Recent advances in neuroimaging have produced intricate maps revealing the complexity of the brain's microscale circuits, with whole-brain coverage. Analyzing and integrating these data have uncovered new patterns of brain organization, including the systematic spatial variation of gene expression (Burt et al., 2018; Fulcher et al., 2019), cytoarchitecture (Goulas et al., 2016), neuron densities (Erö et al., 2018), cortical thickness (Wagstyl et al., 2015), axonal connectivity (Oh et al., 2014), cognitive function (Margulies et al., 2016), and local dynamical properties (Shafiei et al., 2020). Existing evidence suggests that, to a good first approximation, these properties vary together along a dominant hierarchical axis in mouse and human (Burt et al., 2018; Fulcher et al., 2019; Wang, 2020).

To understand the functional role of observed physiological patterns, like systematic spatial variations in brain architecture, we need a way of simulating their effect on whole-brain dynamics. Physiologically based brain models achieve this, using methods from statistical physics to capture the dynamics of large populations of neurons and their interactions (Deco et al., 2008; Breakspear, 2017). Neural population models can capture the complex spatiotemporal dynamics in modern neuroimaging datasets, including persistent activity, intermittent oscillations, and multi-stability (Robinson et al., 2016; Noori et al., 2020; Froudist-Walsh et al., 2021; Mejías and Wang, 2022), and have successfully reproduced a wide range of experimental phenomena, from the alpha rhythm to seizure dynamics (Mejias et al., 2016; Breakspear, 2017; Schneider et al., 2021; Sip et al., 2022). The physiological formulation of these models means that their variables and parameters encode interpretable and biologically measurable properties of neural circuits, like the strengths and timescales of interactions between neuronal populations. This allows them to provide a unique mechanistic account of whole-brain dynamics that can be validated against both physiological experiments and the dynamical patterns observed in neuroimaging experiments.

While most existing brain models involve dynamical rules that are spatially uniform (e.g., the same model parameters in all brain areas), recent work has begun to investigate the effect of non-uniform dynamical rules, constrained by emerging brain-atlas datasets. An early example is the work of Chaudhuri et al. (2015), which incorporated a variation in recurrent excitation corresponding to that of measured spine count in the macaque. More recent work in human has incorporated spatial heterogeneity in model parameters with: the MRI-derived T1w:T2w map (Demirtas et al., 2019); T1w:T2w, the first principal component of gene transcription, and an inferred excitation:inhibition ratio (Deco et al., 2021); a linear combination of T1w:T2w and the principal resting-state functional connectivity (FC) gradient (Kong et al., 2021); a fitted parametric variation that recapitulated an interpretable hierarchical variation (Wang et al., 2019); and a spatial variation in excitability with a spatial map of epileptogenicity in modeling seizure dynamics and spread (Jirsa et al., 2017; Courtiol et al., 2020). These papers have reported improved out-of-sample model fits to empirical data, evaluated according to a range of summary statistics of the resulting dynamics (most typically FC), and provided insights into how spatial variation in biological mechanisms (like recurrent excitation) may underpin whole-brain dynamical regimes. While these studies demonstrate the promise of producing more accurate predictions of measured brain dynamics by incorporating regional heterogeneity—constraining to physiological data, or through large-scale parameter fitting (Wang et al., 2019)—the resulting models are correspondingly complex and challenging to interpret in terms of the mechanisms which underpin their dynamics. The tools of dynamical systems have the potential to reveal the dynamical features that improve model fits to data, including the bifurcation structure that defines the accessible dynamical regimes and the range of such regimes that different brain areas can access, including their vicinity to critical points (Deco and Jirsa, 2012; Deco et al., 2013; Cocchi et al., 2017; Demirtas et al., 2019; Wang et al., 2019). In this work, we show that analyzing the dynamical response of individual brain regions to inputs using bifurcation diagrams provides an understanding of model behavior in terms of accessible dynamical regimes, an approach that is particularly valuable for understanding the increased complexity of spatially non-uniform models.

The mouse is an ideal organism to develop comprehensively constrained physiologically based models of brain dynamics, but models of the mouse brain have been relative few compared to the large number of studies of human cortex. Existing models of mouse-brain dynamics on the macroscale have taken a variety of approaches, from phenomenological—connectome-coupled Kuramoto oscillators (Choi and Mihalas, 2019; Allegra Mascaro et al., 2020) and network diffusion models (Shadi et al., 2020)—through to neural mass models (Lin et al., 2020) coupled via a connectome (Melozzi et al., 2017, 2019) and interacting populations of spiking neural networks (Nunes et al., 2021). Compared to human, there is an abundance of high-resolution, whole-brain physiological data in mouse (Fulcher et al., 2019), including directed tract-tracing axonal connectivity data (Oh et al., 2014; Harris et al., 2019), high-resolution gene-expression maps (Lein et al., 2007), and cell-density atlases (Kim et al., 2017; Erö et al., 2018). High-quality whole-brain neuroimaging data using fMRI in mouse is also available, allowing us to evaluate model predictions in the resting state (Zerbi et al., 2015; Grandjean et al., 2020) and as a result of targeted manipulations (Zerbi et al., 2019; Markicevic et al., 2020, 2022). Prior work has shown that FC is strongly constrained by direct structural pathways (Grandjean et al., 2017), and prior dynamical models have reported the ability of coupled dynamical models to reproduce FC structure, especially when modeling using matching individual structural connectivity (Melozzi et al., 2019). In this work, we develop a neural-mass model of mouse cortical dynamics, and aim to understand the dynamical regimes in which it best captures resting-state fMRI data in mouse. We also aim to characterize the impact of incorporating spatial variations in excitatory and inhibitory cell densities as spatial variations in model parameters from a dynamical systems perspective.

2. Methods

As illustrated in Figure 1, we developed a neural mass model of the right hemisphere of the mouse cortex, across 37 cortical areas, comprising a simple Wilson–Cowan local dynamical model (Figure 1A) coupled via a directed structural connectome (Figure 1B). These regions are shown on the mouse brain in Figure 1C, colored by their relative excitatory cell densities (which are incorporated into the model in section 3.2). Of the 38 cortical regions reported in Oh et al. (2014), we excluded the frontal pole (FRP) due to its small size (likely contributing to noisy, outlying values of excitatory and inhibitory cell densities Erö et al., 2018). In visualizing our results, we grouped cortical regions according to six functional labels: Somatomotor, Medial, Temporal, Visual, Anterolateral, and Prefrontal (Harris et al., 2019) (see Supplementary Table S1 for full list).

FIGURE 1

Figure 1. Simulating and evaluating a coupled neural-mass model of mouse cortical dynamics. (A) The dynamics of individual brain regions follow the Wilson–Cowan equations (Wilson and Cowan, 1972, 1973) which govern interactions between local excitatory (E) and inhibitory (I) neural populations. (B) Regions are coupled together by connections defined by the AMBCA (Oh et al., 2014), represented as a directed adjacency matrix (connections shown black). A schematic shows how these long-range structural connections couple local cortical regions via excitatory projections (Breakspear, 2017). (C) Heterogeneity in local model parameters can be introduced as a perturbation that follows the measured variation in excitatory and inhibitory neural densities. Here the variation in excitatory cell density is plotted across the 37 mouse cortical areas as deviations relative to the mean level (green), using brainrender (Claudi et al., 2021) and data from Erö et al. (2018). (D) Model simulation yields activity time series for each brain region, from which pairwise linear correlations (functional connectivity, FC) are computed. (E) Model simulations are evaluated against empirical FC, averaged across 100 mice, as the Spearman correlation between all unique pairwise FC values, yielding an FC–FC score, ρ_FCFC.

As shown in Figure 1A, a given brain region consists of both an excitatory (E) and an inhibitory (I) neural population, whose dynamics are governed by the Wilson–Cowan equations (Wilson and Cowan, 1972, 1973). Brain regions are coupled via long-range excitatory projections using a binary, directed connectome from the Allen Mouse Brain Connectivity Atlas (AMBCA) (Oh et al., 2014; Fulcher and Fornito, 2016) (Figure 1B). As these data are the result of right-hemisphere viral tracer injections, yielding estimates of ipsilateral cortical connectivity in the right hemisphere, we modeled just the right hemisphere in this work, but note that model of both hemispheres could be developed in future under the assumption of lateral symmetry [e.g., as Melozzi et al. (2017)]. Simulating the model yields dynamics for the E and I populations; we take the activity time series of the excitatory population to evaluate the similarity of pairwise linear correlation structure as functional connectivity (FC), shown in Figure 1D. To assess the goodness of fit, we compare this simulated FC to an empirical FC calculated on a mouse fMRI dataset (Figure 1E). The goodness of fit is assessed as a Spearman correlation coefficient computed between all pairs of FC values from the empirical data and the model (Figure 1E). Spearman's correlation coefficient was used instead of Pearson's correlation coefficient to capture a potentially nonlinear but monotonic relationship.

As our main aim was to develop tools to understand the distributed dynamics of neural mass models, we favored simplicity in focusing on the Wilson–Cowan (W–C) model relative to alternative models. In addition, its physiological formulation is crucial for mapping to experimental cell-density data, as its parameters encode measurable properties with physical units that can be constrained by such data. The W–C model also exhibits a wide range of dynamical behaviors, including bifurcations, hysteresis, stable fixed-points (attractors), and limit cycles (oscillatory attractors) (Wilson and Cowan, 1972, 1973; Cowan et al., 2016), that are common features of dynamical systems in general, including more complex biophysical neural population models. We use a formulation of the Wilson–Cowan equations based on the mean firing rates of coupled populations of excitatory and inhibitory neurons, as

$\begin{array}{lll}{\tau }_eĖ& =-E+(1-E)S\left[a_e{w}_{ee}E-w_{ei}I-B_e+J_e\right],& (1)\end{array}$

$\begin{array}{lll}{\tau }_iİ& =-I+(1-I)S\left[a_i{w}_{ie}E-w_{ii}I-B_i\right],& (2)\end{array}$

where E and I are the mean firing rates of the excitatory and inhibitory populations, respectively (Hz); S(v) = h/[1+exp(−v)] is the sigmoidal firing-rate function; h (which is set to 1 here) is the upper bound for the sigmoid function representing a maximal population firing rate (Hz); a_e, a_i control the gradient scaling for the sigmoid function (V⁻¹); w_xy are the coupling weights from population y to population x, where x and y correspond to excitatory (e) or inhibitory (i) populations (V s); B_e, B_i are the firing thresholds for excitatory/inhibitory cells (V); J_e is the voltage induced by external current injected into the excitatory cells, defined below as a weighted sum over external inputs (V); and τ_e, τ_i are the time constants of excitation and inhibition respectively (s).

Neural masses, corresponding to cortical areas, were coupled via projections between excitatory populations through the external current term, J_e. For a given region a, $J_e^{(a)}(t)$ is computed as

$\begin{array}{ll}{J}_e^{(a)}(t)=G{\displaystyle \sum _b}{A}_{ab}{E}^{(b)}(t),& (3)\end{array}$

where G is a global coupling constant (V s), A_ab is the adjacency matrix corresponding to the structural connectome (unweighted here), and E^(b) is the excitatory activity of region b (Hz). It is helpful to define the quantity

$\begin{array}{ll}{J}_{\text{tot}}^{(a)}(t)=J_e^{(a)}(t)-B_e=G{\displaystyle \sum _b}{A}_{ab}{E}^{(b)}(t)-B_e,& (4)\end{array}$

as the total input that includes the constant offset B_e. We will use this to understand the dynamical response of a brain area to its net input in section 3.1.

In addition to this “homogeneous” model, in which the parameters are identical for all brain regions, we also analyze a heterogeneous model (in section 3.2), in which the coupling parameters, w_ij, vary across regions. We calibrate this variation to estimated cell-density data (Erö et al., 2018), by making the assumption that local connectivity from excitatory and inhibitory cells is uniform, and thus that coupling strengths from a given population are proportional to the density of cells of that population. Thus, we adjust the coupling parameters corresponding to outputs from the excitatory population, w_ee and w_ie, according to measured variations in excitatory cell density across cortical areas, and adjust w_ii and w_ei according to measured variations in inhibitory cell density. Defining nominal parameter values as ŵ_xy (for x and y taking i and e), we can then define linear parameter perturbations for a given region a as:

$\begin{array}{ll}\begin{array}{c}{w}_{ee}^{(a)}={ŵ}_{ee}(1+R_e^{(a)}),w_{ie}^{(a)}={ŵ}_{ie}(1+R_e^{(a)}),\\ w_{ii}^{(a)}={ŵ}_{ii}(1+R_i^{(a)}),w_{ei}^{(a)}={ŵ}_{ei}(1+R_i^{(a)}),\end{array}& (5)\end{array}$

where rescaling factors, R_e and R_i, represent relative variations in excitatory and inhibitory cell density, respectively (see Figure 1C for a visualization of how excitatory cell density varies across cortical areas). To map cell-density measurements to corresponding R_e and R_i values, we first z-score normalized raw excitatory and inhibitory cell-density data, as e^(a) and i^(a), respectively, across all regions, a. We then defined a simple proportional mapping to model parameters via a single scaling parameter, σ≥0, as

$\begin{array}{ll}{R}_e^{(a)}=\sigma e^{(a)},R_i^{(a)}=\sigma i^{(a)}.& (6)\end{array}$

In this formulation, setting σ = 0 sets all $R_e^{(a)}=R_i^{(a)}=0$ and reproduces the spatially homogeneous model; increasing σ increases the level of variation in coupling parameters across areas. Note that there is much scope for defining more complex mappings involving more new parameters, but defining the mapping from cell densities to model parameters in this simple, single-parameter scheme allows us to more clearly tackle our main aim to investigate how the model's dynamical features are shaped by such variation.

For a given system of coupled ODEs defined above, dynamics were simulated using The Virtual Brain (Sanz-Leon et al., 2015; Melozzi et al., 2017), yielding simulated time series for each region. The system was driven by white noise with a mean μ = 0 and standard deviation s = 1.3 × 10⁻⁵ using the Euler–Maruyama method with a fixed time step, Δt = 0.1 ms, for a total simulation length of 1.2 × 10⁵ ms, (or 2 min at 1,000 Hz). Initial transients of 1 s (1,000 time steps) were removed from all simulations to focus on the model's steady-state dynamics. As our aim was to understand the dynamical properties of the model that enable it to match the statistics of measured fMRI dynamics, we chose not to adjust the model output, E^(a)(t), through a simulation of the hemodynamic response function to match the fMRI measurement [but could be done in future using, e.g., a convolution of a canonical hemodynamic response (Boynton et al., 1996) or a biophysical model (Friston et al., 2000; Kim and Ress, 2016)].

fMRI data for 100 wild-type mice are taken from Zerbi et al. (2021), and consisted of blood-oxygen-level-dependent (BOLD) signals recorded from 100 anesthetized mice measured at rest for a period of 15 min using a Biospec 70/16 small animal MR system operating at 7T, equipped with a cryogenic quadrature surface coil for signal detection (Bruker BioSpin AG, Fällanden, Switzerland). The data were processed (see Supplementary Material for details) and parcellated using the Allen Common Coordinate Framework (CCF v3). Using time-series data from each of the 37 cortical regions analyzed here, we computed a functional connectivity (FC) matrix for each mouse as pairwise Pearson correlations. These matrices were averaged across mice to yield a group-average FC that was used as the basis of comparison for computing FC–FC scores.

Models were assessed on their ability to reproduce the pairwise linear correlation structure (FC) of empirical mouse fMRI data, as the Spearman correlation between predicted and measured FC values: the FC–FC score, ρ_FCFC. While we focused here on reproducing pairwise linear correlations using ρ_FCFC, we note that a more comprehensive evaluation of model fit, incorporating aspects of local dynamics and dynamic FC properties, will be important for future investigations to more fully evaluate the rich patterns contained in the dynamics (Cabral et al., 2017; Aquino et al., 2021; Deco et al., 2021). To account for variability in simulated model dynamics due to a finite simulation time and different random seeds, we computed ρ_FCFC for 40 repeats of each simulation using different random seeds. Code for reproducing the simulations and analysis presented here is available at https://github.com/DynamicsAndNeuralSystems/MouseBrainModelling.

3. Results

Here, we aim to understand the dynamical principles underlying coupled dynamical models using a neural mass model of the mouse cortex. First, in section 3.1, we investigate the spatially uniform case in which all brain regions are governed by identical dynamical rules. Focusing on model behavior in the vicinity of saddle-node and Hopf bifurcations, we characterize the model's dynamical regimes that best capture empirical FC structure. We then investigate the spatially heterogeneous case in section 3.2, in which regional variations in parameters are introduced according to variations in excitatory and inhibitory cell-density maps (Erö et al., 2018), which shape the model's local bifurcation properties and resulting dynamical regimes.

3.1. What Dynamical Features Drive High Model Performance?

We first characterize the model's dynamical regimes that best capture the pairwise correlation structure of experimental mouse fMRI, with the aim to understand how the positioning of individual nodes (brain regions) around specific types of bifurcations affects the model's ability to capture empirical FC. In this section we focus on a homogeneous model, in which all brain areas are governed by the same dynamical rules, but differ in their inputs from other regions (via the connectome). We characterize the model's behavior in each of three regimes: (i) in the vicinity of a single stable equilibrium, which we denote as the “Fixed Point” regime [using parameters adapted from Sanz-Leon et al. (2015)]; (ii) in the vicinity of a bistable region separated by saddle-node bifurcations, which we denote as the “Hysteresis” regime [using parameters from Heitmann et al. (2018)]; and (iii) in the vicinity of a pair of Hopf bifurcations, denoted as the “Limit Cycle” regime [using parameters from Borisyuk and Kirillov (1992)]. Parameter values for each of these three regimes are given in Supplementary Table S2. Bifurcation diagrams of excitatory firing, E, as a function of net external input, J_tot = J_e−B_e, are plotted for the Fixed Point regime (Figure 2D), Hysteresis regime (Figure 2E), and Limit Cycle regime (Figure 2F). These plots show how stable states of E vary with J_tot as solid lines (with unstable states shown for the bistable regime in Figure 2E and lower and upper limits of a limit-cycle oscillation in Figure 2F). They thus capture the rules underlying the dynamical behavior of individual brain regions in response to their aggregate input from other brain regions, J_tot, with each parameter setting defining a qualitatively different set of accessible dynamics, and types of response to inter-regional inputs. Importantly, these basic bifurcation structures, and the insights we gain from them, are not specific to the W–C model but are common features of many dynamical models (Strogatz, 2018).

FIGURE 2

Figure 2. Model performance is highly sensitive to the types of dynamical features available to the coupled dynamical network, with high FC–FC found near bifurcations and where external inputs have strong dynamical responses. (A–C) FC–FC score between model and data is plotted as a heat map in G–B_e space for the three model regimes considered here (see text): (A) “Fixed-point” regime, (B) “Hysteresis” regime, and (C) “Limit-cycle” regime. Corresponding J_tot–E bifurcation diagrams [cf. Equation (4)] for each regime are shown in the right-hand panels (D–F), showing stable E fixed points (solid), unstable E fixed points (dotted), and minima and maxima of limit-cycle oscillations (solid lines with shading). Dashed vertical lines represent the minimum J_tot corresponding to selected B_e values. Gray horizontal lines represent the range of J_tot values across regions and time for a sample simulation from the corresponding point in G–B_e space annotated in (A–C). Parameter values for each regime are in Supplementary Table S2.

We can understand the dynamics of an individual brain region in terms of the variation in its inputs over time, J_tot(t) (recalling that J_tot is high when a region has many inputs from other high-activity, or high-E, regions). To understand this in more detail, we consider two key parameters that control the range of J_tot that can be explored by a given brain region. As per Equation (4), these parameters are: (i) the excitatory firing threshold, B_e, which contributes a constant offset to J_tot; and (ii) the global coupling constant, G, which scales each region's response to excitatory inputs from other connected regions. Increasing B_e decreases J_tot for all cortical regions, shifting the range of J_tot explored by network nodes to the left on the J_tot–E bifurcation diagram. We can see this from the annotated levels of input, J_tot, corresponding to selected B_e values (when J_e = 0) as vertical dashed lines in Figures 2D–F, which denote the minimum J_tot for selected B_e values. Low G≈0 removes the effect of inter-regional coupling altogether (J_e≈0), resulting in a very narrow range of J_tot around B_e, while increasing G allows individual regions to respond more strongly to external inputs, and thus span a greater range of J_tot values. Given fixed values of B_e and G, the key factor controlling how different brain regions differ in J_tot in the homogeneous model is their connected neighbors, with high in-degree regions having more inputs and thus the potential to achieve a higher J_e and J_tot than low in-degree regions.

We are now able to analyze how different values of B_e (which adjusts the baseline of J_tot) and G (which scales the excitatory inputs, J_e, relative to this baseline) shape the dynamics of a given brain region. For example, some combinations of B_e and G confine all nodes in the network to a fixed-point attractor, whereas others allow some nodes to span one (or multiple) bifurcations. Different choices of G and B_e control the diversity of dynamical features supported by the model, but what types of configurations yield high FC–FC scores, ρ_FCFC? Our results, comparing across a range of both G and B_e, are shown as heat maps for the Fixed Point (Figure 2A), Hysteresis (Figure 2B), and Limit Cycle (Figure 2C) regimes. To visualize the correspondence between points in G–B_e space and the resulting range of J_tot(t) (and hence accessible dynamical regimes) they correspond to in the model simulation (range taken across time and nodes), we annotated this range in Figures 2D–F for key selected points in each corresponding heat map—labeled as “i,” “ii,” etc. For example, points toward the left of the G–B_e heat map correspond to low G and thus narrowing the range of J_tot, while points near the top of the heat map correspond to high B_e and hence low baseline inputs; hence points labeled “i” correspond to low and narrow ranges of J_tot, as annotated to the bifurcation diagrams in Figures 2D–F.

We first note a wide range of ρ_FCFC in all cases, indicating that model performance depends strongly on the local response to inputs, G and B_e, and thus the types of dynamical regimes available to the nodes of the coupled network. We also see that each dynamical regime exhibits characteristic regions of G–B_e space in which there is high FC–FC correspondence (colored red in Figures 2A–C), which reaches as high as ρ_FCFC = 0.52 (for the Fixed-Point regime), ρ_FCFC = 0.50 (Hysteresis regime), and ρ_FCFC = 0.56 (Limit-Cycle regime). All three model regimes can capture FC better than the direct correlation between SC and FC, ρ_SCFC = 0.42, indicating a benefit of accounting for distributed dynamics via coupled dynamical equations in capturing FC. Furthermore, model performance is consistent with, or higher than recently reported results for mouse cortex using a reduced Wong–Wang model (Wong and Wang, 2006) in a bistable regime [and using a Balloon–Windkessel BOLD filter (Friston et al., 2000) and a linear correlation, ρ_FCFC]: 0.35≾ρ_FCFC≾0.50 (Melozzi et al., 2019).

To understand how the model can produce high FC–FC, ρ_FCFC = 0.52±0.03, in the Fixed Point regime (Figure 2A), we start by exploring the qualitatively different types of input–output responses in Figure 2D. At high excitatory firing threshold, B_e, and low coupling, G (labeled “i” in Figures 2A,D), nodes can only access the relatively flat, low-E steady-state branch, weakening inter-regional communication across the brain and leading to poor FC–FC. A similar suppression of inter-regional communication, and resulting low ρ_FCFC, occurs when the model is confined to the upper branch at low B_e and high G (labeled “iii” in Figures 2A,D). In the intermediate region, labeled “ii” in Figures 2A,D, we obtain high FC–FC scores, up to a maximum ρ_FCFC = 0.52±0.03 (at B_e = 3.3 mV, G = 0.65 mVs). Here, brain areas can access the sharp gradient of the sigmoid-like stable branch in E, and are thus highly sensitive in their response to variations in the activity of neighboring brain regions. This gives us the somewhat surprising result that this very simple model, in a regime in which regions respond to the aggregate activity of their neighbors (but without any complex local dynamical features like bifurcations or oscillations) can produce high ρ_FCFC = 0.52, consistent with results reported recently using more complex models (Melozzi et al., 2019). This is qualitatively consistent with direct structural connections providing a strong constraint on the resulting FC (Grandjean et al., 2017), with non-direct interactions providing a more minor perturbation (Robinson, 2012).

We next investigated a “Saddle Node” model regime [using parameters from Borisyuk and Kirillov (1992)], that involves a pair of saddle-node bifurcations with an intermediate bistable region, shown in Figures 2B,E. We obtained qualitatively similar results to the Fixed-Point regime analyzed above: ρ_FCFC is low when nodes are confined to relatively flat low-E branch (at low G and high B_e, labeled “i” in Figures 2B,E) or the high-E branch (high G and low B_e, labeled “iii” in Figures 2B,E), where responses to external inputs are weak. Stronger FC–FC scores (e.g., a maximum ρ_FCFC = 0.50±0.14 at B_e = 3.7 mV and G = 0.35 mVs) again arise in the intermediate region, where the local activity response is most sensitive to driving inputs, J_tot (labeled “ii” in Figures 2B,E). The difference now is the increased diversity of supported dynamics: regions coexist between the stable low-E and high-E states, and can switch between them. This bistability leads to a greater dynamical repertoire of regions in the network, including longer-timescale switching (cf. Supplementary Figure S2), but this is not reflected in an improved ρ_FCFC.

Finally, we investigated model dynamics in the neighborhood of a stable limit cycle, separated by two Hopf bifurcations [model parameters from Heitmann et al. (2018)], shown as a J_tot–E bifurcation diagram in Figure 2F. As for the two regimes studied above, when nodes are confined to a relatively flat stable branch, labeled “i” and “v” in Figures 2C,F, FC–FC scores are low. For a similar reason, we also find low FC–FC when nodes are confined to a limit-cycle oscillation (labeled “iii” in Figures 2C,F), where nodes have a restricted ability to respond to their inputs in a way that their neighbors can meaningfully respond to [since nodes are coupled via E, cf. Equation (3)]. But the heat map in Figure 2C reveals two regions of G–B_e space with high ρ_FCFC, labeled “ii” and “iv.” In the region labeled “ii” (e.g., ρ_FCFC = 0.38±0.03 at B_e = 2.8 mV, G = 0.45 mVs), nodes sit on a stable branch which has a small but sufficient curvature to enable local activity to respond, albeit weakly, to inputs from connected regions. But the best fits to data, reaching ρ_FCFC = 0.56±0.04 (at B_e = 1.5 mV, G = 0.7 mVs), are found in the region labeled “iv” in Figures 2C,F. In this region of B_e–G space, nodes can access two distinctive types of dynamics: the limit-cycle regime (at low J_tot) and the high-E fixed-point attractor (at high J_tot). High FC–FC scores are also obtained when nodes can also access the low-E stable branch (at high G and high B_e). Importantly, the high-E branch at high J_tot has a relatively sharp dependence on J_tot, a feature that is common to obtaining high-ρ_FCFC scores in all three model regimes. Together, our analyses in this section demonstrate the importance of a model that allows nodes to respond sensitively to inputs from their network neighbors for reproducing FC.

3.1.1. Interpreting Simulated Dynamics in Terms of Bifurcation Diagrams

Bifurcation diagrams provide an understanding of the dynamical regimes accessed by individual nodes, and the way in which they respond to changes in inputs, information that can guide understanding of the complex distributed dynamics that result from a full model simulation. For the Limit-Cycle regime, simulated multivariate time series and corresponding FC matrices are plotted in Figure 3 for points labeled “ii,” “iii,” and “iv” in Figures 2C,F. In “ii,” nodes are confined to the low-E stable branch and, accordingly, the dynamics consist of noisy deviations from a low-E stable fixed point (Figure 3D). These perturbations can drive changes in structurally connected nodes, yielding weak pairwise correlations shown in Figure 3A. In “iii,” when nodes are mostly confined to the limit cycle and ρ_FCFC is low, most nodes exhibit oscillations (with some longer-timescale deviations, cf. Figure 3E), and a minority of other nodes (situated near the low-J_tot Hopf bifurcation) move between noisy deviations from the low-E stable branch and oscillatory limit-cycle dynamics. This results in very high pairwise correlations between groups of synchronized oscillatory nodes, r>0.8, such that the underlying structural connections play less of a role in shaping the pairwise correlation structure, resulting in a low FC–FC score. In “iv,” with the highest ρ_FCFC, the Hopf bifurcations facilitate complex spatiotemporal dynamics shown in Figure 3F. While many nodes spend most of the simulation near the high-E stable branch (those with high J_tot), we observe periods of time during which groups of nodes (near the Hopf bifurcation) display synchronized oscillations, embedded in globally complex and distributed dynamics on longer timescales. These analyses demonstrate how analyzing the response of local nodes to inputs, as ranges of J_tot in a bifurcation diagram (as in Figures 2D–F), can ground an understanding of the complex distributed dynamics that result from the full coupled model, which can be visualized effectively as heat maps (Figure 3).

FIGURE 3

Figure 3. Different dynamical features of the limit-cycle regime yield very different dynamics, including noisy deviations about a stable fixed point, synchronous oscillations, and a complex distributed dynamics featuring intermittent synchronization with high FC–FC. Here, we investigate simulated time series (lower row) and functional connectivity matrices (upper row) for three regions in B_e–G space annotated “ii,” “iii,” and “iv” in Figure 2. (A–C) Simulated functional connectivity matrices are plotted for “ii,” “iii,” and “iv,” respectively. (D–F) Simulated E time series are plotted as a node × time heat map (or “carpet plot” Aquino et al., 2020) for all brain regions for “ii,” “iii,” and “iv,” respectively. Colored bars label the six cortical divisions listed in Supplementary Table S1. In all plots, nodes are ordered as per Supplementary Table S1.

3.1.2. Resolving Inter-regional Differences in Inputs

The variation in qualitative dynamics across individual brain areas in the multivariate time series plotted in Figures 3D–F indicates that different network elements are accessing different dynamical regimes permitted by the model, resulting from substantial variability in the J_tot(t) experienced by different nodes. Since all nodes are governed by the same dynamical rules, and, hence, the same bifurcation diagrams, we can annotate J_tot(t) ranges onto a common bifurcation diagram to understand how the dynamics of individual regions are governed by different types of inputs from their connected neighbors. That is, rather than plotting just the overall range of J_tot (from the minimum to the maximum across all nodes), as in Figures 2D–F, we can resolve the individual ranges of J_tot experienced by each individual node on the J_tot–E bifurcation diagram. An example is shown for the Limit-Cycle regime at “iv” in Figure 4A [where we have plotted J_e instead of J_tot, equivalently, for a fixed B_e = 1.5 mV, cf. Equation (4)]. We see how, even with fixed dynamical rules, the range of J_e experienced by individual nodes varies markedly. Some regions have low J_e across the simulation, like the dorsal retrosplenial area, RSPd (annotated in Figure 4A), and, therefore, only display oscillations, as plotted in Figure 4B. Other regions with high J_e across the simulation, like the posterior parietal association areas, PTLp (annotated in Figure 4A), are confined to the stable high-E branch across the full simulation and display dynamics consistent with noisy deviations from a fixed point, as shown in Figure 4B. Regions like the ventral retrosplenial area, RSPv (annotated in Figure 4A), span the Hopf bifurcation, and thus exhibit more complex patterns that contain both oscillatory dynamics and noisy excursions about a stable fixed-point, depending on fluctuations in inputs, J_e(t). The short samples of E(t) for six annotated Medial regions in Figure 4B reveal some of these dynamics, including dynamic phase relationships between the oscillatory populations. These findings demonstrate the usefulness of interpreting the dynamics of coupled mass models in terms of time-varying inputs to the constituent populations.

FIGURE 4

Figure 4. Resolving different ranges of inputs, J_e, experienced by different network nodes allows us to understand their variable dynamical behavior in a coupled network model. Here we focus on the point labeled “iv” in the Limit-Cycle regime (Figures 2C,F), B_e = 1.5 mV, G = 0.7 mVs, in which nodes differ substantially in their inputs, J_e, and hence their resulting dynamics. (A) Bifurcation diagram for E and a function of J_e (as Figure 2F), with ranges of net excitatory drive, J_e, across the model simulation annotated for each brain region (colored according to the six labeled divisions). All regions are ordered according to Supplementary Table S1, and are labeled for the six Medial regions, which are plotted in (B). (B) E time series for the six Medial regions—PTLp, VISam, VISpm, RSPd, RSPv, and RSPagl—shown for the final 1 s of the simulation.

3.2. Understanding Heterogeneity in Local Dynamical Rules

Above, we used bifurcation diagrams to show that complex distributed dynamics in a neural-mass model can be understood in terms of the responses of individual regions to inputs from their connected neighbors. Despite equivalent local dynamical rules, and hence identical bifurcation diagrams for all brain regions, we found substantial inter-regional variability in accessible dynamical regimes and resulting activity dynamics, due to differences in structural connectivity and resulting J_e(t). In this section, we aim to understand the effect of varying the local dynamical rules themselves, by incorporating spatial heterogeneity in the properties of local cortical circuits (via a corresponding variation in model parameters). Specifically, we varied excitatory and inhibitory coupling strengths of individual brain areas according to excitatory and inhibitory cell-density data (Erö et al., 2018). We focused on the Limit Cycle regime of the W–C model described above, which displayed the richest dynamical repertoire and highest ρ_FCFC. As described in section 2, we used relative variations in excitatory and inhibitory cell densities across cortical areas to define a corresponding variation in R_e and R_i, which proportionally adjust coupling parameters—w_ii, w_ie, w_ei, w_ee—across brain areas. Setting R_e = R_i = 0 for all areas recovers the homogeneous model studied above [see Equation (5) for details]. This simple formulation allows us to understand how varying the excitatory and inhibitory coupling parameters across areas, in accordance with underlying excitatory and inhibitory cell densities, shape the dynamical responses of individual areas to inputs, and, hence, the resulting model dynamics.

3.2.1. Levels of Excitation and Inhibition Perturb Bifurcation Diagrams

To understand how variations in R_e and R_i affect model dynamics, we first analyze how these parameters shape the J_tot–E bifurcation diagrams for an individual area. The effect of ±10% variations to coupling parameters (corresponding to the ranges −0.1 < R_e < 0.1 and −0.1 < R_i < 0.1), are shown as J_tot–E bifurcation diagrams in Figure 5, varying just R_e (Figure 5A), just R_i (Figure 5B), R_e and R_i together with R_e = R_i (Figure 5C), and R_e and R_i such that R_e = −R_i (Figure 5D). We find that even these relatively small, ≈10%, perturbations have a substantial effect on the dynamical responses of individual areas, affecting: (i) the range of J_tot over which model exhibits stable oscillations; (ii) the oscillation amplitudes themselves; and (iii) steady-state activity levels. As shown in Figure 5, cortical areas with a higher excitatory cell density, R_e, have higher-amplitude oscillations, a wider range of J_tot over which stable oscillations are exhibited, and, for the same J_tot, increased activity, E, in the upper branch. Different changes result from modifying the inhibitory cell density, shown in Figure 5B: increasing R_i shifts the same bifurcation and fixed-point structure to higher J_tot (equivalent to raising the firing threshold, B_e). That is, regions with higher inhibitory cell density, R_i, require a greater aggregate input, J_e, to produce the same dynamics. Varying both R_e and R_i, shown in Figures 5C,D, yields combinations of the individual perturbations from R_e and R_i individually. These results demonstrate how relatively small variations in excitatory and inhibitory coupling parameters can have large effects on the bifurcation structure and dynamical regimes exhibited by local cortical regions. The effects are more dramatic for R_e and R_i values in the range from −0.5 to 0.5, where the Hopf bifurcations can be removed altogether from the Limit Cycle regime (Supplementary Figure S3), or additional stable states can be added via saddle-node bifurcations in the Hysteresis regime (Supplementary Figure S4).

FIGURE 5

Figure 5. Variations in excitatory and inhibitory cell density modify the dynamical regimes accessible to cortical regions. We model the effect of variations in excitatory and inhibitory cell density via perturbation parameters R_e and R_i, respectively, as defined in Equation (5). Relative to the nominal bifurcation diagram, R_e = R_i = 0 (black), we investigate variations in −0.1 ≤ R_e ≤ 0.1 and −0.1 ≤ R_i ≤ 0.1. Four types of variation were investigated: (A) R_e only (R_i = 0); (B) R_i only (R_e = 0); (C) R_e and R_i, such that R_e = R_i; and (D) R_e and R_i, such that R_e = −R_i. The legend indicates values of R_e.

3.2.2. Understanding Mouse Cortical Model Dynamics Constrained by Excitatory and Inhibitory Cell Densities

In the heterogeneous model, individual different brain areas differ both in their J_tot values that they receive from their coupled neighbors (due to differences in their structural connections), but also have different dynamical rules, due to different individual combinations of R_e and R_i values. As demonstrated above for the homogeneous model, this understanding of the dynamical responses of individual brain areas to inputs from across the network is crucial to guiding understanding of the complex, distributed dynamics of the full coupled model. In this section, we explore how the impact of local variations in R_e and R_i can be visualized and used to understand the dynamics of the full coupled model. Recall that our heterogeneous model is formulated with a single new parameter, σ, that defines how strongly relative differences in excitatory and inhibitory cell densities are mapped to corresponding changes in the model's coupling parameters (Equation 6). For the Limit-Cycle regime, we investigated how FC–FC scores change as we introduce a greater degree of inter-areal heterogeneity, σ. The variation in ρ_FCFC as a function of σ across the range 0 ≤ σ ≤ 1 is shown in Figure 6E. We did not find a substantial increase in ρ_FCFC when incorporating heterogeneity, σ>0, although there was a modest improvement relative to the homogeneous model (σ = 0) for σ = 0.2, yielding ρ_FCFC = 0.60±0.05. Testing this result against a null distribution (obtained by repeating the procedure but with randomly permuted excitatory and inhibitory cell-density data) using a permutation test yielded p≈0.15, indicating that ρ_FCFC = 0.60 does not constitute a significant improvement relative to the homogeneous model (see section 2 for details). As we discuss later, this result may be contributed to the small number of regions in the model, the simplicity of the dynamical equations, or the dominance of SC in constraining FC in the anesthetized mouse (Grandjean et al., 2017).

FIGURE 6

Figure 6. Modeling spatial variation in local excitatory and inhibitory cell densities produces complex distributed dynamics. (A) Bifurcation diagrams are plotted for all cortical areas according to their excitatory and inhibitory cell densities. Regions are colored according to their labeled anatomical grouping and the homogeneous case (R_e = R_i = 0) is shown in black for comparison. (B) The type of equilibrium dynamics displayed by a given cortical region, limit cycle (blue) or fixed point (red), is plotted as a function of J_e−B_e for all cortical regions for the range of J_e−B_e they experience across the model simulation. Nodes are ordered as per Supplementary Table S1 and shading reflects the six anatomical groupings labeled in A. Dashed lines shown at the top correspond to the uniform case (R_e = R_i = 0) for comparison. (C) Simulated time series for all brain areas are plotted as a heat map. Colors annotated to the right label the six anatomical groupings listed in A. (D) Simulated functional connectivity matrix. (E) FC–FC score as a function of the scaling parameter, σ, Equation (6). Results are shown for the model constrained by excitatory and inhibitory cell-density data (blue) and the permutation-based null distribution shown as mean ± standard deviation (red).

While we did not find evidence of a significant improvement in FC–FC score from a simple incorporation of heterogeneity, our main aim was to demonstrate how tools from dynamical systems can help to understand the complex coupled dynamics of such a spatially heterogenous model. We used the point, σ = 0.2, inferred above as a suitable example for this purpose. With σ = 0.2, we plotted J_e–E bifurcation diagrams for all brain regions on the same plot in Figure 6A. The plot shows how differences in excitatory and inhibitory cell densities results in different bifurcation diagrams, that correspond to similar qualitative changes as analyzed in Figure 5 above. Specifically, brain regions now differ substantially in their: critical values, J_e, that separate limit cycle from fixed-point dynamics; ranges of J_e in which oscillations are stable; oscillation amplitudes; and fixed-point activity levels, E, in the upper branch (for a given J_e). Compared to the homogeneous model, two regions with the same input, J_e, no longer indicates that they will be subject to the same dynamical rules.

To understand how these changes in local dynamical rules affect the resulting cortical dynamics, we next plotted the range of J_e that each node experiences across the simulation. As shown in Figure 6B, this can be represented as a horizontal line, distinguishing J_e values corresponding to what stable dynamical feature—“limit cycle” or “fixed point”—according to each region's individual bifurcation structure. This results in a richer dynamical landscape for the model: some brain regions can access both stable limit cycle and fixed-point dynamics, others can only access the high-E fixed-point equilibrium, while others can access just the limit-cycle attractor. It is useful to connect the range of dynamical regimes each region accesses across the simulation, shown in Figure 6B, with the E dynamics themselves, shown in Figure 6C. We can clearly see the high-E regions on the upper stable branch, as well as the more complex intermittent oscillations of regions that can access limit-cycle dynamics. The functional connectivity matrix from this simulation is shown in Figure 6D. This representation of pairwise correlations in the model dynamics hides much of the richness of the individual time series themselves (Figure 6C), and the dynamical rules that underlie them (Figures 6A,B). The ability to represent qualitative dynamical regimes of individual regions in a coupled network model—as J_e–E bifurcation diagrams with individual ranges of J_e explored for each region—provides a powerful illustration of the dynamics supported by the coupled components of a complex networked dynamical model.

4. Discussion

In this article, we developed a neural-mass model of the mouse cortex. We showed how bifurcation diagrams can be used to understand how regional differences in dynamics result from differences in inputs, J_tot, and delineated the types of dynamical regimes that yield good fits to experimental functional connectivity. We first analyzed a homogeneous model in which all regions are governed by identical dynamical rules to show how regional variations in dynamics result from differences in inputs (driven by differences in structural connectivity). We then extended this treatment to a heterogeneous model in which the bifurcation structures themselves vary across regions due to variation in local excitatory and inhibitory cell densities. Our results provide a useful framework for understanding the mechanisms that underlie complex simulated model dynamics, using a combination of local bifurcation diagrams (annotated with ranges of inputs for different regions) and visualizations of the multivariate time-series dynamics [as “carpet plots” (Aquino et al., 2020)]. These analyses will be particularly important for understanding how the brain's microscale circuits give rise to the complex distributed dynamics observed in brain-imaging experiments. A common scientific goal of modeling a system is to accurately reproduce important properties of it, while also gaining an understanding of how it does so. While successful approaches have been demonstrated for maximizing goodness of fit [sometimes optimizing large numbers of parameters (Wang et al., 2019; Kong et al., 2021; Wischnewski et al., 2021)], obtaining understanding is a key challenge for complex nonlinear models of brain dynamics. The analyses and visualizations demonstrated in this work aim to provide an understanding of the model dynamics in terms of the dynamical regimes that individual regions can access, shaped by their inputs from coupled neighbors. Key analyses include: (i) assessing the role of input parameters B_e and G in shaping empirical FC fits in terms of corresponding ranges spanned across J_e–E bifurcation diagrams (Figure 2); (ii) annotating of J_e for all cortical regions onto a common bifurcation diagram (Figure 4A); (iii) analyzing perturbations to bifurcation diagrams due to variations in local circuit properties (Figure 5); and (iv) annotating region-specific qualitative equilibrium dynamics across ranges of J_e for all regions in a single plot (Figure 6B). As whole-brain models develop to incorporate whole-brain datasets—including whole-brain maps of gene-expression and cell types (Fulcher et al., 2019; Yao et al., 2021)—these types of analyses will be crucial to understanding how this complexity shapes the underlying dynamical mechanisms, both at the level of individual brain regions, and their distributed interactions.

While many studies focus on determining an optimal working point, i.e., structural connectome scaling G, we find that the offset (B_e in the present model) is also critical in determining how local regions respond to inputs, and hence the resulting ρ_FCFC. We also found strong fits to empirical FC whenever brain regions were able to respond to inputs with sufficient gain, likely reflecting the strong role of direct structural connections in shaping FC in the anesthetized mouse (Grandjean et al., 2017). In particular, even in the Fixed-Point regime, in which the model exhibits the most constrained dynamical repertoire, we report ρ≈0.52, consistent with other published results in the literature [FC–FC scores up to ≈0.5 (Melozzi et al., 2019) using a Wong–Wang model (Wong and Wang, 2006)]. Only a small improvement was found when the model operated near a Hopf bifurcation, ρ_FCFC = 0.56. This highlights the ability of simple dynamical features to capture aspects of measured dynamics, consistent with prior comparisons demonstrating high performance of simple models (Messé et al., 2014, 2015; Nozari et al., 2021). The results also demonstrate the importance of comparing model performance against simpler benchmarks, and justifying increased model complexity only if it accompanies enhanced explanatory power.

Incorporating spatial variations in local dynamical rules according to whole-brain maps has immense potential in allowing us to connect new physiological understanding of neural circuits to the whole-brain dynamics that they enable. In this work, we incorporated spatial variations in excitatory and inhibitory cell densities as a corresponding perturbation to coupling parameters between E and I populations, with a single scaling parameter, σ. However, there are alternative ways in which this heterogeneity could be implemented and constrained in future, for example, by allowing σ to differ for excitatory (σ_e) and inhibitory (σ_i) populations. Incorporating more detailed physiological data into correspondingly more complex biophysical models (e.g., incorporating multiple inhibitory cell types), brings further parametric freedom that needs to be properly constrained from a combination of physiological and neuroimaging data. Our approach for assessing the improvement in ρ_FCFC after incorporating cell-density data involved a permutation approach against randomized assignment of the data to regions (preserving the match between e and i, but permuting their assignment to brain regions), and did not reveal a significant improvement relative to null gradients (p≈0.15). This may be due to the relatively small number of brain regions included, and the focus on FC–FC as an evaluation metric rather than a more comprehensive set of evaluations. Other ways of assessing the improvement of the spatially heterogeneous model could also be explored, such as testing the e:i ratio against alternative spatial gradients [as Deco et al. (2021)], or taking an optimization approach to estimate the optimal e and i gradients, and then assess their similarity to the measured excitatory and inhibitory cell-density data [as Wang et al. (2019)].

As our aim here was to demonstrate methods for understanding the dynamics of coupled neural models with heterogeneity using a simple modeling approach, many aspects of the model could be improved in future work. First, we have focused here on a specific simple biophysical model, the Wilson–Cowan model (Wilson and Cowan, 1972, 1973; Cowan et al., 2016), that allowed us to incorporate variations in excitatory and inhibitory cell-density data. We have focused on the behavior of the model in three specific dynamical regimes (a fixed point with gain, hysteresis, and limit cycle), but the results should be qualitatively applicable to those same dynamical regimes of other models. However, we note that other models with different dynamical features may display different behavior, such as the Wong–Wang model (Wong and Wang, 2006; Deco et al., 2013, 2021; Murray et al., 2017; Demirtas et al., 2019; Wang et al., 2019), or models that incorporate cortico–thalamic interactions (Wilson and Cowan, 1973; Robinson et al., 2015; Lin et al., 2020; Müller et al., 2020). We also note that while our aim here was to understand the model dynamics directly, it is common practice to simulate a hemodynamic response, such as the Balloon–Windkessel model (Friston et al., 2000) or a more sophisticated hemodynamic response function (Aquino et al., 2014). Simulating a slower hemodynamic response would introduce challenges in mapping bifurcation diagrams in E to the corresponding BOLD dynamics, and could lead to substantial qualitative differences between the dynamics of the neural model and the HRF-filtered dynamics. As a result, our specific conclusions about model performance in different dynamical regimes may not generalize to different choices of hemodynamic responses, but this could be achieved in future work by attempting to construct an effective bifurcation diagram of the dynamics of the BOLD forward solution as a function of the model parameters. We note, however, the body of evidence showing improved performance of linear models over nonlinear, biophysically informed models, in capturing the dynamical properties of fMRI data (Messé et al., 2014, 2015), and a recent finding that the performance of nonlinear neural mass models can drop when including HRF (Nozari et al., 2021). This suggests that, in the absence of thoroughly validated neural-mass models at the level of population neural activity (Lin et al., 2020), and a clearly demonstrated improvement of a BOLD forward model, neural mass models may be more conservatively viewed as a phenomenological means of capturing different types of dynamics and dynamical interactions, for which our simple approach, here, is valid and useful.

We also highlight our relatively simple treatment of structural connectivity, as a binary adjacency matrix, even though estimates of axonal connectivity strengths vary over at least four orders of magnitude (Oh et al., 2014). It remains an open question what greater structural connection “strengths” (approximated by the number of axons connecting two brain areas), corresponds to dynamically, e.g., faster connection speeds, a stronger effect on local population mean dynamics, or some alternative type of response. While the model here does not include time delays (assuming fast inter-regional interactions on the timescale of neural dynamics), they are likely to be crucial in shaping the brain's distributed dynamics (Petkoski and Jirsa, 2019) and should be explored in future work. We next note a major simplifying assumption in using a neural-mass model, which involves representing the spatially continuous cortical sheet as a set of 37 discrete cortical areas, abstracted away from their physical embedding (Robinson, 2019). Given the spatial resolution of modern mouse-brain maps, and the often continuous spatial variation they reveal, it will be important to develop models that accurately capture this physical continuity, e.g., using a neural field approach (Robinson et al., 2003).

Finally, we highlight the limitation of evaluating our model with respect to its ability to match only the linear correlation structure, FC, of the empirical dataset. fMRI data have a much richer dynamical structure than is captured by the static FC, including the dynamics of FC across a recording (Cabral et al., 2017; Demirtas et al., 2019; Aquino et al., 2021; Deco et al., 2021) and the organization of regional timescales (Sethi et al., 2017; Shafiei et al., 2020). For example, despite producing very different patterns in simulated time series, we found similar fits, ρ_FCFC, across the Fixed-Point, Hysteresis, and Limit-Cycle regimes of our homogeneous model, and when incorporating heterogeneity. The more complex distributed dynamics, including intermittent synchronization seen in carpet plots from the Limit Cycle regime (Figure 3F) and when incorporating regional heterogeneity (Figure 6C), qualitatively match the types of patterns seen in empirical fMRI dynamics better than in the fixed-point regime. This highlights the simplicity of the FC–FC score, ρ_FCFC, in capturing only the pairwise linear correlation structure in the data, and indicates the need for future work to perform a more comprehensive evaluation. This should include similar visualizations of model performance across B_e–G space (as Figures 2A–C), where the most distinctive models features for reproducing a greater range of characteristics of fMRI dynamics may be more clearly distinguished.

With the increasing availability of high-resolution neuroscience data, in space and time, the need for tools to provide interpretable accounts of their dynamics is pressing. Our work demonstrates a range of useful tools to analyze the behavior of coupled dynamical models of brain dynamics, helping them to provide understanding of the dynamical mechanisms that underpin their predictions. Our results emphasize the importance of benchmark comparison (e.g., a simple fixed-point model yields high FC–FC), visualization (e.g., very different dynamical patterns exhibited in carpet plots can yield similar correlation structures in FC), and proper statistical testing (e.g., while the heterogeneous model yields improved FC–FC, it is not significantly better than repeating the process on randomized data), practices that may help guide progress in the field.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics Statement

The animal study was reviewed and approved by the Ethical Committee of the Canton Zurich, Switzerland.

Author Contributions

BF and EM contributed to conception and design of the study. PS performed all simulations and analysis, supervised by BF and EM. Mouse fMRI data were processed by VZ. PS wrote an initial draft of the manuscript, which was refined for submission by BF and EM. All authors contributed to manuscript revision.

Funding

This work was supported by the Physics Foundation, School of Physics, The University of Sydney.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2022.847336/full#supplementary-material

References

Allegra Mascaro, A. L., Falotico, E., Petkoski, S., Pasquini, M., Vannucci, L., Tort-Colet, N., et al. (2020). Experimental and computational study on motor control and recovery after stroke: toward a constructive loop between experimental and virtual embodied neuroscience. Front. Syst. Neurosci. 14, 31. doi: 10.3389/fnsys.2020.00031

PubMed Abstract | CrossRef Full Text | Google Scholar

Aquino, K. M., Fulcher, B., Oldham, S., Parkes, L., Gollo, L., Deco, G., et al. (2021). On the intersection between data quality and dynamical modelling of large-scale fMRI signals. NeuroImage. 119051. doi: 10.1016/j.neuroimage.2022.119051

PubMed Abstract | CrossRef Full Text | Google Scholar

Aquino, K. M., Fulcher, B. D., Parkes, L., Sabaroedin, K., and Fornito, A. (2020). Identifying and removing widespread signal deflections from fMRI data: Rethinking the global signal regression problem. NeuroImage 212, 116614. doi: 10.1016/j.neuroimage.2020.116614

PubMed Abstract | CrossRef Full Text | Google Scholar

Aquino, K. M., Robinson, P. A., Schira, M. M., and Breakspear, M. J. (2014). Deconvolution of neural dynamics from fMRI data using a spatiotemporal hemodynamic response function. NeuroImage 94, 203–215. doi: 10.1016/j.neuroimage.2014.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Borisyuk, R. M., and Kirillov, A. B. (1992). Bifurcation analysis of a neural network model. Biol. Cybern. 66, 319–325.

PubMed Abstract | Google Scholar

Boynton, G. M., Engel, S. A., Glover, G. H., and Heeger, D. J. (1996). Linear systems analysis of functional magnetic resonance imaging in human V1. J. Neurosci. 16, 4207–4221.

PubMed Abstract | Google Scholar

Breakspear, M. J. (2017). Dynamic models of large-scale brain activity. Nat. Neurosci. 20, 340–352. doi: 10.1038/nn.4497

PubMed Abstract | CrossRef Full Text | Google Scholar

Burt, J. B., Demirtas, M., Eckner, W. J., Navejar, N. M., Ji, J. L., Martin, W. J., et al. (2018). Hierarchy of transcriptomic specialization across human cortex captured by structural neuroimaging topography. Nat. Neurosci. 27, 889. doi: 10.1038/s41593-018-0195-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Kringelbach, M. L., and Deco, G. (2017). Functional connectivity dynamically evolves on multiple time-scales over a static structural connectome: models and mechanisms. NeuroImage 160, 84–96. doi: 10.1016/j.neuroimage.2017.03.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaudhuri, R., Knoblauch, K., Gariel, M.-A., Kennedy, H., and Wang, X.-J. (2015). A large-scale circuit mechanism for hierarchical dynamical processing in the primate cortex. Neuron 88, 419–431. doi: 10.1016/j.neuron.2015.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, H., and Mihalas, S. (2019). Synchronization dependent on spatial structures of a mesoscopic whole-brain network. PLoS Comput. Biol. 15, e1006978. doi: 10.1371/journal.pcbi.1006978

PubMed Abstract | CrossRef Full Text | Google Scholar

Claudi, F., Tyson, A. L., Petrucco, L., Margrie, T. W., Portugues, R., and Branco, T. (2021). Visualizing anatomically registered data with brainrender. eLife 10, e65751. doi: 10.7554/eLife.65751

PubMed Abstract | CrossRef Full Text | Google Scholar

Cocchi, L., Gollo, L. L., Zalesky, A., and Breakspear, M. J. (2017). Criticality in the brain: a synthesis of neurobiology, models and cognition. Progr. Neurobiol. 158, 132–152. doi: 10.1016/j.pneurobio.2017.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Courtiol, J., Guye, M., Bartolomei, F., Petkoski, S., and Jirsa, V. K. (2020). Dynamical mechanisms of interictal resting-state functional connectivity in epilepsy. J. Neurosci. 40, 5572–5588. doi: 10.1523/JNEUROSCI.0905-19.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Cowan, J. D., Neuman, J., and van Drongelen, W. (2016). Wilson–cowan equations for neocortical dynamics. J. Math. Neurosci. 6, 1. doi: 10.1186/s13408-015-0034-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., and Jirsa, V. K. (2012). Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J. Neurosci. 32, 3366–3375. doi: 10.1523/JNEUROSCI.2523-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Kringelbach, M. L., Arnatkeviciute, A., Oldham, S., Sabaroedin, K., Rogasch, N. C., et al. (2021). Dynamical consequences of regional heterogeneity in the brain's transcriptional landscape. Sci. Adv. 7, eabf4752. doi: 10.1126/sciadv.abf4752

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Ponce-Alvarez, A., Mantini, D., Romani, G. L., Hagmann, P., and Corbetta, M. (2013). Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J. Neurosci. 33, 11239–11252. doi: 10.1523/JNEUROSCI.1091-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Robinson, P. A., Jirsa, V. K., Breakspear, M. J., and Friston, K. J. (2008). The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4, e1000092. doi: 10.1371/journal.pcbi.1000092

PubMed Abstract | CrossRef Full Text | Google Scholar

Demirtas, M., Burt, J. B., Helmer, M., Ji, J. L., Adkinson, B. D., Glasser, M. F., Van Essen, D. C., et al. (2019). Hierarchical heterogeneity across human cortex shapes large-scale neural dynamics. Neuron 101, 1181–1194.e13. doi: 10.1016/j.neuron.2019.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Erö, C., Gewaltig, M.-O., Keller, D., and Markram, H. (2018). A cell atlas for the mouse brain. Front. Neuroinf. 12, e17727. doi: 10.3389/fninf.2018.00084

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Mechelli, A., Turner, R., and Price, C. J. (2000). Nonlinear responses in fMRI: the balloon model, volterra kernels, and other hemodynamics. NeuroImage 42, 649–662. doi: 10.1016/j.neuroimage.2008.04.262

PubMed Abstract | CrossRef Full Text | Google Scholar

Froudist-Walsh, S., Bliss, D. P., Ding, X., Rapan, L., Niu, M., Knoblauch, K., et al. (2021). A dopamine gradient controls access to distributed working memory in the large-scale monkey cortex. Neuron 109, 3500–3520.e13. doi: 10.1016/j.neuron.2021.08.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Fulcher, B. D., and Fornito, A. (2016). A transcriptional signature of hub connectivity in the mouse connectome. Proc. Natl. Acad. Sci. U.S.A. 113, 1435–1440. doi: 10.1073/pnas.1513302113

PubMed Abstract | CrossRef Full Text | Google Scholar

Fulcher, B. D., Murray, J. D., Zerbi, V., and Wang, X.-J. (2019). Multimodal gradients across mouse cortex. Proc. Natl. Acad. Sci. U.S.A. 116, 4689–4695. doi: 10.1073/pnas.1814144116

PubMed Abstract | CrossRef Full Text | Google Scholar

Goulas, A., Uylings, H. B. M., and Hilgetag, C. C. (2016). Principles of ipsilateral and contralateral cortico-cortical connectivity in the mouse. Brain Struct. Funct. 252, 1–15. doi: 10.1007/s00429-016-1277-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandjean, J., Canella, C., Anckaerts, C., Ayrancı, G., Bougacha, S., Bienert, T., et al. (2020). Common functional networks in the mouse brain revealed by multi-centre resting-state fMRI analysis. NeuroImage 205, 116278. doi: 10.1016/j.neuroimage.2019.116278

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandjean, J., Zerbi, V., Balsters, J., Wenderoth, N., and Rudina, M. (2017). The structural basis of large-scale functional connectivity in the mouse. J. Neurosci. 37, 0438–17–8101. doi: 10.1523/JNEUROSCI.0438-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, J. A., Mihalas, S., Hirokawa, K. E., Whitesell, J. D., Choi, H., Bernard, A., et al. (2019). Hierarchical organization of cortical and thalamic connectivity. Nature 575, 195–202. doi: 10.1038/s41586-019-1716-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Heitmann, S., Aburn, M. J., and Breakspear, M. (2018). The brain dynamics toolbox for matlab. Neurocomputing 315, 82–88. doi: 10.1016/j.neucom.2018.06.026

CrossRef Full Text | Google Scholar

Jirsa, V. K., Proix, T., Perdikis, D., Woodman, M. M., Wang, H., Gonzalez-Martinez, J., et al. (2017). The virtual epileptic patient: individualized whole-brain models of epilepsy spread. NeuroImage 145, 377–388. doi: 10.1016/j.neuroimage.2016.04.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J. H., and Ress, D. (2016). Arterial impulse model for the BOLD response to brief neural activation. NeuroImage 124, 394–408. doi: 10.1016/j.neuroimage.2015.08.068

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, Y., Yang, G. R., Pradhan, K., Venkataraju, K. U., Bota, M., García del Molino, L. C., et al. (2017). Brain-wide maps reveal stereotyped cell-type-based cortical architecture and subcortical sexual dimorphism. Cell 171, 456–469.e22. doi: 10.1016/j.cell.2017.09.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, X., Kong, R., Orban, C., Wang, P., Zhang, S., Anderson, K., et al. (2021). Sensory-motor cortices shape functional connectivity dynamics in the human brain. Nat. Commun. 12, 6373. doi: 10.1038/s41467-021-26704-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Lein, E., Hawrylycz, M. J., Ao, N., Ayres, M., Bensinger, A., Bernard, A., et al. (2007). Genome-wide atlas of gene expression in the adult mouse brain. Nature 445, 168–176. doi: 10.1038/nature05453

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, I.-C., Okun, M., Carandini, M., and Harris, K. D. (2020). Equations governing dynamics of excitation and inhibition in the mouse corticothalamic network. bioRxiv Preprint. 2020.06.03.132688. doi: 10.1101/2020.06.03.132688

CrossRef Full Text | Google Scholar

Margulies, D. S., Ghosh, S. S., Goulas, A., Falkiewicz, M., Huntenburg, J. M., Langs, G., et al. (2016). Situating the default-mode network along a principal gradient of macroscale cortical organization. Proc. Natl. Acad. Sci. U.S.A. 113, 12574–12579. doi: 10.1073/pnas.1608282113

PubMed Abstract | CrossRef Full Text | Google Scholar

Markicevic, M., Fulcher, B. D., Lewis, C., Helmchen, F., Rudin, M., Zerbi, V., et al. (2020). Cortical excitation:inhibition imbalance causes abnormal brain network dynamics as observed in neurodevelopmental disorders. Cereb. Cortex 30, 4922–4937. doi: 10.1093/cercor/bhaa084

PubMed Abstract | CrossRef Full Text | Google Scholar

Markicevic, M., Sturman, O., Bohacek, J., Rudin, M., Zerbi, V., Fulcher, B. D., et al. (2022). Neuromodulation of striatal D1 cells shapes BOLD fluctuations in anatomically connected thalamic and cortical regions. bioRxiv Prepint 2022.03.11.483972. doi: 10.1101/2022.03.11.483972

CrossRef Full Text | Google Scholar

Mejias, J. F., Murray, J. D., Kennedy, H., and Wang, X.-J. (2016). Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex. Sci. Adv. 2, e1601335. doi: 10.1126/sciadv.1601335

PubMed Abstract | CrossRef Full Text | Google Scholar

Mejías, J. F., and Wang, X.-J. (2022). Mechanisms of distributed working memory in a large-scale network of macaque neocortex. eLife 11, e72136. doi: 10.7554/eLife.72136

PubMed Abstract | CrossRef Full Text | Google Scholar

Melozzi, F., Bergmann, E., Harris, J. A., Kahn, I., Jirsa, V., and Bernard, C. (2019). Individual structural features constrain the mouse functional connectome. Proc. Natl. Acad. Sci. U.S.A. 116, 26961–26969. doi: 10.1073/pnas.1906694116

PubMed Abstract | CrossRef Full Text | Google Scholar

Melozzi, F., Woodman, M. M., Jirsa, V. K., and Bernard, C. (2017). The virtual mouse brain: a computational neuroinformatics platform to study whole mouse brain dynamics. eNeuro 4, ENEURO.0111–17.2017. doi: 10.1523/ENEURO.0111-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Messé, A., Rudrauf, D., Benali, H., and Marrelec, G. (2014). Relating structure and function in the human brain: relative contributions of anatomy, stationary dynamics, and non-stationarities. PLoS Comput. Biol. 10, e1003530. doi: 10.1371/journal.pcbi.1003530

PubMed Abstract | CrossRef Full Text | Google Scholar

Messé, A., Rudrauf, D., Giron, A., and Marrelec, G. (2015). Predicting functional connectivity from structural connectivity via computational models using MRI: an extensive comparison study. NeuroImage 111, 65–75. doi: 10.1016/j.neuroimage.2015.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, E. J., Munn, B. R., and Shine, J. M. (2020). Diffuse neural coupling mediates complex network dynamics through the formation of quasi-critical brain states. Nat. Commun. 11, 6337. doi: 10.1038/s41467-020-19716-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Murray, J. D., Jaramillo, J., and Wang, X.-J. (2017). Working memory and decision-making in a frontoparietal circuit model. J. Neurosci. 37, 12167–12186. doi: 10.1523/JNEUROSCI.0343-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Noori, R., Park, D., Griffiths, J. D., Bells, S., Frankland, P. W., Mabbott, D., et al. (2020). Activity-dependent myelination: a glial mechanism of oscillatory self-organization in large-scale brain networks. Proc. Natl. Acad. Sci. U.S.A. 117, 13227–13237. doi: 10.1073/pnas.1916646117

PubMed Abstract | CrossRef Full Text | Google Scholar

Nozari, E., Bertolero, M. A., Stiso, J., Caciagli, L., Cornblath, E. J., He, X., et al. (2021). Is the brain macroscopically linear? A system identification of resting state dynamics. arXiv [Preprint]. arXiv:2012.12351. doi: 10.48550/arXiv.2012.12351

CrossRef Full Text | Google Scholar

Nunes, R. V., Reyes, M. B., Mejias, J. F., and de Camargo, R. Y. (2021). Directed functional and structural connectivity in a large-scale model for the mouse cortex. Netw. Neurosci. 5, 874–889. doi: 10.1162/netna00206

PubMed Abstract | CrossRef Full Text | Google Scholar

Oh, S. W., Harris, J. A., Ng, L., Winslow, B., Cain, N., Mihalas, S., et al. (2014). A mesoscale connectome of the mouse brain. Nature 508, 207–214. doi: 10.1038/nature13186

PubMed Abstract | CrossRef Full Text | Google Scholar

Petkoski, S., and Jirsa, V. K. (2019). Transmission time delays organize the brain network synchronization. Philosoph. Trans. R. Soc. Math. Phys. Eng. Sci. 377, 20180132. doi: 10.1098/rsta.2018.0132

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A. (2012). Interrelating anatomical, effective, and functional brain connectivity using propagators and neural field theory. Phys. Rev. E 85, 011912. doi: 10.1103/PhysRevE.85.011912

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A. (2019). Physical brain connectomics. Phys. Rev. E 99, 012421. doi: 10.1103/PhysRevE.99.012421

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A., Postnova, S., Abeysuriya, R. G., Kim, J. W., Roberts, J. A., McKenzie-Sell, L., et al. (2015). “A multiscale working brain model,” in Validating Neuro-Computational Models of Neurological and Psychiatric Disorders. (Cham: Springer), 107–140.

Robinson, P. A., Rennie, C. J., Rowe, D. L., O'Connor, S. C., Wright, J. J., Gordon, E., et al. (2003). Neurophysical modeling of brain dynamics. Neuropsychopharmacology 28, S74–S79. doi: 10.1038/sj.npp.1300143

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A., Zhao, X., Aquino, K. M., Griffiths, J. D., Sarkar, S., and Mehta-Pandejee, G. (2016). Eigenmodes of brain activity: neural field theory predictions and comparison with experiment. NeuroImage 142, 79. doi: 10.1016/j.neuroimage.2016.04.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanz-Leon, P., Knock, S. A., Spiegler, A., and Jirsa, V. K. (2015). Mathematical framework for large-scale brain network modeling in The Virtual Brain. NeuroImage 111, 385–430. doi: 10.1016/j.neuroimage.2015.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider, M., Broggini, A. C., Dann, B., Tzanou, A., Uran, C., Sheshadri, S., et al. (2021). A mechanism for inter-areal coherence through communication based on connectivity and oscillatory power. Neuron 109, 4050–4067.e12. doi: 10.1016/j.neuron.2021.09.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Sethi, S. S., Zerbi, V., Wenderoth, N., Fornito, A., and Fulcher, B. D. (2017). Structural connectome topology relates to regional BOLD signal dynamics in the mouse brain. Chaos Interdiscipl. J. Nonlin. Sci. 27, 047405. doi: 10.1063/1.4979281

PubMed Abstract | CrossRef Full Text | Google Scholar

Shadi, K., Dyer, E., and Dovrolis, C. (2020). Multisensory integration in the mouse cortical connectome using a network diffusion model. Netw. Neurosci. 4, 1030–1054. doi: 10.1162/netna00164

PubMed Abstract | CrossRef Full Text | Google Scholar

Shafiei, G., Markello, R. D., Vos de Wael, R., Bernhardt, B. C., Fulcher, B. D., and Misic, B. (2020). Topographic gradients of intrinsic dynamics across neocortex. eLife 9, e62116. doi: 10.7554/eLife.62116

PubMed Abstract | CrossRef Full Text | Google Scholar

Sip, V., Guye, M., Bartolomei, F., and Jirsa, V. (2022). Computational modeling of seizure spread on a cortical surface. J. Comput. Neurosci. 50, 17–31. doi: 10.1007/s10827-021-00802-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Strogatz, S. H. (2018). Nonlinear Dynamics and Chaos with Student Solutions Manual: With Applications to Physics, Biology, Chemistry, and Engineering. CRC Press.

Wagstyl, K., Ronan, L., Goodyer, I. M., and Fletcher, P. C. (2015). Cortical thickness gradients in structural hierarchies. NeuroImage 111, 241–250. doi: 10.1016/j.neuroimage.2015.02.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, P., Kong, R., Kong, X., Liégeois, R., Orban, C., Deco, G., et al. (2019). Inversion of a large-scale circuit model reveals a cortical hierarchy in the dynamic resting human brain. Sci. Adv. 5, eaat7854. doi: 10.1126/sciadv.aat7854

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X.-J. (2020). Macroscopic gradients of synaptic excitation and inhibition in the neocortex. Nat. Rev. Neurosci. 21, 169–178. doi: 10.1038/s41583-020-0262-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24.

PubMed Abstract | Google Scholar

Wilson, H. R., and Cowan, J. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik 13, 55–80.

PubMed Abstract | Google Scholar

Wischnewski, K. J., Eickhoff, S. B., Jirsa, V. K., and Popovych, O. V. (2021). Towards an efficient validation of dynamical whole-brain models. Sci Rep. 12, 4331. doi: 10.1038/s41598-022-07860-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, K.-F., and Wang, X.-J. (2006). A recurrent network mechanism of time integration in perceptual decisions. J. Neurosci. 26, 1314–1328. doi: 10.1523/JNEUROSCI.3733-05.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, Z., van Velthoven, C. T. J., Nguyen, T. N., Goldy, J., Sedeno-Cortes, A. E., Baftizadeh, F., et al. (2021). A taxonomy of transcriptomic cell types across the isocortex and hippocampal formation. Cell 184, 3222–3241.e26. doi: 10.1016/j.cell.2021.04.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Zerbi, V., Floriou-Servou, A., Markicevic, M., Vermeiren, Y., Sturman, O., Privitera, M., et al. (2019). Rapid reconfiguration of the functional connectome after chemogenetic locus coeruleus activation. Neuron 103, 702–718.e5. doi: 10.1016/j.neuron.2019.05.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Zerbi, V., Grandjean, J., Rudin, M., and Wenderoth, N. (2015). Mapping the mouse brain with rs-fMRI: an optimized pipeline for functional network identification. NeuroImage 123, 11–21. doi: 10.1016/j.neuroimage.2015.07.090

PubMed Abstract | CrossRef Full Text | Google Scholar

Zerbi, V., Pagani, M., Markicevic, M., Matteoli, M., Pozzi, D., Fagiolini, M., et al. (2021). Brain mapping across 16 autism mouse models reveals a spectrum of functional connectivity subtypes. Mol. Psychiatry 26, 1–11. doi: 10.1038/s41380-021-01245-4

PubMed Abstract | CrossRef Full Text | Google Scholar