1 Introduction

Neuroscience aims to characterize the dynamics of the brain and understand how they emerge as a product of anatomy and physiology. This is challenging however, due to the complexity of the brain’s multi-scale organization and its study in situ. Large-scale, data-driven, biophysically-detailed models (Markram et al., 2015; Billeh et al., 2020) offer a tool which can complement laboratory investigation. This “bottom-up” approach combines detailed constituent models into a 3D model of a given brain volume, thus consolidating disparate, multi-scale data sources. Such models aim to provide a general tool for multi-scale investigation; access and manipulation of any model component enables predictions about how these components combine to shape emergent dynamics. With increasing biological data, model refinement and validation must be continuous. This necessitates 1) open sourcing of models with high quality software tools for iterative, community-driven refinement and validation, and 2) rigorous methodology for model parameterization and validation, such that model iterations can be meaningfully compared. Here and in our companion paper (Reimann et al., 2022a) we present our methodology for parameterizing the anatomy and physiology of a cortical model, and validating its emergent dynamics. This enabled predictions about how cortical activity is shaped by high-dimensional connectivity motifs in local and mid-range connectivity, and spatial targeting rules by inhibitory subpopulations. Elsewhere, this has allowed the presented model to be used to study the formation of cell assemblies (Ecker et al., 2024), functional synaptic plasticity (Ecker et al., 2023), the role of non-random connectivity motifs on reliability (Egas Santander et al., 2024), the composition of high-level electrical signals such as the EEG (Tharayil et al., 2024), and how spike sorting biases population codes.

Specifically, we built and validated a model of the entire non-barrel primary somatosensory cortex (nbS1) comprising eight subregions. Whilst our previous data-driven, biophysically detailed model (Markram et al., 2015) provided insights at the scale of a single cortical column (Reimann et al., 2013, 2017, 2022b; Nolte et al., 2019, 2020; Newton et al., 2021), the new model is ∼140 times larger and to our knowledge offers the first simulations of in vivo-like spontaneous and stimulus-evoked activity in a biophysically detailed cortical model with interregion connectivity. In the companion paper, we introduce the anatomical model (Fig. 1, Step 1), describing how neuron morphologies were placed within an atlas-based geometry and connected through local and mid-range synapses. Here, we describe our improved techniques to model and validate the electrical properties of neurons and synapses (Figure 1, Steps 2 & 3 (Reva et al., 2023; Barros-Zulaica et al., 2019)), and to compensate for input from missing brain areas (Figure 1, Steps 4). These improvements enabled enhanced validation of emerging in vivo-like activity (Figure 1, Steps 5), including the reproduction and extension of five published studies in rodent sensory cortex under a single in vivo-like regime (Figure 1, Steps 6).

Figure 1:

Validating that network activity emerges from the same interactions driving in vivo dynamics, requires 1) an approach to tackle the large parameter space without overfitting and 2) comparison of emerging dynamics with laboratory experiments. With respect to 1), we strongly prefer parameterization with directly measured quantities over fitting parameters to yield the correct emerging activity. Additionally, where fitting is applied, we adhere to a principle of compartmentalization of parameters. That is, once a parameter has been parameterized at one biological level, it is no longer a free parameter at a higher level. For example, the maximal conductance of a synapse is fit to biological amplitudes of postsynaptic potentials, but is then never updated in the process of reaching in-vivo like population activity. If a connection is valid on the single-cell level, its contribution at the population level should be equally valid. As a result, the emerging in-vivo like activity is the consequence of only 10 free parameters representing the extrinsic input into the individual layers from other brain regions. These parameters were fitted using a novel methodology presented here and were in accordance with the mean number of missing synapses for each population (Felleman and Van Essen, 1991; Harris et al., 2019; Gao et al., 2022). With respect to 2), we conducted the following validations: Spontaneous activity reproduced layer-wise in vivo firing rates (or a specified proportion of in vivo firing rates to account for in vivo recording bias (Wohrer et al., 2013)), varied along a spectrum of asynchronous to synchronous activity, exhibited spatially structured fluctuations and produced long-tailed firing rate distributions with sub 1 Hz peaks as in vivo (Wohrer et al., 2013; Buzśaki and Mizuseki, 2014). Under the same parameterization, the model also reproduced precise millisecond dynamics of layer-wise populations in response to simple stimuli and demonstrated selective propagation of stimulus-evoked activity to downstream areas. Again under the same parameterization, we reproduced and extended more complex experiments through accurate modeling of targeted optogenetic stimulation and lesions. Importantly, we highlight where emergent activity shows discrepancies with in vivo activity, to guide future data-driven model refinement.

The model generated a number of predictions (Table S1), including about the role of different layers in driving layer 2/3 stimulus responses and how inhibitory interneuron types encode contrast, synchronous and rate-coded information. Additionally, with access to the full structural connectome and tools for precisely editing it (Pokorny et al., in prep), we were able to make predictions about the relationship between structure and function. For example, we predict that an increase in the prevalence of non-random connectivity motifs towards deeper layers leads to a matched increase in spiking correlations, and that subregions more strongly innervated by mid-range connectivity have higher correlated activity locally. Additionally, we utilized new analysis of inhibitory connectivity (Schneider-Mizell et al., 2023) in the MICrONS electron microspcopy dataset (MICrONSConsortium et al., 2021) to make predictions about the role of specific targeting rules by different inhibitory neuron types. Particularly, we compared activity in the original connectome to a rewired connectome with increased perisomatic targeting by PV+ neurons, and increased targeting of inhibitory populations by VIP+ neurons. Inhibitory populations were more strongly inhibited (increasingly towards central layers) and required more non-local drive to reach the firing rates observed in vivo. Evoked responses increased and decreased in superficial and deeper excitatory populations respectively, suggesting layerspecific roles of the more specific inhibitory targeting.

To provide a framework for further studies and integration of experimental data, a sizable part of the model and simulation tools are made available (Figure 1, Step 7). The detailed modeling approach provides a one-to-one correspondence with most types of experimental data, allowing different datasets to be readily integrated. Due to the incredible speed of discovery in neuroscience an integrative model will always be lagging behind the latest available data. We believe the solution is to provide a scientifically solid, validated model with clearly characterized strengths and weaknesses, along with the tools to advance or customize it for individual projects. We have therefore also made our tools for building and improving the model openly available: www.github.com/BlueBrain.

2 Results

The companion paper describes the full anatomical nbS1 model containing 4.2M morphological neuron models, connected through 9.1B local synapses and 4.1B mid-range synapses (Reimann et al., 2022a). Each neuron is modeled as a multi-compartmental model belonging to one of 60 morphological types (m-types; Figure 2A1), and is either an instance or statistical variant of an exemplar in a pool of 1,017 morphological reconstructions. Neurons were placed (orientated towards the surface) within an atlas-based geometry (Figure 2A2) based on estimated layer-wise density profiles of different morphological types. Local connectivity is based on axo-dendritic overlap with neighbouring neurons, whilst mid-range connectivity combines data on interregion connectivity and laminar innervation profiles. Thalamocortical afferents were also modeled based on laminar innervation to the barrel cortex. To simulate emergent activity it is necessary to model and validate the electrical properties of these neurons and synapses (Figure 1 Steps 2 & 3). These modeling steps are based on published methods and data sources (Table 1), and are summarized first. The method for compensating for missing synapses and the remaining results are then described (Figure 1 Steps 4-6).

Figure 2:

Table 1:

2.1 Modeling and validation of neuron physiology

Electrical properties of single neurons were modeled by optimising ion channel densities in specific compartment-types (soma, axon initial segment (AIS), basal dendrite, and apical dendrite) (Figure 2B) using an evolutionary algorithm (IBEA; Van Geit et al., 2016) so that each neuron recreates electrical features of its corresponding electrical type (e-type) under multiple standardized protocols. The methodology and resulting electrical models are described in Reva et al. (2023) (see also Methods). Electrical features included firing properties (e.g. spike frequency, interspike interval), action potential waveforms (e.g. fall and rise time, width) and passive properties (e.g. input resistance). The optimisation was performed for a subset of neuron models. The resulting ion channel densities were generalised to other neuron models of the same e-type. For excitatory neurons such generalization was only made to neurons within the same layer. The resulting electrical activity of each neuron was validated against the corresponding electrical features, including the characteristic firing properties of the 11 e-types (Figure 2C). Model neurons with a mean zscore (over the electrical features) further than two standard deviations from the experimental mean were discarded.

For each of the 60 morphological types (m-types), the corresponding fractions of etypes were determined from experimental data as in Markram et al. (2015), resulting in 208 morpho-electrical types (me-types; Figure 2D). Correspondence with predominant expression of biological markers PV (parvalbumin), Sst (somatostatin) or 5HT3aR (serotonin receptor 3A) was determined based on me-type, as previously done in Markram et al. (2015) (Figure 2D). Finally, for layer five thick-tufted pyramidal cell (L5 TTPC) morphologies, we found that dendritic electrical features, namely the attenuation of back propagating action potentials and excitatory postsynaptic potentials (EPSPs), reproduced experimental measurements (Figure 2E, Reva et al., 2023).

The required data to constrain region-specific me-type distributions was not available when building the model; to our knowledge Yao et al. (2023) is the first study that explored the question systematically (in mouse cortex). The atlas-based geometry does however impose heterogeneity in morphologies that e-types are paired ith: only morphologies that are in line with the region-specific local curvature of the region are used throughout the volume. We demonstrate the effect of geometry on morphological composition in the companion paper.

2.2 Modeling and validation of synaptic physiology

Synapses are modeled with a stochastic version of the Tsodyks-Markram model (Tsodyks and Markram, 1997; Markram et al., 1998; Fuhrmann et al., 2002; Loebel et al., 2009), upgraded to feature multi-vesicular release as described in Barros-Zulaica et al. (2019) and Ecker et al. (2020). The model assumes a pool of available vesicles that is utilized by a presynaptic action potential, with a release probability dependent on the extracellular calcium concentration ([Ca²⁺]_o; Ohana and Sakmann, 1998; Rozov et al., 2001; Borst, 2010). Additionally, single vesicles spontaneously release with a low frequency as an additional source of variability. The utilization of vesicles leads to a postsynaptic conductance with bi-exponential kinetics. Short-term plasticity (STP) dynamics in response to sustained presynaptic activation are either facilitating (E1/I1), depressing (E2/I2), or pseudo-linear (I3). E synaptic currents consist of both AMPA and NMDA components, whilst I currents consist of a single GABA_A component, except for neurogliaform cells, whose synapses also feature a slow GABA_B component. The NMDA component of E synaptic currents depends on the state of the Mg²⁺ block (Jahr and Stevens, 1990), with parameters fit to cortical recordings from Vargas-Caballero and Robinson (2003) by Chindemi et al. (2022).

