Tools and Resources

Neuroscience

Modeling and simulation of neocortical micro- and mesocircuitry (Part I, anatomy)

Blue Brain Project, École polytechnique fédérale de Lausanne (EPFL), Campus Biotech, Switzerland
Riga Business School, Riga Technical University, Latvia
Nottingham Trent University, United Kingdom
ELKH-University of Debrecen, Neuroscience Research Group, Hungary
University of Aberdeen, United Kingdom
Laboratory for Topology and Neuroscience (UPHESS), Brain Mind Institute, School of Life Sciences, École polytechnique fédérale de Lausanne (EPFL), Switzerland
Neural Circuits Laboratory, Newcastle University, United Kingdom

Jan 20, 2026

https://doi.org/10.7554/eLife.99688.3

Open access
Copyright information

eLife Assessment

This manuscript reports a detailed model of juvenile rat somatosensory cortex, consisting of 4.2 million morphologically and biophysically detailed neuron models, arranged in space and connected according to diverse experimental data - a valuable tool for the field. The construction of the model is based on a methodology with solid supporting evidence. It should be noted that, by necessity, such a large-scale model development involves many assumptions, interpolations, and decisions that could have compounding downstream effects on further analyses that may be difficult to disambiguate.

https://doi.org/10.7554/eLife.99688.3.sa0

Significance of the findings:

Valuable: Findings that have theoretical or practical implications for a subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Solid: Methods, data and analyses broadly support the claims with only minor weaknesses

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

The function of the neocortex is fundamentally determined by its repeating microcircuit motif, but also by its rich, interregional connectivity. We present a data-driven computational model of the anatomy of non-barrel primary somatosensory cortex of juvenile rat, integrating whole-brain scale data while providing cellular and subcellular specificity. The model consists of 4.2 million morphologically detailed neurons, placed in a digital brain atlas. They are connected by 14.2 billion synapses, comprising local, mid-range and extrinsic connectivity. We delineated the limits of determining connectivity from neuron morphology and placement, finding that it reproduces targeting by Sst+ neurons, but requires additional specificity to reproduce targeting by PV+ and VIP+ interneurons. Globally, connectivity was characterized by local clusters tied together through hub neurons in layer 5, demonstrating how local and interegional connectivity are complicit, inseparable networks. The model is suitable for simulation-based studies, and the model is made openly available to the community.

Introduction

Cortical dynamics underlie many cognitive processes and emerge from complex multi-scale interactions. These emerging dynamics can be explored in large-scale, data-driven, biophysically-detailed models (Markram et al., 2015; Billeh et al., 2020), which integrate different levels of organization. The strict biological and spatial context enables the integration of knowledge and theories, the testing and generation of precise hypotheses, and the opportunity to recreate and extend diverse laboratory experiments based on a single model. This approach differs from more abstract models in that it emphasizes anatomical completeness of a chosen brain volume rather than implementing a specific hypothesis. Using a ‘bottom-up’ modelling approach, many detailed constituent models are combined to produce a larger, multi-scale model. To the best possible approximation, such models should explicitly include different cell and synapse types with the same quantities, geometric configuration and connectivity patterns as the biological tissue it represents.

Investigating the multi-scale interactions that shape perception requires a model of multiple cortical subregions with inter-region connectivity, and for certain aspects, the subcellular resolution provided by a morphologically detailed model is also required. In particular, Barabási et al., 2023 argued that the function of the healthy or diseased brain can only be understood when the true physical nature of neurons is taken into account and no longer simplified into point-neuron networks. Also, Einevoll et al., 2019 pointed out that simulations of large-scale models are essential for bridging the scales between the neuron and system levels in the brain. In that regard, modern electron-microscopic datasets have reached a scale that allows the reconstruction of a ground truth wiring diagram of local connectivity between several thousand neurons (Bae et al., 2021). However, this only covers a small fraction of the inputs a cortical neuron receives. While afferents from outside the reconstructed volume are detected, one can only speculate about the identity of their source neurons and connections between them. The scale required to understand inter-regional interactions is only available at lower resolutions in the form of region-to-region or voxel-to-voxel connectivity data.

To help better understand cortical structure and function, we present a general approach to create morphologically detailed models of multiple interconnected cortical regions based on the geometry of a digitized volumetric brain atlas, with synaptic connectivity predicted from anatomy and biological constraints (Figure 1). We used it to build a model of the juvenile rat non-barrel somatosensory (nbS1) regions (Figure 1, center). These regions were selected for the wealth of available experimental data from various labs, and to build upon our previous modeling work (Table 1). The workflow is based on the work described in Markram et al., 2015, with several additions, refinements and new data sources that have been independently described and validated in separate publications (Table 1). The model captures the morphological diversity of neurons and their placement in the actual geometry of the modeled regions through the use of voxelized atlas information. We calculated at each point represented in the atlas the distance to and direction towards the cortical surface (Figure 1; step 1). We used that information to select from a pool of morphological reconstructions anatomically fitting ones and orient their dendrites and axons appropriately (Figure 1; step 2). As a result, the model was anatomically complete in terms of the volume occupied by dendrites in individual layers. We then combined established algorithms for the prediction of local (Reimann et al., 2015) and mid-range (Reimann et al., 2019) connectivity (Figure 1; step 3) as well as extrinsic connectivity from thalamic sources (Markram et al., 2015, Figure 1; step 4) to generate a connectome at subcellular resolution that combines those scales.

Figure 1 with 1 supplement see all

Download asset Open asset

Overview of the model building and analysis workflow.

Step 1: Building was based on a volumetric atlas of the modeled regions: (1) S1J; (2) S1FL; (3) S1Tr; (4) S1HL; (5) S1Sh; (6) S1DZ; (7) S1DZO; (8) S1ULp. Additional atlases of biological cell densities and local orientation towards the surface were built. Step 2: Neuron morphologies were reconstructed and classified into 60 morphological types (*m-types*). They were placed in the volume according to the densities and orientations from step 1. Step 3: The anatomy of intrinsic synaptic connectivity was derived as the union of one algorithm for *local* connectivity and one for *mid-range* connectivity. Step 4: Extrinsic inputs from two thalamic sources were placed on modeled dendrites according to published methods. Step 5: Taken together, these steps allowed us to predict the topology of connectivity at scale with (sub-)cellular resolution. Step 6: The anatomical model served as the basis of a physiological model, ready to be simulated. This is presented in an accompanying manuscript. Step 7: The model, simulation and analysis tools have been made publicly available. Left: During modeling, three types of generalization had to be made to fulfill input data requirements: from mouse to rat, from adult to juvenile, and from one cortical region to another. Generalizations used are indicated in each step.

We characterized several emerging aspects of connectivity (Figure 1; step 5). First, we found that brain geometry, that is differences in cortical thickness and curvature have surprisingly large effects on how much individual layers contribute to the connections a neuron partakes in. Second, we characterized the predicted structure of connectivity at an unprecedented scale and determine its implications for neuronal function. In particular, we analyzed how the widths of thalamo-cortical axons constrains the types of cortical maps emerging. Furthermore, we characterized the global topology of interacting local and mid-range connectivity, finding highly complex topology of local and mid-range connectivity that specifically requires neuronal morphologies. Finally, we systematically analyzed the higher-order structure of connectivity beyond the level of pairwise statistics, such as connection probabilities. Doing so, we characterized highly connected clusters of neurons, distributed throughout the volume that are tied together by mid-range synaptic paths mediated by neurons in layer 5, which act as ‘highway hubs’ interconnecting spatially distant neurons in the model. The highly non-random structure of higher-order interactions in the model’s connectivity was further validated in a range of follow-up publications using the model (Ecker et al., 2024b; Egas Santander et al., 2024a; Reimann et al., 2024).

Finally, we present an accompanying manuscript that details neuronal and synaptic physiology modeled on top of these results, describes the emergence of an in vivo-like state of simulated activity, and delivers a number of in silico experiments generating insights about the neuronal mechanisms underlying published in vivo and in vitro experiments (Isbister et al., 2024; Figure 1; step 6).

The anatomically detailed modeling approach provides a one-to-one correspondence between most types of experimental data and the model, allowing the data to be readily integrated. However, this also leads to the difficult challenge to curate the data and decide which anatomical trend should be integrated next. Due to the incredible speed of discovery in the field of neuroscience an integrative model will always be lagging behind the latest results, and due to its immense breadth, there is no clear answer to which feature is most important. We believe the solution to this is to provide a validated model with clearly characterized strengths and weaknesses, along with the computational tools to customize it to fit individual projects. We have therefore made not only the model, but also most of our tool chain openly available to the public (Figure 1; step 7).

