References

\pdfbookmark[1]TitleTitlePage

Fundamental activity constraints lead to specific interpretations of the connectome

Jannis Schuecker^{1,\Yinyang}, Maximilian Schmidt ^{1,\Yinyang},
Sacha J. van Albada^{1}, Markus Diesmann^{1,2,3}, and Moritz
Helias^{1,3}

Institute of Neuroscience and Medicine
(INM-6) and Institute for Advanced Simulation (IAS-6) and JARA BRAIN
Institute I, Jülich Research Centre, 52428 Jülich, Germany

Department of Psychiatry, Psychotherapy
and Psychosomatics, Medical Faculty, RWTH Aachen University, 52062
Aachen, Germany

Department of Physics, Faculty
1, RWTH Aachen University, 52062 Aachen, Germany

^{\Yinyang}
Both authors contributed equally
to this work.

Correspondence to: Jannis Schuecker, Wilhelm-Johnen-Straße, 52425 Jülich j.schuecker@fz-juelich.de

## Abstract

The continuous integration of experimental data into coherent models of the brain is an increasing challenge of modern neuroscience. Such models provide a bridge between structure and activity, and identify the mechanisms giving rise to experimental observations. Nevertheless, structurally realistic network models of spiking neurons are necessarily underconstrained even if experimental data on brain connectivity are incorporated to the best of our knowledge. Guided by physiological observations, any model must therefore explore the parameter ranges within the uncertainty of the data. Based on simulation results alone, however, the mechanisms underlying stable and physiologically realistic activity often remain obscure. We here employ a mean-field reduction of the dynamics, which allows us to include activity constraints into the process of model construction. We shape the phase space of a multi-scale network model of the vision-related areas of macaque cortex by systematically refining its connectivity. Fundamental constraints on the activity, i.e., prohibiting quiescence and requiring global stability, prove sufficient to obtain realistic layer- and area-specific activity. Only small adaptations of the structure are required, showing that the network operates close to an instability. The procedure identifies components of the network critical to its collective dynamics and creates hypotheses for structural data and future experiments. The method can be applied to networks involving any neuron model with a known gain function.

## Author summary

The connectome describes the wiring patterns of the neurons in the brain, which form the substrate guiding activity through the network. The influence of its constituents on the dynamics is a central topic in systems neuroscience. We here investigate the critical role of specific structural links between neuronal populations for the global stability of cortex and elucidate the relation between anatomical structure and experimentally observed activity. Our novel framework enables the evaluation of the rapidly growing body of connectivity data on the basis of fundamental constraints on brain activity and the combination of anatomical and physiological data to form a consistent picture of cortical networks.

## Introduction

The neural wiring diagram, the connectome, is gradually being filled through classical techniques combined with innovative quantitative analyses (Markov et al., 2014a, b) and new technologies (Wedeen et al., 2008; Axer et al., 2011). The connectivity between neurons is considered to shape resting-state and task-related collective activity (Cole et al., 2014; Shen et al., 2015). For simple networks, clear relationships with activity are known analytically, e.g., a dynamic balance between excitatory and inhibitory inputs in inhibition-dominated random networks leads to an asynchronous and irregular state (van Vreeswijk & Sompolinsky, 1996; Amit & Brunel, 1997a; Tetzlaff et al., 2012). Structures and mechanisms underlying large-scale interactions have been identified by means of dynamical models (Deco & Jirsa, 2012; Cabral et al., 2014) and graph theory (Markov et al., 2013; Goulas et al., 2014). Furthermore, the impact of local network structure on spike-time correlations is known in some detail (Ostojic et al., 2009; Pernice et al., 2011; Trousdale et al., 2012). Under certain conditions, there is a one-to-one mapping between correlations in neuronal network activity and effective connectivity, a measure that depends on the network structure and on its activity (van Albada et al., 2015). In conclusion, anatomy does not uniquely determine dynamics, but dynamical observations help constrain the underlying anatomy. Despite advances in understanding simple networks, a complete picture of the relationship between structure and dynamics remains elusive.

Fig. 1 visualizes the integrative loop between experiment, model, and theory that guides the investigation of structure and dynamics. In the traditional forward modeling approach, structural data from experimental studies determine the connectivity between model neurons. Combined with the specification of the single-neuron dynamics and synaptic interactions, simulations yield a certain network dynamics. There is a fundamental problem with this approach.

Despite the impressive experimental progress in determining network parameters, any neural network model is underdetermined, because of the sheer complexity of brain tissue and inevitable uncertainties in the data. For instance, counting synapses on a large scale presently takes a prohibitive amount of time, and no available technique allows determining precise synaptic weights for entire neural populations. Although it may be possible to quantify the full connectome of an individual in future, inter-individual variability will require modelers to express connectivity as connection probabilities or rules of self-organization if they want to learn about general principles of the brain. In modeling studies, parameters are usually tuned manually to achieve a satisfactory state of activity, which becomes unfeasible for high-dimensional models due to the size of the parameter space. In particular, it is a priori unclear how parameters of the model influence its activity. In consequence, modifications cannot be performed in a targeted fashion, and it is difficult to find a minimal set of modifications necessary for aligning a model with experimental activity data.

Overcoming this problem requires a shift of perspective. Instead of regarding the model exclusively in a forward manner, generating predicted activity from structure, we in addition consider the system in a reverse manner, predicting the structure necessary to explain the observed activity. Our theory, starting from a mean-field description, provides a direct link between structure and activity. In contrast to simulations, the theory is invertible, which we exploit to identify connections critical for the dynamics and to find a minimal set of modifications to the structure yielding a realistic set point of activity. The predicted alterations generate hypotheses on brain structure, thus feeding back to anatomical experiments.

Applying the method to a multi-scale network model of the vision-related areas of macaque cortex, we derive targeted modifications of a set of critical connections, bringing the model closer to experimental observations of cortical activity. Based on the global properties of the bistable phase space of the model, the method refines the model’s construction principles within experimental uncertainties and identifies the connections that decisively shape the dynamics. Preliminary results have been presented in abstract form (Schmidt et al., 2015).

## Results

### Global stability in a simple network

We consider networks with neurons structured into interconnected populations. A neuron in population receives incoming connections from neurons in population , each with synaptic efficacy . Additionally, each neuron is driven by Poisson sources with rate and synaptic efficacy . All neurons in one population have identical parameters, so that we can describe the network activity in terms of population-averaged firing rates .

Our method is applicable if the employed neuron model has a known gain function, either analytically or as a function obtained by interpolating numerical results from simulations. In this study, we model single cells as leaky integrate-and-fire model (LIF) neurons (see Methods). The possible stationary states of these networks are characterized by the firing rates that are equilibria of

(1) |

where is a pseudo-time. The gain function is known analytically and indicates its dependence on the model parameters (see Methods).

The input-output relationship typically features a non-linearity which, in combination with feedback connections, can cause multi-stability in the network. In particular excitatory-excitatory loops cause the system defined by (1) to exhibit multi-stable behavior in the stationary firing rates. A necessary condition for the bistability is that the transfer function has an inflection point. The LIF neuron model can have such an inflection point, originating from the interplay of its leak term and the threshold. Less realistic neuron models, such as the perfect integrate-and-fire model, do not have such an inflection point. To illustrate its origin, we first consider the noiseless case (Dayan & Abbott, 2001) without absolute refractoriness (). The transfer function initially grows from zero with infinite slope due to the threshold and crosses the identity line (Fig. 2C). For larger input the leak term can be neglected and the transfer function approaches a linear function with finite slope (see, e.g., (Kriener et al., 2014), eq. 11), equivalent to a perfect integrator. This is only possible with a negative curvature at intermediate rates, i.e., a reduction in the slope, which makes the transfer function cross the identity line another time, causing the bistability. Network-generated noise only affects the low-rate regime where it smears out the kink causing the transfer function to grow from zero with positive curvature (see, e.g., (Brunel, 2000), eq. 22). Importantly, the qualitative picture, i.e., the bistable behavior, is not altered. A finite refractory period only has an effect for very high input values where the transfer function saturates at for the given parameters.