The workflow for determining a dense parameter set for all synaptic pathways, starting with sparse data from the literature, is described for the use case of hippocampal CA1 in Ecker et al. (2020) and briefly in the Methods. We combined data sources used in Markram et al. (2015), with a large number of recent data sources (Qi and Feldmeyer, 2016; Barros-Zulaica et al., 2019; Yang et al., 2020, 2022). The resulting pathway-specific parameters are listed in Tables S2, S3, S4, the most common short-term dynamics are depicted in Figure 3A1-2, and the assignment of STP profiles to different pathways are shown in Figure 3A3. Postsynaptic potential (PSP) amplitudes and their coefficient of variation (CV; std/mean) closely matched their biological counterparts (r = 0.99, n = 27; Figure 3B1; Table S5 and r = 0.63, n = 10; Figure 3C1; Table S6, respectively). The dense parameter set also allowed prediction of PSP amplitudes and CVs for all cortical pathways (Figure 3B2, C2). The frequency of miniature postsynaptic currents (mPSCs) were also in line with in vitro measurements (r = 0.92, n = 5; Figure 3D; Table S7).

Figure 3:

The anatomical model includes fibers from the thalamus, based on fibres projecting to the barrel cortex from the ventral posteriomedial (VPM) and posteriormedial (POm) thalamic nuclei. These fibres make synaptic contacts within a radius of the fiber probabilistically based on laminar innervation profiles (see Reimann et al., 2022a; Figure 3E). Compared to the previous model (Markram et al., 2015), POm projections are added and the physiology of synapses from VPM improved (see Methods; Figure 3E, F). Latencies of layer-wise EPSPs increase with distance from the thalamus (Figure 3G). Additionally, thalamocortical EPSP amplitudes normalized relative to a single population were compared to normalized EPSPs in response to optogenetic stimulation targeting bundles of thalamic fibers (Sermet et al., 2019). This provided contrasting insights, however. For example, whilst VPM to L6 I EPSPs match the initial validation data (Figure 3F), VPM to L6 PV+ responses appear too strong relative to other populations (Figure 3H1). The results suggest that the model’s POm to L5 E pathway is too weak, when compared to other POm to E and all POm to PV+ pathways (Figure 3H2, right).

2.3 Defining subvolumes and populations

To enable smaller simulations and targeted analyses, we defined standardized partitions comprising all neurons (and their connections) contained in a subvolume of the full nbS1 model. We decomposed the model into full-depth hexagonal prisms of a certain diameter that we call hexagonal subvolumes. These are slightly curved, following the geometry of the cortex, and have an intact layer structure. When the diameter of the hexagons is 520 µm, comparable to the size of the single cortical column model from (Markram et al., 2015), we call these subvolumes columns. Taken together, a central column and the six columns surrounding it define a seven column subvolume (Figure 4A). Additionally, when we separately analyze the E and I neurons in different layers, we refer to them as neuron populations of the model or a subvolume. Groups of I neurons that predominantly express either PV, Sst or 5HT3aR markers, are referred to as inhibitory subpopulations.

Figure 4:

2.4 Compensating for missing brain regions

Initial simulations produced no activity, and it was necessary to compensate for missing excitatory input from neurons external to the model. We estimated that the number of missing synapses for the full nbS1 and seven hexagon subvolume as ∼ two and seven times the number of internal synapses, based on Oh et al. (2014). We model the effect of these missing synaptic inputs explicitly as time-varying and statistically independent somatic conductance injections, using Ornstein-Uhlenbeck (OU) processes that mimic aggregated random background synaptic inputs (Destexhe et al., 2001; Figure 4B1; see also Discussion). The mean (OU_µ) and standard deviation (OU_σ) are defined as a percentage of an individual cell’s input conductance at its resting membrane potential. The choice of OU_µ and OU_σ for different populations and [Ca²⁺]_o (Figure 4B2), determine the model’s emergent dynamics and mean firing rates (Figure 4B3). Initial simulations showed that using the same value of OU_µ and OU_σ for all neurons makes some populations highly active whilst others are silent, or produces network-wide bursts. Finding population-specific OU parameters which produce stable in vivo-like activity is challenging, however, due to the non-linearity and computational cost of simulations.

We developed an algorithm that efficiently calibrates population-specific OU parameters, for different values of [Ca²⁺]_o, which affects the strength and reliability of synaptic transmission (see “Synaptic Physiology”). Additionally, it fixes the ratio OU_σ/OU_µ, denoted R_OU, across all neurons to specific values representing the amount of noise in the extrinsic inputs (Figure 4B2). As extracellularly-derived firing rates are known to be severely overestimated to an unknown degree (Olshausen and Field, 2006; Wohrer et al., 2013; Buzśaki and Mizuseki, 2014), the algorithm finds parameters which produce realistic inter-population firing rate ratios at different global levels of activity (Figure 4B3). Specifically, it targets mean spontaneous firing rates of neuron populations that are a constant fraction P_FR of in vivo reference values, for 10 P_FR values between 0.1 and 1. We refer to [Ca²⁺]_o, R_OU and P_FR as meta-parameters. The initial calibration was made for the seven hexagonal subvolume, and later generalised to the full nbS1 model.

The algorithm was first run for [Ca²⁺]_o = 1.1 mM, R_OU = 0.4. Target firing rates for the 10 P_FR values were reached (Figure 4C) after five iterations of 10 simulations of 6.5 seconds (Figure 4D, left upper). Only three iterations were needed when calibration began using previously calibrated parameters for different R_OU and [Ca²⁺]_o (Figure 4D, right upper). This allowed us to parameterize extrinsic input for a range of combinations of the three meta-parameters (Figure 4D, lower & 4E, left) within biologically relevant ranges ([Ca²⁺]_o: 1.05 mM - 1.1 mM (Jones and Keep, 1988; Massimini and Amzica, 2001; Amzica et al., 2002; Gonzalez et al., 2022); P_FR: 0.1 - 1.0; R_OU : 0.2 - 0.4 (Destexhe et al., 2001)). We found that the resulting mean conductance injections were highly population-specific (Figure 4E, centre). The range of OU_µ values between populations was on average 5.2 times higher than the range of OU_µ values for a single population across meta-parameter combinations (Figure 4E, right). This was expected, as the missing mid-range innervation is known to have specific and heterogeneous laminar profiles (Felleman and Van Essen, 1991; Harris et al., 2019; Gao et al., 2022). We estimated the number of missing synapses per neuron assuming a total density of E synapses of 1.1 synapses/µm, based on mean spine densities (Larkman, 1991; Datwani et al., 2002; Kawaguchi et al., 2006), and sub-tracting the number of synapses present in the model. We confirmed a strong correlation between this measure and required conductance injection (Figure 4F).

2.5 In vivo-like spontaneous activity

Histograms of spontaneous activity for two meta-parameter combinations are shown in Figure 4G1 (and for all combinations in Video 1). The two combinations have the same population firing rates (P_FR = 0.9) but the first simulation uses a higher R_OU and lower synaptic transmission reliability ([Ca²⁺]_o). In the second simulation, we found larger fluctuations, correlated between E and I populations across layers. Overall, increasing [Ca²⁺]_o and P_FR, and decreasing R_OU, increases the amplitude of correlated fluctuations (Video 1; Figure 4G2) and increases power at lower frequencies (Figure 4G3). This supports the notion that decreasing R_OU or increasing [Ca²⁺]_o shifts the model along a spectrum from externally (noise) driven to internally driven correlated activity.

We confirmed that spontaneous firing rate distributions were long-tailed with peaks below 1 Hz for all of the meta-parameter combinations, in line with in vivo activity (Figure S3). In the horizontal dimensions, under higher correlation regimes, fluctuations are spatially coordinated and global within the central hexagon of the simulated subvolume (Figure 4H; Video 2). In depth, the size of fluctuations and level of correlated activity increases from supragranular to infragranular layers (Figure 4G1, I, left; Video 1), suggesting that activity in deeper layers of the model is more internally driven. We characterized the effect of recurrent connectivity by measuring the layer-specific decrease in firing rates when neurons in the model are disconnected from each other. As expected, we found that the effect (quantified by the ratio between connected and disconnected firing rates: R_C/U) increased from suprato infragranular layers (Figure 4I, right).

2.6 High-dimensional connectivity motifs predict layer-wise spontaneous activity

By analysing the underlying network structure, we observed a corresponding gradient in the topology of intra-layer synaptic connectivity from supra- to infragranular layers. We calculated the mean node participation of neuron populations in various dimensions (see Methods; Figure 4J, top). For dimension one, this is simply the degree of a neuron; for dimensions above one, this generalizes the notion to counting participation in dense, directed motifs of increasing size (directed simplices; Reimann et al., 2017). We found that correlations of this measure with the ratio of connected and disconnected firing rates increased with dimension (Figure 4J, bottom, 100% of data), indicating the importance of large directed simplices in shaping activity. Curiously, the correlations for higher dimensions were driven by a small number of neurons with a very high node participation. This was indicated by a sharp drop of correlation when neurons above a given value were excluded (Figure 4J, bottom). This demonstrates high variability between neurons, in line with biology, both structurally (Towlson et al., 2013; Nigam et al., 2016) and functionally (Wohrer et al., 2013; Buzśaki and Mizuseki, 2014).

2.7 Stimulus-responses reproduce in vivo dynamics with millisecond-scale precision

We compared stimulus-responses with in vivo barrel cortex responses to both single whisker deflections and active whisker touches under anaesthetized (Reyes-Puerta et al., 2015) and awake states (Yu et al., 2019) respectively. While the model is of non-barrel somatosensory regions, this nevertheless provides validations of overall excitability and the laminar structure of responses reflecting general trends of cortical processing. We activated a percentage (F_P) of the VPM fibers projecting to the central column of the seven column subvolume (Figure 5A1). For each selected thalamic fiber, spike times were drawn from VPM PSTHs (peristimulus time histogram) recorded in vivo for the two stimulus types (Figure 5A2; Methods). For both stimulus types, a single thalamic stimulus was presented 10 times at 1 s intervals at three intensities (F_P : 5%, 10%, 15%).

Figure 5:

For the parameter combination [Ca²⁺]_o = 1.05 mM, R_OU = 0.4, P_FR = 0.3 and F_P = 10%, E and I populations in each layer show clear responses on each trial (Figure 5B; three trials shown), except for L1. Videos 3 and 4 show all simulated combinations of the four meta-parameters. Out of 72 parameter combinations, 21 passed an initial assessment of in vivo similarity to the in vivo data based on latencies and decays (see Methods, Figure 5C, S4). Responses remained localized to within 600 µm of the stimulus location (Figure 5D, Video 5). In vivo, layer-wise I populations precede the corresponding E populations and response latencies increase with distance from thalamus, except for the fast response of L4 I. The latencies and decays are matched by criteria passing simulations (Figure 5C, E, left), except for a longer L4 I latency, and shorter and longer sustained responses of L2/3 E and L4 E. PV+ and Sst+ subpopulations respond reliably in layers 2-6, but the 5HT3aR+ subpopulation responses are weaker and less reliable (Figure S4). Under the active whisker touch paradigm, PV+ latencies precede Sst+ latencies in vivo, and either precede or are simultaneous with Sst+ latencies in silico (Figure 5E, right).