Already the process of adding a new data source to drive a refinement of the model serves to provide important insights. We demonstrate this by comparing the model to connectivity characterized through electron microscopy (Bae et al., 2021), finding mismatches, and describing the changes required to fix them. This allowed us to determine which rules are required to predict connectivity from the locations and densities of neuronal processes. Previously, simple overlap of distributions of axonal and dendritic segments has been proposed (Peters’ rule; Peters and Feldman, 1976; Garey, 1999), and contrasted with findings of preference for specific cell types or subcellular domains (White and Keller, 1987; Mishchenko et al., 2010). Our approach to local connectivity combines overlap with the principle of cooperative synapse formation (Fares and Stepanyants, 2009; Reimann et al., 2015), additionally Schneider-Mizell et al., 2024 proposed a combination of overlap and targeting preferences. Our comparison to electron microscopy let us uncover the strength and nature of the targeting preferences shaping connectivity beyond neuronal and regional anatomy. We found that cooperative synapse formation explains some forms of apparent targeting. Additionally, we found that the distribution of postsynaptic compartments targeted by connections from somatostatin (Sst)-positive neurons is readily predicted from overlap only, while for parvalbumin (PV)-positive and vasoactive intestinal peptide (VIP)-positive neurons additional specificity plays a role. We found no indication of additional specificity for excitatory neurons. The model is available both with and without the characterized inhibitory specificity.

Results

Atlas-based neuron placement

The workflow for modeling the anatomy of juvenile rat nbS1 is based on the work described in Markram et al., 2015, with several additions and refinements. Most of the individual steps and data sources have already been independently described and validated in separate publications (Table 1). The basis of the work was a digital atlas of juvenile rat somatosensory cortex, based on the classic work of Paxinos and Watson, 2007 (Figure 1, step 1.). The atlas provides region annotations, that is each voxel is labeled by the somatosensory subregion it belongs to. It also provided a spatial context for the model, with the non-modeled barrel region being surrounded on three sides by modeled regions (specifically: S1DZ and S1ULp). The atlas was enhanced with spatial data on cortical depth and local orientation towards the cortical surface (Bolaños-Puchet et al., 2024; Bolaños-Puchet et al., 2023). We further added annotations for cortical layers, placing layer boundaries at fixed normalized depths (Supplementary file 5). Using laminar profiles of cell densities (Figure 1—figure supplement 1; Keller et al., 2019), we generated additional voxelized datasets with cell densities for each morphological type (see Materials and methods).

Table 1

References to publications of input data and methods employed for individual modeling steps.

An asterisk next to a reference indicates that substantial adaptations or refinements of the data or methods have been performed that will be explained in this manuscript. In the other cases, a basic summary will be provided and an exhaustive description in the Methods.

Data
Stage	Topic	Reference
Atlasing	Region annotation atlas	Paxinos and Watson, 2007
	Orientation, depth and flat map	Bolaños-Puchet et al., 2023
	Cell density profiles	Keller et al., 2019
Volume filling	Neuron reconstructions	Markram et al., 2015
	in vivo neuron reconstructions	*New, original data
Synaptic connectivity	Bouton densities and numbers of synapses per connection	Reimann et al., 2015
	Pathway strengths, synapse density profiles and topographical mapping	Reimann et al., 2019
Thalamic inputs	Bouton density profiles	Meyer et al., 2010
	Projection axon lengths and widths	Economo et al., 2016
Modeling methods
Stage	Topic	Reference
Atlasing	Cell density volume generation	Keller et al., 2019
Volume filling	Neuron classification	Kanari et al., 2019
	Neuron placement	*Markram et al., 2015
Synaptic connectivity	Local connectivity	Reimann et al., 2015
	Inter-region connectivity	Reimann et al., 2019
Thalamic inputs	Input generation	*Markram et al., 2015

Next, we filled the modeled volume of the atlas with neurons (Figure 1, step 2.). Neuronal morphologies were reconstructed in slices, and repaired algorithmically (Anwar et al., 2009; Markram et al., 2015). Out of 1017 morphologies, 58 were new in vivo stained reconstructions used for the first time in this work (see Materials and methods). The morphologies were classified into 60 morphological types (m-types, see ; 18 excitatory Figure 2A, 42 inhibitory), based on expert knowledge and objectively confirmed by topological classification (Kanari et al., 2019).

Figure 2 with 5 supplements see all

Download asset Open asset

Anatomy of the model.

(A) Exemplar 3D reconstructions of the 18 excitatory m-types. Rendering and visualization was done in NeuroMorphoVis (Abdellah et al., 2018). Dendritic diameters are scaled (x3) for better resolution. (B) Modeled cortical layers, with exemplars of excitatory morphological types placed in the model. (C) The placement of each morphological type recreates the biological laminar anatomy of dendrites and axons, which then serves as the basis of local connectivity. Inset shows modeled brain regions (solid colors) in the context of the non-modeled regions (transparent).

As only 58 morphologies were in vivo stained reconstructions, the rest potentially suffered from slicing artifacts (see Figure 2—figure supplement 2 for examples), despite applying a repair algorithm (Markram et al., 2015). A topological comparison (Kanari et al., 2018; Kanari et al., 2019) between axons and dendrites of neurons in all layers (Figure 2—figure supplements 1 and 3) revealed that in vitro reconstructions from slices could not capture detailed axonal properties beyond 1000 µm, but could faithfully reproduce dendritic arborization.

Cell bodies for all m-types were placed in atlas space according to their prescribed cell densities. At each soma location, a reconstruction of the corresponding m-type was chosen based on the size and shape of its dendritic and axonal trees (Figure 2—figure supplement 4). Additionally, it was rotated according to the orientation towards the cortical surface at that point. These steps ensured that manually identified features of the morphologies (see Supplementary file 6) landed in the correct layers (Figure 2C, Figure 2—figure supplement 5). Additionally, a random rotation around an axis orthogonal to layer boundaries was applied.

Biological dendrite volume fractions emerge in the model

The distribution over 60 m-types in the model captured the great morphological diversity of cortical neurons and matched the reference data used (Figure 3A, see Materials and methods). The placement of morphological reconstructions matched expectation, showing an appropriately layered structure with only small parts of neurites leaving the modeled volume (Figure 2C). For a more quantitative validation, we calculated the fraction of the volume occupied by neurites in 100 depth bins of a cylindrical volume spanning all layers and with a radius of 100 µm (Figure 3B). This included an estimated volume for axons forming the mid-range connectivity between modeled regions, based on its synapse count (see below and Materials and methods).

Figure 3

Download asset Open asset

Volumetric anatomy of the model.

(A) M-type composition per layer: For each layer, stacked histograms of the relative fractions of each m-type, comparing the model (right bar for each layer) to the input data (left bar). For simplicity, the layer designation is stripped from each m-type. Left: For excitatory types; right: Inhibitory types. (B) Stacked histograms of the fraction of space filled by neurites at various depths in the model. The y-axis indicates the distance in µm from the bottom of layer 6. (B1) For neurites of neurons in different layers indicated in different colors. Grey: estimated lower bound for the volume of axons supporting the mid-range connectivity. (B2) For neurites of different m-types. Colors as in A, but inhibitory types are grouped together. (B3) For different types of neurites. (C) Comparing fractions for axons and dendrites to the literature (Santuy et al., 2018). X-marks indicate overall means, other marks means in individual layers; teal: reference, black: model.

We found a clearly layered structure also for the neurite volumes, with apical dendrites contributing substantially to the volume in layers above their somas. The volume was clearly dominated by dendrites, filling between 23% and 47% of the space, compared to 2% to 11% for axons (Figure 3B3). The range of values for dendrites matched literature closely (between 26% and 46%, Santuy et al., 2018). Conversely, values for axons were lower than literature (between 17% and 24%); however, this is explained by the vast amount of external input axons from non-somatosensory and extracortical sources that are not part of the model.

Local, mid-range and extrinsic connectivity modeled separately

Synaptic connectivity in comparable models is often modeled by imposing experimentally characterized connection probabilities, potentially distance-dependent, onto the various neuronal pathways (Potjans and Diesmann, 2014; Pronold et al., 2024). As this approach directly parameterizes structural connections strengths at the population level it allows easy recreation of meso-scale architectures, as characterized in, e.g., Jiang et al., 2015. On the other hand, it has been shown that top-down imposed connection probabilities, even distance-dependent ones, cannot fully capture the non-random higher order structure of connectivity at cellular resolution (Gal et al., 2021). Such non-random connectivity has been repeatedly demonstrated in cortical circuitry (Song et al., 2005; Perin et al., 2011; Reimann et al., 2024; Egas Santander et al., 2024a) and has functional impact (Reimann et al., 2017b; Egas Santander et al., 2024a; Nolte et al., 2020). Instead, we used a previously published approach (Reimann et al., 2015) that selects synapses as a subset of axo-dendritic, axo-somatic and axo-axonic appositions, based on biologically motivated rules (Figure 1, step 3; 4 A). Note that this algorithm is not equivalent to the so-called ‘Peters’ rule’ (Peters and Feldman, 1976; Garey, 1999; Rees et al., 2017), as the subset is not selected completely randomly. Crucially, the concept of cooperative synapse formation (Fares and Stepanyants, 2009) is used, which introduces statistical dependence for potential synapses between the same pair of neurons: If one is selected the probability that others are selected increases. The algorithm has been demonstrated to result in a nonrandom higher-order structure of connectivity matching biological references (Reimann et al., 2015; Gal et al., 2017; Reimann et al., 2017a; Reimann et al., 2017b). Conversely, connection probabilities are not directly controlled by the modeler. We measured the emerging connection probabilities by emulating a multi-patch procedure that is often used for that purpose in vitro (Figure 4B, Figure 4—figure supplement 2A, see Materials and methods). We compared them to connection probabilities reported in 124 literature sources gathered by Zhang et al., 2019 (Figure 4—figure supplement 2B–D). We note that reported connection probabilities vary drastically, for example, ranging from zero to 88% for L2/3 EXC to L2/3 PV. Significant variability persisted when only data sources for rat or the age of the model (P14) were considered, with ranges of values stretching up to 51% and 54%, respectively. While values for the model were sometimes outside the reported range, this was mostly restricted to pathways with only one or two sources.