The basic problem is that there is a trade-off between excitation at the fixed points and their stability. In particular, exciting the model to bring a fixed point closer to experimental observations requires a method to preserve its stability. We achieve this by controlling the influence of the excitatory-excitatory loops on the phase space of the network.

As an illustration we first study the mechanism using the simple network architecture depicted in Fig. 2A: a population of excitatory neurons is coupled to itself with indegree and is driven by external Poisson sources with the same indegree and rate . All connections have identical synaptic weights . An increase in the external drive shifts the input-output relationship of this one-dimensional system (Fig. 2C) to the left. The bifurcation diagram is shown in Fig. 2E: for low there is only one fixed point with low activity (LA). When increasing , an additional pair of fixed points of which one is stable and the other is unstable emerges via a saddle-node bifurcation, leading to a bistable system. The second stable fixed point exhibits high firing rates, denoted as the high-activity (HA) state. For even higher values of , the LA state loses stability.

The equilibria, given by the zeros of the velocity in the bistable case, are shown in Fig. 2F. An increase of the drive on the one hand increases the firing rate of the LA fixed point (inset) but on the other hand reduces its basin of attraction, indicated by the colored bars in the top left corner. For illustrative purposes, we extend the problem to two dimensions by splitting the excitatory population into two subpopulations of equal size (Fig. 2B), mimicking a loop between excitatory populations in the model of the vision-related areas of macaque cortex. The corresponding (symmetric) two-dimensional phase space is shown in Fig. 2D. The basin of attraction for the LA fixed point, limited by the separatrix (Strogatz, 1994), is reduced with increasing external drive.

Since we have a bistable system, there must be at least one unstable fixed point on the separatrix at the intersection of the nullclines, i.e., the subspace for which the velocity in direction vanishes. We use the unstable fixed point to preserve the basin of attraction when the external drive is increased. For this purpose, we modify the connectivity to reverse the shift of the unstable fixed point due to the parameter change (see Methods for a detailed derivation). Since the separatrix follows the unstable fixed point, this approximately restores the original basin of attraction.

The resulting velocity of the system ( defines the system (1)) is shown in Fig. 2F. The firing rate in the LA fixed point is increased as desired (inset), and the unstable fixed point coincides with that obtained in the original system . This pattern of fixed points is also indicated by the zero vectors of the velocity field (Fig. 2D). The separatrix follows the unstable fixed point, and the basin of attraction in the system ) is restored to that in the original system ). Fig. 2G shows the behavior of the LA fixed point in more detail. The modification of does not noticeably change the location of the LA fixed point. In conclusion, the method allows us to increase the firing rates in the LA fixed point without modifying its basin of attraction.

The purely excitatory network is the simplest model to explain a phase space configuration with a LA and a HA fixed point. Inhibitory feedback is not necessary for this bistability, but it would certainly alter the input-output relationship. For example the classical excitatory-inhibitory network (Brunel, 2000) in the balanced regime has an input-output relationship with a negative slope and thus only one fixed point exists. However, if a pair of such balanced E-I networks is coupled by sufficiently strong mutually excitatory connections, these connections cause a bistability in a similar manner. Thus the mechanism shown in the purely excitatory network can also lead to the emergence of a HA attractor in more complex networks with inhibition.

### Bistability in the multi-scale network model

We investigate a multi-scale network model of cortical areas to understand the structural features essential for a realistic state of baseline activity. The model extends and adapts the microcircuit model presented in (Potjans & Diesmann, 2014), which covers surface area of early sensory cortex (Fig. 3A), to all vision-related areas of macaque cortex (Fig. 3B). Based on the microcircuit model, an area is composed of 4 layers (2/3, 4, 5, and 6) each having an excitatory (E) and an inhibitory (I) neural population, except parahippocampal area TH, which consists of only 3 layers (2/3, 5, and 6). A detailed description of the data integration is given in (Schmidt et al., 2016).

Simulations of the model (Fig. 3C) reveal that, though realistic levels of activity can be achieved for populations in layers 2/3 and 4, populations 5E and 6E of the majority of areas show vanishingly low or zero activity in contrast to empirical data (Swadlow, 1988; de Kock & Sakmann, 2009). Inputs from subcortical and non-visual cortical areas are modeled as Poissonian spike trains, whose rate is a free, global parameter. To elevate the firing rates in the excitatory populations in layers 5 and 6, we enhance the external Poisson drive onto these populations parametrized by (see Methods). However, already a perturbation of a few percent leads to a state with unrealistically high rates (Fig. 3D), caused by the reduced basin of attraction of the low-activity state similar to Fig. 2D. Our aim is to improve the working point of the model such that all populations exhibit spiking activity while preventing the model from entering a state with unrealistically high rates of (figure 13 of (Swadlow, 1988), (de Kock & Sakmann, 2009)). The employed technique exposes the mechanism giving rise to the observed instability and identifies the circuitry responsible for this dynamical feature.

### Targeted modifications preserve global stability

We apply the procedure derived in Methods and find targeted modifications to the connectivity that preserve the global stability of the low-activity fixed point for increased values of the external drive .

In the following we choose the inactive state as the initial condition. The exact choice is not essential since we are only interested in the fixed points of the system. Fig. 4B shows the integration of (1) over pseudo-time for different levels of the external drive to populations 5E and 6E parametrized by . For low values of , the integration converges to the LA fixed point shown in Fig. 4D, and is in agreement with the activity emerging in the simulation (Fig. 3C). For increased values of , the system settles in the HA fixed point (Fig. 4E), again in agreement with the simulation. The population-specific firing rates in the HA state found in the mean-field predictions (Fig. 4E) are also close to those in the simulation (Fig. 3D), but minor deviations occur due to the violation of the assumptions made in the diffusion approximation. In particular, at these pathologically high rates, the neurons fire regularly and close to the reciprocal of their refractory period, while in the mean-field theory we assume Poisson spike statistics. Still, the mean-field theory predicts the bistability found in the simulation. Since the theory yields reliable predictions in both stable fixed points, we assume that also the location of the unstable fixed point in between these two extremes is accurately predicted by the theory.

To control the separatrix we need to find the unstable fixed point of the system. This is nontrivial since the numerical integration of (1) for finding equilibria by construction only converges to stable fixed points. If the unstable fixed point has only one repelling direction (Fig. 5A), it constitutes a stable attractor on the dimensional separatrix. The separatrix is a stable manifold (Strogatz, 1994), and therefore a trajectory originating in its vicinity but not near an unstable fixed point initially stays in the neighborhood. If an initial condition just outside the separatrix is close to the basin of attraction of a particular unstable fixed point, the trajectory initially approaches the latter. Close to the fixed point the velocity is small. Ultimately trajectories diverge from the separatrix in the fixed point’s unstable direction, as illustrated in Fig. 4A. In conclusion, we expect a local minimum in the velocity along the trajectories close to the unstable fixed point. To estimate the location of the unstable fixed point in this manner, we need to find initial conditions close to the separatrix. Naively, we would just fix the value of and vary the initial condition. However, due to the high dimensionality of our system this is not feasible in practice. Instead, we vary for a fixed initial condition. Fig. 4B shows the firing rate averaged across populations for two trajectories starting close to the separatrix, where the first one converges to the LA fixed point and the second one to the HA state. The trajectories diverge near the unstable fixed point and thus we define the last local minimum of the Euclidean norm of the velocity vector as the critical time at which we assume the system to be close to the unstable fixed point (Fig. 4C). We find four relevant and distinct unstable fixed points, of which two are shown in Fig. 6.

