Understanding Neocortical Micro- and Mesocircuitry: A Deep Dive into Anatomy

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. 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. 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 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.