Figure 4 with 4 supplements see all

Download asset Open asset

Intrinsic connectivity as union of local and mid-range connectivity.

(A) Local connectivity was derived as in Markram et al., 2015; Reimann et al., 2015, that is based on neuron placement and morphologies. Right: Axo-dendritic appositions were considered as potential synapses. They were then filtered to yield among other biological constraints biologically realistic bouton densities. (B) Range of resulting connection probabilities within.100 µm Indicated are the minimum (inner circle) and maximum (outer circle) found over 20 repetitions of sampling 1000 pairs in a simulated multi-patch procedure (Figure 4—figure supplement 2A). See Materials and methods for mapping from morphological types to inhibitory subclasses (PV/SST/INH). (C) An exemplary layer 5 PC axon reconstruction from Janelia MouseLight (mouselight.janelia.org, Gerfen et al., 2018) that includes mid-range collaterals, shown in the context of mouse somatosensory regions. The blue circle highlights local branches all around the soma. Red highlights depict more targeted collaterals into neighboring regions. (D1) Predicted pathway strengths as indicated in Figure 4—figure supplement 1. (D2) Pathway strength emerging from the application of the apposition-based connectivity algorithm described in Reimann et al., 2015. (E) Values of the diagonal of (D1) compared to the diagonal of (D2). (F) Schematic of the strategy for connectivity derivation: Within a region only apposition-based connectivity is used; across regions the union of apposition-based and mid-range connectivity. (G) Connection strength constraints for the mid-range connectivity derived by subtracting D2 from D1 and setting elements $< 0$ to 0 (colors as in D). (H) Resulting total density from both types of modeled synaptic connections in individual regions compared to the data in (D1).

To create biologically realistic overall neuron out-degrees, our algorithm requires axon reconstructions to be complete. While we try to ensure the completeness of axonal arborizations through the use of in vivo stained reconstructions and repair algorithms, this becomes more challenging, the more an axon stretches away from its soma, due to lack of dye penetration and potentially slicing artifacts. According to our analysis of the morphologies used, most axons were only accurately reconstructed up to 1000 µm from the soma (see above). This contrast between local axon around the soma and more mid-range collaterals stretching into neighboring regions can also be seen in in vivo stained reconstructions (Figure 4C). Therefore, to model connections at larger distances, we used a second, previously published algorithm (Reimann et al., 2019). It places mid-range synapses according to three biological principles that all need to be separately parameterized based on experimental data (as detailed in the subsequent paragraphs and Figure 4—figure supplement 1): First, connection strength: ensuring that the total number of synapses in a region-to-region pathway matches biological estimates. Second, layer profiles: ensuring that the relative number of synapses in different layers matches biological estimates. Third, topographical mapping: Ensuring that the specific locations within a region targeted by mid-range connections of neurons describe a biologically parameterized, topographical mapping. We note that in the literature on macro-scale connectomics, the pathways between somatosensory regions would be considered ‘short-range, inter-regional connections’. We will still refer to connections from the two algorithms as the local and mid-range connectomes respectively, for simplicity and because of the intuitive contrast it invokes.

To parameterize connection strengths, we used data from axon tracing experiments provided by the Allen Institute (Harris et al., 2019), adapted from mouse to rat, yielding expected densities of projection synapses between pairs of regions (Figure 4D1). First, we asked to what degree the local connectome algorithm suffices to model the connectivity within a region. To that end, we compared the mean excitatory synapse densities of local connectivity to the target values adapted from Reimann et al., 2019, diagonal Figure 4D1 vs. Figure 4D2; Figure 4E. These target values were derived from relative region-to-region connectivity strengths, estimated by axon tracing, scaled to match overall volumetric synapse densities estimated from electron microscopy. We found that the overall average matches the data fairly well, however the variability across regions was lower in the model ( $0.123 \pm 0.017$ $μ m^{- 3}$ ; mean ± std in the model vs. $0.097 \pm 0.06$ $μ m^{- 3}$ ). Based on these results, we decided that the local connectome sufficed to model connectivity within a region. It also created a number of connections across region borders (Figure 4D2). Consequently, we parameterized the strengths of additional mid-range connections to be placed as the difference between the total strength from the data and the strengths resulting from local connectivity, with connection strengths within a region set to zero (Figure 4G). As a result synaptic connections between neighboring regions will be placed by both algorithms (Figure 4F), with a split ranging from 20% local to 70% local. The lower spread of apposition-based synapse density within a region (see above) will also reduce the variability of combined synapse density from both algorithms (Figure 4H). While this will halve the coefficient of variation of density across regions from 0.34 to 0.17, the overall mean density over all regions is largely preserved (0.23 data vs. 0.24 for the combined algorithms; Figure 4H, red dot).

Finally, the dendritic locations of synaptic inputs from thalamic sources were modeled as in Markram et al., 2015, using experimental data on layer profiles of bouton densities of thalamo-cortical axons, and morphological reconstructions of these types of axons (Figure 1, step 4; Figure 4—figure supplement 4A; Materials and methods). Based on these data, each thalamic input fiber was assigned an innervated domain in the model, recreating the layer profile along an axis orthogonal to layer boundaries and spreading equally in the other dimensions (Figure 4—figure supplement 4B–E). We modeled two types of thalamic inputs, based on the inputs into barrel cortex from the ventral posteromedial nucleus (VPM-based) and from the posterior medial nucleus (POm-based), respectively (Harris et al., 2019; Shepherd and Yamawaki, 2021). While barrel cortex was not a part of the model, we used these projections as examples of a core-type projection, providing feed-forward sensory input (VPM-based) and a matrix-type projection, providing higher-order information (POm-based).

The numbers of thalamic input fibers innervating the model were estimated as follows. Laminar synapse density profiles were summed over the volume to estimate the total number of thalamo-cortical synapses, and the number of synapses per neuron was estimated from the lengths of thalamo-cortical axons. The ratio of these numbers resulted in 72,950 fibers for the POm-based matrix-type projection, and 100,000 fibers for the VPM-based core-type projection. These numbers are consistent with the volume ratio of the two thalamic nuclei ( $1.25 m m^{3}$ for POm to $1.64 m m^{3}$ for VPM; ratio: 0.76).

Specificity of axonal targeting

For the axons of various inhibitory neuron types, certain aspects of connection specificity have been characterized (Tremblay et al., 2016), such as a preference for peri-somatic innervation for PV+ basket cells. The local connectivity in the model is based axo-dendritic overlap, combined with a pruning rule that prefers multi-synaptic connections, but does not take the post-synaptic compartment type into account (Reimann et al., 2015). As such, it captures targeting preferences resulting from specific axonal morphologies, but not potential effects due to potential molecular mechanisms during synapse formation or pruning. This is based on the idea that the developmental mechanisms underlying connectivity are manifold and complicated, but their cumulative effect is at least partially visible in the characteristic shapes of the neurites. It has been demonstrated by us and other groups that anatomy-based approaches suffice to create highly nonrandom network topologies that match many biological trends (Gal et al., 2017; Reimann et al., 2017a; Udvary et al., 2022), but the analyses have so far been focused mainly on excitatory sub-networks. To investigate to what degree the approach suffices to explain observed targeting trends, we compared the model to characterizations of targeting specificity from the literature, recreating the electron microscopy studies of Motta et al., 2019 and Schneider-Mizell et al., 2024 in silico.

Motta et al., 2019 analyzed an approximately 500,000 µm³ volume of tissue in layer 4, considering the postsynaptic targets of the contained axons. As axons in the volume were fragments, specificity was assessed by comparing the data to a binomial control fit against observations of axons forming at least one synapse onto a postsynaptic compartment type (see Materials and methods). Axons with unexpectedly high fractions onto a compartment type were then considered to target that type. We calculated targeting fractions of axons forming at least 10 synapses in a comparable layer 4 volume of our model (Figure 5A). Compared to Motta et al., 2019, fractions of synapses onto apical dendrites were elevated, and more inhibitory synapses were considered (Figure 5B). As in the original study, for almost all compartment types, observed distributions were more long tailed than expected with a number of axons showing a significantly high targeting fraction (Figure 5C). The fractions of axons with such specificities also matched the reference, except for an even higher fraction being specific for apical dendrites in the model (Figure 5D). The match resulted from local connectivity in our model using the principle of cooperative synapse formation (Fares and Stepanyants, 2009), where all synapses forming a connection are kept or pruned together. This creates a statistical dependence between the synapses that results in a significant deviation from the binomial control models. We demonstrate this with a simple stochastic model of cooperative synapse formation, fitted to the data of Motta et al., 2019. For each axon in the data, we generated a random list of compartments types for potential synapses matching the overall frequencies observed in the data. For each, the number generated was 2.5 times the number of synapses of the corresponding axon. Next, we pruned synapses to the original counts in two ways: First, we grouped potential synapses onto the same type into groups of three that are kept or discarded together, thereby introducing statistical dependence. The results provide the same amount of apparent targeting as the data (Figure 5—figure supplement 1A-C). Second, we pruned without the grouping where no targeting was found (Figure 5—figure supplement 1D-F). We do not claim that the stochastic model and its parameters provide an accurate description of the underlying biological processes, only that it demonstrates the effect of cooperative synapse formation. We conclude that the non-random trends observed by Motta et al., 2019 are an indication of cooperative synapse formation captured by our model.