The passing parameter combinations form a contiguous region in the parameter space (Figure 5F), characterized by combinations of low P_FR, [Ca²⁺]_o and F_P with higher R_OU, which displayed weakly correlated spontaneous activity (compare Figure 5F with Figure 4G2). This offers a prediction that true in vivo firing rates are within the range of these lower FRs. The remaining results apply to the criteria passing simulations. We characterized the magnitude of responses by R_E, the ratio between the peak of the evoked response and the pre-stimulus baseline activity level. Assuming that extracellular sampling biases are similar during spontaneous and evoked activity, R_E provides a bias-free metric for comparing response magnitude with in vivo. Decreasing R_OU and P_FR, and increasing [Ca²⁺]_o and F_P led to a global increase of R_E, and results matched in vivo closely for a large part of the parameter space (Figure 5F).

Interestingly, the parameter combination closest to in vivo was central in the parameter space of criteria passing simulations (Figure 5G, left). Values of R_E, normalized by the mean R_E across populations, were remarkably constant (Figure 5G, right), indicating that variation of the parameters provides a global scaling of response magnitude. Moreover, I populations matched the in vivo reference closely, while the responses of L2/3 E were slightly greater, and for L4 E and L5 E slightly weaker. The region of the parameter space where responses match in vivo for the active whisker touch stimulus was the same as for single whisker deflections (Figure S5). Latencies were also similar to in vivo, with the main discrepancy being slower responses in L4 (E and I). With respect to the metaparameters, we found that the best match for the anesthetized state was characterized by lower firing rates (P_FR) and reduced external noise levels (lower R_OU) compared to the awake active whisker touch paradigm (white areas in Figures 5F and S5B). These results fit expectation for anaesthetized and awake states respectively, and provide predictions of their spontaneous firing rates, noise levels and corresponding in silico regimes.

Response sparsity (i.e. the proportion of neurons spiking at least once following a stimulus) varied across populations and meta-parameter combinations, but was around the range of 10-20% reported in vivo in the barrel cortex (Barth and Poulet (2012); Figure S4D). Interestingly, excitatory population responses were sparser than the corresponding inhibitory layer-wise populations across layers. The observed heterogeneous layer-wise population sparsities offer a prediction for in vivo. Moreover, the majority of responding cells spike only once to single whisker whisker deflections (Figure S4E), as reported in vivo (Isbister et al., 2021).

2.8 Validation and prediction through reproduction and extension of complex laboratory experiments using a single model parameterization

Building upon the in vivo-like activity states, we combined various simulation techniques to recreate laboratory experiments from the literature (see Methods). For these experiments, we use the parameter combination: [Ca²⁺]_o = 1.05 mM, R_OU = 0.4 and P_FR = 0.3, as this combination supports realistic in vivo-like responses to the two stimulus types over a range of fiber percentages (F_P).

2.8.1 Exploring the canonical model: layer-wise contributions to L2/3 stimulus-responses

In the canonical model of the cortex (reviewed e.g., in Lübke and Feldmeyer, 2007; Feldmeyer, 2012) information from the thalamus arrives to L4, then propagates to L2/3, from there to L5 (which serves as the main cortico-cortical output layer) and lastly to L6 (which projects back to the thalamus). The coordinated action of all layers of S1 is required for high-level behavioral tasks (Park et al., 2020). As the canonical model is based on the highest layer-wise density of axons, it cannot describe all interactions in the cortex. For example VPM does not only innervate L4, but also the border of L5 and L6, as well as the bottom of L2/3 (Meyer et al., 2010; Constantinople and Bruno, 2013; Sermet et al., 2019).

To study how L4 contributes to the stimulus preference of L2/3 PCs, Varani et al. (2022) used optogenetic inactivation of L4 PCs during whisker stimulation and quantified the changes in the subthreshold response of L2/3 PCs. They found, that the early phase of the subthreshold response significantly differed from the control condition, if the whisker was deflected in either the most or the least preferred direction (see the top and bottom rows of their Figure 5B, C). From this they concluded that both L4 and VPM contribute to the direction tuning of L2/3 PCs. After reproducing their experimental conditions in silico (Figure 6A,B; Figure S6; see Methods) we confirmed that we can reproduce their results, i.e., subthreshold responses of L2/3 PCs decreased, when L4 PCs were inhibited (Figure 6C for preferred direction whisker stimulation; see Methods).

Figure 6:

We then leveraged our in silico setup to study what Varani et al. (2022) could not, because of methodological limitations. In our reading, the authors aimed to test how direct E connections from L4 PCs to L2/3 PCs influence the stimulus representation in L2/3. This connection can not specifically be blocked in vivo, instead (95% of) the L4 PC population is inhibited (as well as some lower L3 PCs). In our setup we could selectively block the connection and found almost the same result (compare Figure 6C and D, left). This extends the conclusion of Varani et al. (2022): L4 PCs contribute to the stimulus preference of L2/3 PCs via direct E connections, and not via disynaptic inhibition.

The authors also discussed studying L5 PCs’ contribution to L2/3 responses (as a large fraction of L5 PC axons terminate in L2/3), but this is infeasible with current mouse lines. Leveraging our model, we found that L5 contributes much less to subthreshold L2/3 traces than L4 (Figure 6D, right). Extending to other presynaptic layers, we found that the contribution of L2/3 is similar to that of L4, whereas inputs from L6 are negligible (Figure 6D, right). Whilst mouse lines targeting L5 PCs might arrive soon (which could validate our prediction), blocking L2/3 connections between L2/3 cells without hyperpolarizing the same L2/3 population seems only achievable in silico.

2.8.2 Contrast tuning by inhibitory interneuron sub-types under optogenetic activation

Shapiro et al. (2022) compared the modulatory effects of optogenetic activation of interneuron subtypes on contrast tuning of neurons in V1. Whilst visual regions differ from somatosensory regions, the study provided an opportunity to explore the generality of interactions between E and I subpopulations. We therefore created a version of the study in silico, presenting spatio-temporally modulated patterns with various contrast levels through the VPM inputs (Figure 6E,F; see Methods).

Overall, we found within the single simulated column 228 PCs and 259 PV+ interneurons with robust contrast tuning of their responses (see Methods). Within spike detection range of a vertically penetrating extracellular electrode (50 µm; Henze et al. (2000); Buzśaki (2004)), we found 15 PCs and 13 PV+ interneurons with robust tuning, which is roughly comparable to the original study (average of 8.2 tuned PCs and 3.0 PV+ neurons per mouse Shapiro et al., 2022). Unlike the original study we found no Sst+ neurons with robust tuning.

Additionally, we found the following results, all in line with (Shapiro et al., 2022): The firing rates of both PCs and PV+ neurons were affected by optogenetic activation of the PV+ subpopulation (Figure 6G), resulting in changes to the contrast tuning curves of these populations (Figure 6H). The main effect of the optogenetic manipulation was an increase in contrast tuning for PV+ and a decrease for PCs. The tuning curves could be accurately fit by sigmoidal functions (r_PV² = 0.995 ± 0.008; r_PC² = 0.990 ± 0.012), and consequently the effect of optogenetic activation could be quantified in terms of changes to the fitted parameters (Figure 6H1 inset, S7; Methods). We found great variability between neurons with increasing trends for parameter m in PV+ and c₅₀ in PCs, and decreasing trends for R_max in PV+ and both m and R_max in PCs. Unlike in (Shapiro et al., 2022), we found in addition an increasing trend of c₅₀ also in PV+ interneurons (Figure S7).

Comparing three abstract mathematical models (Shapiro et al., 2022; Methods), we found that the saturation additive model best captured the effect of the optogenetic manipulation for PV+ neurons (Figure S10A, B1), as in the reference study. The relative contributions of saturation and addition varied strongly between individual neurons (Figure S10C1). At low intensities of the optogenetic manipulation, this distribution matched qualitatively the results of Shapiro et al. (2022), but at higher intensities we observed a shift towards stronger saturation. Combining the saturation additive model with a simplified, population-level model of PC firing (Shapiro et al., 2022; Methods) provided a good description of the observations (Figure S10B2) with relative contributions of saturation and addition that qualitatively match the reference (Figure S10C2).

2.8.3 Encoding of synchronous and rate-coded signals by inhibitory interneuron sub-types

Prince et al. (2021) studied how different I subpopulations represent signals that are encoded as changes in the firing rate or synchronicity of PCs (Figure 6I). We repeated their study in silico (Figure 6J), but using in vivo-like conditions as the authors have suggested in their Discussion. The original study pointed out the importance of membrane properties and short-term dynamics of synaptic inputs onto the subpopulations in shaping the results. Note, that the lower synaptic reliability in vivo weakens synaptic depression (Borst, 2010), and the ongoing background input leads to higher membrane conductance (Destexhe et al., 2001); therefore our results are expected to differ substantially, and provide an independent prediction for a different dynamic state. While in the original in vitro study stimulation of 10 PCs was sufficient to find a decodable signal, we expected a larger number to be required under our in vivo-like conditions, due to the lower [Ca²⁺]_o, and the higher level of noise (see Methods).

Indeed, we found little discernible mutual information in I spiking for the activation of 10 PCs, but that mutual information increased with the number of activated neurons (Figure 6K, right). For 1000 activated neurons each of the three I subpopulations in L2/3 showed clear but qualitatively unique modulation for both encoding schemes (Figure 6J). As in (Prince et al., 2021) we found differences in the encoding capabilities of the I subpopulations (Figure 6K). We found strongest encoding in PV+ neurons for both schemes, and overall low encoding strength in Sst+ neurons. Further, we found overall low amounts of mutual information for the synchronous stimulus, which is in line with Prince et al. (2021), as they linked encoding of that stimulus to depressing synaptic dynamics, which would be weakened in vivo.

2.9 Exploring the effect of inhibitory targeting trends observed in electron microscopy

We next explored how certain trends in inhibitory targeting observed in a recent mouse electron microscopic connectome (Schneider-Mizell et al., 2023; Figure 7A) affect spontaneous and evoked activity. Consequently, an alternative connectome which recreates these trends has been built (Reimann et al., 2022a). We refer to this as the Schneider-Mizell compatible connectome (SM-connectome). It ensures that PV+ neurons target the perisomatic region, VIP+ neurons target other inhibitory neurons, and cells in L1 connect predominantly monosynaptically. Consequently, inhibition onto E neurons is sparser and closer to the soma (Figure 7A2). After validating the presence of these trends, we characterized how this affects the dynamics of the spontaneous and evoked states (Figure 7A3). Note, that individual synapses were recalibrated to preserve the pathway specific reference PSP amplitudes. As expected, the largest change was for PV+ neurons where the average number of vesicles in the release-ready pool N_RRP was increased by 2.4 times.

Figure 7:

