Fundamental activity constraints lead to specific interpretations of the connectome

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 primate visual 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 minimal 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.


Introduction
The neural wiring diagram, the connectome, is gradually being filled through classical techniques combined with innovative quantitative analyses [1,2] and new technologies [3,4]. The connectivity between neurons is considered to shape resting-state and task-related collective activity [5,6]. 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 [7][8][9]. Structures and mechanisms underlying large-scale interactions have been identified by means of dynamical models [10,11] and graph theory [12,13]. Furthermore, the impact of local network structure on spiketime correlations is known in some detail [14][15][16]. Under certain conditions, there is a one-toone mapping between correlations in neuronal network activity and effective connectivity, a measure that depends on the network structure and on its activity [17]. 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. The integrative loop between experiment, model, and theory. Black arrows represent the classical forward modeling approach: Experimental anatomical data are integrated into a model structure, which gives rise to the activity via simulation. The model activity is compared with experimental functional data. The usual case of disagreement leads to the need to change the model definition. By experience and intuition the researcher identifies the parameters to be changed, proceeding in an iterative fashion. Once the model agrees well with a number of experimental findings, it may also generate predictions on the outcome of new experiments. Red arrows symbolize our new approach: informed by functional data, an analytical theory of the model identifies critical connections to be modified, refining the integration of data into a consistent model structure and generating predictions for anatomical studies.
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 [18].