Figure 5 with 1 supplement see all

Download asset Open asset

Recreating an electron-microscopic specificity analysis In the top part we analyzed synapses, axons and dendrites in an approximately cubic volume in L4, emulating the techniques of Motta et al., 2019.

In the bottom part we analyzed a 100 x 100 µm volume of the MICrONS dataset defined by Schneider-Mizell et al., 2024. (A) Fraction of synapses placed on somata (SOM), proximal dendrites (PD), smooth dendrites (SD), apical dendrites (AD) and axon initial segments (AIS) for axon fragments in the sampled volume. Black bars indicate mean values over axons, arrows indicate binomial probabilities fit to observations of axons forming at least a single synapse; pink: for excitatory axons, black: inhibitory. Fractions missing from 100% are onto other compartment types. (B) Overall count on synapses on different postsynaptic compartment classes inside the studies volumes, comparing the data of Motta et al., 2019 to the model. Colors as indicated in the legend. (C) Distributions of the fractions of synapses onto different compartment types over excitatory (left) and inhibitory (right) axons. Expected from the binomial model in A (grey) against the observations in our anatomical model (black outlines). (D) Fraction of axons with significantly increased synapse counts onto different compartment types compared to the binomial control. Indicated for two values of the false detection rate criterion (q=0.05, 0.3, Storey and Tibshirani, 2003). Comparing the reference data to our anatomical model. (E) Top numbers of neurons in four connectivity-derived classes in the 100x100 µm volume of the MICrONS dataset, defined by Schneider-Mizell et al., 2024. Bottom: Numbers in a comparable volume of the model. For assignment of m-types into the four classes see the main text. (F) Derivation of alternative connectivity with targeting specificity, using the example of the ‘PeriTC’ class. The axo-dendritic appositions of an axon (top) are classified as matching the targeting (green box; here: appositions with somata or proximal dendrites) or not. Non-matching appositions are removed with probability $p_{n t}$ (middle). This is followed by a non-specific removal of connections (formed by single or multiple synapses) until the biological density of synapses on the axon is met (bottom). (G) Bouton densities resulting from the process in F. Left to right: For PeriTC neurons, targeting somata and proximal dendrites; for InhTC neurons, targeting inhibitory neurons; for SparTC neurons, targeting the first synapse of a connections; for SparTC neurons, without targeting preference. Respective optimized $p_{n t}$ values reported in the main text. (H) Resulting targeting of postsynaptic compartments, fraction of synapses in multisynaptic connections and clumped synapses (see Schneider-Mizell et al., 2024). Red: Default model; blue: optimized alternative model with specific targeting; grey: optimized alternative model without specific targeting; orange: data of Schneider-Mizell et al., 2024 Boxes outline the characteristic feature of each class. (I) Inhibitory targeting specificities combining the best performing models.

Recently, the MICrONS dataset (Bae et al., 2021) has been analyzed with respect to the axonal targeting of inhibitory subtypes in a subvolume of 100 x 100 µm² surface area spanning all layers of the cortex (Schneider-Mizell et al., 2024). Similar to Motta et al., 2019 they considered their distributions of types of postsynaptic compartments. But due to the larger reconstructed volume, they were able to analyze complete or almost complete axons, allowing for a quantitative comparison rather than relying on a fitted binomial control model. Additionally, they calculated the fraction of synapses that are part of a multisynaptic connections and the fraction thereof that is within 15 µm of another synapse of the same connection.

A comparable volume of the model (see Materials and methods) contained 173 interneurons vs. 163 in the original study. Their distribution into four connectivity classes according to morphological and molecular determinants hypothesized by Schneider-Mizell et al., 2024 approximately matched as well (Figure 5E; perisomatic targeting: Basket Cells; distal targeting: Sst+; sparsely targeting: Neurogliaform Cells and L1; inhibitory targeting: bipolar and VIP+). The postsynaptically targeted compartments matched for the distal targeting group (Figure 5H, top right; Figure 5I), indicating that their preference can be explained by their axonal morphology, specifically its trend to ascend to and then branch in superficial layers. For the other three classes, their respective eponymous trends were not recreated (Figure 5H, red vs. orange).

We then explored the strengths of targeting mechanisms required to explain postsynaptic targets in a modified version of our connectivity algorithm: Beginning with all axo-dendritic appositions as potential synapses, we first remove appositions that are not placed on the preferred postsynaptic compartment with a probability $p_{n t}$ (Figure 5F; top to middle). This replaces a non-specific, but otherwise identical first pruning steps in the regular version of our algorithm. This is followed by removing connections non-specifically, until the biological density of synapses on the axon is matched (Figure 5F; middle to bottom). As reference for biological densities of synapses on axons, we use the sources listed in Reimann et al., 2015.

For the perisomatic targeting class, probability $p_{n t}$ to remove non-proximal, non-soma synapses was optimized against the data of Schneider-Mizell et al., 2024 to 97%. The remaining synapses had a density on the axons of perisomatic targeting cells that was three times higher than biology, which could be reduced in the second, non-specific pruning step removing two thirds of the connections (Figure 5G, left). This indicates substantial room for rewiring through structural plasticity while preserving the targeting specificity of perisomatic targeting cells. The resulting specificity of postsynaptic compartments and multi-synaptic connections then match the reference data (Figure 5H, top left, blue vs orange; Figure 5I).

For the inhibitory targeting class, probability $p_{n t}$ to remove synapses on non-inhibitory neurons was optimized to a similar value of 96.5%. Curiously, this first step already reduced the resulting axonal density of synapses to the biological value, indicating that this class of interneurons cannot perform substantial rewiring without losing its targeting specificity (Figure 5G, second from left).

For sparsely targeting cells, we evaluated two hypotheses: First, we note that this targeting class is associated with Neurogliaform Cells, which are known to have volumetrically transmitting synapses. It is possible that the sparseness of their targeting can be explained by very few of their synapes having an anatomical postsynaptic partner, rather than by a targeting mechanism. Indeed a non-specific removal of 96% of all synapses recreated the sparsity of these connections found in Schneider-Mizell et al., 2024 (Figure 5H, bottom left, grey vs orange; Figure 5I). This reduced the axonal density of synapses to 30% of the biological value, implying that the remaining 70% may be volumetrically transmitting (Figure 5G, right). Second, we randomly picked a ‘first’ synapse from each connection formed by this class that we considered to be targeted. Of the remaining synapses, we removed $p_{n t}$ = 95%, which recreated the characteristic sparsity of the connections equally well (Figure 5H, bottom left, blue vs orange; Figure 5I). In this case, the second, non-specific pruning step was required, indicating substantial room for rewiring (Figure 5G, second from right).

Using a recently developed computational tool (‘Connectome-Manipulator’; Pokorny et al., 2024a), we created a rewired connectome instance that combines the existing excitatory connectivity with the additional rules for the inhibitory connectivity in a single connectome. This rewired connectome is openly available in SONATA format (Dai et al., 2020), see key resources table in the Materials and methods.

Structure of thalamic inputs

Though we have found that the anatomy-based prediction of connectivity underestimates the specificity of some inhibitory connection types, it remains a powerful tool that has been demonstrated to recreate non-random trends of excitatory connections that make up the majority of synapses (Reimann et al., 2015; Reimann et al., 2017b; Gal et al., 2017). We therefore set out to characterize the anatomy and topology of connections at all scales considered in the model. We began with the thalamic input connections, whose axons were generated based on the axon density profiles and the statistics of their horizontal spread (Figure 6—figure supplement 1). Synapses were established based on their volumetric overlap with cortical dendrites to match the laminar synapse density profiles of the presynaptic thalamic axons (see Materials and methods and Figure 4—figure supplement 4). The individual postsynaptic dendritic morphologies, cell placement, and orientation were considered in this process. The process did not prefer the dendrites of one neuron type over another. In contrast, Cruikshank et al., 2007; Cruikshank et al., 2010; Sermet et al., 2019 found physiologically stronger thalamo-cortical inputs onto excitatory and PV-positive neurons than onto SST-positive neurons. While the authors of the studies make compelling arguments that this at least partially reflects the anatomical strength of the pathways, others have cautioned against mixing of the anatomy and physiology of connectivity (Sporns, 2013). We therefore believe that further validation of the thalamo-cortical pathways is best performed in the accompanying manuscript (Isbister et al., 2024), where the experimental methods of the references are closely recreated.