Inspite of this, the rewiring affected the activity sufficiently that recalibration of missing synaptic compensation was needed to attain the target firing rates for the original meta-parameter combinations (Figure S11). Mean conductances onto I neurons had to be increased, indicating the efficacy of the newly introduced I-I targeting by VIP+ neurons (Figure 7B1). On the other hand, conductances for E neurons remained roughly the same, except for L6. Curiously, we found that the firing rates of I populations were reduced when the internal connectivity was activated compared to when it was inactive, further highlighting the efficacy of the I-I targeting (Figure 7B2, S11). For E neurons in L2-5, activating internal connectivity increased firing rates, but less than in the original connectome. As the number of inhibitory synapses innervating them decreased (Figure 7A2), this indicates that the placement of PV+ synapses closer to the soma increased the efficacy of inhibitory as expected. Taken together, this indicates that in the SM-connectome, I populations are more externally driven and more internally inhibited (as a function of proximity to L4), and L6 E more internally driven during spontaneous activity. Over the different meta-parameter combinations, correlation in the network decreased relative to the original circuit with some exceptions (Figure 7B3, S11). This indicates, that activity in the SM-connectome is shifted slightly towards a more asynchronous state.

By correlating the structural and functional changes we concluded that the increased conductance injection into the I populations can be explained by the increased inhibitory synapse count (Figure 7C). Conversely, the mostly decreased synapse counts for E populations did not lead to a decrease of conductance injection, further indicating that the more perisomatic positioning of PV+ synapses compensates for the reduced synapse counts. The lack of effect onto most of the E populations disappeared in the evoked state (Figure 7D). Overall, responses increase strongly in L23 and slightly decreased in L5 and L6 (Figure 7D1), while latencies of responses increased in every layer (Figure 7D2), suggesting different roles of the more specific inhibitory targeting in different layers.

2.10 Mid-range connectivity determines independent functional units and selective propagation of stimulus-responses at larger spatial scales

After calibrating the model of extrinsic synaptic input for the seven column subvolume, we tested to what degree the calibration generalizes to the entire nbS1. Notably, this included addition of mid-range connectivity (Reimann et al., 2022a). The total number of local and mid-range synapses in the model was 9138 billion and 4075 billion, i.e., on average full model simulations increased the number of intrinsic synapses onto a neuron by 45%. Particularly, the final ϕ_1.05,0.4 calibrated for the seven column subvolume was used in conjunction with P_FR [0.1, 0.15, …, 0.3]. Each of the full nbS1 simulations (using [Ca²⁺]_o = 1.05 mM and R_OU = 0.4) produced stable non-bursting activity, except for the simulation for P_FR = 0.3, which produced network-wide bursting activity.

Activity levels in the simulations of spontaneous activity were heterogeneous, ranging from a mean firing rate of 0.15 Hz to 0.3 Hz (when activity was averaged in flat space for P_FR = 0.15; Figure 8A1; Video 6). The subregion with consistently higher firing rates was the upper-limb representation (S1ULp), separated from subregions of lower activity by the dysgranular zones (S1DZ, S1DZO). This demonstrates that the larger spatial scale of the model supports several independently acting functional units. Locations with higher firing rates were accompanied by higher spiking correlations of the local population (Figure 8A2), indicating activity that was more driven by intrinsic excitation. Unlike in the simulations of smaller volumes, the correlation emerged mostly from short, transient packets of activity in L6 (Figure 8B). Consequently, deviation from the target P_FR in those cases was largest in L6 (Figure 8C). Taken together, it appears as if the additional synapses from mid-range connectivity moved the model along the spectrum from extrinsically vs. intrinsically driven, but further for some regions than for others. We confirmed a divergence into two populations: one with correlation slightly increased from the simulations of smaller volumes, and one with almost maximal correlations (Figure 8D). This distinction was determined mostly in L6 with increases of correlations in other layers remaining small.

Figure 8:

We repeated the previous single whisker deflection evoked activity experiment in the full model, providing a synchronous thalamic input into the forelimb sub-region (S1FL; Figure 8E; Video 7). Responses in S1FL were remarkably similar to the ones in the seven column subvolume, including the delays and decays of activity (Figure 8F). However, in addition to a localized primary response in S1FL within 350 µm of the stimulus, we found several secondary responses at distal locations (Figure 8E; Video 7), which was suggestive of selective propagation of the stimulus-evoked signal to downstream areas efferently connected by mid-range connectivity. The response of the main activated downstream area (visible in Figure 8E) was confined to L6 (Figure 8G).

As the model is homogeneous across the different subregions in terms of neuronal morphologies and physiology, we expected the differences in spontaneous and evoked activity across the model to be explained by the topology of synaptic connectivity (Reimann et al., 2022a). While the extrinsic input from the OU injections was uncorrelated by design, synaptic input from local connectivity would feature at least some degree of correlation. This led to asking, how correlated was the input from mid-range connections? We estimated for each pair of mid-range connections innervating a 460 µm hexagonal subvolume its spiking correlation and compared the resulting distribution to the correlation within the subvolume itself (see Methods; Figure 8H, upper). We found low correlations for most pairs, but with a significant number of high-correlation outliers for subvolumes with high internal correlations. Moreover, when we considered the mean of the measure, we found that it explains most of the variability in internal correlations, for both spontaneous (Figure 8H, lower) and evoked states (Figure S12).

This led to the following tentative explanation: Local connectivity provides inputs from within 500 µm, which is virtually guaranteed to be correlated to some degree (Figure 8I upper vs. lower, orange). This is a result of connections being locally dense, as they are concentrated within a small neighborhood (Figure 8 lower, green vs. orange). mid-range connectivity originates from further away, sampling from different functional units that may or may not be correlated, and being more diluted (Figure 8I lower, green vs. blue). Consequently, mid-range inputs are likely to act not much different from the uncorrelated OU process, unless this is overcome by a non-random topology of mid-range connectivity.

Specifically, we considered the mid-range connectivity at meso-scale in terms of the graph of synapse counts between 460 µm subvolumes. When the graph was thresholded at eight million synapses, most subvolumes lacked connections at that strength completely, but others formed a single connected component, indicating the formation of a rich club (Colizza et al., 2006; Figure 8J, upper, red and orange). We found that subvolumes with high internal correlation were part of the rich club, or directly adjacent to one of its members (Figure 8J). In the presence of strong thalamic inputs, this could be overcome, as the increased activity level spread along sufficiently strong mid-range connections (Figure 8K), albeit with substantial delays between 60 and 100 ms. This could induce high correlations in strongly connected parts of the model, outside the rich club as well (Figure S12) and shows that the model supports selective propagation of activity to areas efferently connected with mid-range connectivity.

3 Discussion

We presented novel methods to build, simulate and validate models of cortical tissue that correspond directly with 3D digitised brain atlases. As demonstrated, this enables laboratory experiments to have a simulatable in silico counterpart, and vice versa: predictions made by the model automatically have a precise correspondence in biology. By recreating and extending five laboratory experiments under a single model parameterization, we provided strong model validation and demonstrated the model’s natural versatility. To our knowledge the simulations of the full nbS1 represent the first simulations of stable in vivo-like spontaneous and stimulus-evoked activity in a large-scale biophysically detailed model of multiple cortical subregions connected through local and mid-range connectivity. The model generated an initial set of predictions about the relationship between cortical structure and function, and provides a basis for future studies and predictions. Table S1 describes each quantitative prediction made in this manuscript, and how it might be tested in vivo. We made the model and simulation tools available for use, enabling community driven testing, exploration and refinement.

Excitation from non-modeled brain regions is modeled by injecting somatic conductances and is calibrated for specific populations. For the purpose of calibration, we presented a novel technique, which allows rapid and methodical characterisation and validation of emergent model activity. The efficiency of the technique allows rapid recalibration after changes to the model such as adding new brain regions or changing anatomical details. Connected divided by unconnected firing rates measured during the process are also a prediction of the degree to which the activity of a population is determined by incoming connections from local vs extrinsic sources. We found an increase of the importance of local connectivity from superficial to deeper layers, in line with the canonical view of the cortex (Lübke and Feldmeyer, 2007; Feldmeyer, 2012), which places layers 5 and 6 at the end of local information flow.

Although abstract relative to the model, the injection amplitudes could be predicted from the number of missing synapses, suggesting it provides a functional representation of their anatomical counterparts. The modular model of extrinsic connectivity can be refined or replaced in the future, or when larger fractions of the brain are modeled. Our approach to obtaining in vivo-like activity contrasts with a recent approach in a largescale biophysical model of primary visual cortex of similar scale to our smaller subvolume (Billeh et al., 2020). There, extrinsic inputs were delivered through dendritic synapses instead of somatic injection, which is more anatomically realistic, and intrinsic recurrent weights were adjusted to match extracellularly recorded firing rates. They find a single activity state, instead of the in vivo-compatible spectrum here that allowed us to contrast anesthetized and awake states. The model in general differs from hybrid models, which jointly use biologically-detailed models and point neuron models (Billeh et al., 2020; DuraBernal et al., 2023), and also from the work of Egger et al. (2020) who modeled activation of L5 pyramidal cells by constraining patterns of synaptic input based on receptive fields and predicted anatomical innervation.

Given the complexity of the model, how can we be confident that the its activity is in vivo-like? First, we demonstrated correlated spontaneous dynamics in the form of global and local fluctuations and dynamic E-I balance, as found experimentally (Renart et al., 2010). Second, population firing rate distributions were long-tailed with sub 1 Hz peaks, and were similar for spontaneous and stimulus-evoked activity, as reported (Wohrer et al., 2013). Moreover, mean firing rates below 1 Hz are required by metabolic constraints (Attwell and Laughlin, 2001; Lennie, 2003). In response to simple whisker stimuli, response sparsity, spike counts, and the temporal profile and amplitudes of layer-wise populations were similar to in vivo under anesthetized and awake conditions, for different regions of the meta-parameter space. The effect of the spontaneous meta-parameters on evoked responses, also shows that the networks spontaneous state affects stimulus-responses, as observed in vivo (Isbister et al., 2021). Although from the barrel system, these stimuli offer some of the simplest stimulus-response paradigms, enabling principles of neural dynamics and information processing to be studied; in particular, the correspondence between atomic units of sensory information and neural representations. The model also predicts the number of thalamic fibers stimulated for whisker-flick stimuli. Such validations are crucial for complex models, as they provide context for more complex validations. Moreover, in a rich and complex nonlinear system, the source of discrepancies with in vivo activity are simpler to ascertain under simpler protocols.

For more complex validations, we reproduced protocols from three laboratory experiments. These experiments used a wide array of techniques, which could be accurately recreated due to the model’s close correspondence with biological tissue. Additionally, we could go beyond the original experiments. For example, we predicted an increased number of neurons would be required to have a decodable effect on postsynaptic activity under in vivo conditions over the in vitro conditions of the original experiment. This can be explained by lower synaptic reliability under in vivo conditions. This demonstrates that the approach of modeling entire brain regions can be used to further probe the topics of the original articles and cortical processing. This could be continued in the future using paradigms that take advantage of the multiple regions offered by the model to study processing in cortical hierarchies.