To counteract the shift of the separatrix caused by the increase in , we follow the procedure described in Methods. We subject the modifications of connectivity to the additional following constraints. In line with the anatomical literature, we do not allow for changes of the connectivity that would lead to cortico-cortical connections originating in the granular layer 4 (Felleman & Van Essen, 1991), and we also disallow inhibitory cortico-cortical connections, as the vast majority of long-range connections are known to be excitatory (Salin & Bullier, 1995; Tomioka & Rockland, 2007). In addition, we naturally restrict indegrees to positive values. We find that four iterations (numbered by index ) corresponding to the four distinct unstable fixed points suffice to preserve the basin of attraction of the LA state with respect to an increase of the external drive up to . In the following we concentrate on iterations 1 and 2, where the second one is also representative for iterations 3 and 4, which are qualitatively alike. To derive the required modifications of the indegree matrix, we decompose into its eigenmodes and quantify the contribution of each eigenmode to the shift of the unstable fixed point (see Methods). This allows us to identify the most effective eigendirection: in each iteration there is exactly one unstable eigendirection with an eigenvalue (Fig. 5A). The associated critical eigenvector is approximately anti-parallel to the shift of the fixed point, (inset of Fig. 5B), and of similar length. The critical eigendirection (red dot in Fig. 5B) constitutes the most effective modification, giving the largest contribution to the desired shift while requiring only a small change of in average total indegrees. In the chosen space of eigenmodes, the modifications are minimal in the sense that only this most effective eigenmode is changed.

The associated eigenvector predominately points into the direction of populations 4E and 5E of areas FEF and 46 (Fig. 5C), while has large entries in the 5E populations of two areas (Fig. 5D). The high rates of these populations at the unstable fixed points (cf. Fig. 6A,B with Fig. 5C,D) reflect that the instability is caused by increased rates in excitatory populations, particularly in population 5E. Each iteration shifts the transition to the HA state (the value of for which the separatrix crosses the initial condition) to higher values of and increases the attainable rates of populations 5E and 6E in the LA state (Fig. 7A). After all four iterations, the average total indegrees (summed over source populations) of the system are changed by . The first iteration mainly affects connections within and between areas 46 and FEF (Fig. 7B). In particular, the excitatory loops between the two areas are reduced in strength, especially those involving layer 5 (Fig. 7C). We thus identify two areas forming a critical loop. In the remaining iterations, the changes are spread across areas and especially connections originating in layer 5 are weakened (Fig. 7D). In conclusion, the method identifies critical structures in the model both on the level of areas and on the level of layers and populations, and leads to a small but specific structural change of the model.

### Analysis of the modifications

In the following we analyze the modifications of the connectivity with respect to the internal and inter-area connections in detail. The intrinsic circuits of the areas are modified in different directions, as shown for two exemplary areas V4 and CITv in Fig. 8A. Despite this heterogeneity, significant changes affect mostly excitatory-excitatory connections (Fig. 8A, bottom panel) with connections from population 5E experiencing the most significant changes (top panel of Fig. 8A). In fact, the anatomical data (Binzegger et al., 2004) underlying the microcircuit model (Potjans & Diesmann, 2014) contain only two reconstructed excitatory cells from layer 5, but considerably more for other cell types, indicating a higher uncertainty for layer 5 connections. Fig. 8B shows the correlation between intrinsic connectivity changes for all pairs of areas, with areas ordered according to a hierarchical clustering using a farthest point algorithm (Voorhees, 1986) on the correlation matrix. We find four clusters each indicating a group of areas which undergo changes with similar patterns. The groups are displayed in different colors in the histogram in Fig. 8B. The areas of the model are categorized into architectural types based on cell densities and laminar thicknesses (see Methods). Areas with architectural type 4, 5 and 6 are distributed over several clusters. We can interpret this as a differentiation of these types into further subtypes. The resulting changes of the intra-areal connectivity are small (Fig. 8A), but still significant for network stability.