Regarding the laminar structure, we found for both projections that the peaks of the mean number of thalamic inputs per neuron occur at lower depths than the peaks of the synaptic density profiles (Figure 6A). This is consistent with synapses on apical dendrites of PCs being higher than their somas, but the fact that most peaks occur at places where the synapse density is close to zero gives a clear indication that synapse density profiles alone can be misleading about the location of innervated neurons. At the level of individual neurons, the number of thalamic inputs varied greatly, even within the same layer (Figure 6B). Overall, the matrix-type projection innervated neurons in layers 1 and 2 more strongly than the core-type projection, while in layers 3, 4, and 6, the roles were reversed. Neurons in layer 5 were innervated on average equally strongly by both projections, although layer 5a preferred the matrix-type and layer 5b the core-type projection.

Figure 6 with 2 supplements see all

Download asset Open asset

Anatomy of thalamic innervation.

(A) Depth profiles of synapse densities (dashed lines) and mean number of thalamic inputs per neuron (solid lines) for core- (red) and matrix-type (blue) thalamo-cortical projections. Shaded area indicates the standard error of mean. For a region-specific validation of synapse densities (dashed lines) against experimental data, see Figure 6—figure supplement 2. (B) Mean and standard deviation of the number of thalamic inputs for neurons in individual layers or all neurons. Colors as in A. (C) Common thalamic innervation (CTI) of an exemplary neuron (black dot) and neurons surrounding it, calculated as the intersection over union of the sets of thalamic fibers innervating each of them. Scale bar: 200 µm. (D) CTI of pairs of neurons at various horizontal distances. Dots indicate values for 125 randomly picked pairs; lines indicate a sliding average with a window size of 40 µm. We perform a Gaussian fit to the data, extracting the amplitude at 0 µm (A) and the standard deviation ( $σ$ ). (E) Values of A (unitless) and $σ$ (in flatmap coordinates) for pairs in the individual layers or all pairs. Colors as in A.

To characterize the horizontal structure, we introduced the common thalamic innervation (CTI, see Materials and methods) as a measure of the overlap in the thalamic inputs of pairs of neurons. As pairs with many common inputs are likely to have similar stimulus preferences, the magnitude and range of this effect has consequences for the emergence of functional assemblies of neurons. As expected from the horizontal extent of individual fibers, the CTI was distance-dependent (Figure 6C,D), and showed strong variability. Even directly neighboring pairs might not share a single thalamic afferent, leading to a sparse spatial distribution of pairs with strong overlap. A Gaussian fit of the distance dependence of CTI revealed roughly equally strong overlapping innervation for both core- and matrix-type projections (Figure 6E1). The strength of the overlap increased for lower layers in the case of core-type and decreased for matrix-type projections, while being equally strong in layer 5. The horizontal range of common innervation was larger for matrix-type projections in all layers (Figure 6E2). In summary, while the anatomy of thalamo-cortical projections introduces a spatial bias into the emergence of cortical maps, it is relatively weak on its own and supports different stimulus preferences even for neighboring pairs of neurons.

Local brain geometry affects connectivity

Cortex is often seen as a homogeneous structure with parallel layers, but at the larger spatial scales considered in this model, significant variability in its height and curvature can be observed. Our approach to connectivity allowed us to predict the impact of local brain geometry on connectivity (see also Ronan et al., 2011). We partitioned the model into columnar subvolumes with radii of approximately 230 µm (Figure 7A) that we then analyzed separately. We do not claim that these columns have a biological meaning on their own, they only serve to discretize the model into individual data points each representing localized circuitry affected by local geometry.

Figure 7 with 2 supplements see all

Download asset Open asset

A1, A2: Parcellation of the modeled volume into 230 µm radius columns.

Exemplary slice of columns highlighted in green. (B) Geometrical metrics of column subvolumes in the flat view. Peripheral columns masked out (grey); green outline: highlighted columns in A. Left: Conicality, defined as the slope of a linear fit of depth against column radius. Negative values indicate narrowing towards L6. Right: Column height, i.e. cortical thickness at the location of the column. (C1) Modularity of the networks of connections within each column. (C2) Column conicality and height colored by modularity; (D1) Conicality of columns against their laminar neuronal composition, normalized against the overall composition of the model. Colored lines indicate linear fits. (D2) Conicality against the density of connections in subnetworks given by the intersections of columns with individual layers. (E) Counts of afferents formed onto neurons in individual columns from neurons in the entire model. (E1) Normalized in-degrees originating from neurons in individual layers, plotted against conicality. (E2) In-degree of neurons in individual layers normalized by the overall in-degree in the model into each layer plotted against column height. (F) r-values of linear fits against generalized, n-dimensional in-degree as in E2. (F1) Of generalized in-degree against height; (F2) Of generalized in-degree against conicality.

We quantified the two main variable factors of geometry, measuring the height and conicality (a feature that measures the convexity of a column, see Materials and methods) of each column (Figure 7B). We found that differences in these factors lead to differences in relevant topological parameters, such as the modularity of the local network (Figure 7C). This effect was mostly mediated by differences in the neuronal composition, both in terms of total neuron count (not shown) and relative counts for individual layers (Figure 7D1). However, conicality also affected the density of connections in local, layer-specific subnetworks (Figure 7D2), that is, a measure that is normalized against neuron counts.

Modularity in particular, is strongly correlated with column volume and neuronal count (Figure 7—figure supplement 2A,B), which are in turn driven by the column’s height and conicality. This effect is shaped by the network structure and not just by its size, since it is not present for controls where the number of neurons and their connections are maintained but their pairing is assigned at random (i.e. ER-controls). In fact, for ER-controls, modularity is anti-correlated with volume and neuron count (Figure 7—figure supplement 2B). One way in which the results could be artificial is as follows: Modularity is first defined for a partition of neurons into modules; next the partition maximizing the measure is considered. This can only be done heuristically, and it is possible that optimal partitions are more readily found in certain connectomes than others (see also Materials and methods). To rule this out, we focused on the most modular columns (i.e. those with modularity higher than 0.4). For these, we showed that on average random subnetworks on 50% of the neurons maintained the modularity values of the full columns and that this is not an artifact of the algorithm computing the modularity values (Figure 7—figure supplement 2C).

Going beyond the local, purely internal networks, we considered how much each column was innervated by the individual layers of the entire model. We predict a surprisingly strong impact of geometry, for example the ratio of inputs from layer 3 to inputs from layer 6 shifts from almost 2–1 in convex regions to 1–2 in concave regions (Figure 7E1). On the other hand, we consider the total amount of inputs received, in each layer of a column. In layers 4, 5, and 6, neurons had a higher in-degree if they were members of tall columns (i.e. placed at locations of large cortical thickness). For layers 1, 2, and 3 the trend was weakened or nonexistent (Figure 7E2). The notion of in-degree can be generalized to the n-dimensional in-degree measuring participation in specific, directed motifs of n+1 neurons (see Materials and methods). While trends differed between individual layers, overall the dependence of generalized in-degree on geometrical measures increased with dimension (Figure 7F1 and F2 ‘full’). This was particularly driven by neurons in layer 6. Curiously, in that layer the sign of the r-values, with respect to conicality, inverted from dimensions 1 and 2 to dimensions above 2, indicating the overall innervation of layer 6 is stronger in convex regions, but the participation in higher-order motifs is stronger in concave regions.

The effect of conicality can be explained by convex / concave regions allocating more or less relative space to individual layers, thereby weighing the contributions of the corresponding subnetworks differently. The effect of column height is twofold: First, taller columns will contain more neurons. Second, a different selection of neuron morphologies will be placed, depending on the column height. Briefly, taller morphologies will be selected for taller columns (see also Figure 2—figure supplement 4). Indeed, the selected morphologies in layers 2–5 were significantly different from the overall composition in columns taller than 2 mm or shorter than 1.4 mm; in layer 6 only the short columns led to significant deviance (Figure 7B, Figure 7—figure supplement 1A). Due to the finite number of morphology reconstructions used there is the risk of an artificial effect: It is possible that in geometrically outlying columns only a small number of non-representative morphologies were placed. Testing this, we found that the morphological diversity overall increased from superficial to deeper layers, as more morphological reconstructions were available for the thicker, deeper layers (Figure 7—figure supplement 1B). Over columns, diversity was relatively uniform. The largest decrease in diversity was observed for short columns and layer 5. In that layer, the least diverse column used around a third of the available reconstructed morphologies (Figure 7—figure supplement 1C). As this still amounted to 909 different morphological reconstruction with significant spread between them, we do not believe that dominance by non-representative morphologies explains the qualitative differences in connectivity.

The complexity of local and mid-range connectivity requires neuronal morphologies

So far, we have analyzed single (thalamo-cortical) pathways making up less than 5% of the synapses in the model, and connectivity at scales that are already readily achievable in electron-microscopic reconstructions. The large size of the model also allowed us to characterize predicted connectivity at cellular resolution, but scales that are experimentally only accessible with regional or voxelized resolution (Oh et al., 2014; Bota et al., 2015; Scannell et al., 1995; Scholtens et al., 2014), thereby bridging the scales as outlined in the introduction. Note that the topological methods in the following sections represent and analyze the connectome as a graph, with neurons as nodes and connections between neurons as directed edges, irrespective of the number of synapses in the connection.