Responses to whisker flick stimuli were closer to biology for values of the P_FR metaparameter between 0.1 and 0.5, according with the known inflation of mean firing rates stemming from the bias of extracellular spike sorting techniques towards larger and more active neurons (Olshausen and Field, 2006). Supporting our estimates, patch-clamp experiments show ubiquitous neuronal silence (Crochet and Petersen, 2006), and spontaneous and stimulus-evoked barrel cortex firing rates as much as 10 times lower than for extracellular recordings (0.05-0.15 Hz vs. 0.8-1.5 Hz; Olshausen and Field, 2006). Each recording technique comes with challenges and biases (Shoham et al., 2006; Olshausen and Field, 2006; Barth and Poulet, 2012; Wohrer et al., 2013). Calcium imaging and patchclamp experiments are less biased towards frequently spiking neurons (Wohrer et al., 2013), yet the former can only infer spiking activity and may be biased towards neurons that favor marker expression. Patch-clamp techniques can reduce firing rates, be biased towards larger or more active neurons (Olshausen and Field, 2006); and are limited to recording a few neurons simultaneously, preventing characterizations of single trial population dynamics. An earlier version of the model has been used to simulate cortical extracellular potentials (Reimann et al., 2013). We can now simulate modern extracellular electrodes and compare multi-unit activity and spike sorting results to the in vivo-like ground truth to gain further insights into the potential over-estimation of mean firing rates.

During simulations of the entire nbS1, we found emerging spatial inhomogeneities in network activity. Specifically, we observed sharp transitions in firing rates and pairwise correlations at disgranular zones. This highlights that the model supports several, independently acting functional units. Analysis of spatial correlations suggests a radius of approximately 400 µm for a single unit. This property was not present in the previous intra-regional scale model (Markram et al., 2015), but is rather an emergent property of the interregional scale, demonstrating that we achieved our goal of building a model that can be used to study mid-range interactions between regions. Analyzing the source of the inhomogeneities suggested a prominent role of a non-random rich-club structure of mid-range connections at the meso-scale. Within functional units, we found that inhomogeneous increases of spiking correlations were confined to spontaneous packets of L6 activity.

It is important to highlight where the model deviates from biological data. First, the model omits several known anatomical elements, such as glia and the presence of gap junctions between certain neuron-types, but the model’s spatial context provides a natural scaffold for them. Second, we model only a single electrical type for pyramidal cells, while some literature sources propose several (Gouwens et al., 2019). However, we are not aware of any method which can recreate the electrical features of different types of pyramidal cells using only generic ion channel models. To achieve the firing pattern behavior of more specific electrical types, usually ion channel kinetics are tweaked, and this would validate the compartmentalization of parameters. Our modeling techniques assume generality of rules and parameters, unless indicated by the data. Therefore, violations indicate that biological rules and data may be more specific. For example, model thalamic inputs innervate any dendrite placed in the targeted layers equally, but this did not reproduce the timing of inhibitory subpopulation responses. Stronger innervation of PV+ neurons than other inhibitory neurons may rectify this. However, it is unclear whether innervation should be anatomically or physiologically stronger, as the best available data mixes structural and physiological strength (Sermet et al., 2019).

Finally, when required, we generalized data from different animals and brain regions to build our rat nbS1 model. This is the accepted state-of-the-art in computational neuroscience, e.g., 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 advanced mouse V1 model of the Allen Institute (Billeh et al., 2020). While, a first truly single species model would be a great advance, it is not required to be the basis of in silico research. Our model (and previous ones) can be used to study the many cortical mechanisms that are common to the closely related organisms rat and mouse. This is illustrated by the follow-up research using the model in the fields of spiking codes (Ecker et al., 2024; Egas Santander et al., 2024), plasticity (Ecker et al., 2023) and electrical signals (Tharayil et al., 2024).

However, we cannot always predict the effect caused by mixing data sources. An interesting example of this is described in the companion paper (Reimann et al., 2022a), stating that certain rules of inhibitory connectivity found in a mouse electron-microscopic (EM) connectome (Schneider-Mizell et al., 2023) cannot be explained by an non-selective pruning axo-dendritic appositions. We demonstrated the utility of the model to explore the effect of any connectivity rule, using the example of a mouse EM-compatible “SMconnectome”. After rewiring, and refitting synaptic, and conductance injection parameters, we found that the purely perisomatic targeting by PV+ neurons and the increased inhibitory targeting by VIP+ neurons shifted how much individual layers were driven by intrinsic and extrinsic excitation. This could not be explored with current techniques in vivo and highlights a strength of in silico experimentation. On the other hand, this data source contradicts some of our rat nbS1 parameters, e.g., the number of synapses per connection from light microscopy (Markram et al., 2015) do not match with the ones obtained from mouse V1 EM. Thus we used the rat data sources for building our baseline model, but make both models available in a simulatable format (see Data Availability).

Acknowledgements

The authors would like to thank Fabien Delalondre, Adrien Devresse, Hugo Dictus, Juan B. Hernando, Nicolas Cornu, Daniel Fernandez, Jeremy Fouriaux and Ioannis Magnakaris for helpful discussions and contributions to the model; Eva Kenny for help with project management; Cyrille Favreau, Marwan Abdellah and Nadir Roman for support with graphics and figure design and Karin Holm and Akiko Sato for support of manuscript development and helpful discussions. The authors would also like to thank Heiko Luhmann, Vicente Reyes Puerta and Jhy-Jung Sun for supporting access to data used in the primary validation of the model.

Funding

This study was supported by funding to the Blue Brain Project, a research center of the École polytechnique fédérale de Lausanne (EPFL), from the Swiss government’s ETH Board of the Swiss Federal Institutes of Technology.

Author contributions

In alphabetical order under each point.

• Conceptualization: H.M., E.B.M., S.R., M.W.R.

• Data curation: S.B.P., A.E., W.V.G, J.B.I, D.M., S.R., R.R., M.R., A.T.

• Formal analysis: A.A., S.B.P., J.B.I., A.E., D.E.S., P.L., C.P., M.W.R., M.R, C.R., V.S., A.T., W.V.G.

• Investigation: S.B.P., A.E., D.E.S., J.B.I., P.L., C.P., M.W.R.

• Methodology: S.B.P., N.B.Z., D.E.S., A.E., J.B.I., P.L., D.M., C.P., M.W.R., M.R., A.R., A.T., W.V.G.

• Project administration: J.D.C., J.K., S.L., H.M., J.P., A.R., S.R., M.W.R., F.S.

• Resources: J.D.C., M.G., W.J., P.K., J.K., P.K., F.P., F.S.

• Software: A.A., O.A., J.B.A., J.D.C., T.Da., T.De., G.F., M.G., J.H., G.I., W.J., J.K., P.K., F.P., M.W.R., V.S., A.T., M.W.

• Supervision: H.M., S.R., M.W.R.

• Validation: A.A., N.B.Z., S.B.P., G.C., A.D., A.E., M.G., J.B.I., D.K., C.P., M.W.R., M.R., V.S., W.V.G.

• Visualization: E.B., S.B.P, A.E., D.E.S., J.B.I., J.P., C.P., M.W.R.

• Writing - original draft: S.B.P, D.E.S., A.E., J.B.I., C.P., M.W.R., M.R.

• Writing - review & editing: S.B.P., D.E.S., A.E., J.B.I., D.K., C.P., S.R., M.W.R., V.S.

Declaration of interests

The authors declare no competing interests.

4 Methods

Resource availability

Lead contact

Further information and requests for data and code should be directed to and will be fulfilled by the lead contact: Michael W. Reimann (michael.reimann@epfl.ch)

Materials availability

No materials were used in this computational work.

Data and code availability

• Neuron morphology reconstructions, ion channel models, neuron electrical models, volumetric atlases, and the description of the model in SONATA format (Dai et al., 2020) have been deposited at Zenodo and are publicly available as of the data of the publication. The DOI is listed in the key resources table. Raw single cell recordings and paired synaptic recordings have been deposited in the Blue Brain portal and are publicly available as of the date of publication. The links are listed in the key resources table.

• All original code has been deposited at Zenodo and is publicly available as of the date of publication. DOIs are listed in the key resources table.

• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

4.1 Method details

4.1.1 Optimisation of ion channel conductance densities

Neuron models from Reva et al. (2023) were made as follows. Ion channel conductance densities were optimised for a pool of reconstructed neurons, for which in vitro recordings had previously been made. Ion channel densities were then generalised to other neuron models of the same e-type. Resulting neuron models were only used in the construction of the model if the value of the cost function was within 5F standard deviations of the experimental mean.

4.1.2 AIS Adjustment

For some electrical models, further adjustments were made to AIS sections to handle the large input currents present under conditions of high network activity. In short, for certain pyramidal cell morphologies, the AIS was replaced by a cylinder with a length of 60 µm. Its diameter was manually calibrated such that the ratio of the AIS input resistance over the input resistance of somato-dendritic compartments was equal to 2. Afterwards, it was confirmed that the cost function of their electrical models remained below the threshold.

4.1.3 Parametrization of synaptic physiology

The Tsodyks-Markram synaptic model (Tsodyks and Markram, 1997; Markram et al., 1998; Fuhrmann et al., 2002; Loebel et al., 2009), upgraded to feature multi-vesicular release (Barros-Zulaica et al., 2019; Ecker et al., 2020) comprises the following parameters: peak conductance g^, depression and facilitation time constants D, F, the decay time constant of the PSC τ_decay, the release probability U_SE, the average number of vesicles in the release-ready pool N_RRP, the Hill coefficient of the nonlinear [Ca²⁺]_o dependent scaling of release probability U_Hill, and finally g^_ratio: the NMDA/AMPA ratio for excitatory synapses, and the GABA_B/GABA_A ratio for inhibitory synapses. Where possible, pathway and synapse-type specific parameter values were taken from the literature, and after eventual corrections (for differences in the temperature and solutions used in the experiments), were put directly into the model using the technique described for the use case of hippocampal CA1 in Ecker et al. (2020). For example, parameters U_SE, D, F, and U_Hill have previously been fitted for a number of pathways using data acquired from paired-recording experiments. Furthermore, g^_ratio parameter has been previously determined from outside-out patch recording experiments.

The remaining parameters (g^ and N_RRP) were iteratively calibrated using an in silico setup reproducing the paired-recordings experiments detailed below. For this parameterization, all paired-recording data sources from Markram et al. (2015) were re-used. To further constrain the variance of EPSP amplitudes in L5 TTPCs, fitted parameters from Barros-Zulaica et al. (2019) were used. The dataset was also enriched with recent recordings from L6 (Qi and Feldmeyer, 2016; Yang et al., 2020, 2022). Instead of generalizing L4

E peak synaptic conductances for all thalamic fibers (as in Markram et al., 2015), synaptic conductances were constrained using EPSP amplitude measurements from thalamocortical slices (Beierlein and Connors, 2002; Beierlein et al., 2003). Fibers originating from POm were not present in Markram et al. (2015). Due to the lack of in vitro measurements, physiological parameters from VPM synapses were generalized to describe POm synapses.