Connections between areas can be characterized by their and (see Methods). The reflects the overall strength of an inter-areal connection and is only weakly affected across connections (Fig. 8C), with a correlation between original and modified logarithms of of . Significant variations in the occur mostly for very weak connections that are likely to have substantial relative uncertainties in the experimental data. The two overlapping red dots in Fig. 8C represent the connections between areas 46 and FEF, which are modified in the first iteration (Fig. 7C). The determines the laminar pattern of the location of source neurons for cortico-cortical connections. Overall, data are available for of the inter-areal connections in the parcellation of Felleman & van Essen (Felleman & Van Essen, 1991), while the for the rest are derived from the sigmoidal law. The majority of connections undergo small changes in their laminar source pattern (Fig. 8D) and connections with large modifications () are weak (average compared to in the model as a whole). Because weak connections are represented by low counts of labeled neurons, they have a relatively large uncertainty in their laminar patterns, justifying larger adjustments. Spearman’s rank correlation between the of the original model that were directly taken from experiments and the logarithmic ratios of cell densities is (, p-value of a two-sided test for uncorrelated data). For the modified model, we take the of all connections into account and obtain (, indicating a reduced, but still significant, monotonic dependence between and the logarithmic ratios of cell densities.

To judge the size of the modification to the connectivity, we compare it to the variability of measured cortico-cortical connection densities (Markov et al., 2014a). We quantify the latter as the average inter-individual standard deviation of the logarithmic , i.e., , where the overbar denotes the average over injections and the average over connections. This variability equals while the average modification of the logarithmic is . The main experimental connection probabilities used to construct the intra-areal connectivity of the model have an average relative standard deviation of across electrophysiological experiments (cf. Table 1 of (Potjans & Diesmann, 2014)) while the intra-areal connection probabilities of the model are modified by on average. The authors of (Scannell et al., 2000) report even greater variability in their review on cortico-cortical and thalamocortical connectivity. These considerations show that on average, the changes applied to the connectivity are well within the uncertainties of the data. Overall, 35 out of 603 connections were removed from the network. In the CoCoMac database, of these are indicated by only a single tracer injection, while the overall proportion of connections measured by a single injection is .

For the modified connectivity and , which we choose to avoid being too close to the transition (Fig. 7A), the theory predicts average rates in populations 5E and 6E of and , which is closer to experimentally observed rates compared to the original model. Furthermore we find that the modified connectivity allows us to decrease the inhibition in the network to . Simulating the full spiking network model then results in reasonable rates across populations and areas (Fig. 9B, D). The average rates in populations 5E and 6E are increased compared to a simulation of the original model from and to and , respectively. All populations exhibit firing rates within a reasonable range of to (Fig. 9D), as opposed to the original state in which a considerable fraction of excitatory neurons is silent (Fig. 3E). The theoretical prediction is in excellent agreement with the rates obtained in the simulation (Fig. 9A, C). Small discrepancies are caused by violations of the employed assumptions, i.e., Poissonian spiking statistics (Fourcaud & Brunel, 2002). Differences between theory and simulation are small, and negligible for the central aim of the study: the integration of activity constraints into the data-driven construction of multi-scale neuronal networks.

## Discussion

This study investigates the link between experimentally measured structural connectivity and neuronal activity in a multi-scale spiking network model of the vision-related areas of macaque cortex (Schmidt et al., 2016). We devise a theoretical method that systematically combines anatomical connectivity with physiological activity constraints. Already weak constraints, demanding the activity to neither vanish nor be pathologically high, yield a set of specific but small structural modifications necessary to increase the model’s excitation to a realistic level. We do not fit the model parameters to experimental activity data in the sense of minimizing an error function, since the sparse relevant experimental data do not allow for defining such a function in a meaningful way. Nevertheless, we considerably improve the network activity with a change in only a small set of parameters. The procedure constrains the experimentally obtained connectivity maps to a realization that is compatible with physiological experiments. This establishes a path from experimentally observed activity to specific hypotheses about the anatomy and demands for further experiments.

Connections are modified both inside and between cortical areas, on average well within the known uncertainties of the underlying data. The model areas are based on the early sensory cortex model presented in (Potjans & Diesmann, 2014). This circuit is adapted to individual areas by taking into account neuronal densities and laminar thicknesses. The model definition renders areas with equal architectural type similar in their internal connectivity, a drastic but inevitable simplification due to the lack of more detailed experimental data. The proposed method softens this assumption: small adaptations of internal connectivity distinguish the architectural types into further subtypes. These modifications are significant for the global stability of the network. Thus, our approach enables purely anatomy-based area categorizations to be refined with dynamical information.

Connections between areas are changed in terms of total strength and laminar patterns. Overall, the changes are small, but significant for specific connections. The loop formed by areas 46 and FEF is critical to the global stability of the network. Both areas have been investigated in (Markov et al., 2014a), albeit in a different parcellation scheme than the scheme used here (Felleman & Van Essen, 1991). Our method suggests a weaker coupling of these two areas than found in the anatomical data set. Uncertainties, partly due to the mapping between parcellations, leave room for this interpretation. Areas 46 and FEF belong to prefrontal cortex and are multimodal, indicating that the influence of other parts of cortex could stabilize them, a mechanism outside the scope of the present model of vision-related areas. Both explanations can be tested either in an experimental study or with an extended model.

A few weak connections undergo large changes in their laminar patterns. With the present activity constraints, the method hereby weakens hierarchical relations in the model structure, indicating that preserving these relations requires additional dynamical constraints such as layer-specific coherence between areas (Kerkoerle14_14332; Bastos et al., 2015). Conversely, it is possible to achieve satisfactory population rates with a less pronounced hierarchical structure.

Our analysis reveals that layer 5 excitatory cells play a critical role in the model’s dynamics, in line with the observed ongoing activity in mammalian neocortex (Sanchez-Vives & McCormick, 2000; Beltramo et al., 2013). This critical role is often attributed to single-neuron properties, with a subset of layer 5 neurons displaying pacemaker activity (Le Bon-Jego & Yuste, 2007; Lőrincz et al., 2015; Neske et al., 2015). In addition, we here find that the network architecture itself already explains the strong impact of layer 5 on the phase space of the network, suggesting that single-neuron properties and network structure jointly enable layer 5 to exert its dominant influence.

The cortico-cortical connectivity of the model is compiled from the extensive dataset of (Markov et al., 2014a) combined with the CoCoMac database (Stephan et al., 2001; Bakker et al., 2012), which collects data from hundreds of tracing studies. One could consider alternative methods for combining this information into one connectivity graph, for instance taking into account how consistently a given connection is reported across studies (Schmitt et al., 2014), and compare different methods by analyzing the resulting network dynamics. The presented mean-field theory could then be used to estimate the firing rates of each network instance without performing time-consuming simulations. However, here we first choose the integrative approach that accumulates all available experimental evidence into a single model and afterwards modify the resulting connectivity with our analytical procedure, thereby effectively discarding uncertain connections.

We restrict this study to networks of leaky integrate-and-fire model neurons, consistent with the key concept of the models we consider: individual cells are modeled in a simple manner to expose the impact of structural connectivity on the network dynamics. Moreover, the current-based leaky integrate-and-fire neuron can reproduce in-vivo like activity (Rauch et al., 2003; Jolivet et al., 2006) and is analytically tractable, which enables the identification of mechanisms underlying specific network effects. More complex neuron models can be incorporated into the method by replacing the gain function of each neuronal population with an analytical expression or an interpolated function obtained from spiking single-neuron simulations. For example, one could use a conductance-based point-neuron model for which the network dynamics can be described by population rate models (Shriki et al., 2003), featuring a non-monotonic gain function: the gain is reduced if excitatory and inhibitory inputs are increased in a balanced manner (Kuhn et al., 2004). Generally, this renders a system more stable. However, the bistability considered in our work is caused by excitatory inputs. Since conductance-based models also have a monotonically increasing gain function in dependence of the excitatory conductance alone, we expect the bistability to occur for such models as well.

Importantly, our method only requires a description of the system’s fixed points and their dependence on model parameters. The employed theory uses the diffusion approximation to derive a self-consistency equation for the stationary population rates valid for high indegrees and small synaptic efficacies. These requirements are fulfilled in the multi-area model and therefore the theoretical prediction agrees with the stationary activity of the simulation. Moreover, the theory predicts the bistability of the model, which is non-trivial, as the mean-field assumption of Poisson statistics of the activity is generally violated in the high-activity state. Nevertheless, since the activity in this state is mostly driven by strong mean inputs and the theory converges to the noiseless solution in the mean-driven limit, its predictions still provide viable approximations. The firing rates in the unstable fixed points are predominantly low, while the exceptions with very high rates are again mean-dominated. Presumably this explains the accuracy by which the theory predicts the locations of these fixed points and the resulting global stability properties.

Since the high-activity attractor of the model under consideration is unrealistic, we aim to prevent a transition to this state. However, in the high-dimensional system it is not a priori clear in which direction the separatrix has to be shifted to ensure stable dynamics. We therefore choose a pragmatic approach and shift the separatrix back to its initial location, inverting the shift which reduced the global stability of the low-activity fixed point. We achieve this to a good approximation by preserving the location of unstable fixed points on the separatrix. To this end, we use a linearization around these locations and an eigenmode decomposition to identify the set of connections to adjust. In the multi-area model, this linearization is justified because the system operates close to an instability so that only minor modifications are required. The method can be generalized to larger modifications by changing parameters iteratively in small steps.

Biological networks have various stabilization mechanisms not considered here, which render them less critical. For instance, during growth, homeostatic mechanisms guide the system toward the right structure. Furthermore, short-term synaptic plasticity (reviewed in (Morrison et al., 2008)), homeostatic synaptic scaling (Turrigiano et al., 1998) and spike-frequency adaptation (e.g. summarized in (Benda & Herz, 2003)) may prevent the system from entering the high-activity state. However, introducing these self-organizing mechanisms increases model complexity, causing a more intricate relation between structure and activity. Therefore, we start from a mean-field description on the level of neuronal populations, ignoring details of synaptic dynamics. Mild constraints on the activity lead to a network structure within the anatomical range of parameters. This network yields globally stable activity, suggesting that additional stabilization mechanisms are not required to achieve this. Nonetheless, they can potentially render the network more robust against external stimulation.

The observed inter-individual variability may reflect that mechanisms of self-organization and homeostasis find structurally different implementations of the same function. Thus, in studies across individuals we cannot expect that progress in experimental technology narrows down the variability of parameters indefinitely. The combined ranges of parameters rather specify the solution space, and our method provides a way to find a particular solution.

In principle, the method applies to any model parameter. It would be possible to modify, for example, synaptic weights. Since experimental data on synaptic weights are sparse, this is another natural choice. Moreover, such an analysis may provide hints about suitable synaptic plasticity rules that dynamically stabilize the model. The method can be applied to networks with more complex sets of attractors compared to the bistable case considered here. Though in high-dimensional systems such as our multi-area model, a larger number of attractors would make it more challenging to find all relevant unstable fixed points, the underlying idea of preserving the location of a separatrix is general. In contrast to the model considered here, transitions between fixed points can have a functional meaning in certain neuronal networks with multiple attractors. The specific location of the separatrix is then functionally relevant. Our method exposes the sensitivity of the location of the separatrix to certain model parameters and allows controlling its location in a specific manner.

In this work, we analyzed the global stability properties of the neuronal network on the population level. In contrast, Ostojic (Ostojic, 2014) performs a local stability analysis on the level of single neurons of an initially stable fixed point in a system with only one attractor. The author investigates the point at which the real parts of the eigenvalues of the Jacobian matrix evaluated at this fixed point become positive, i.e., the fixed point turns spectrally unstable and the system undergoes a transition to a heterogeneous asynchronous state. Analyzing the spectral stablilty on the single-neuron level does not reveal the global stability properties required in the current work: While a local stability analysis only considers infinitesimal perturbations, studying the basin of attraction gives information about the size of fluctuations against which the fixed point is stable. We expect both attractors to be spectrally stable because they do not show strong rate fluctuations and the mean-field theory predicts the activity accurately in both cases. Nevertheless, a heterogeneous state could occur if the synaptic weights were increased or the external drive was stronger. However, (Goedeke et al., 2016) show that the transition in stochastic systems that quantitatively resemble the spiking network (Grytskyy et al., 2013), does not coincide with the loss of spectral stability.

One striking feature of the heterogeneous state is bursty spiking behavior of individual cells. Bursty spiking is also observed in the multi-area model for increased synaptic weights of inter-area connections (Schmidt et al., 2016). The fixed points cannot be accurately described in this case because the fluctuations need to be taken into account (Mastroguiseppe & Ostojic, 2016). Simulations show (figures 4 and 5 of (Schmidt et al., 2016)) that the modifications obtained in this study are still able to prevent the system from a transition to a HA attractor also in the presence of bursting neurons. This indicates that the phase space does not change qualitatively and our results are robust against such bursting behavior.

Experimental data on stationary activity in cortex are sparse. We therefore restrict ourselves to fundamental constraints and increase the drive to the model in an area-unspecific way to fulfill them. The resulting heterogeneity of the firing rates across areas is thus not imposed by the method, but rather arises from the connectivity that remains strongly informed by anatomical data. Alternatively, one could predefine a desired state and investigate the parameter changes necessary to achieve it.

The presented analytical method that combines anatomy and activity data into a consistent model is restricted to stationary firing rates. In future studies, also higher-order statistical measures of activity can be used as constraints. Resting-state fMRI, for example, provides information on the functional connectivity between areas as a second-order measure. When combined with analytical predictions of functional connectivity, the method may shed light on the anatomical connection patterns underlying inter-area communication.

## Methods

In this study, we model single cells as leaky integrate-and-fire model neurons with exponentially decaying postsynaptic currents. Table 1 specifies the model parameters.

Synapse parameters | ||
---|---|---|

Name | Value | Description |

excitatory synaptic strength | ||

(Fig. 3) (Fig. 9) | relative inhibitory synaptic strength | |

local excitatory transmission delay | ||

local inhibitory transmission delay | ||

inter-areal transmission delay, with the distance between areas | ||

transmission speed |

Neuron model | ||
---|---|---|

Name | Value | Description |

membrane time constant | ||

absolute refractory period | ||

postsynaptic current time constant | ||

membrane capacity | ||

reset potential | ||

fixed firing threshold | ||

leak potential |

In the diffusion approximation, the dynamics of the membrane potential and synaptic current is (Fourcaud & Brunel, 2002)

(2) |

where is the membrane time constant and the synaptic time constant, respectively. The membrane resistance has been absorbed into the definition of the current. The input spike trains are approximated by a white noise current with fluctuations and mean value . Here is a centered Gaussian white process satisfying and . Whenever the membrane potential crosses the threshold the neuron emits a spike and is reset to the potential , where it is clamped during . All neurons in one population have identical parameters, so that we can describe the network activity in terms of population-averaged firing rates that depend on population-dependent input determined by the connectivity. Using the Fokker-Planck formalism, the stationary firing rates for each population are given by (Fourcaud & Brunel, 2002)

(3) | |||||

(4) | |||||

(5) |

which is correct up to linear order in and where with denoting the Riemann zeta function (Abramowitz & Stegun, 1974). Here, is chosen from the set of model parameters . If is a matrix, we vectorize it by concatenating its rows and indicate this by lower case, i.e., following (Magnus & Neudecker, 1995). If the chosen parameter is a scalar we denote it with .

We find the fixed points of (3) by solving the first-order differential equation (1) (Wong & Wang, 2006)

using the classical fourth-order Runge-Kutta method (RK4) with step size , where denotes a dimensionless pseudo-time. The same approach can be used to solve the activity on a single neuron level (Sadeh & Rotter, 2015). Note that (1) does not reflect the real time evolution of the population rates, but rather is a mathematical method to obtain the system’s fixed points. In contrast to (Wong & Wang, 2006) we do not only search for stable fixed points, but also use (1) to obtain unstable attractors (cf. Results), an idea originating from the study of simple attractor networks ((Amit & Brunel, 1997b) esp. their figure 2 and eq. 12).

In a bistable situation, the initial condition of (1) determines which fixed point the system settles in. However, studying the behavior for a particular initial condition is of minor interest, since the actual spiking network is a stochastic system which fluctuates around the fixed points of the deterministic system defined by (1). Even if we knew that (1) would relax to the LA fixed point for one particular initial condition, this would not necessarily imply that this state is indefinitely stable. Global stability is determined by the size of the basin of attraction of the LA fixed point.

In the following, we derive the equations leading us to targeted modifications of a parameter necessary to compensate for the changes in the global stability induced by the change of another parameter . To this end, we study the behavior of the fixed points with respect to an infinitesimal change in the chosen model parameter. Let and be the corresponding locations of the fixed points and their separation. We can then expand into a Taylor series up to first order in and obtain

(6) |

with

(7) | |||||

where we notice that and only depend on the target population . We accordingly define two diagonal matrices and with and . We further define the connectivity matrix , where denotes element-wise multiplication, also called the Hadamard product (see (Cichocki et al., 2009), for a consistent set of symbols for operations on matrices). The derivatives with respect to have the compact expressions

with the Jacobian of some vector-valued function and

where we use the Hadamard product again to define the matrix . Inserting the total derivatives into (7), we derive the final expression for , reading

(8) |

where we use for the identity matrix and define the effective connectivity matrix and the matrix . The latter has dimensionality , where is the dimension of (for example, for the vector of indegrees). With the aid of (8), evaluating (6) at the unstable fixed point predicts the shift of the separatrix (Fig. 2D) to linear order. We now consider an additional parameter which is modified to counteract the shift of the unstable fixed point caused by the change in parameter , i.e.,

(9) |

where we used that the inverse of appears on both sides of the equation and hence drops out. Note that the tuple may represent any combination of model parameters, for example external input, indegrees, synaptic weights, etc. For a particular choice of we solve (9) for . For the illustrative example shown in Fig. 2, where and , (9) simplifies to

since and appearing in the respective s cancel each other.

To determine critical connections in a more complex model, we choose , i.e. the vector of indegrees, and solve (9) with a decomposition into eigenmodes. We can write the right-hand side as

(10) |

The equation holds because are only affected by connections to population , and therefore their derivatives and hence , vanish for . We now make the ansatz

(11) |

which decomposes the changes into eigenmodes of the effective connectivity. The and are the -th right and left eigenvectors of as defined in (8), fulfilling the bi-orthogonality and the completeness relation Inserting (10) with (11) into (9) yields

(12) |

Thus we can solve for by multiplying from the left with

Our goal is to find a set of connections which dominate the size of the basin of attraction of the LA fixed point. Inserting (12) into (6) leads to

where are the eigenvalues of , which are either real or complex conjugate pairs since . To determine the influence of each eigenmode on the shift of the fixed point, we project the eigenvectors onto the fixed-point shift by multiplying each side with and solve again for to obtain

(13) |

where we define the (possibly complex-valued) coefficients We aim at a decomposition of into real components. If , is real, so we can work directly with . Complex eigenvalues and corresponding eigenvectors come in conjugate pairs, so in this case we combine the corresponding coefficients , to have all contributions and by construction (13). Each quantifies how much of the total fixed-point shift can be attributed to the -th eigenmode, which allows identification of the most effective eigendirection (see Results), where we apply the ansatz (11) to a multi-area, multi-layer model of the vision-related areas of macaque cortex.

The spiking simulations of the network model were carried out on the JUQUEEN supercomputer (Jülich Supercomputing Centre, 2015). All simulations were performed with NEST version 2.8.0 (Eppler et al., 2015) with optimizations for the use on the supercomputer which will be included in a future release. The simulations use a time step of and exact integration for the subthreshold dynamics of the leaky integrate-and-fire neuron model (reviewed in (Plesser & Diesmann, 2009)).

### Multi-area model

The multi-area model of the vision-related areas of macaque cortex uses the microcircuit model of (Potjans & Diesmann, 2014) as a prototype for all 32 areas in the FV91 parcellation (Felleman & Van Essen, 1991) and customizes it based on experimental findings on cortical structure. From anatomical studies, it is known that cortical areas in the macaque monkey are heterogeneous in their laminar structure and can be roughly categorized into 8 different architectural types based on cell densities and laminar thicknesses. This distinction was originally developed for prefrontal areas (Barbas & Rempel-Clower, 1997), and then extended to the entire cortex (Hilgetag et al., 2016). The visual cortex, and thus the model, comprises areas of categories 2, 4, 5, 6, 7 and 8. Precise layer-specific neuron densities are available for a number of areas, while for other areas, the neuron density is estimated based on their architectural type (see (Schmidt et al., 2016) for details).

The inter-areal connectivity is based on binary data from the CoCoMac database (Stephan et al., 2001; Bakker et al., 2012; Felleman & Van Essen, 1991; Rockland & Pandya, 1979; Barnes & Pandya, 1992) indicating the existence of connections, and quantitative data from (Markov et al., 2014a). The latter are retrograde tracing data where connection strengths are quantified by the fraction of labeled neurons () in each source area. The original analysis of the experimental data was performed in the M132 atlas (Markov et al., 2014a). Both the FV91 and the M132 parcellations have been registered to F99 space (Van Essen, 2002), a standard macaque cortical surface included with the software tool Caret (Van Essen et al., 2001). This enables mapping between the two parcellations.

On the target side, we use the exact coordinates of the injections to identify the equivalent area in the FV91 parcellation. To map the data on the source side from the M132 atlas to the FV91 parcellation, we count the number of overlapping triangles on the F99 surface between any given pair of regions and distribute data proportionally to the amount of overlap using the F99 region overlap tool on http://cocomac.g-node.org. In the model, this is mapped to the indegree the target area receives from source area divided by its total indegree, i.e., . Here, is defined as the total number of synapses between and divided by the total number of neurons in . On the source side, laminar connection patterns are based on CoCoMac (Felleman & Van Essen, 1991; Barnes & Pandya, 1992; Suzuki & Amaral, 1994; Morel & Bullier, 1990; Perkel et al., 1986; Seltzer & Pandya, 1994) and on fractions of labeled neurons in the supragranular layers () (Markov et al., 2014b). Gaps in these data are bridged exploiting a sigmoidal relation between and the logarithmized ratio of overall cell densities of the two areas, similar to (Beul et al., 2015). We map the to the ratio between the number of synapses originating in layer 2/3 and the total number of synapses between the two areas, assuming that only excitatory populations send inter-area connections, i.e., , where the indices and go over the different populations within area and , respectively. In the context of the model, we use the terms and to refer to the relevant relative indegrees given here. On the target side, the CoCoMac database provides data from anterograde tracing studies (Felleman & Van Essen, 1991; Rockland & Pandya, 1979; Jones et al., 1978; Seltzer & Pandya, 1991).

Missing inputs in the model, i.e., from subcortical and non-visual cortical areas, are replaced by Poissonian spike trains, whose rate is a free, global parameter. In the original model all populations of a particular area receive the same indegree of external inputs . The only exception to this rule is area TH where the absence of granular layer 4 is compensated by an increase of the external input to populations 2/3E and 5E by . To elevate the firing rates in the excitatory populations in layers 5 and 6, we increase the external drive onto these populations. The possibility of a higher drive onto these populations is left open by the sparseness of the corresponding experimental data. We enhance the external Poisson drive of the 5E population, parametrized by the incoming connections per target neuron (indegree), in all areas by increasing . The simultaneous increase in the drive of 6E needs to be stronger, since the firing rates in population 6E of the original model (Fig. 3C) are even lower than the rates of 5E (averaged across areas: for 5E compared to for 6E). We thus scale up linearly with such that results in .

## Acknowledgement

The authors gratefully acknowledge the computing time granted (jinb33) by the JARA-HPC Vergabegremium and provided on the JARA-HPC Partition part of the supercomputer JUQUEEN at Forschungszentrum Jülich. Partly supported by Helmholtz Portfolio Supercomputing and Modeling for the Human Brain (SMHB), the Helmholtz young investigator group VH-NG-1028, EU Grant 269921 (BrainScaleS). This project received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 720270. All network simulations carried out with NEST (http://www.nest-simulator.org).

## References

- Abramowitz & Stegun (1974) Abramowitz, M & Stegun, IA (1974). Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables. (New York: Dover Publications).
- Amit & Brunel (1997a) Amit, DJ & Brunel, N (1997a). Dynamics of a recurrent network of spiking neurons before and following learning. Network: Comput Neural Systems 8:373–404.
- Amit & Brunel (1997b) Amit, DJ & Brunel, N (1997b). Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex. Cereb Cortex 7:237–252.
- Axer et al. (2011) Axer, M, Grässel, D, Kleiner, M, Dammers, J, Dickscheid, T, Reckfort, J, Hütz, T, Eiben, B, Pietrzyk, U, Zilles, K, et al. (2011). High-resolution fiber tract reconstruction in the human brain by means of three-dimensional polarized light imaging. Front Neuroinformatics 5.
- Bakker et al. (2012) Bakker, R, Thomas, W, & Diesmann, M (2012). CoCoMac 2.0 and the future of tract-tracing databases. Front Neuroinformatics 6.
- Barbas & Rempel-Clower (1997) Barbas, H & Rempel-Clower, N (1997). Cortical structure predicts the pattern of corticocortical connections. Cereb Cortex 7:635–646.
- Barnes & Pandya (1992) Barnes, CL & Pandya, DN (1992). Efferent cortical connections of multimodal cortex of the superior temporal sulcus in the rhesus monkey. J Compar Neurol 318:222–244.
- Bastos et al. (2015) Bastos, AM, Vezoli, J, Bosman, CA, Schoffelen, JM, Oostenveld, R, Dowdall, JR, De Weerd, P, Kennedy, H, & Fries, P (2015). Visual areas exert feedforward and feedback influences through distinct frequency channels. Neuron 85:390–401.
- Beltramo et al. (2013) Beltramo, R, D’Urso, G, Maschio, MD, Farisello, P, Bovetti, S, Clovis, Y, Lassi, G, Tucci, V, Tonelli, DDP, & Fellin, T (2013). Layer-specific excitatory circuits differentially control recurrent network dynamics in the neocortex. Nat Neurosci 16:227–234.
- Benda & Herz (2003) Benda, J & Herz, V (2003). A universal model for spike-frequency adaptation. Neural Comput 15:2523–2546.
- Beul et al. (2015) Beul, SF, Barbas, H, & Hilgetag, CC (2015). A predictive structural model of the primate connectome. arXiv preprint arXiv:1511.07222.
- Binzegger et al. (2004) Binzegger, T, Douglas, RJ, & Martin, KAC (2004). A quantitative map of the circuit of cat primary visual cortex. J Neurosci 39:8441–8453.
- Brunel (2000) Brunel, N (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J Comput Neurosci 8:183–208.
- Cabral et al. (2014) Cabral, J, Kringelbach, ML, & Deco, G (2014). Exploring the network dynamics underlying brain activity during rest. Prog Neurobiol 114:102–131.
- Cichocki et al. (2009) Cichocki, A, Zdunek, R, Phan, AH, & Amari, Si (2009). Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation. (John Wiley & Sons).
- Cole et al. (2014) Cole, MW, Bassett, DS, Power, JD, Braver, TS, & Petersen, SE (2014). Intrinsic and task-evoked network architectures of the human brain. Neuron 83:238–251.
- Dayan & Abbott (2001) Dayan, P & Abbott, LF (2001). Theoretical Neuroscience. (Cambridge: MIT Press).
- de Kock & Sakmann (2009) de Kock, CPJ & Sakmann, B (2009). Spiking in primary somatosensory cortex during natural whisking in awake head-restrained rats is cell-type specific. Proc Natl Acad Sci USA 106:16446–16450.
- Deco & Jirsa (2012) Deco, G & Jirsa, VK (2012). Ongoing cortical activity at rest: Criticality, multistability, and ghost attractors. J Neurosci 32:3366–3375.
- Eppler et al. (2015) Eppler, JM, Pauli, R, Peyser, A, Ippen, T, Morrison, A, Senk, J, Schenck, W, Bos, H, Helias, M, Schmidt, M, et al. (2015). NEST 2.8.0.
- Felleman & Van Essen (1991) Felleman, DJ & Van Essen, DC (1991). Distributed hierarchical processing in the primate cerebral cortex. Cereb Cortex 1:1–47.
- Fourcaud & Brunel (2002) Fourcaud, N & Brunel, N (2002). Dynamics of the firing probability of noisy integrate-and-fire neurons. Neural Comput 14:2057–2110.
- Goedeke et al. (2016) Goedeke, S, Schuecker, J, & Helias, M (2016). Noise dynamically suppresses chaos in neural networks. arXiv. 1603.01880v1 [q-bio.NC].
- Goulas et al. (2014) Goulas, A, Bastiani, M, Bezgin, G, Uylings, HB, Roebroeck, A, & Stiers, P (2014). Comparative analysis of the macroscale structural connectivity in the macaque and human brain. PLoS Comput Biol 10:e1003529.
- Grytskyy et al. (2013) Grytskyy, D, Tetzlaff, T, Diesmann, M, & Helias, M (2013). A unified view on weakly correlated recurrent networks. Front Comput Neurosci 7:131.
- Hilgetag et al. (2016) Hilgetag, CC, Medalla, M, Beul, S, & Barbas, H (2016). The primate connectome in context: principles of connections of the cortical visual system. NeuroImage p. (in press).
- Jolivet et al. (2006) Jolivet, R, Rauch, A, Lüscher, HR, & Gerstner, W (2006). Predicting spike timing of neocortical pyramidal neurons by simple threshold models. J Comput Neurosci 21:35–49.
- Jones et al. (1978) Jones, E, Coulter, J, & Hendry, S (1978). Intracortical connectivity of architectonic fields in the somatic sensory, motor and parietal cortex of monkeys. J Compar Neurol 181:291–347.
- Jülich Supercomputing Centre (2015) Jülich Supercomputing Centre (2015). JUQUEEN: IBM Blue Gene/Q supercomputer system at the Jülich Supercomputing Centre. Journal of large-scale research facilities 1.
- Kriener et al. (2014) Kriener, B, Helias, M, Rotter, S, Diesmann, M, & Einevoll, G (2014). How pattern formation in ring networks of excitatory and inhibitory spiking neurons depends on the input current regime. Front Comput Neurosci 7:1–21.
- Kuhn et al. (2004) Kuhn, A, Aertsen, A, & Rotter, S (2004). Neuronal integration of synaptic input in the fluctuation-driven regime. J Neurosci 24:2345–2356.
- Le Bon-Jego & Yuste (2007) Le Bon-Jego, M & Yuste, R (2007). Persistently active, pacemaker-like neurons in neocortex. Front Neurosci 1:123–129.
- Lőrincz et al. (2015) Lőrincz, ML, Gunner, D, Bao, Y, Connelly, WM, Isaac, JT, Hughes, SW, & Crunelli, V (2015). A distinct class of slow ( 0.2–2 Hz) intrinsically bursting layer 5 pyramidal neurons determines UP/DOWN state dynamics in the neocortex. J Neurosci 35:5442–5458.
- Magnus & Neudecker (1995) Magnus, JR & Neudecker, H (1995). Matrix differential calculus with applications in statistics and econometrics. (John Wiley & Sons).
- Markov et al. (2013) Markov, NT, Ercsey-Ravasz, M, Van Essen, DC, Knoblauch, K, Toroczkai, Z, & Kennedy, H (2013). Cortical high-density counterstream architectures. Science 342.
- Markov et al. (2014a) Markov, NT, Ercsey-Ravasz, MM, Ribeiro Gomes, AR, Lamy, C, Magrou, L, Vezoli, J, Misery, P, Falchier, A, Quilodran, R, Gariel, MA, et al. (2014a). A weighted and directed interareal connectivity matrix for macaque cerebral cortex. Cereb Cortex 24:17–36.
- Markov et al. (2014b) Markov, NT, Vezoli, J, Chameau, P, Falchier, A, Quilodran, R, Huissoud, C, Lamy, C, Misery, P, Giroud, P, Ullman, S, et al. (2014b). Anatomy of hierarchy: Feedforward and feedback pathways in macaque visual cortex. J Compar Neurol 522:225–259.
- Mastroguiseppe & Ostojic (2016) Mastroguiseppe, F & Ostojic, S (2016). Intrinsically-generated fluctuating activity in excitatory-inhibitory networks. arXiv p. 1605.04221.
- Morel & Bullier (1990) Morel, A & Bullier, J (1990). Anatomical segregation of two cortical visual pathways in the macaque monkey. Visual neuroscience 4:555–578.
- Morrison et al. (2008) Morrison, A, Diesmann, M, & Gerstner, W (2008). Phenomenological models of synaptic plasticity based on spike-timing. Biol Cybern 98:459–478.
- Neske et al. (2015) Neske, GT, Patrick, SL, & Connors, BW (2015). Contributions of diverse excitatory and inhibitory neurons to recurrent network activity in cerebral cortex. J Neurosci 35:1089–1105.
- Ostojic (2014) Ostojic, S (2014). Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons. Nat Neurosci 17:594–600.
- Ostojic et al. (2009) Ostojic, S, Brunel, N, & Hakim, V (2009). How connectivity, background activity, and synaptic properties shape the cross-correlation between spike trains. J Neurosci 29:10234–10253.
- Perkel et al. (1986) Perkel, DJ, Bullier, J, & Kennedy, H (1986). Topography of the afferent connectivity of area 17 in the macaque monkey: A double-labelling study. J Compar Neurol 253:374–402.
- Pernice et al. (2011) Pernice, V, Staude, B, Cardanobile, S, & Rotter, S (2011). How structure determines correlations in neuronal networks. PLoS Comput Biol 7:e1002059.
- Plesser & Diesmann (2009) Plesser, HE & Diesmann, M (2009). Simplicity and efficiency of integrate-and-fire neuron models. Neural Comput 21:353–359.
- Potjans & Diesmann (2014) Potjans, TC & Diesmann, M (2014). The cell-type specific cortical microcircuit: Relating structure and activity in a full-scale spiking network model. Cereb Cortex 24:785–806.
- Rauch et al. (2003) Rauch, A, La Camera, G, Lüscher, H, Senn, W, & Fusi, S (2003). Neocortical pyramidal cells respond as integrate-and-fire neurons to in vivo like input currents. J Neurophysiol 90:1598–1612.
- Rockland & Pandya (1979) Rockland, KS & Pandya, DN (1979). Laminar origins and terminations of cortical connections of the occipital lobe in the rhesus monkey. Brain Res 179:3–20.
- Sadeh & Rotter (2015) Sadeh, S & Rotter, S (2015). Orientation selectivity in inhibition-dominated networks of spiking neurons: effect of single neuron properties and network dynamics. PLoS Comput Biol 11:e1004045.
- Salin & Bullier (1995) Salin, PA & Bullier, J (1995). Corticocortical connections in the visual system: structure and function. Physiol Rev 75:107–154.
- Sanchez-Vives & McCormick (2000) Sanchez-Vives, MV & McCormick, D (2000). Cellular and network mechanisms of rhythmic recurrent activity in neocortex. Nat Neurosci 3:1027–1034.
- Scannell et al. (2000) Scannell, JW, Grant, S, Payne, BR, & Baddeley, R (2000). On variability in the density of corticocortical and thalamocortical connections. Phil Trans R Soc B 355:21–35.
- Schmidt et al. (2016) Schmidt, M, Bakker, R, Diesmann, M, & van Albada, S (2016). Full-density multi-scale account of structure and dynamics of macaque visual cortex. arXiv preprint arXiv:1511.09364v4.
- Schmidt et al. (2015) Schmidt, M, Schuecker, J, Diesmann, M, & Helias, M (2015). Shaping phase space of neural networks via connectivity. In Proceedings of the 11th Göttingen Meeting of the German Neuroscience Society p. T26 7C.
- Schmitt et al. (2014) Schmitt, O, Eipert, P, Kettlitz, R, Leßmann, F, & Wree, A (2014). The connectome of the basal ganglia. Brain Structure and Function pp. 1–62.
- Seltzer & Pandya (1991) Seltzer, B & Pandya, DN (1991). Post-rolandic cortical projections of the superior temporal sulcus in the rhesus monkey. J Compar Neurol 312:625–640.
- Seltzer & Pandya (1994) Seltzer, B & Pandya, DN (1994). Parietal, temporal, and occipita projections to cortex of the superior temporal sulcus in the rhesus monkey: A retrograde tracer study. J Compar Neurol 343:445–463.
- Shen et al. (2015) Shen, K, Hutchison, RM, Bezgin, G, Everling, S, & McIntosh, AR (2015). Network structure shapes spontaneous functional connectivity dynamics. J Neurosci 35:5579–5588.
- Shriki et al. (2003) Shriki, O, Hansel, D, & Sompolinsky, H (2003). Rate models for conductance-based cortical neuronal networks. Neural Comput 15:1809–1841.
- Stephan et al. (2001) Stephan, K, Kamper, L, Bozkurt, A, Burns, G, Young, M, & Kötter, R (2001). Advanced database methodology for the collation of connectivity data on the macaque brain (CoCoMac). Phil Trans R Soc B 356:1159–1186.
- Strogatz (1994) Strogatz, SH (1994). Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering. (Reading, Massachusetts: Perseus Books).
- Suzuki & Amaral (1994) Suzuki, WL & Amaral, DG (1994). Perirhinal and parahippocampal cortices of the macaque monkey: cortical afferents. J Compar Neurol 350:497–533.
- Swadlow (1988) Swadlow, HA (1988). Efferent neurons and suspected interneurons in binocular visual cortex of the awake rabbit: Receptive fields and binocular properties. J Neurophysiol 59:1162–1187.
- Tetzlaff et al. (2012) Tetzlaff, T, Helias, M, Einevoll, G, & Diesmann, M (2012). Decorrelation of neural-network activity by inhibitory feedback. PLoS Comput Biol 8:e1002596.
- Tomioka & Rockland (2007) Tomioka, R & Rockland, KS (2007). Long-distance corticocortical GABAergic neurons in the adult monkey white and gray matter. J Compar Neurol 505:526–538.
- Trousdale et al. (2012) Trousdale, J, Hu, Y, Shea-Brown, E, & Josic, K (2012). Impact of network structure and cellular response on spike time correlations. PLoS Comput Biol 8:e1002408.
- Turrigiano et al. (1998) Turrigiano, GG, Leslie, KR, Desai, NS, Rutherford, LC, & Nelson, SB (1998). Activity-dependent scaling of quantal amplitude in neocortical neurons. Nature 391:892–896.
- van Albada et al. (2015) van Albada, SJ, Helias, M, & Diesmann, M (2015). Scalability of asynchronous networks is limited by one-to-one mapping between effective connectivity and correlations. PLoS Comput Biol 11:e1004490.
- Van Essen (2002) Van Essen, DC (2002). Windows on the brain: the emerging role of atlases and databases in neuroscience. Curr Opin Neurobiol 12:574–579.
- Van Essen et al. (2001) Van Essen, DC, Drury, HA, Dickson, J, Harwell, J, Hanlon, D, & Anderson, CH (2001). An integrated software suite for surface-based analyses of cerebral cortex. Journal of the American Medical Informatics Association 8:443–459.
- van Vreeswijk & Sompolinsky (1996) van Vreeswijk, C & Sompolinsky, H (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science 274:1724–1726.
- Voorhees (1986) Voorhees, EM (1986). Implementing agglomerative hierarchic clustering algorithms for use in document retrieval. Information Processing & Management 22:465–476.
- Wedeen et al. (2008) Wedeen, VJ, Wang, R, Schmahmann, JD, Benner, T, Tseng, W, Dai, G, Pandya, D, Hagmann, P, D’Arceuil, H, & de Crespigny, AJ (2008). Diffusion spectrum magnetic resonance imaging (dsi) tractography of crossing fibers. NeuroImage 41:1267–1277.
- Wong & Wang (2006) Wong, KF & Wang, XJ (2006). A recurrent network mechanism of time integration in perceptual decisions. J Neurosci 26:1314–1328.