We began by analyzing the global structure of neuron-to-neuron connectivity in the entire model, considering local and mid-range connectivity separately. The topology of synaptic connectivity at single neuron resolution has previously been described in terms of the over-expression of directed simplices (Reimann et al., 2017b). A directed simplex of dimension $n$ (or $n$ -simplex, plural $n$ -simplices) is a neuron motif of $n + 1$ neurons that are connected in a purely feed-forward fashion, with a single source neuron sending connections to all others, a single sink neuron receiving connections from all others, and the connections between non-sink neurons forming an $n - 1$ -simplex (Figure 8A1 inset; Figure 8—figure supplement 2A). In particular, 0-simplices correspond to single nodes, and 1-simplices to directed edges. Simplex counts of different dimensions in a network provide a metric of network complexity and can be used to discern their underlying structure (Kahle, 2009; Curto et al., 2013; Giusti et al., 2015). Regarding function, high-dimensional simplices have been demonstrated to shape the structure of spiking correlation between neurons and membership in functional cell assemblies (Reimann et al., 2017b; Ecker et al., 2024a).

Figure 8 with 2 supplements see all

Download asset Open asset

Global connectivity structure local vs mid-range.

(A1) Simplex counts of the local connectivity network and several types of random controls (see text). Examples of $d$ -simplices for $d =$ 1, 2, and 3. Intuitively, an $n$ -simplex is formed by taking an $n - 1$ -simplex and adding a new node to which all neurons connect (sink). (A2) In orange shades, simplex counts of the mid-range connectivity network and several types of random controls. In gray, simplex counts of the combined network and in blue the sum of the local and the mid-range simplex counts. Inset: local, mid-range and the sum of simplex counts on a linear scale. (B1) Node participation per layer, normalized to add up to one in each dimension. Top/bottom row: Local/mid-range connectivity network. Left to right: Participation as source, in any position, and as a sink of a simplex. (B2) From left to right: Spatial location of the cells in the simplicial $n$ -cores in flat coordinates for all $n$ . Depth coordinates for each connected component in the simplicial cores (cells participating in simplices of maximal dimension) each dot marks the depth of a neuron in that component. Adjacency matrix of the simplicial cores with respect to local (green dots) and mid-range (orange dots) connections. (B3) Number of unique neurons in the simplicial cores present in each position of a simplex (i.e., position along the order from source to target. numbers: total counts, bars: counts relative to the total). (C1) Right: Distribution of path distances between pairs of neurons in the local simplicial core; green: along only local edges; grey: along all edges. Left: Euclidean distances between pairs at a given path distance. Arrows: See C2. (C2) Total Number of neurons in all paths of length 3 between neurons of the local core which are at distance 3 in the combined circuit (black arrow in C1) split by location. Red/purple: at positions 2 and 3 respectively; grey: expected from randomly assigned m-types while maintaining the global distribution. Cross-/stripe-patterned: For pairs at path distances greater/less-or-equal than 3 along local edges only (green arrows in C1).

In line with previous results (Reimann et al., 2017b), we found simplices up to dimension seven in the local connectivity (Figure 8A1, green). The maximal dimension did not increase compared to Markram et al., 2015 even with the larger scale of the present model, in accordance with using the same algorithm for local connectivity. However, the addition of mid-range connectivity did produce a major change. In mid-range connectivity alone, simplices of dimension up to 15 were observed (Figure 8A2, orange). This holds true even though local and mid-range connectivity have roughly the same number of edges (2.1 billion local, 2.5 billion mid-range), indicating that the higher simplex counts are not simply due to a larger number of connections. The simplex counts in the combined network of local and mid-range connections are not simply the sum of the local and mid-range simplex counts (Figure 8A2). They are consistently higher and also attain a higher dimensionality generating motifs of up to dimension 18, indicating a strong structural link between the two systems.

We compared the simplex counts of the model to a range of relevant controls that capture simple anatomical properties, such as the density of connections or cortical layers, but ignore the impact of neuronal morphologies (Figure 8A1,A2). This allowed us to assess the degree to which the neuronal geometry generates the complexity of the network. The control models and the parameters on which they capture were the following: the Erdős–Rényi (ER) model used the overall connection density, the stochastic block model (SBM) used density in m-type-specific pathways, the configuration model (CM) used sequences of in- and out-degrees, and finally the distance block model (DBM) used distance-dependence and neuron locations for individual m-type-specific pathways. See Materials and methods for details.

We found that the structure of connectivity is not only determined by the parameters captured by the controls. The CM control was the closest control for mid-range connectivity, suggesting that the effect of degree is important, as expected from the very long-tailed degree distributions of the mid-range connectome (Figure 8—figure supplement 1A). Nonetheless, this control model still largely underestimates the complexity of the mid-range network. At the same time, the DBM control was the closest control for local connectivity, showing that much of its structure is indeed determined by spatial distance under certain morphological constraints.

Simplicial cores define central subnetworks, tied together by mid-range connections

Next, we investigated in which layers the neurons forming these high-dimensional structures resided by measuring node participation, that is, the number of simplices to which a node belongs, and a measure of the node’s centrality (Sizemore et al., 2018; see Materials and methods for details). We found that in local connectivity, most simplices of dimensions 2, 3, and 4 have their source in layer 3 and their sink in layer 5 Figure 8B1 top. Yet, for dimensions 5 and above, we found a shift towards layer 6, with both sources and sinks found mostly in that layer. In mid-range connectivity, the structures are much more concentrated on layer 5, which contains both source and sink Figure 8B1 bottom. Only for simplices of dimension greater than 8 does layer 4 provide more sources. Taken together, this indicates a robust local flow of information from superficial to deeper layers and within deep layers, with layer 5 forming a backbone of structurally strong mid-range connectivity. Neurons in layer 6 form numerous simplices among themselves with no apparent output outside of layer 6, but they are known to be the source of many cortico-thalamic connections (Shepherd and Yamawaki, 2021) that are not a part of this model.

The model made an explicit distinction between m-types forming outgoing mid-range connections (i.e. projecting m-types) and those that do not (Figure 8—figure supplement 1C, ‘proper sinks’) leading to bimodal distributions of out-and total degrees in the corresponding DBM and SBM controls (Figure 8—figure supplement 1A). In contrast, the actual mid-range network has a unimodal, long-tailed degree distribution, similar to biological neuronal networks (Giacopelli et al., 2021). This further demonstrates that the m-types alone can not capture the specific targeting between neuron groups without their actual morphologies.

One type of connectivity specificity that has been found in biological neuronal networks (Towlson et al., 2013) is the formation of a rich club (Zhou and Mondragon, 2004; van den Heuvel and Sporns, 2011). This is characterized by a rich-club curve, which measures whether high degree nodes are more likely to connect together, compared to the CM control (see Materials and methods). We observe that the local network has a rich-club effect, but it is not stronger than expected in a spatially finite model with distance-dependent connectivity (Figure 8—figure supplement 2B1/2). On the other hand, the mid-range network does not exhibit a rich-club effect, see Figure 8—figure supplement 2B2 bottom, showing that high degree nodes are not central to the network.

When we generalized the rich club analysis to higher dimensions by considering node participation instead of degree (Opsahl et al., 2008; see Materials and methods for details) a different picture emerged (see Figure 8—figure supplement 2B3). We found that the rich club coefficients in the mid-range network increased with dimension, while the opposite happens for the local network. We speculate that this effect persists after normalization with respect to relevant controls. Unfortunately, this cannot be currently verified since generating appropriate random controls for these curves is an open problem currently investigated in the field of random topology (Unger and Krebs, 2024). Nonetheless, this indicates that central nodes strongly connect to each other forming a structural backbone of the network that is not determined by degree alone.

We studied these structural backbones, by focusing on the simplicial cores of the networks, which is a higher dimensional generalization of the notion of network core, which is determined by degree alone (see Materials and methods). The general trend observed is that the local network is distributed, while the mid-range network is highly localized. First in terms of the locations of neurons participating in the cores (Figure 8B2 left; green vs. orange), and their laminar locations (Figure 8B2 right). Second, even though the numbers of neurons in the local and mid-range cores are in the same order of magnitude, the local core has 26 disconnected component while the mid-range core is fully connected (Figure 8B3 bottom). Finally, the mid-range core was successively more strongly localized towards the source position, where only seven unique neurons were the sources of all of the 15,108 highest-dimensional mid-range simplices (Figure 8B3 top).

Finally, we studied the interactions between the local and mid-range networks, specifically, the ability of mid-range connnectivity to form short-cuts between neurons of the local core. Even though the sub-network on the nodes in the local core has multiple connected components (Figure 8B2), when paths through nodes outside the core are allowed, it is fully connected. When only local connectivity was considered, the path distances between core neurons were widely distributed between one and seven with a median of four and a strong dependence on their Euclidean distance (Figure 8C1, green). On the other hand, when mid-range connections were added the maximum path distance drops to five, with a negligible number of pairs at that path distance. The median path distance drops to three, and the dependence of the path distance on the Euclidean distance of the pairs nearly disappears (Figure 8C1, grey). Crucially, even though the number of edges almost doubled with the addition of mid-range connections, the number of direct connections remained negligible; instead paths of length three seem to be dominant for information exchange between neurons.