Not all synaptic peak conductances can be measured with paired recordings due to the space-clamp artifact (Markram et al., 2015; Ecker et al., 2020), and estimating the number of release-ready vesicles is also challenging experimentally (Loebel et al., 2009; Barros-Zulaica et al., 2019). Thus, for synaptic pathways for which there was in vitro reference data, these two parameters were determined in a simulation-driven iterative process. Firstly, 50 pairs of in silico neurons were randomly sampled. Second, the postsynaptic cell was current-clamped, and the PSP at the soma of the postsynaptic cell in response to a presynaptic action potential was measured. This was repeated 35 times for each pair as in Ecker et al. (2020). Thridly, the mean and the CV of these amplitudes were extracted and compared with reference in vitro data. Finally, peak conductance (g^), and the average number of vesicles in the release-ready pool (N_RRP) were adjusted to match the amplitude and CV of the in vitro reference (Ecker et al., 2020). In some cases, this resulted in reaching the lower bound of univesicular release (Barros-Zulaica et al., 2019).

Where experimental data was not available, mean synaptic physiology parameters from similar pathways were used (Markram et al., 2015; Ecker et al., 2020). For the VPM to L5 E pathway, the maximum of VPM to L4 E and VPM to L6 E conductances (instead of the mean) were used, based on resulting evoked firing rates found in a later modeling step (Figure S4C). To calibrate the mPSC frequencies of different pathways, single cell simulations with different values of the spontaneous release frequency for all synapses of a set of 1000 cells in a given pathway were run. In these simulations, in silico voltage-clamp recordings were performed to measure the resulting mPSC frequency at the soma. This data was then fitted with a logarithmic function and the value of the spontaneous release frequency matching the in vitro reference value for the mPSC frequency was interpolated. As in vitro paired recording data is sparse, all available sources to determine synaptic parameters were re-used for validation.

4.1.4 Datasets used for validation of in vivo-like dynamics

Activity in the model is compared with spike-sorted extracellular recordings from the barrel cortex of anaesthetised rats (Reyes-Puerta et al., 2015). The dataset contains single unit spike times for both spontaneous activity and responses to single whisker deflections, across 6 layer-wise neuron populations (L2/3 E, L2/3 I, L4 E, L4 I, L5 E & L5 I). We extend the dataset using the average of L2/3 I, L4 I and L5 I as a reference for L1 I and L6 I, and an alternative reference for L6 E (De Kock et al., 2007). The spontaneous firing rates of the 9 neuron populations are collected in the vector of reference firing rates, V_FR.

As a secondary comparison, we use a dataset of in vivo juxtacellular and whole-cell recordings of barrel cortex responses in awake mice during an active whisker touch paradigm (Yu et al., 2019). While the animal model is different, using this data enables comparison to awake responses of layer-wise E populations and I subpopulations (VIP+ (vasoactive intestinal peptide, a part of our 5HT3aR+ subpopulation), Sst+, FS (fast spiking, corresponding to our PV+ subpopulation)). Moreover, this allows validating the model against properties of cortical responses across similar but different paradigms, reducing the risk of over fitting.

4.1.5 Somatic conductance injection

Random background synaptic inputs, representing uncorrelated inputs from other brain regions not present in the model, are modelled as a conductance-based somatic injection (Figure 4B1). Specifically, we inject conductance signals described by an Ornstein-Uhlenbeck (OU) process (similar to Destexhe et al. (2001)), given by the following equation:

where g(t) is the conductance signal, g₀ its mean value, τ is the decay time constant, D is the diffusion coefficient, and χ(t) is a Gaussian white noise process. The diffusion co-efficient can also be expressed in terms of the standard deviation σ of the signal: . Note that we only consider an excitatory conductance signal representing mid-range excitatory drive, whereas local inhibitory input is deemed sufficient within the circuit model. As we are modeling a non-specific random background noise, we use a different, uncorrelated OU process for each cell.

Given the diversity of morphologies and electrical behaviors of cells in our circuit, we scale the mean and standard deviation of the injected signal based on individual cell properties. Specifically, we compute the input resistance R_in of each cell at its resting potential, then take its reciprocal value G_in = 1/R_in (units of conductance). The mean g₀ and standard deviation σ of the injected conductance signal are then expressed as percentages OU_µ and OU_σ of this input conductance G_in:

The OU_µ and OU_σ percentage values are used in the calibration technique described in the following sections. To reiterate, OU_µ and OU_σ scale the mean of the injected conductance and the size its fluctuations respectively.

4.1.6 Details of spontaneous activity calibration

OU_µ and OU_σ are parameterised separately for each of 9 populations (layer-wise E and I combinations; Figure 4B2), leading to 18 parameters, grouped into two 9-dimensional vectors OU_µ and OU_σ. Spontaneous firing rates across populations, denoted by the vector C_FR, are then determined by the vectors OU_µ and OU_σ, and [Ca²⁺]_o. The objective of the calibration procedure is to determine a mapping that outputs the OU_µ and OU_σ required to achieve a particular set of firing rates C_FR, for a given [Ca²⁺]_o.

This requires two simplifications: First, we reduce the space of possible OU_µ and OU_σ values by fixing OU_σ = OU_µ ∗ R_OU, where R_OU is a pre-specified constant that controls the level of noise assumed to be present in the external synaptic input. Second, we learn the mapping only for dynamical states where the firing rates of the 9 populations are a proportion of reference in vivo firing rates, i.e. C_FR = P_FR ∗ V_FR, with P_FR ∈ [0, 1] (Figure 4B2). This is based on the assumption that extracellular spike detection sampling biases (Olshausen and Field, 2006; Wohrer et al., 2013) affect all populations equally.

The mapping can be written as Ψ_Ca2+_,ROU : C_FR 1→ OU_µ and is simplified into two separate mappings which are determined separately and then combined (Figure S1). The first mapping finds values of the 9 output parameters OU_µ, for a given R_OU, which produce particular unconnected firing rates U_FR (i.e., where the neurons in the network are not connected to each other). The second then maps between U_FR and C_FR.

The first mapping χ : (R_OU, U_FR) 1→ OU_µ, is attained by running 1 second simulations over various combinations of OU_µ and OU_σ (Figure S1). R_OU constrains OU_µ and OU_σ to a line of possible combinations in the 2D space. We generate an approximation of χ based on linear interpolation between the results of simulations along the line. Note in particular, that χ is independent of [Ca²⁺]_o.

Given the mapping χ, the remaining problem is to determine the values of U_FR that achieve certain target firing rates given by P_FR ∗ V_FR, for P_FR ∈ [0, P_max] and each combination of [Ca²⁺]_o and R_OU, i.e., to determine the mapping ϕ_Ca2+_,R: U_FR 1→ C_FR.

This is done by measuring values of C_FR produced for given values of U_FR in a number of simulations. However, instead of exploring the whole parameter space of possible U_FR values, we only explore a small number of values of U_FR which we expect to produce the target C_FR values. This is done by iteratively learning the mapping ϕ_Ca2+_,ROU : U_FR 1→ C_FR and in each iteration using ϕ to predict U_FR values that bring the dynamics closer to the target firing rate, i.e., by sampling values given by U_FR = ϕ⁻¹(C_FR) = ϕ⁻¹(P_FR ∗ V_FR) for 10 equidistant values of P_FR ∈ [0, P_max]. The initial guess of ϕ is simply the identity, i.e., assuming no effect of connecting the network, with P_max = 0.5. All other iterations use P_max = 1.0.

We update ϕ by fitting an exponential: C_FR = α ∗ e^{βUF R} + κ to the values of U_FR, C_FR measured on the latest iteration (Figure S1). In a recurrently connected network such as our model, the connected firing rate of a population depends on the firing rates of all other populations in complex and nonlinear ways. This is addressed by our simplification to restrict the compatible dynamic states, where by definition the target firing rate of one population determines the target firing rates of all others.

The firing rates of the L2/3 I, and L4 I and L6 I populations could not be sufficiently lowered for P_FR < 0.2 and P_FR < 0.1 respectively, as they spontaneously fired even at their resting potentials (Figure S1D).

4.1.7 Generalization to other network states

We first approximated ϕ for [Ca²⁺]_o = 1.1 mM and R_OU = 0.4. To reduce the number of iterations steps needed for other combinations of [Ca²⁺]_o and R_OU, we initialized the optimization process with the mapping ϕ obtained for this parameter pair instead of basing it on the identity (Figure S1D, lower; Figure S1). Only for P_FR above 0.8 and [Ca²⁺]_o = 1.1 mM did we have to slightly relax our acceptance criteria (Figure 4D, lower). Additionally [Ca²⁺]_o = 1.1 mM, P_FR = 1.0, and R_OU = 0.2 led to bursting activity uncharacteristic of in vivo activity and was excluded from further study.

4.1.8 Evoked activity comparison and validation

The spike times of single whisker deflection stimuli were drawn from a PSTH of VPM neurons in response to 3 ms mechanical single whisker deflections made at 1 Hz in urethane anesthetized rats (Diamond et al., 1992; Figure 5A2). Since the in vivo data for cortical responses to single whisker deflections was from responses to 2 ms mechanical whisker deflections, we compressed the tail of the VPM PSTH to allow a direct comparison with in silico responses (Figure 5A2). A VPM PSTH for the active whisker touch stimulus was used from the same dataset as the corresponding cortical responses (Yu et al., 2019).

Firstly, in silico responses were compared to data from the C2 barrel of the anaesthetized rat in response to single whisker deflections made at < 1Hz. For the active whisker touch stimulus, in vivo activity was used from aggregated over principal barrels corresponding to the untrimmed whisker.

As an initial test of similarity with the corresponding in vivo dataset, trial-averaged PSTHs were first calculated for each parameter combination, both for all neurons, and for each of the layer-wise E and I populations. We then tested whether the latencies, 50% and 25%-decay points of the trial-averaged PSTHs were respectively no more than 10 ms, 10 ms and 40 ms later than those of the corresponding in vivo populations, and that there was no secondary rise in the PSTHs after the initial decay (by testing that none of the PSTHs for individual populations were more than 35% above baseline after 75 ms). L1 I was excluded from these tests, as it showed little or no response to the stimulus and we had no in vivo reference. These tests allowed us to split the 4D parameter space into two regions, one where in silico responses are consistent with in vivo responses, and another where they are not (see Results).

4.1.9 In-silico experimental methods

Simulation of a morphologically detailed model first requires specifying the conditions (parameters and inputs) under which to perform the simulation. Given the correspondence between variables in the model and properties of biological systems, we are able to mimic to a certain extent the conditions and protocols of experimental studies. These are combinations of existing techniques, or even techniques that cannot or have not yet been performed experimentally. On top of the compensation for missing external input and different [Ca²⁺]_o, we implemented further mechanisms to simulate specific experimental conditions. On a technical level, these comprise somatic current or conductance injections, adjustment of synaptic connectivity parameters, and selective activation of thalamic fibers. These are static conditions, however, in that their time course must be determined before the simulation is run. Somatic injections and connectivity adjustments can be targeted with single neuron or single pair resolution, respectively. To each thalamic fiber or neuron in the model, we can assign an arbitrary spike train, triggering synaptic release from all of its synapses with anatomically determined delays (from axonal conduction speed and path length). These mechanisms can target specific groups of neurons, e.g. based on neuron properties, such as location, layer, m-type or e-type, or inhibitory subpopulation (determined based on me-type as in Markram et al. (2015); Figure 2B).