Global stability in a simple network
We consider networks with neurons structured into N interconnected populations. A neuron in population i receives K ij incoming connections from neurons in population j, each with synaptic efficacy J ij . Additionally, each neuron is driven by K ext Poisson sources with rate ν ext and synaptic efficacy J ext . All neurons in one population have identical parameters, so that we can describe the network activity in terms of population-averaged firing rates ν i .
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 where s is a pseudo-time. The gain function Φ is known analytically and A indicates its dependence on the model parameters {K, J, ν ext , . . .} (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 excitatoryexcitatory loops cause the system defined by Eq (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 [19] without absolute refractoriness (τ r = 0). 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 1 t m R V y À V r (see, e.g., [20], 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. [21], 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 1/τ r = 500 spikes/s 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 K and is driven by external Poisson sources with the same indegree K ext = K and rate ν ext . All connections have identical synaptic weights J. An increase in the external drive shifts the input-output relationship Φ(ν, ν ext ) of this one-dimensional system (Fig 2C) to the left. The bifurcation diagram is shown in Fig 2E: for low ν ext there is only one fixed point with low activity (LA). When increasing ν ext , 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 ν ext , the LA state loses stability.
The equilibria, given by the zeros of the velocity _ n 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 [22], 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 _ n i in direction i vanishes. We use the unstable fixed point to preserve the basin of attraction when the external drive ν ext is increased. For this purpose, we modify the connectivity K ! K 0 to reverse the shift of the unstable fixed point due to the parameter change n ext ! n 0 ext (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 Fðn 0 ext ; K 0 Þ (F defines the system Eq (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 F(ν ext , K). 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 Φðn 0 ext ; K 0 ) is restored to that in the original system Φ(ν ext , K). Fig 2G shows the behavior of the LA fixed point in more detail. The modification of K 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 [21] 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 [23], which covers 1 mm 2 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 [24].
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 [25,26]. Inputs from subcortical and non-visual cortical areas are modeled as Poissonian spike trains, whose rate ν ext 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 ≳ 0.05 spikes/s while preventing the model from entering a state with unrealistically high rates of ≳ 30 spikes/s (figure 13 of [25], [26]). 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 K 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 ν(0) = (0, . . ., 0) T 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 Eq (1) over pseudo-time s 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 Eq (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 N − 1 dimensional separatrix. The separatrix is a  [23], with permission). B Sketch of the most common laminar patterns of cortico-cortical connectivity of the multi-area model. C Population-averaged firing rates encoded in color for a spiking network simulation of the multi-area model with low external drive ðk ¼ 1:0Þ. D As C but for increased external drive k ¼ 1:125. The color bar refers to both panels. Areas are ordered according to their architectural type along the horizontal axis from V1 (type 8) to TH (type 2) and populations are stacked along the vertical axis. The two missing populations 4E and 4I of area TH are marked in black and firing rates < 10 À 2 Hz in gray. E Histogram of population-averaged firing rates shown in C for excitatory (blue) and inhibitory (red) populations. The horizontal axis is split into linear-(left) and logscaled (right) ranges. F as E corresponding to state shown in D.
doi:10.1371/journal.pcbi.1005179.g003 stable manifold [22], 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 s c 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 [27], and we also disallow inhibitory cortico-cortical connections, as the vast majority of longrange connections are known to be excitatory [28,29]. In addition, we naturally restrict indegrees to positive values. We find that four iterations (numbered by index j) 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 κ = 1.15. 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 K into its N 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 j there is exactly one unstable eigendirection with an eigenvalue Re l ðjÞ c À Á > 1 (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 2.3% 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 u ð1Þ c predominately points into the direction of populations 4E and 5E of areas FEF and 46 (Fig 5C), while u ð2Þ c has large entries in the 5E populations of two areas (Fig 5D). The high rates of these populations at the unstable fixed points (cf. Figs 6A and 6B with 5C and 5D) 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 11.3%. 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 [30] underlying the microcircuit model [23] 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 [31] 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 FLN and SLN (see "Methods"). The FLN 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 FLN of r Pearson = 0.79. Significant variations in the FLN 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 SLN determines the laminar pattern of the location of source neurons for cortico-cortical connections. Overall, data are available for 24% of the inter-areal connections in the parcellation of Felleman & van Essen [27], while the SLN 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 (|δSLN| > 0.5) are weak (average FLN ¼ 6 Á 10 À 4 compared to FLN ¼ 10 À 2 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 SLN of the original model that were directly taken from experiments and the logarithmic ratios of cell densities is ρ = −0.63 (p = 3 Á 10 −11 , p-value of a two-sided test for uncorrelated data). For the modified model, we take the SLN of all connections into account and obtain ρ = −0.40 (p = 6 Á 10 −20 ), indicating a reduced, but still significant, monotonic dependence between SLN 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 [1]. We quantify the latter as the average interindividual standard deviation of the logarithmic FLN, i.e., s ¼ h ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi logFLN À logFLN À Á 2 q i, where the overbar " Á denotes the average over injections and hÁi the average over connections. This variability equals 2.17 while the average modification of the logarithmic FLN is 1.34. The main experimental connection probabilities used to construct the intra-areal connectivity of the model have an average relative standard deviation of 30% across electrophysiological experiments (cf. Table 1 of [23]) while the intra-areal connection probabilities of the model are modified by 9% on average. The authors of [32] 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, 83% of these are indicated by only a single tracer injection, while the overall proportion of connections measured by a single injection is 59%.
For the modified connectivity and κ = 1.125, which we choose to avoid being too close to the transition (Fig 7A), the theory predicts average rates in populations 5E and 6E of 1.3 and 0.18 spikes/s, 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 g = −11. Simulating the full spiking network model then results in reasonable rates across populations and areas (Fig 9B and 9D). The average rates in populations 5E and 6E are increased compared to a simulation of the original model from 0.09 and 2 Á 10 −5 spikes/s to 3.0 and 0.4 spikes/s, respectively. All populations exhibit firing rates within a reasonable range of 0.05 to 30 spikes/s (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 and 9C). Small discrepancies are caused by violations of the employed assumptions, i.e., Poissonian spiking statistics [33]. 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 [24]. 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 [23]. 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 [1], albeit in a different parcellation scheme than the scheme used here [27]. 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 [34,35]. 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 [36,37]. This critical role is often attributed to single-neuron properties, with a subset of layer 5 neurons displaying pacemaker activity [38][39][40]. 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 singleneuron 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 [1] combined with the CoCoMac database [41,42], 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 [43], 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 [44,45] 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 [46], featuring a non-monotonic gain function: the gain is reduced if excitatory and inhibitory inputs are increased in a balanced manner [47]. 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 [48]), homeostatic synaptic scaling [49] and spike-frequency adaptation (e.g. summarized in [50]) may prevent the system from entering the high-activity state. However, introducing these selforganizing 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 [51] 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, [52] show that the transition in stochastic systems that quantitatively resemble the spiking network [53], 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 interarea connections [24]. The fixed points cannot be accurately described in this case because the fluctuations need to be taken into account [54]. Simulations show (figures 4 and 5 of [24]) 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.
In the diffusion approximation, the dynamics of the membrane potential V and synaptic current I s is [33] where τ m is the membrane time constant and τ s the synaptic time constant, respectively. The membrane resistance τ m /C m has been absorbed into the definition of the current. The input spike trains are approximated by a white noise current with fluctuations / σ 2 and mean value μ. Here ξ is a centered Gaussian white process satisfying hξ(t)i = 0 and hξ(t)ξ(t 0 )i = δ(t − t 0 ).
Whenever the membrane potential V crosses the threshold θ, the neuron emits a spike and V is reset to the potential V r , where it is clamped during τ r . All neurons in one population have identical parameters, so that we can describe the network activity in terms of population-averaged firing rates ν i that depend on population-dependent input μ i , σ i determined by the connectivity. Using the Fokker-Planck formalism, the stationary firing rates for each population i are given by [33] 1 which is correct up to linear order in ffiffiffiffiffiffiffiffiffiffiffi t s =t m p and where g ¼ jzð1=2Þj= ffiffi ffi 2 p ; with z denoting the Riemann zeta function [55]. Here, A is chosen from the set of model parameters {K, J, ν ext , . . .}. If A is a matrix, we vectorize it by concatenating its rows and indicate this by lower case, i.e., a = (a 00 , a 01 , . . ., a 0N , a 10 , . . ., a 1N , a N0 . . ., a NN ) = vec(A T ) T following [56]. If the chosen parameter is a scalar we denote it with a.
We find the fixed points of Eq (3) by solving the first-order differential equation Eq (1) [57] _ ν ≔ dν ds ¼ Φðν; AÞ À ν; using the classical fourth-order Runge-Kutta method (RK4) with step size h = 0.01, where s denotes a dimensionless pseudo-time. The same approach can be used to solve the activity on a single neuron level [58]. Note that Eq (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 [57] we do not only search for stable fixed points, but also use Eq (1) to obtain unstable attractors (cf. "Results"), an idea originating from the study of simple attractor networks ( [59] esp. their figure 2 and eq 12). In a bistable situation, the initial condition of Eq (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 Eq (1). Even if we knew that Eq (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 b necessary to compensate for the changes in the global stability induced by the change of another parameter a. To this end, we study the behavior of the fixed points with respect to an infinitesimal change δa = a 0 − a in the chosen model parameter. Let ν Ã (a) and ν Ã (a 0 ) be the corresponding locations of the fixed points and δν Ã = ν Ã (a 0 ) − ν Ã (a) their separation. We can then expand ν Ã (a 0 ) into a Taylor series up to first order in δa and obtain where we notice that S i and T i only depend on the target population i. We accordingly define two diagonal matrices S and T with S ii = S i and T ii = T i . We further define the connectivity matrix W = K ⊛ J, where ⊛ denotes element-wise multiplication, also called the Hadamard product (see [60], for a consistent set of symbols for operations on matrices). The derivatives with respect to a j have the compact expressions with the Jacobian ðD a f Þ ij ≔ @f i @a j of some vector-valued function f and ds 2 where we use the Hadamard product again to define the matrix W 2 ≔ K ⊛ J ⊛ J. Inserting the total derivatives into Eq (7), we derive the final expression for Δ a , reading where we use 1 for the identity matrix and define the effective connectivity matrix M and the matrix " Δ a . The latter has dimensionality N × P, where P is the dimension of a (for example, P = N 2 for a = k the vector of indegrees). With the aid of Eq (8), evaluating Eq (6) at the unstable fixed point predicts the shift of the separatrix (Fig 2D) to linear order. We now consider an additional parameter b which is modified to counteract the shift of the unstable fixed point caused by the change in parameter a, i.e., of areas, while for other areas, the neuron density is estimated based on their architectural type (see [24] for details).
The inter-areal connectivity is based on binary data from the CoCoMac database [27,41,42,66,67] indicating the existence of connections, and quantitative data from [1]. The latter are retrograde tracing data where connection strengths are quantified by the fraction of labeled neurons (FLN) in each source area. The original analysis of the experimental data was performed in the M132 atlas [1]. Both the FV91 and the M132 parcellations have been registered to F99 space [68], a standard macaque cortical surface included with the software tool Caret [69]. 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 FLN is mapped to the indegree K AB the target area A receives from source area B divided by its total indegree, i.e., FLN AB ¼ K AB = P B 0 K AB 0 . Here, K AB is defined as the total number of synapses between A and B divided by the total number of neurons in A. On the source side, laminar connection patterns are based on CoCoMac [27,67,[70][71][72][73] and on fractions of labeled neurons in the supragranular layers (SLN) [2]. Gaps in these data are bridged exploiting a sigmoidal relation between SLN and the logarithmized ratio of overall cell densities of the two areas, similar to [74]. We map the SLN 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., SLN AB ¼ where the indices i and j go over the different populations within area A and B, respectively. In the context of the model, we use the terms FLN and SLN to refer to the relevant relative indegrees given here. On the target side, the CoCoMac database provides data from anterograde tracing studies [27,66,75,76].
Missing inputs in the model, i.e., from subcortical and non-visual cortical areas, are replaced by Poissonian spike trains, whose rate ν ext is a free, global parameter. In the original model all populations of a particular area receive the same indegree of external inputs K ext . 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 20%. 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 K 5E,ext incoming connections per target neuron (indegree), in all areas by increasing κ = K 5E,ext /K ext . 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: 0.09 spikes/s for 5E compared to 2 Á 10 −5 spikes/s for 6E). We thus scale up K 6E,ext linearly with κ such that κ = 1.15 results in K 6E,ext /K ext = 1.5.