We therefore studied these paths of length three in the combined circuit in more detail, labeling the neurons along them from position 1 (start-node) to position 4 (end-node). We found that layer 5 (the layer with the highest outgoing edge probability, see Figure 8—figure supplement 1B) is over-represented in nodes in positions 2 and 3 compared to a random assignment of m-types (Figure 8C2). Moreover, the over-representation depends on the path distance between the start-node and end-node within the local network: when the value is larger than three (bars hatched by crosses), the effect is stronger than for pairs at a distance three or less (bars hatched by horizontal bars). This demonstrates and quantifies how neurons in layer 5 act as ‘highway hubs’ providing shortcuts between neurons that are far away from each other in the local circuit.

Discussion

We have presented a model of the non-barrel primary somatosensory cortex of the juvenile rat that represents its neuronal and - particularly - synaptic anatomy in high detail. The model comprises a spatial scale that allows for the study of cortical circuits not only as isolated functional units, but also their interactions along inter-regional connections. It also demonstrates how novel insights that are not readily apparent in disparate data and individual models can be gained when they are combined in a way that creates a coherent whole. Specifically, we were able to make multiple predictions about the structure of cortical connectivity that required integration of all anatomical aspects potentially affecting connectivity. At the scale of this model, anatomical aspects affecting connectivity went beyond individual neuronal morphologies and their placement, and included intrinsic cortical curvature and other anatomical variability. This was taken into account during the modeling of the anatomical composition, for example by using three-dimensional, layer-specific neuron density profiles that match biological measurements, and by ensuring the biologically correct orientation of model neurons with respect to the orientation towards the cortical surface. As local connectivity was derived from axo-dendritic appositions in the anatomical model, it was strongly affected by these aspects.

One type of literature resource our approach to connectivity did not use is experimental measurements of local connection probabilities. While these resources are valuable data for our understanding of cortical processing, they are difficult to use in a biophysically detailed model covering large spatial scales. First, connection probabilities are often identified through a PSP response at the soma (Song et al., 2005; Perin et al., 2011; Thomson et al., 2002). However, for connections formed by low numbers of synapses on distal dendrites the PSP may be attenuated too much to be measurable at the soma (Neher, 1992). The connection may still be relevant, e.g., by collaboratively triggering non-linear dendritic events or affecting plasticity of nearby synapses (Farinella et al., 2014; Iacaruso et al., 2017; Tazerart et al., 2020). These effects are not considered in simplified models, leading to a focus on somatically visible connections only. However, in our case, we want to be able to study these effects with our model. Second, some studies report connection probabilities, but only incidentally, being primarily interested in the physiology of a synaptic pathway (Thomson et al., 2002; Markram et al., 1997; Mason et al., 1991). Hence, they employ sampling strategies that aim to maximize the number of connected pairs encountered. Third, even studies aiming to characterize the anatomy of circuitry can provide contradictory estimates: Jiang et al., 2015 report 0 connections for 150 sampled pairs of PCs in layer 5 of adult mouse V1, while in the proofread portion of the data of Bae et al., 2021 13 out of 154 pairs (8.4%) within 100 µm are connected. Fourth, connection probabilities are often measured by sampling pairs of neurons at inter-soma distances below 100 µm or even only 50 µm (see Table 4.1 of Zhang et al., 2019 for an overview). While connection probabilities drop with distance, in a model covering large spatial dimensions, the number of potentially connected pairs grows with the square of distance. Consequently, a connection probability of, for example 0.15% at 500 µm would represent as many connections as a probability of 15% at 50 µm, and it cannot be rounded to zero. Recently, unbiased connection sampling techniques have been developed that provide estimates at larger distances (Chou et al., 2023), that are likely to solve some of these issues in the future.

Connectivity derived from axo-dendritic appositions alone was insufficient at the large spatial scale of the model, as it was limited to connections at distances below 1000 µm. While we found that it generated the right amount of connectivity within a somatosensory subregion, we combined it with a second algorithm for inter-regional connectivity. The algorithm is parameterized by a combination of overall pathway strength, topographical mapping and layer profiles, which together describe a probability distribution for the segments of mid-range axons, that is, an average axonal morphology, using established concepts. On the dendritic side, the algorithm takes individual neuronal morphologies and their placement into account. We have demonstrated that the sum of local and non-local synapse densities matches the reference. It is important to consider whether such a mixture of two independent algorithms captures potentially important statistical interactions between the local and mid-range connectivity of individual neurons. We note that local connectivity is fully constrained by morphology, and the non-local connectivity is parameterized independently for individual morphological types. The combined connectome therefore captures important correlations at that level, such as stronger and weaker non-local cortico-cortical connections from slender-tufted and thick-tufted layer 5 PCs, respectively. When the model was created, interactions beyond this level had not yet been described because local axon reconstructions are reconstructed from slices, which prevents the possibility of obtaining long range information. Even in vivo staining still only reconstructs axons within a region (Buzás et al., 2006). On the other hand, long-range axon reconstructions are obtained through whole brain staining, which has low accuracy for local connectivity (Winnubst et al., 2019). Analysis of new EM datasets, such as a characterization of distance-dependent targeting of excitatory/inhibitory neurons by the axons of L5 thick-tufted pyramidal cells (Bodor et al., 2023) using the data of Bae et al., 2021 could be incorporated in future models. Finding non-random correlations between local and non-local connections would require a strong null model to compare measurements to and our model can serve as that.

Recreating electron-microscopic analyses of connectivity in silico, we could then predict limitations of predicting connectivity from neuron placement and morphology and characterize the additional mechanisms shaping connectivity. Conceptually, a number of mechanisms determine the structure of synaptic connectivity, which we will list from general and large-scale to specific and micro-scale: First, large-scale anatomical trends over hundreds of µm, given by non-homogeneous (e.g. layered) soma placement and broad morphological trends (e.g. ascending axons). These are trends that can be captured by simple distance or offset-dependent connectivity models (Gal et al., 2020). Second, small-scale morphological trends captured by axonal and dendritic morphometrics such as branching angles and tortuosity. These are trends that require the consideration of individual morphologies and their variability instead of average ones. Third, the principle of cooperative synapse formation (Fares and Stepanyants, 2009) formalizing an avoidance of structurally weak connections. Fourth, any type-specific trends not captured by neuronal morphologies, such as local molecular mechanisms and type-specific synaptic pruning. Fifth, non-type specific synaptic rewiring, for example through structural plasticity. We will refer to these aspects as (L)arge-scale, (S)mall-scale, (C)ooperativity, (T)ype-specificity and (P)lasticity respectively, and we argue that for the explanation of the connectome, the more general explanations should be exhausted first, before moving on to more specific ones. Also note that this classification is not considering the underlying developmental causes, but instead associates anatomical and non-anatomical predictors with the structure of the connectome. (L) has been shown to accurately predict large-scale connectivity trends, giving rise to Peters’ rule (Peters and Feldman, 1976; Garey, 1999, although the exact meaning of Peters’ rule is debated, see Rees et al., 2017). It has been demonstrated that (L) alone does not suffice to explain non-random higher order trends in excitatory connectivity, while the combination of (L,S,C) does (Gal et al., 2017; Reimann et al., 2017a; Gal et al., 2020; Udvary et al., 2022). More specifically, Reimann et al., 2017b found that (L,C) explains overexpression of reciprocal connectivity in cortical circuits, but (L,S,C) is required to match the biological trend for clustered connectivity, and Udvary et al., 2022 found (L,S) recreates non-random triplet motifs. Billeh et al., 2020 combined (L) with highly refined (P)-type rules, demonstrating great functional match of the resulting model.

Our comparison to the results of Motta et al., 2019 cannot capture the role of (L), as the scale of the considered volume is too small. But it demonstrates that some non-random targeting trends can be explained by (S,C) and highlights the importance of (C). Further, the overall lower specificity of excitatory axon fragments indicates that (T) may not play a role for them. Similarly, the comparison to Schneider-Mizell et al., 2024 shows that (T) is not required to explain the connectivity of ‘distal-targeting’ (i.e., Sst+) neuron types. Conversely, for ‘perisomatic-targeting’ types (i.e., PV+ basket cells) (T) was required to match the distribution of postsynaptic compartment types. Additionally, the number of potential synapses remaining after applying an optimized (T)-type mechanism was too large to be sustained by the axons, implying a crucial role of (P) in reducing it further. This is in contrast to ‘inhibitory-targeting’ neurons, where an optimized (T)-type pruning lowered the synapse count so much that no space for (P) remained. For the ‘sparsely-targeting’ neurons of Schneider-Mizell et al., 2024 we developed two competing hypotheses, one predicting no role for (T) and 70% of their synapses volumetrically transmitting, that is, without clear postsynaptic partner, and the other predicting a role for (T) similar to the perisomatic-targeting neurons. We summarize our predictions in Table 2.

Table 2

Anatomical, morphological and other aspects affecting connectivity and our predictions for their relevance for efferents of different neuron types.

See main text for an explanation of the individual aspects. The signs in the table indicate whether we predict a certain aspect to be not relevant (-), relevant (+), or highly relevant (++) for a given neuron type.

Axon type	(L)arge-scalemorphological	(S)mall-scalemorphological	(C)ooperativity	(T)ype-specificity	(P)lasticity
EXC.	+	++	++	-	+
BCs	+	+	+	+	+
Sst+	++	+	+	-	+
VIP+	+	-	+	++	-
NGC and L1	++	+	-	+/-	+