We combined mechanisms to simulate various in vitro or in vivo experimental paradigms. Pathway lesions, for example, are simulated by selecting sets of pre- and postsynaptic neurons and removing all synaptic connections between them (Figure S2A). Simulated optogenetic stimulation, in addition to being targeted at specific groups of neurons, took into account the attenuation of light with depth for the wavelength being used. Neurons were then inhibited or excited according to the intensity of light reaching them (Figure S2B). Finally, sensory stimuli were simulated by generating instances of stimulus-specific stochastic spike trains activating the thalamic input fibers (Figure S2C).

4.1.10 Optogenetic inhibition or activation

We model optogenetic inhibition or activation of a neuron population through a current injection at the soma of each cell in the population, with an intensity proportional to the cell’s threshold current (see Reva et al. (2023); but with the technical caveat that we did not apply a hyperpolarizing holding current at the same time). To mimic the conditions of surface illumination, we considered the dependence of effective depolarization strength on cortical depth using a modified Beer-Lambert law approximation for the exponential attenuation of light intensity through scattering tissue (Al-Juboori et al., 2013; Azimipour et al., 2014, Figure S2B):

where I(d) describes the light intensity at depth d (in mm), with a maximum light intensity I₀ (on the surface of the cortex) and an effective attenuation coefficient given by:

Values for the absorption coefficient µ_a and the reduced scattering coefficient µ^′ were interpolated for the chosen wavelength from Mesradi et al. (2013). For implementation reasons, the targeted cells were grouped into depth bins and for each group the intensity at the center of the bin was used. Bins were equally distributed in terms of the injected current, and not in terms of depth. The value of I₀ was calibrated independently for each experimental paradigm.

4.1.11 Modeling sensory inputs

Sensory inputs were simulated using a three step procedure (Figure S2C2): First, a time series representing sensory stimulation is assigned to each selected thalamic fiber; second, the time series is optionally transformed with a fiber-dependent transfer function; third, a spike train is generated from the time series through a stochastic spiking process. Thalamic input fibers were associated with roughly columnar, overlapping volumes of the model that they formed synaptic connections in (Reimann et al., 2022a). Their centers were projected into the plane using a flat map of the model (Bolaños-Puchet and Reimann, prep), yielding [x_ζ, y_ζ], the flat locations of each fiber ζ. A sensory stimulus was defined based on fiber location as ρ_ζ(t) = ρ(x_ζ, y_ζ, t). The stimulus function could be partially stochastic and its value was evaluated for each time bin of the simulation (Figure S2C1). Next, a fiber could be associated with a transfer function υ_ζ that transforms the results of ρ. If none is mentioned, identity was used as transfer function. Finally, a spiking process ψ was used to instantiate a spike train for each fiber, based on its transformed time series. That is, the spike train associated with ζ was Γ_ζ = ψ(υ_ζ(ρ_ζ)). The processes used were stochastic, leading to different spike trains for fibers associated with the same time series. Different, specific ρ, υ and ψ were used for different experimental paradigms that will be further described below.

4.1.12 Recreating Varani et al., 2022

To study how input from L4 contributes to L2/3 subthreshold responses Varani et al. (2022) used a 500 ms long whisker hold paradigm, while patch-clamping PCs in L2/3 in anesthetized and awake mice. We encoded the whisker hold stimulus as a step function (ρ_ζ(t) = 1, if 2000 ms ≤ t < 2500 ms, 0 otherwise) in 10% of the VPM fibers within the seven column subvolume. One of four transfer functions were assigned to each fiber, based on the types of kinetic response properties of thalamic neurons identified in (Petersen et al., 2008). The types were selective for whisker position (υ_pos), velocity (υ_vel), acceleration (υ_acc), or direction (υ_dir), and were implemented as:

where r_max = 150 Hz denotes the firing rate of a thalamic fiber when its associated feature property is at the fiber’s preferred value. Transfer functions were randomly assigned to fibers with the fractions identified in Petersen et al. (2008) (Figure ??A, 11% coding for position and acceleration, 58% for velocity and 20% for direction). The spiking process ψ was an adapting Markov process (Muller et al., 2007) with an adaptation time constant of 100 ms. As the whisker movements were short lasting, rates were evaluated at submillisecond (0.1) resolution.

Similarly to before, only the seven column subvolume (210k neurons) was used for the simulations, and the in vivo-like state was realized as Ca²⁺ = 1.05 mM, R_OU = 0.4, and P_FR = 0.3.

The optogenetic inhibition in this experiment targeted 95% (in line with Varani et al. (2022)) of excitatory cells in layer 4. The authors found a few cells which also tested positive for Halo at the bottom of L3 as well, but as they did not quantify it, we decided not to target any lower L3 PCs in the in silico version of the experiment. Based on the 595 nm wavelength (yellow light) parameters were set to µ_a ≈ 0.49 mm⁻¹ and µ^′ ≈ 4.12 mm. For the depth-based spatial binning of cells, 5 bins were used in L4. After scanning several values, I₀ was set to -200% as that reproduced the ≈ 10 mV hyperpolarization of L4 PCs observed in vivo. In line with the in vivo experiment, the optogenetic stimulus ended in a (100 ms long) ramp to avoid rebound spikes (Figure ??B).

When going beyond reproducing the same experimental conditions and instead lever-aging the in silico nature of our setup, synaptic pathways were lesioned by selecting the excitatory population in a given layer as the pre-synaptic population and the excitatory population in L2/3 as the post-synaptic population and not instantiating the connecting synapses during the simulation.

L2/3 PCs had to meet three criteria to be included in the subsequent analysis. Firstly, their activity was required to remain subthreshold during the 500 ms long whisker hold stimulus and in 200 ms long time windows before and after the stimulus, both in control and in silico optogenetic runs. Second, they had to be innervated by at least one (active) VPM fiber. Third, the derivative of their voltage trace had to cross the 1 mV/ms threshold in a 20 ms time window after stimulus onset in the control simulation. The last two were motivated by comparing subthreshold voltages to voltage traces from Varani et al. (2022) that showed large, stimulus evoked EPSPs. Around 8% of L2/3 PCs in the central column met all the above criteria and their voltages were averaged to arrive to the traces shown in Figure ??E-F. Thus, unlike in the original analysis, cells rather than trials were averaged. The motivation for this approach is that while in vivo it is easier to repeat the same paradigm after establishing stable recording conditions in a given cell, in silico it is quicker to record from all cells in a single simulation, instead of repeating the stimulus several times.

4.1.13 Recreating Shapiro et al., 2022

Since our model comprises the nbS1 and not the V1 brain area, we modelled visual drifting gratings in a more abstract way, without taking specifics of the visual system into account. In particular, we defined ρ(x_ζ, y_ζ, t) as a spatio-temporal rate pattern corresponding to linear sinusoidal drifting gratings for a fraction F_P of 937 VPM fibers projecting to a single simulated column (30k neurons, 520 µm diameter). For all other VPM fibers ρ was zero. The gratings had a temporal frequency of f_temp = 2 Hz and a spatial frequency of 0.03 cycles/degree which we translated to f_spat = 0.001 cycles/µm by assuming a cortical magnification factor in rat of 30 µm/degree (Gias et al., 2005). In Figure ??A, the spatial grating patterns are illustrated at different points in time.

As in Shapiro et al. (2022), we used five contrast levels C ∈ [0.06, 0.12, 0.24, 0.5, 1.0] which were defined as the Michelson contrast given the minimum and maximum luminance values L_min and L_max:

Since the physical quantity of luminance does not have a clear correspondence in our model, we used normalized luminance values L_min and L_max between 0.0 and 1.0 by computing the inverse Michelson contrast centered around a mean normalized luminance of 0.5. The resulting values of L_min and L_max where then scaled and shifted to the minimum and maximum rates of the sinusoidal modulation R_min and R_max respectively such that the peak firing rate at contrast 1.0 was given by R_peak and the mean of the modulation corresponding to the background firing rate at contrast 0.0 was given by R_bk. The resulting sinusoidal input rate signal can be written as:

where t is the time in seconds, l_ζ the linear position of ζ within the grating, and |₊ denotes that firing rates are truncated to positive values only. Figure ??B illustrates a random series of rate signals corresponding to different contrasts. No transfer function was used, and spiking process ψ was an adapting Markov process (Muller et al., 2007) with an adaptation time constant of 100 ms. Values for F_P, R_peak, and R_bk were determined beforehand by a parameter optimization so that the resulting grating responses were qualitatively comparable with the ones reported by Shapiro et al. (2022), see next subsection for details.

Calibration of drifting grating stimulus We ran a parameter scan to determine optimal values for F_P (fraction of 937 VPM fibers projecting to the central column of our model to apply grating stimulus to), R_peak (peak firing rate at contrast 1.0), and R_bk (background firing rate at contrast 0.0). This calibration was done by running 45 simulations using all combinations of parameter values in the ranges F_P ∈ [0.5, 0.75, 1.0], R_peak ∈ [5.0, 10.0, 15.0, 20.0, 25.0] Hz and R_bk ∈ [0.05, 0.1, 0.2] Hz, and selecting the optimal combination amongst them. Each simulation lasted 40 s during which 20 contrast stimuli were presented for 1 s, followed by a 1 s (blank) inter-stimulus interval. Four contrast levels C ∈ [0.06, 0.12, 0.5, 1.0] were presented five times each in random order, which was sufficient to fit sigmoidal tuning functions with four parameters (see Eq. 9).

We observed that especially under strong stimulus conditions the response rates of PCs to the first and second cycle of the sinusoidal grating pattern (1 s stimulus at f_temp = 2 Hz) were quite different, with the second response largely attenuated due to synaptic depletion. Therefore, we extracted first and second peak firing rates r₁ and r₂ respectively from the peristimulus time histograms (PSTHs) computed with 1 ms resolution and 20 ms Gaussian smoothing (Figure S8). We defined a measure of the normalized peak difference as the Michelson contrast of the peak rates r₁ and r₂, given by

Additionally, we extracted average firing rates for each contrast level of the whole population of PCs within the full (1 s) as well as the first and second halves (0.5 s each) of the stimulation intervals. We then fitted sigmoidal tuning functions with parameters c₅₀, m, n, and R_max (see Eq. 9) to these average tuning responses (Figure S9). Finally, as summarized in Figure S8 and S9, we selected the best combination of parameters based on the following selection criteria:

• The peak firing rates in response to the grating stimulus should be sufficiently strong, covering a range of values including the one reported in (Shapiro et al., 2022, cf. Figure 1D). So we imposed the constraint that the maximum peak rates over the whole population of PCs of the first and second peaks r₁ and r₂ respectively should be at least r_th = 30 Hz at maximum contrast.

• The peak responses to the first and second cycle of the grating stimulus should not be too different. So we aimed for a low peak difference r^_diff at maximum contrast.