In addition to targeting specificity, we predict the following: First, we were able to predict the effect of cortical anatomical variability on neuronal composition by increasing or decreasing the space available for individual layers. If we assume that each layer has a given computational purpose (Felleman and Van Essen, 1991), then this may have functional consequences. It is possible that cortical circuits compensate for this effect, either anatomically (e.g. using different axon or dendrite morphologies), or non-anatomically. Either case would imply the existence of an active mechanism with the possibility of malfunction. Alternatively, function of cortical circuits is robust against the differences in wiring we characterized. This can be studied either in vivo, or in silico based on this model.

Second, we predicted constraints on the emergence of cortical maps from the anatomy of thalamo-cortical innervation, that is, from the combination of: shape and placement of cortical dendrites, the specific layer pattern formed by thalamo-cortical axons, and their horizontal reach. We demonstrated how differences in these parameters may affect the distances between neuron pairs with similar stimulus preferences. At the lower end, very different preferences are supported even between neighboring neurons. At the higher end, neurons further than ∼350 µm apart are likely to sample from non-overlapping sets of thalamic inputs. Other mechanisms will ultimately affect this - chief among them structural synaptic plasticity - this can be thought of as the ground state that plasticity is operating on.

Third, we predicted the topology of synaptic connectivity with neuronal resolution at an unprecedented scale, that is, combining local and mid-range connectivity. We found that the structure of neither can be explained by connection probabilities, degree-distributions, or distance-dependence alone, not even when individual pathways formed between morphological types are taken into account. We expect the mid-range network to have a small-world topology, but computing small-world coefficients for a network of this size is infeasible (see Materials and methods).

Fourth, we predict that the mid-range connectivity forms a strong structural backbone distributing information between a small number of highly connected clusters. The paths within the mid-range network strongly rely on neurons in layer 5 and form short-cut paths for neurons further away than ∼2 mm; for smaller distances, local connectivity provides equivalent or shorter paths.

All these insights required the construction of an anatomically detailed model in a three-dimensional brain atlas; additionally, most of them required the large spatial scale we used. Their functional and computational implications are difficult to predict without also considering data on neuron activity. To that end, an accompanying manuscript describes our modeling of the physiology of neurons and synapses, along with a number of in silico simulation campaigns and their results. Additionally, more insights were gained from using the model in several publications: In Ecker et al., 2024b; Ecker et al., 2024a; Egas Santander et al., 2024a it has been demonstrated that non-random higher order structure affects functional plasticity, assembly formation and the reliability and efficiency of coding of subpopulations. Moreover, two of these insights were validated against electron microscopy data. Furthermore, in Reimann et al., 2024 the model served as an important null model to compare an electron microscopic connectome to, enabling the discovery of specific targeting of inhibition at the cellular level.

This iteration of the model is incomplete in several ways. First, several compromises and generalizations were made to be able to parameterize the process with biological data, most importantly, generalizing from mouse to rat. We note that such mixing is the accepted state-of-the-art in the field, and also neuroscience in general. All 19 data-driven models of rodent microcircuitry listed in Figure 2 of the recent review of Ramaswamy, 2024 conduct some sort of mixing, including the very advanced mouse V1 model of the Allen Institute (Billeh et al., 2020). In this iteration we focused on investigating the general features of the (multi-region) mammalian cortex, e.g., high-order motifs, connected by L5 neurons across subregions or the effect of curvature on connectivity. In the future, more specific aspects of different cortical regions could be investigated using more specific data sources. This would lead to different versions of the model for different regions that can be compared and contrasted using the various structural metrics we developed in this work. Further, we made a number of assumptions about the biological systems that we explicitly list in Supplementary file 2. Should improved and more specific data sources, such as EM-based connectomes or whole-brain neuron reconstructions, become available in the future, our modeling pipepline is designed to readily utilize them for refinement. There are also known limitations to the actual model building steps. While we enforce cell density profiles along the depth-dimension, densities are assumed to be otherwise uniform within a given subregion. For additional structure, the algorithms will need to be updated; this will be required for example, for the modeling of the barrel field. Similarly, during volume filling, orientation and placement of neuron morphologies is considered only along the depth-axis. Modeling of the barrel system will also require the consideration of orientation and location with respect to nearby barrel centers. Furthermore, we assumed uniform relative layer thickness throughout the modeled region, due to sparsity of the required data. However, given that Yusufoğulları et al., 2015 has demonstrated there are substantial differences in layer thickness between rat hindlimb and barrel field regions and Narayanan et al., 2017 and Wagstyl et al., 2020 found lower but still significant differences within developed rat barrel cortex and human somatosensory cortex respectively, the model pipeline should be updated once sufficient data becomes available.

This demonstrates that detailed modeling requires constant iteration, as a model can never be proven to be right, only to be wrong. Thus, we made the tools to improve our model also openly available (see Data and Code availability section). Figure 1 provides a technical overview of individual modeling steps, the software tools involved and the duration required to run them, to serve as a guide to interested contributors.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Additional information
Software, algorithm	Model loading and interaction	https://doi.org/10.5281/zenodo.8026852
Software, algorithm	Model analysis	https://doi.org/110.5281/zenodo.8016989
Software, algorithm	Topological analysis of connectivity	Egas Santander et al., 2024b; https://github.com/danielaegassan/connectome_analysis	see also Egas Santander et al., 2024a; version 0.0.1
Software, algorithm	atlas-splitter	Gevaert et al., 2024a; https://github.com/BlueBrain/atlas-splitter/tree/v0.1.5	Release 0.1.5
Software, algorithm	atlas-densities	Gevaert et al., 2024b; https://github.com/BlueBrain/atlas-densities/tree/v0.2.5	Release 0.2.5
Software, algorithm	Brainbuilder	Povolotsky et al., 2024; https://github.com/BlueBrain/brainbuilder/tree/v0.20.0	Release 0.20.
Software, algorithm	NeuroR	Arnaudon et al., 2024; https://github.com/BlueBrain/NeuroR/tree/v1.7.0	Release 1.7.0
Software, algorithm	NeuroC	Gevaert et al., 2024c; https://github.com/BlueBrain/neuroc/tree/v0.3.0	Release 0.3.0
Software, algorithm	placement-algorithm	Povolotsky et al., 2023; https://github.com/BlueBrain/placement-algorithm/releases/tag/placement-algorithm-v2.4.0	Release 2.4.0
Software, algorithm	Appositionizer	Wolf and Berchet, 2024; https://github.com/BlueBrain/appositionizer/releases/tag/v1.0.0	see also Kozloski et al., 2008; release 1.0.0
Software, algorithm	Functionalizer	Wolf et al., 2024; https://github.com/BlueBrain/functionalizer/releases/tag/v1.0.0	Release 1.0.0
Software, algorithm	Projectionizer	Gevaert and Herttuainen, 2022; https://github.com/BlueBrain/projectionizer/releases/tag/projectionizer-v3.0.0.dev0	Release 3.0.0.dev0
Software, algorithm	Connectome-Manipulator	Pokorny et al., 2024b; https://github.com/BlueBrain/connectome-manipulator/tree/v1.0.0	see also Pokorny et al., 2024a; release 1.0.0

Method details

An overview of the model building steps, tools used and approximate durations can be found in Figure 9.

Figure 9

Download asset Open asset

Share this article

Cite this article

Overview of the model building and analysis workflow.

References to publications of input data and methods employed for individual modeling steps.

Anatomy of the model.

Volumetric anatomy of the model.

Intrinsic connectivity as union of local and mid-range connectivity.

Recreating an electron-microscopic specificity analysis In the top part we analyzed synapses, axons and dendrites in an approximately cubic volume in L4, emulating the techniques of Motta et al., 2019.

Anatomy of thalamic innervation.

A1, A2: Parcellation of the modeled volume into 230 µm radius columns.

Global connectivity structure local vs mid-range.

Anatomical, morphological and other aspects affecting connectivity and our predictions for their relevance for efferents of different neuron types.

Overview of the building workflow and tools used.

Author details

Michael W Reimann

Contribution

Contributed equally with

For correspondence

Competing interests

Additional information

Sirio Bolaños-Puchet

Contribution

Contributed equally with

For correspondence

Competing interests

Additional information

Jean-Denis Courcol

Contribution

Contributed equally with

Competing interests

Additional information

Daniela Egas Santander

Contribution

Contributed equally with

Competing interests

Additional information

Alexis Arnaudon

Contribution

Competing interests

Benoît Coste

Contribution

Competing interests

Fabien Delalondre

Contribution

Competing interests

Thomas Delemontex

Contribution

Competing interests

Adrien Devresse

Contribution

Competing interests

Hugo Dictus

Contribution

Competing interests

Alexander Dietz

Contribution

Competing interests

András Ecker

Contribution

Competing interests

Cyrille Favreau

Contribution

Competing interests

Gianluca Ficarelli

Contribution

Competing interests

Mike Gevaert

Contribution

Competing interests

Joni Herttuainen

Contribution

Competing interests

James B Isbister

Contribution

Competing interests

Lida Kanari

Contribution

Competing interests

Daniel Keller

Contribution