• The overall tuning response of the PCs should be in a regime where the tuning curve would have a sigmoidal shape, meaning that it should be increasing but saturating with increasing contrasts. To fulfill this constraint, we aimed for

  • – Low c₅₀, i.e., the inflection point of the tuning curve would be at low contrasts

  • – High n, i.e., a steep non-linear increase of the tuning curve until saturation

In order to combine these criteria independent of the actual scales of r^_diff, c₅₀ and n, we computed their individual rank scores in decreasing (r^_diff, c₅₀: lower is better) or increasing (n: higher is better) order. We then selected the parameter combination with the highest product of rank scores, excluding the ones with peak firing rates r₁ and r₂ below r_th and the ones with c₅₀ values close to 1.0 (border cases). We found an optimal parameter combination of F_P = 1.0, R_peak = 10.0 Hz, and R_bk = 0.20 Hz which we used throughout this in-silico experiment.

Optogenetic stimulation Optogenetic stimulation was targeted at either 1654 PV+ or 822 Sst+ interneurons in a single column. We used parameters µ_a ≈ 0.46 mm⁻¹ and _s ≈ 5.38 mm, based on the wavelength of 470 nm (blue light; see Figure S2). I₀ was increased from 0% to 300% in steps of 50%.

Quantification of contrast tuning responses by sigmoidal functions We quantified contrast tuning responses in the same way as described by Shapiro et al. (2022), by least-squares fitting sigmoidal functions to the normalized tuning curves. Normalized tuning curves were obtained by computing the time-averaged firing rates of all 1 s stimulus intervals and dividing them by the mean baseline firing rate (i.e., w/o optogenetic stimulation) at the highest contrast level. The sigmoidal function was given by

where R(c) describes the response amplitude at contrast c, m is the baseline response at minimum contrast, R_max is the maximum increase above baseline, n defines the steepness of the curve, and c₅₀ is the contrast at half R_max. We used the coefficient of determination (r² score) as a measure of the goodness of fit.

Detection of neurons with robust contrast tuning We identified PV+, Sst+, and pyramidal neurons with robust contrast tuning behavior under all conditions. In the experimental study of Shapiro et al. (2022), tetrode recordings were used together with spike sorting. Correspondingly, we only considered neurons firing at rates above 0.5 Hz under all stimulus conditions, meaning they could potentially be detected by spike sorting (Pedreira et al., 2012). In addition, we considered neurons as being robustly tuned if they had strictly monotonically increasing tuning curves.

Modelling the effects of optogenetic stimulation of interneurons Direct photostimulation effects on interneurons were modelled by a divisive scaling model R_div(c), a subtractive shifting model R_sub(c), or a saturation additive model R_sat(c) (Shapiro et al., 2022). For fitting these models, the parameters c₅₀, m, n, and R_max of the underlying contrast tuning function R(c) were kept constant at values obtained from baseline fits (i.e., w/o optogenetic stimulation). The divisive scaling model was defined as

with a scaling term g. The subtractive shifting model was given by

with a shifting term h and rectification to rates equal or above zero. The saturation additive model was defined as

with a saturation term S and an additive term A. We used the r² score to measure the goodness of model fits.

Indirect effects on PCs receiving inhibitory input from optogenetically activated interneurons were modelled by a conductance-based model R_cond(c), assuming a saturating additive model description of interneurons (Shapiro et al., 2022). Response rates under this conductance-based model were given by

with a spike threshold V_th = 3.4 mV and rectification to rates equal or above zero. The membrane potential ΔV (c) as a function of contrast was given by

with values for leak conductance g_L = 6 nS, leak reversal potential R_L = −50 mV, excitatory reversal potential R_E = 0 mV, inhibitory reversal potential R_I = −65 mV, and resting potential V_r = −50 mV as reported by Shapiro et al. (2022). The excitatory synaptic conductance was given by

with the excitatory conductances at low/high contrast given by g_{E min} and g_{E max} respectively. The inhibitory synaptic conductance was given by

with an inhibitory conductance offset Δg_{IE min} = 2 nS at low contrast relative to g_{E min}.

In a first step, parameters S and A were set to zero and the model parameters c50, g_{E min}, g_{E max}, and n were fit to baseline responses of PCs. In a second step, those parameter values were kept constant and the model parameters S and A were fit to PC responses under photostimulation conditions. Again, we used the r² score to measure the goodness of model fits.

4.1.14 Recreating Prince et al., 2021

We also compared the in vivo state of our model with the results of a recent in vitro study (Prince et al., 2021), which explored how neurons with different biophysical properties encode different types of signals. The in vitro study optically activated groups of 10 PCs in slices of L2/3 barrel cortex from mouse lines labelling fast-spiking (FS) and regular-spiking (RS) interneurons. Each PC was targeted individually and a binary signal was encoded either through changes in the rate of optical pulses applied to the 10 stimulus neurons, or through changes in the synchronicity of the pulses. The mutual information shared between the binary signal (when encoded either as changes in firing rate or syn-chronicity) and the firing rate of different inhibitory sub-types (recorded using whole-cell patch clamping in vitro) was analysed.

While the in vitro stimuli used optical pulses targeted at single neurons with timings drawn from Poisson processes, it was not necessary to explicitly model such an optical stimulus in silico, as we can instead elicit spikes directly in the model. In the original study, the 10 stimulus neurons were uniformly separated in a grid-like pattern at horizontal and vertical distances of ∼50 µm. The binary signal alternated between two states at random intervals of 2-7 s. For the rate coding paradigm, the timings of the optical pulses for the 10 stimulus neurons were drawn from 10 independent inhomogeneous Poisson processes, with the rates of optical activation varying between 5 Hz and 0.5 Hz for the up and down states respectively (average firing rate: 2.7 Hz). For the synchronous case, the timings of the optical pulses for the 10 stimulus neurons were drawn from a single inhomogeneous Poisson process during the up state, or 10 independent inhomogeneous Poisson processes during the down state. For the synchronous case, the firing rates of the up and down state were both 2.7 Hz. We therefore tested stimulus encoding using a range of stimulus neuron counts: from 10 to 1,000.

In our interpretation of the study, mutual information was measured between the binary signal and neural activity during a 0-5 ms and 5-50 ms window following a change in the signal. These represented ‘early’ and ‘late’ stimulus encoding windows following the change in the signal, respectively. As the synchronous up state activated synchronous patterns at a rate of 2.7 Hz, we calculated that a synchronous pattern would only be activated during the early and late windows with probability 0.0135 (= 2.7 ∗ 5/1000) and 0.1125 (= 2.7 ∗ 45/1000) respectively. We therefore chose to use higher FR of 20 Hz for the synchronous up and down states, and to analyse mutual information between the binary signal and all 50 ms bins. To afford comparison with the rate code experiment type, we used 30 Hz and 10 Hz as the FRs of the up and down states respectively. We also compare stimulus coding during the first 50 ms following the stimulus change for both the synchronous and rate stimulus types.

4.1.15 Analyzing mid-range connectivity & correlations

Neurons were split into hexagonal subvolumes of a specified diameter d considering their locations in a flattened view (Bolaños-Puchet and Reimann, in prep). Then, the number of mid-range connections within and between the hexagonal subvolumes were counted, yielding in S, the adjacency matrix of the resulting directed graph. Distances between hexagonal subvolumes were calculated by considering the centers of the hexagons in the flattened view. Furthermore, correlations of spiking activity of pairs of hexagonal sub-volumes were also calculated. To that end, all spikes within a hexagonal subvolume were pooled and their number in 5 ms time bins; and lastly the Pearson correlations of the resulting time series were calculated. The correlation within a hexagonal subvolume was calculated similarly, but using the separate time series of the E and I populations. Together, this yielded F, the matrix of correlations of subvolumes. Based on S and F, the expected correlations of mid-range inputs was calculated as follows: Let i be a hexagonal subvolume and S_i the column of S associated with input counts into i. Then P = S_i · S^⊤ is in the matrix of counts of pairs of inputs for all pairs of hexagonal subvolumes. Combining P and F allows one to estimate the distribution of correlations of mid-range inputs into i, albeit at the population rather than single-cell level. For the correlation in Figure 8A2, H, J1 a hexagon size of 400µm was used, while in Figure 8I2 a smaller size of only 50µm was used.

4.1.16 Node participation

Given a connectivity graph G, a directed n-simplex in G is a set of n + 1 nodes which are all all-to-all connected in G in a feed-forward fashion, i.e., such that any subset of these has a unique source and a unique sink (see Reimann et al. (2017) for more details). The n-node participation of a node v in G is the number of directed n-simplices this node is part of. In particular, for n = 1, this is the total degree of the node in G. Thus, this can be though of a generalization of the notion of degree that takes into account account higher order interactions and has been shown to strongly correlate with other node centrality metrics Sizemore et al. (2018).

4.2 Quantification and statistical analysis

Details of all statistical analyses can be found in figures and figure legends.

4.3 Additional Resources

No other additional resources.

Videos

The following supplementary videos are available at: https://www.biorxiv.org/content/10.1101/2023.05.17.541168v5.supplementary-material.

• 1. Layer-wise E and I population rasters and max-normalised histograms of spontaneous activity for the central column of the seven column subvolume for each of the 60 meta-parameter combinations. Meta-parameter values are shown in brackets in the following order: Ca²⁺, P_FR, R_OU.

• 2. Visualisation of spontaneous activity for the seven column subvolume for the parameter combination Ca²⁺ = 1.1 mM, P_FR = 0.9, R_OU = 0.2 after collapsing activity to flatspace (activity binned and smoothed).

• 3. Layer-wise E and I population rasters and max-normalised histograms of evoked activity during the 10 single whisker deflection protocol for each of the simulated meta-parameter combinations. Meta-parameter values are shown in brackets in the following order: Ca²⁺, P_FR, R_OU, F_P.

• 4. Trial-averaged max-normalised histograms of evoked activity during the 10 single whisker deflection protocol for each of the simulated meta-parameter combinations. Rasters of spiking activity shown for the first of 10 trials activity. Dashed lines show in vivo single whisker deflection reference data. Meta-parameter values are shown in brackets in the following order: Ca²⁺, P_FR, R_OU, F_P.

• 5. Visualisation of trial-averaged response of the seven column subvolume over the 10 single whisker deflection protocol for the parameter combination Ca²⁺ = 1.1 mM, P_FR = 0.3, R_OU = 0.4, F_P = 20% after collapsing activity to flatspace (activity binned and smoothed).

• 6. Visualisation of spontaneous activity for the full nbS1 model for the parameter combination Ca²⁺ = 1.05 mM, P_FR = 0.15, R_OU = 0.4 after collapsing activity to flatspace (activity binned and smoothed).

• 7. Visualisation of trial-averaged response of the full nbS1 model over the 10 single whisker deflection protocol for the parameter combination Ca²⁺ = 1.1 mM, P_FR = 0.15, R_OU = 0.4, F_P = 20% after collapsing activity to flatspace (activity binned and smoothed).

Note

This reviewed preprint has been updated to include a sentence to the start of the section titled “Videos” to provide a direct URL to the supplementary videos.