Large-scale biophysically detailed model of somatosensory thalamocortical circuits in NetPyNE

The primary somatosensory cortex (S1) of mammals is critically important in the perception of touch and related sensorimotor behaviors. In 2015, the Blue Brain Project (BBP) developed a groundbreaking rat S1 microcircuit simulation with over 31,000 neurons with 207 morpho-electrical neuron types, and 37 million synapses, incorporating anatomical and physiological information from a wide range of experimental studies. We have implemented this highly detailed and complex S1 model in NetPyNE, using the data available in the Neocortical Microcircuit Collaboration Portal. NetPyNE provides a Python high-level interface to NEURON and allows defining complicated multiscale models using an intuitive declarative standardized language. It also facilitates running parallel simulations, automates the optimization and exploration of parameters using supercomputers, and provides a wide range of built-in analysis functions. This will make the S1 model more accessible and simpler to scale, modify and extend in order to explore research questions or interconnect to other existing models. Despite some implementation differences, the NetPyNE model preserved the original cell morphologies, electrophysiological responses and spatial distribution for all 207 cell types; and the connectivity properties of all 1941 pathways, including synaptic dynamics and short-term plasticity (STP). The NetPyNE S1 simulations produced reasonable physiological firing rates and activity patterns across all populations. When STP was included, the network generated a 1 Hz oscillation comparable to the original model in vitro-like state. By then reducing the extracellular calcium concentration, the model reproduced the original S1 in vivo-like states with asynchronous activity. These results validate the original study using a new modeling tool. Simulated local field potentials (LFPs) exhibited realistic oscillatory patterns and features, including distance- and frequency-dependent attenuation. The model was extended by adding thalamic circuits, including 6 distinct thalamic populations with intrathalamic, thalamocortical (TC) and corticothalamic connectivity derived from experimental data. The thalamic model reproduced single known cell and circuit-level dynamics, including burst and tonic firing modes and oscillatory patterns, providing a more realistic input to cortex and enabling study of TC interactions. Overall, our work provides a widely accessible, data-driven and biophysically-detailed model of the somatosensory TC circuits that can be employed as a community tool for researchers to study neural dynamics, function and disease.

Introduction

The primary somatosensory cortex (S1) of mammals is critically important in the perception of touch and works closely with other sensory and motor cortical regions in permitting coordinated activity with tasks involving grasp (Bosman et al., 2011; Petrof et al., 2015; Barthas and Kwan, 2017). Moreover, the communication of these cortical areas with the thalamus is crucial for maintaining functions, such as sleep and wakefulness, considering that the thalamocortical (TC) circuit is essential for cerebral rhythmic activity (O’Reilly et al., 2021). A greater understanding of S1 cortical circuits will help us gain insights into neural coding and be of assistance in determining how disease states such as schizophrenia, epilepsy and Parkinson’s disease lead to sensory deficits or uncoordinated movement (Vázquez et al., 2013; Petrof et al., 2015; Azarfar et al., 2018; Peña-Rangel et al., 2021).

There exists an impressive, highly detailed model of rat S1 developed by the Blue Brain Project (BBP) (Markram et al., 2015), incorporating anatomical and physiological information from a wide range of experimental studies. This groundbreaking model includes over 31,000 neurons of 55 layer-specific morphological and 207 morpho-electrical neuron subtypes, and 37 million synapses capturing layer- and cell type-specific connectivity patterns and synaptic dynamics. Simulation results matched in vitro and in vivo experimental findings, and the model has been used over the years to reproduce additional experimental results and generate predictions of the dynamics and function of cortical microcircuits (Reimann et al., 2015, 2017a,b; Gal et al., 2017; Hagen et al., 2018; Amsalem et al., 2020). Although the BBP S1 model is state-of-the-art, certain constraints limit its reproducibility and use by the community, as well as its extension or modification to connect to other regions or update model features. The size and complexity of any model of this scope is daunting. Due to its scale and complexity, the original model must be run and analyzed on large High Performance Computing platforms (HPCs), which are not available to many users. Although the model is simulated using NEURON (Carnevale and Hines, 2006; Migliore et al., 2006) a widely used platform within the computational neuroscience community, it also requires other custom libraries specifically designed to facilitate this workflow. These libraries are used to build, manage simulations and analyze the model. However, not all of these libraries and workflows are publicly available (Markram et al., 2015), making it somewhat difficult to modify the code, and scale or simplify the model for simulation on smaller computers, overall reducing its accessibility and reproducibility (McDougal et al., 2016).

Here we implemented the original BBP S1 model in NetPyNE (Dura-Bernal et al., 2019) in order to make it more accessible and simpler to scale, modify and extend. NetPyNE is a python package that provides a high-level interface to the NEURON simulator, and allows the definition of complex multiscale models using an intuitive declarative standardized language. NetPyNE translates these specifications into a NEURON model, facilitates running parallel simulations, automates the optimization and exploration of parameters using supercomputers, and provides a wide range of built-in analysis functions.

Conversion to NetPyNE also makes it easier to connect to previous models developed within the platform, such as our primary motor cortex model (Sivagnanam et al., 2020; Dura-Bernal et al., 2022b), and models implemented in other tools (e.g., NEST) by exporting to the NeuroML or SONATA standard formats. In prior work, we ported a classic model of generic sensory cortical circuits (Potjans and Diesmann, 2014) to our NetPyNE platform (Romaro et al., 2021) in order to make it both more scalable and facilitate modification of cell models and network parameters. The original model used integrate-and-fire neurons and we replaced these with more complex multi-compartment neuron models.

Although we have primarily focused on simplifying the network description, we have also made the model more complex, and more complete, by adding the associated somatosensory thalamic circuits and bidirectional connectivity with cortex to allow interplay of these two highly coordinated areas (Meyer et al., 2010). The deepening of knowledge about the cortico-thalamo-cortical loop (Shepherd and Yamawaki, 2021) should contribute to investigations on rhythmic dysfunctions, such as epilepsy and schizophrenia. But in contrast to cortical microcircuitry, few detailed models exist for the thalamus (Hill and Tononi, 2005; Izhikevich and Edelman, 2008; Murray and Anticevic, 2017; Iavarone et al., 2019).

In this study we present a NetPyNE implementation of the BBP S1 model, capturing most of the original single-cell physiology and morphology, synaptic mechanisms, connectivity and basic simulation results. With the addition of detailed thalamic circuits, we extend the results to show synchronous activity across cortical and thalamic populations, and open the door to new investigations on corticothalamic dynamics. The model is able to port readily across machines and can utilize a fast and efficient implementation on CPUs and GPUs using CoreNEURON. This extension allows the original BBP S1 model to be readily available to be used by the wider community to study a wide range of research questions.

Materials and methods

Individual neuron models

Cell reconstructions were based on the compartmental model Hodgkin-Huxley formalism, with membrane properties represented as components of an electric circuit, and ionic channels modeled as variable conductances. To port the somatosensory microcircuit model in NetPyNE (Dura-Bernal et al., 2019), we recreated the single neuron models using cell files from the Neocortical Microcircuit Collaboration (NMCP)¹ (Ramaswamy et al., 2015). The full dataset comprises 207 morpho-electrical (me) cell types, with 5 examples for each, totaling 1,035 cell models, each stored with morphology file, descriptions of ion channels, and a NEURON HOC template to instantiate the cell, which can be imported directly to NetPyNE (Figure 1). Neuron morphologies from the BBP S1 model (specifically, L1_DLAC, L4_DBC, L23_PC, and L6_TPC_L4) imported into NetPyNE were visualized using the NetPyNE GUI (Figure 1A). The full name of the 207 cell types as well as the corresponding acronym can be found in Supplementary Table 1 in the Supplementary material.

FIGURE 1

Figure 1. Reproduction and validation of BBP S1 cell types in NetPyNE. (A) 3D reconstructions of 4 pairs of m-type example neurons visualized using the NetPyNE graphical user interface: inhibitory cells L1_DLAC and L4_DBC (red), and excitatory cells L23_PC and L6_TPC_L4 (blue). (B–E) Somatic membrane potential of the neurons in (A) under current clamp with amplitude 120% of the neuron firing threshold. NetPyNE results (red) compared to the original BBP model results (blue). For L4_DBC_bIR cells (D) we used a deterministic version of the BBP stochastic potassium channel (StochKv) resulting in divergent results; using same deterministic channel in BBP (BBPdet, blue dotted line) restores the match to NetPyNE results. (F,G) Comparison of BBP and NetPyNE firing rate and time to first spike in response to current-clamp with amplitude 0.1 nA during 2 s for each cell type. Only some cell types with the StochKv show differences. (H) Backpropagating action potential in a L5_TTPC cell. (I) Post synaptic potentials (PSPs) of one dendritic connection with 5 excitatory (blue circles), 10 excitatory (red squares), and 20 inhibitory (cyan diamonds) synapses. (J) Comparison BBP and NetPyNE mean firing rate for all m-type populations. Due to the small number of cells with StochKv (3.6%), NetPyNE population firing rates closely match those of BBP.

Benchmark testing validated physiological responses (Figures 1B–E) at 3 current clamp amplitudes (120%, 130%, and 140% of threshold; only 120% shown). Slight differences were observed in the cell types with a stochastic version of the K⁺ channel mechanism (StochKv; Figure 1D) where we used a deterministic version of the channel from OpenSourceBrain (Gleeson et al., 2019). The StochKv NMODL (.mod) mechanism required additional code outside of NetPyNE in order to update its state, and the inclusion of stochastic variables in each section of the cells significantly increased the simulation time. In order to understand the effect of StochKv on cell response, we applied a current clamp (0.1 nA, 2s) to the soma of each of the 1,035 cells, and used the Electrophys Feature Extraction Library (eFEL)² to compare BBP and NetPyNE mean firing rate (Figure 1F) and time to first spike (Figure 1G) for those with and without the StochKv channel. As expected, the variability for cells with the StochKv channel in the original model was pronounced. Although present in 54/207 me-types, the StochKv channels are only in 3.63% of all cells. Within each m-type (morphology-type) those with StochKv also correspond to a minority of e-types (electrical) types; for example, only 32% of L4_DBC cells have e-type bIR (with StochKv channels). Given the small proportion of cells with StochKv channels (3.63%), the NetPyNE mean firing rates per m-type population closely matched those of the original BBP model (Figure 1J). Furthermore, the deterministic version of StochKv preserved irregular cell spiking patterns (CV_BBP = 0.25 ± 0.16; CV_NetPyNE = 0.16 ± 0.13; where CV is the inter spike interval coefficient of variation; see Supplementary Figure 1) as well as the neural firing rate (FR_BBP = 30.24 ± 24.33, FR_NetPyNE = 28.67 ± 21.20) in the current-clamp simulation with amplitude 0.1 nA during 2 s. For stimulation amplitude 0.8 nA, the CV (BBP = 0.14 ± 0.17; NetPyNE = 0.15 ± 0.30) and FR (BBP = 115.58 ± 70.20, NetPyNE = 110.15 ± 66.58) were similar in both model implementations (Supplementary Figure 1).

The NetPyNE implementation perfectly reproduced the original neuronal intrinsic dynamics since all model parameters were directly imported from the original HOC files, the same NMDOL files were used (except StochKv), and the underlying simulation engine was NEURON in both cases (see Figures 1B–E). To validate this, we simulated somatodendritic backpropagating action potentials (Figure 1H) and dendrosomatic postsynaptic potentials (Figure 1I) in an example L5_TTPC cell. Results were identical in the NetPyNE implementation and the original BBP cell models. To model dendrosomatic postsynaptic potentials (PSPs), we added excitatory connections with 5 and 10 synapses, and an inhibitory connection with 20 synapses, to the L5_TTPC neuron. Additionally, we provided the same three subthreshold inputs within a short time interval, to demonstrate temporal integration of PSPs (Figure 1I).

Distribution and connectivity of cortical populations

Rather than instantiating the connectivity from a list of individual synapses based on anatomical overlap of neuronal arbors (Reimann et al., 2015), we created our S1 port using probability rules for both neuron distribution and connections. The network consisted of 31,346 cells in a cylindrical volume 2,082 μm height and radius of 210 μm as in the original model (Figure 2). Each population was randomly distributed within its specific layer (L1, L2/3, L4, L5, or L6). The number of cells in each one of 207 me-types was taken from the NMCP (Ramaswamy et al., 2015) the minicolumn data available was not used to distribute cells. A 2D representation of the cell distribution within the cylindrical volume is shown in Figure 2A, with layer thicknesses (in μm) for L1, L23, L4, L5, and L6 set to 165, 502, 190, 525, and 700, respectively. We used the S1 connectome (Gal et al., 2017) from NMCP, following the approach in Reimann et al. (2017b): 7 stochastic instances of a model microcircuit based on averaged measurements of neuron densities were used to calculate distance-dependent probabilities of connection. In each microcircuit instance, we calculated the connection probability for each pair of neurons based on the 2D somatic distance (horizontal XZ-plane) for each of the 1,941 pathways. To estimate the distance-dependent probability, we calculated the probability in evenly spaced intervals starting at 25 ± 25 μm, in 50 μm intervals, up to 375 ± 25 μm. Next, we calculated the mean probability across the 7 microcircuits in evenly spaced intervals and used the mean values to fit the connection probability rules. We evaluated multiple functions for each pathway and selected the one that provided the best fit to the data. Figures 2B–D shows how this approach was used to calculate the connection probability of 3 example projections: data from the 7 microcircuit instances (mcs; cyan circles) was averaged across microcircuits (green diamonds) and fitted to either a single exponential (Figure 2B, red line); an exponential with a linear saturation rule (Figure 2C); or a single gaussian (Figure 2D, dashed line). Because the original S1 model shows high variability in the number of synapses per connection, we calculated the mean values for each pathway and used it as a parameter in our model. The result is a representative reconstruction of the S1 column connectivity in NetPyNE, with approximately 27.6 million excitatory synapses and 9.6 million inhibitory synapses.

FIGURE 2

Figure 2. Reproduction and validation of BBP S1 neuron distribution and connectivity in NetPyNE. (A) 2D representation of the location of 31,346 cells in a cylinder with 2,082 μm height and 210 μm radius, with each subtype (different colors) randomly distributed within its layer (L1, L2/3, L4, L5, or L6). (B–D) Probability of connection as a function of neuron pairwise 2D distance for three example pathways, each with a different best fit function: a single exponential (B, red line), exponential with a linear saturation rule (C, red line), and single gaussian fit (D, black dashed line). Cyan circles represent data from 7 microcircuit instances, and green diamonds represent the mean across the 7 instances. (E,F) Comparison of the number of connections between NetPyNE and BBP for each of the 1,941 pathways (E) and 4 projection types (p-types) (F). (G) Postsynaptic potential (PSP) generated by connection between L23_PC neurons (Table 1, #18) simulated in NetPyNE; mean PSP trace (black line) across 20 PSP random instances (gray lines).

Using the fitted connectivity rules, we reconstructed an entire S1 column in NetPyNE and compared the two versions using the mean number of connections. To avoid overfitting, we generated 7 different instances using different connectivity seeds for both the NetPyNE and BBP models. The number of connections was similar in both models for each of the 1,941 pathways (Figure 2E) and for each of the four projection types (p-type) (EE, EI, IE, II) (Figure 2F).

Synaptic physiology

The original BBP S1 model included detailed synaptic properties (conductances, post-synaptic potentials (PSP), latencies, rise and decay times, failures, release probabilities, etc.) recapitulating published experimental data. Short-term dynamics were used to classify synapses into the following types (s-types): inhibitory facilitating (I1), inhibitory depressing (I2), inhibitory pseudo-linear (I3), excitatory facilitating (E1), excitatory depressing (E2), and excitatory pseudo-linear (E3). A set of rules were then derived from experimental data to assign an s-type to each broad class of connections. Based on the NMCP data, there were 29 classes of connections as determined by the combination of pre- and post-synaptic me-types. The synaptic properties, s-type and p-type for each class of connections are summarized in Table 1.

TABLE 1

Table 1. Synaptic properties, s-type, p-type, and rules for each class of connections implemented in NetPyNE.

The dual-exponential conductance model with rise time (τ_rise) 0.2 ms was used for all synapses. Moreover, synaptic properties included the kinetic parameters: peak conductance (g_syn; in nS) and decay time (τ_decay; in ms); and dynamic parameters: neurotransmitter release probability (U), time constant for recovery from depression (D; in ms) and time constant for recovery from facilitation (F; in ms). The NetPyNE implementation reproduces the original PSP amplitudes from Markram et al. (2015). An example of PSPs for a connection between L23_PC neurons (Table 1, #18) simulated in NetPyNE is shown in Figure 2G. The mean PSP peak amplitude across 20 PSPs (with different randomization seeds) was 1.0 mV, which matches the value obtained in Markram et al. (2015). We also included a compact description of the rules to determine what connections belong to each class, based on the pre- and postsynaptic cell types (Table 1). For clarity, we rearranged the classes of connections by s-types in the sequence I1, I2, I3, E1, E2, and E3 (from 0 to 28), and included the original BBP class label for reference. The parameters D and F correspond to the synapses with short term plasticity (STP), which could be optionally added to recurrent S1 connections, and connections from thalamus to S1.

The s-types for each class of connections and for each of the 1,941 pathways are color-coded and illustrated in Figure 3A. Since pathways depend on m-types but connection classes depend on me-types (each m-type includes multiple me-types), it is possible to have multiple s-types for the same pathway; in those cases we simply labeled it as either I2 or E2. To implement the dynamics of each s-type in NetPyNE we used a deterministic version of the dual-exponential synaptic model (Fuhrmann et al., 2002; Hennig, 2013). Example simulations of the PSP for the different s-types are shown in Figure 3B. For each example, we ran 20 simulations with 5 different post-synaptic cells of the same me-type and 4 random synaptic distributions. Pre- and post-synaptic neurons of specific me-types were selected to illustrate each of the six s-types (see Figure 3).

FIGURE 3

Figure 3. Matrix of s-types for each of the 1,941 pathways and simulated PSPs examples for each s-type in NetPyNE. (A) Color-coded s-types for each class of connections (top) and for each of the 1,941 pathways (bottom). Note that for pathways with multiple s-types only either I2 or E2 is shown. (B) Example simulations of post-synaptic potentials to illustrate each of the six s-types. Each example shows the results of 20 simulations with five different post-synaptic cells (different colors) of the same me-type, and 4 random synaptic distributions. Inhibitory s-types were simulated using pathway L23_SBC:L23_PC, which included different s-types depending on pre-synaptic e-type: I1 for e-type cAC, I2 for e-type dNAC, and I3 for e-type bNAC (as shown in Table 1). Excitatory s-types E1 and E2 were simulated using pathway L23_PC:L23_LBC, e-types cAC, and dNAC, respectively; and s-type E3 was simulated using pathway L6_TPC_L4:L6_TPC_L4.

Extending the model to include thalamic populations and connectivity

We extended the model to include somatosensory thalamic populations with cell type-specific dynamics, intra-thalamic connectivity and bidirectional projections with cortex. In the original model, thalamic inputs were modeled as spike generators that only provided feedforward inputs to S1. Our somatosensory thalamus model is composed of the excitatory ventral posterolateral (VPL), ventral posteromedial (VPM) and the posteromedial (POm) nuclei, and the inhibitory reticular nucleus (RTN). We used single compartment cell models with dynamics tuned to reproduce previous studies on the interaction between the thalamic relay and reticular cells (Destexhe et al., 1996a), but adjusted to work in large-scale networks (Moreira et al., 2021). The thalamic circuit architecture consisted of six stacked populations as a rough approximation of the thalamic anatomical layout (Figure 4A). The top three were inhibitory populations comprising the outer, middle and inner sectors of the RTN, and spanning a height of 78, 78, and 156 μm, respectively. Below these were the three excitatory populations, VPL, VPM, and POm, with heights of 156, 156, and 312 μm, respectively. The horizontal dimensions (XZ-plane) for all populations were 420 μm × 420 μm. Cells were randomly distributed across each nuclei with the number of cells in each population based on cellular density obtained from the Cell Atlas for the Mouse Brain³ (Erö et al., 2018). Although POm was larger than VPL and VPM, we reduced its cell density by 50%, resulting in approximately the same population size. This lower density accounts for the proportion of coexisting, but functionally isolated, M1-projecting TC cells present in POm with no projections to S1 (Guo et al., 2020).

FIGURE 4

Figure 4. NetPyNE model of somatosensory thalamic populations and connectivity extending the original S1 model. (A) Distribution of neurons across the six different thalamic populations, roughly mimicking the thalamus anatomy (x and y axes in μm). (B) Schematic of bidirectional connectivity between thalamic regions and cortex. Bidirectional connections between S1/RTN and first order (FO) regions VPL and VPM are topological, whereas those with high order (HO) region POm are non-topological and implemented using divergence rules. (C) Convergence connectivity between thalamic regions and S1 (D) RTN cells have sector-specific distance-dependent connectivity (E) convergence connectivity matrix across all thalamic populations.

The intrathalamic connectivity was based on data of axonal and dendritic footprints for each nucleus. The VPL and VPM are considered first-order nuclei (FO), which means they receive afferent information from peripheral sensory organs (not modeled here) and are interconnected with cortex (Sugitani et al., 1990; Ma, 1991; Luczyńska et al., 2003) and RTN (Lam and Sherman, 2011) in a topological fashion. On the other hand, POm is considered a higher-order (HO) nucleus, so input arrives mainly from the cortex, in this case, from S1 L5 and L6 (Ohno et al., 2012; O’Reilly et al., 2021). The connectivity pattern of HO nuclei has not been properly characterized, but literature reports a decreased level of organization of HO nuclei inputs to RTN (Lam and Sherman, 2011), as it sends projections to S1.

We therefore adopted three connectivity strategies. In the first, neurons from FO nuclei projected to RTN with a column-like topological organization. We implemented this by combining a probability of connection that decreased exponentially with the horizontal distance between the pre- and post-synaptic cells with a decay constant proportional to the footprint radius, and which was truncated to 0 outside of the footprint radius (this denotes the maximum distance of connection in the XZ-plane). The following footprint diameters were derived from experimental data (or estimated in the case of no literature reports) for the different axonal footprints of each thalamic projection: RTN→VPL and RTN→VPM: 64.33 μm (Lam and Sherman, 2007); VPL→RTN: 97.67 μm; and VPM→RTN: 103.57 μm (Lam and Sherman, 2011). The second strategy applies to RTN→RTN connectivity and implements a sector-specific distance-dependent connectivity. More specifically, within each RTN sector, the probability of connection decayed exponentially and was truncated to 0 based on a footprint radius of 264.63 μm (Lam et al., 2006). In strategies one and two, the maximum distance in the Y-plane was set to 10% of the footprint radius, following the disc-like morphology from the axonal projections of the relay cells and the dendritic trees of reticular cells (Murray Sherman and Guillery, 2001; Lam et al., 2006). The third strategy was a divergence rule, with the number of projections from and to HO nuclei having a fixed value and being distributed without spatial constraints. This divergence value was adjusted so that the HO dynamics resembled that of the FO nuclei. This allowed us to replicate a column-like topological organization in FO nuclei using single-compartment cells (Lam and Sherman, 2007), and distribute the HO connections to behave in a functionally similar fashion (Figure 4B).

All excitatory thalamic nuclei were indirectly interconnected through their RTN projections, which was divided into three sectors, in line with reports of preferred innervation zones by each of the thalamic nuclei (Lam and Sherman, 2011). Synapses within RTN were mediated by GABA_a, those from RTN to the excitatory nuclei by a combination of GABAa and GABAb with equal weight, and those from the excitatory nuclei to RTN and cortex by AMPA (Destexhe et al., 1996a). The probability and weight of connections were the targets of parameter optimization. The matrix with the convergence of intra-thalamic connections is shown in Figure 4E.

Feedback corticothalamic connectivity originated from S1 m-types L5_TTPC2 and L6_TPC_L4 (O’Reilly et al., 2021). Similar to the topological rules described above, we implemented connectivity with convergence of 30 (i.e., number of pre-synaptic cells projecting to each post-synaptic cell), but only if the horizontal distance between the pre- and post-synaptic neurons was lower than 50.0 μm (Figure 4).

TC connectivity from VPL, POm and VPM to S1 was implemented using convergence values estimated from previous studies (Meyer et al., 2010; Figure 4C). Convergence values for each of the 55 m-types were calculated based on the weighted average of the populations in each layer. The convergence values for inhibitory populations were multiplied by a scaling factor derived from the original mode (∼0.595). This thalamic convergence factor for inhibitory cells was estimated by dividing the IE ratio of VPM thalamic innervation (83/775 = 0.107) by the average IE population ratio (4,779/26,567 = 0.18). The resulting S1 column received approximately 4.95 M synapses from VPM, 4.95 M from VPL, and 3.1 M from POm. This is consistent with values that can be derived from experimental studies (Meyer et al., 2010), with 4.27 M synapses from VPM and, and 2.66 M from POm. We approximated TC synaptic physiology using the parameters of model #25 in Table 1, and using 9 synapses per connection, following the Markram et al. (2015) characterization of TC synapses as excitatory depressing (E2).

Background inputs

Each cell in the S1 circuit received 10 synaptic inputs from Poisson-distributed spike generators (NetStims) to represent the global effect of spontaneous synapses, background, and other noise sources from non-modeled brain regions projecting to S1. These stimuli were randomly distributed remove across all sections. The quantal synaptic conductance was calculated based on the average quantal conductance for excitatory and inhibitory synapses. We tuned the excitatory and inhibitory stimuli rates using grid search parameter exploration to obtain average excitatory firing rates of ∼1 Hz and physiological firing rates for most S1 populations.

Model building

We used the NetPyNE modeling tool (Dura-Bernal et al., 2019) to build, manage simulations, and analyze results of the S1 and thalamic circuit model. NetPyNE employs NEURON (Carnevale and Hines, 2006; Migliore et al., 2006; Lytton et al., 2016) as backend simulation engine, with either the standard or CoreNEURON libraries (Kumbhar et al., 2019). The high-level Python-based declarative language provided by NetPyNE facilitated the development of this highly complex and extensive circuit model. This language enabled us to easily import existing morphological and biophysical parameters of different cell types, and define complex connectivity and stimulation rules. We used NetPyNE to explore and optimize the model parameters through automated submission and managing of simulations on supercomputers. We also employed NetPyNE’s built-in analysis functions to plot 2D representations of cell locations, connectivity matrices, voltage traces, raster plots, local field potentials (LFPs), 3D synapses representations, and firing rate statistics. NetPyNE can also be used to export the model into the NeuroML (Gleeson et al., 2010) and SONATA (Dai et al., 2020) standard formats.

Model parameters are based on experimental data and the original model (Markram et al., 2015). Nonetheless, parameter optimization was necessary to ensure the model reproduces experimental measures such as population firing rates and PSP. The parameters optimized for the S1 TC circuit were the background rate for excitatory and inhibitory connections. For the intrathalamic projections, we optimized connection weight (range 0–2 mV), connection probability (range 0–1), y-axis connection radius (1%, 2%, 5%, or 10%) and connectivity divergence of the HO populations (5, 10, 20, or 40 cells). For the TC and corticothalamic projections we optimized connection weight (range 0–2 mV) and connection probability (range 0–1). More details about model parameter optimization/exploration are described in section 2 of the Supplementary material.

Simulation of local field potentials

We simulated LFP extracellular recordings using the “line source approximation” (Buzsáki et al., 2012; Łęski et al., 2013; Parasuram et al., 2016), which is based on the sum of the transmembrane currents generated by each segment of each neuron, divided by the distance between the segment and the electrode. This method assumes that the electric conductivity (sigma = 0.3 mS/mm) and permittivity of the extracellular medium are constant everywhere and do not depend on frequency. LFP calculation, analysis and visualization was performed using NetPyNE.

Given the computational cost and memory requirements of simulating the full S1 model with morphologically-detailed neurons while recording LFPs, we calculated transmembrane currents only for the most central cells within an 84 μm (20% of 420 μm) diameter cylinder. This means that only 4.4% of the neurons were simulated in full detail, i.e., using full morphological reconstructions and with all synapses from the full model. The dynamics of the remaining cells (∼96% S1 and thalamus) were simulated using spike generators (VecStims in NEURON) using the spiking activity previously recorded in full simulations. That is, the inputs and activity of the 4% of fully detailed neurons from which the LFPs were calculated were identical to those of the full network simulations (when 100% of the neurons are simulated in detail).

Results

Reproduction of cell morphologies, physiological responses, spatial distribution and connectivity

Cells imported into NePyNE using the files from The Neocortical Microcircuit Collaboration NMCP (Ramaswamy et al., 2015), reproduced the morphological and electrophysiological characteristics of the original model (Figure 1): Mean firing rate and time to the first spike after a current clamp stimulation were fitted for all the 1,035 cell types. Firing dynamic differences were observed in cells with the stochastic K channel (StochKv), but their firing irregularity was partly preserved (Supplementary Figure 1) and, given their low proportion (3.63%), the average firing rates of all m-type populations closely matched those in the original model (Figure 1H).

We were able to recreate the general characteristics across the 7 BBP S1 microcircuit instances: the 31,346 cells were distributed randomly by layer, and probabilistic connections were generated for each of the 1,941 pathways (Figure 2). Here, we replaced the original connectivity method, based on the overlap between axonal and dendritic fields, with one based on connection probability based on cell type, layer, inter-cell distance, and dendritic pattern of post-synaptic locations. This network parameterization allowed us to rescale the microcolumn and generate different instances by changing the random number generator seed. Our probabilistic rules best reproduced the original number of connections using a Gaussian fit in most projection pathways (1,303 of 1,941) and an exponential fit plus a linear saturation in the remaining 638 cases (Figure 2).

Extension to include detailed thalamic circuits

We extended the model to include the somatosensory thalamic populations with projections to S1: RTN, POm, VPL, and VPM. The number of thalamic cells was adapted to fit a cylindrical column with the same radius as the S1 column. This facilitated the inclusion of topological connectivity rules between the two regions. We reproduced the firing dynamics of the different thalamic cell types using a single compartment neuron model (Moreira et al., 2021). The connections from TC cells to S1 were based on convergence rules derived from experimental data (Meyer et al., 2010), and synaptic physiological mechanisms were generalized from the BBP VPM projections to S1 layers 4 and 5 (Markram et al., 2015). Feedback connections originated from S1 cell types L5_TTPC2 and L6_TPC_L4 and targeted VPL and VPM following a topological organization, and POm in a following a non-topological broader distribution (Figure 4). The parameters of the thalamic circuit were adjusted to reproduce a stable self-sustained activity with rhythmic bursting and spindle oscillations (Destexhe and Contreras, 2011), as well as a shift in dynamics following localized excitatory input in the relay cells (Bonjean et al., 2012; Moreira et al., 2021).

Cortical and thalamic circuits independent response to background inputs (no thalamocortical connections)

We first evaluated the response to background inputs of the S1 cortical circuit and the thalamus circuit independently, i.e., without any connections between cortex and thalamus. When driven with background inputs, the S1 model generated spontaneous activity with most populations (48 out of 55) firing within physiological rates (Figure 5). To achieve this, excitatory and inhibitory background inputs were tuned via grid search parameter optimization (see section Materials and methods). Figure 5 illustrates the S1 spontaneous activity results, including a spiking raster plot of all 31,346 cortical cells, examples of voltage traces for each of the 207 me-type population, and the average firing rates for each of the 55 m-type populations. The thalamic populations, disconnected from S1 and driven by background inputs, exhibited stable self-sustained activity with rhythmic bursting at theta ∼6 Hz (Kim and McCormick, 1998; Figures 5C,D). These oscillations, which were most prominent in the RTNi and POm populations, emerged despite the lack of rhythmicity in the background inputs. The thalamic circuit oscillatory dynamics are consistent with the recurrent interactions between thalamic relay and reticular neurons described in previous studies (Destexhe et al., 1996a).

FIGURE 5

Figure 5. NetPyNE S1 and thalamus circuit response to background inputs (spontaneous activity). (A) Spiking raster plot of the 31,346 cells in the S1 column during 1 s (the first second was omitted to allow the network to reach a steady state). (B) Example voltage traces for each of the 207 me-types grouped by rows into their respective 55 m-types (same time period as raster plot). (C) Spiking raster plot of the 7,266 thalamic neurons during 1 s showing intrinsic oscillations. (D) Example voltage traces for each of the 6 thalamic populations. (E) Mean firing rates of each of the 55 m-types for NEURON (red) vs. CoreNEURON (blue). (F) Comparison of the time required to create the network and run the simulation on a 40-core Google Cloud virtual machine using NEURON (red) or CoreNEURON (blue).

Simulations were run using NetPyNE and NEURON on a Google Cloud virtual machine with 40 cores. We compared the S1 results using the standard NEURON simulation engine vs. CoreNEURON, a state-of-the-art solver optimized for large scale parallel simulations on both CPUs and GPUs (Kumbhar et al., 2019). Both simulation engines produced very similar firing rates for each population (Figure 5E), with excitatory and L1-L3 inhibitory cells showing overall lower firing rates than L4-L6 inhibitory cells. The overall average firing rate across the 2 simulated seconds was 0.95 Hz in both cases (NEURON: 59,779 spikes; CoreNEURON: 59,749 spikes). This demonstrates the consistency of results obtained from both simulation engines, making CoreNEURON a viable alternative to study the S1 network. CoreNEURON was 2.4x faster to create the network and 2.2x faster to run the simulation (Figure 5F).

Somatosensory cortex circuit response to background inputs with short term plasticity (no thalamocortical connections)

We simulated the response of the S1 cortical circuit to background inputs but including short term plasticity (STP) in its local synaptic connections (Figure 6A). Adding STP resulted in the emergence of synchronous bursting within the S1 cortical column at approximately 1 Hz frequency (compare S1 raster in Figures 5A, 6A). The spontaneous synchronous bursts first appeared in L5, and then spread to all S1 cells within 100 ms. Figure 6B shows an amplified raster plot of L4-L6 with 70 ms of activity at the time when spontaneous synchronous bursts started. Figure 6C shows example voltage traces of cortical and thalamic neurons, illustrating the spike synchrony of S1 and the thalamic bursts. These results are comparable to the simulations presented in Figures 11B,C of the original publication (Markram et al., 2015).

FIGURE 6

Figure 6. NetPyNE S1 circuit response to background inputs with short-term plasticity (STP). (A) Spiking raster plot of S1 with STP. S1 and thalamus were not interconnected; only intracortical connections were included. (B) Amplified spiking raster plot (A) showing the 70 ms around the time when synchronous bursts first occur in L5 (black) and then propagate to L6 (red) and L4 (blue). (C) Example traces from (A) showing spike synchrony across cortical populations. Rasters in A show 2.5 s after steady state was reached.

Somatosensory cortex and thalamic circuit response with bidirectional thalamic connectivity and cortical short term plasticity

We then simulated the full circuit with bidirectional connections between S1 and thalamus and STP in the thalamus to S1 connections (Figure 7A). The full cortico-thalamo-cortical circuit exhibited overall increased activity with S1 oscillations around 6 Hz frequency, and strong thalamic oscillatory activity at the same frequency. Oscillations were now synchronized across all S1 and thalamic populations. Figure 7B shows the voltage traces of several cortical and thalamic neurons, illustrating the spike synchrony of S1 and thalamic populations. Finally, in Figure 7C we compare the mean firing rate for all S1 and thalamic populations with (red bars) and without (blue bars) bidirectional TC connectivity. All 55 model populations now exhibited physiological firing rates. Adding bidirectional TC connectivity resulted in a modest increase of the overall mean firing rate, from 4.96 to 5.29 Hz, with more pronounced increases in the average firing rates of L1 and L2/3 inhibitory populations. These results do not have a direct correspondence to any in the original BBP publication, since the original model did not include thalamic populations bidirectionally connected to cortex.

FIGURE 7

Figure 7. NetPyNE S1 and thalamic circuit response with bidirectional thalamic connectivity and cortical STP. (A) Spiking raster plot of the fully connected circuit model, including bidirectional connections between S1 and thalamus (shows 2.5 s of simulation after steady state was reached). Oscillations at ∼6 Hz were now synchronized across all S1 and thalamic populations. (B) Example traces from (A) during 800 ms showing spike synchrony across cortical and thalamic populations. (C) Comparison of mean firing rates of each of the 55 S1 and 6 thalamic m-types with (S1-TH connected) and without (S1-TH disconnected) bidirectional thalamocortical connectivity (compare rasters in panel A and Figure 6A, respectively).

Somatosensory cortex and thalamic circuit response after reducing the extracellular calcium concentration to reproduce asynchronous in vivo-like state

Experimental evidence shows that extracellular calcium concentration ([Ca²⁺]_o) in vivo is lower than in vitro, and, as a consequence, PSP amplitudes are also lower (Borst, 2010). Markram et al. (2015) divided the dependency of PSPs on [Ca²⁺]_o into three classes for specific connection types: steep, intermediate, and shallow. Here, the PSP amplitudes were set to have steep dependence for connections between PC-PC and PC-distal targeting cell types (DBC, BTC, MC, BP) and a shallow dependence for connections between PC-proximal targeting (LBCs, NBCs, SBCs, ChC). An intermediate level of dependence was assumed for other connections. To simulate reduced [Ca²⁺]_o in the NetPyNE implementation we decreased the cao parameter from 2.0 to 1.2 in all cells, and modified the use parameter of synaptic transmission (U) adding a factor to multiply its value from 0.25 to 0.75. This resulted in a transition from synchrony (in vitro-like, Figures 6, 7) to asynchrony (in vivo-like, Figure 8) network states, as in Markram et al. (2015). The increased asynchrony happened both for the S1-TH disconnected (Figure 8A) and the S1-TH connected (Figure 8B) cases. Decreasing extracellular calcium concentration resulted in decreased firing rates for most populations (compare Figures 7C, 8C). In both the in vitro and in vivo conditions, cortical firing rates were generally slightly higher for the S1-TH connected case. However, bidirectional thalamic connectivity (S1-TH connected) resulted in increased thalamic population firing rates in vitro, whereas under in vivo conditions (low [Ca²⁺]_o), it lowered thalamic firing rates and decreased synchrony.

FIGURE 8

Figure 8. NetPyNE S1 and thalamic circuit response with low extracellular calcium (in vivo-like asynchronous states). (A) Spiking raster plot of spontaneous activity with only intracortical and intrathalamic connections (S1-TH disconnected). (B) Spiking raster plot of the fully connected circuit model, including bidirectional connections between S1 and thalamus (S1-TH connected). (C) Comparison of mean firing rates of each of the 55 S1 and 6 thalamic m-types without (S1-TH disconnected) and with (S1-TH connected) bidirectional thalamocortical connectivity.

Local field potentials recorded from the in vivo-like somatosensory cortex circuit

We simulated extracellular LFP recordings at multiple depths and horizontal distances in the S1 cortical column during the in vivo-like state (Figure 9). The LFP calculation was based on the transmembrane currents across all segments of neurons. To reduce the computational cost of the calculation, we included only the 1,376 morphologically detailed neurons (4.4% of the total neurons) within a central cylinder of 84 μm diameter (Figures 9A,B). The remaining 29,970 S1 and 7,266 thalamic neurons were simulated using artificial spike generators (VecStims) to ensure the dynamics of the detailed neurons were identical as in the full scale simulation (Figure 9D). The simulated morphologically-detailed neurons therefore included the same 2,702,107 synapses with STP as those in the full in vivo simulation. We inserted recording electrodes at 4 different cortical depths (500, 1000, 1500 and 2000 um) and 2 radial (x-z plane) distances (0 and 297 μm) from the cylinder center (Figure 9C). Recorded LFP amplitudes were in the order 1–1,000 μV consistent with the experimental literature (Reimann et al., 2013; Hagen et al., 2018; Figures 9E,F). The amplitudes of LFPs recorded further away from the cylinder were attenuated ∼10 to 20x compared to those closer to the cylinder center, for example, the peak amplitudes for electrodes 1 and 5 were 401 and 24 μV, respectively. This is consistent with LFP amplitude being inversely proportional to the squared distance between electrode and current sources. When compared to the in vitro recorded LFPs (see Supplementary Figures 2, 3), which exhibited stronger slow frequency oscillations, the attenuation measured at the distant electrodes was only 5x. This is consistent with the observed frequency-dependent attenuation phenomenon, where high frequency signals are attenuated more than low frequency oscillations (Buzsáki et al., 2012; Reimann et al., 2013). The LFP power spectral densities generally depict an inverse relationship between power and frequency, which is typically described in animal LFP recordings. Overall, these preliminary results demonstrate the model can be used to simulate and capture several physiological features of extracellular LFPs.

FIGURE 9

Figure 9. Local field potentials (LFPs) recorded from the in vivo-like S1 circuit. (A,B) Lateral and top-down 2D representation of the location of 1,376 morphological cells within a cylinder with 2,082 μm height and 84 μm radius. Morphologically-detailed neurons are shown in red (L1, L4, and L6) and blue (L23 and L5); whereas the 29,970 simplified cells (spike generators) are shown in orange (L1, L4, and L6) and cyan (L23 and L5) circles. (C) 3D representation of the morphologically-detailed neurons with the location of all synapses (red dots), and the location of the 8 LFP recording electrodes at 4 different depths and 2 radial distances (color triangles). (D) Spiking raster plot of the morphologically-detailed neurons used to calculate the LFP. (E) LFP signals recorded at the 4 electrodes in the center of the cylinder (colors correspond to triangles in panel C). Electrodes numbered 0, 1, 2, and 3 correspond with cortical depths (y) 500, 1,000, 1,500, and 2,000 μm, respectively. (F) Same as E, but for the LFPs recorded at a radial distance of 297 μm; electrodes 4–7. (G–J) Power spectral densities (PSDs) for electrodes 0, 3, 4, and 7 calculated over a 10-s simulation (initial transient period was not included). PSDs exhibit an inverse relationship between power and frequency.

Discussion

We provide here the first large-scale S1 model that is accessible to the wider community, building on the details of the prior state-of-the-art BBP S1 model. The model closely reproduced the original cell morphologies and electrophysiological responses for the 207 morpho-electrical (me) cell types, with 5 examples for each, totaling 1,035 cell models (Figure 1); the spatial distribution of these cells across layers; and the connectivity properties of the 1,941 pathways, including synaptic dynamics and short-term plasticity (Figures 2, 3). After tuning, the simulations produced reasonable dynamics with rates and activity patterns corresponding to in vivo measures of cortical activity (Figures 5, 6). There was no direct comparison to the full network dynamics of the original BBP model since original simulation data was not available. However, firing rates and overall 1 Hz underlying oscillation when using STP was comparable to that seen in the original model version paper (Markram et al., 2015; Figure 11). We also extended the model by adding thalamic circuits, including 6 distinct thalamic populations that reproduced cell and circuit-level dynamics, and with intrathalamic, TC and corticothalamic connectivity derived from experimental data (Figure 4). The addition of the thalamic circuit resulted in distinct activity patterns and synchronous activity across cortical and thalamic populations (Figure 7). Finally, we decreased the extracellular calcium concentration ([Ca²⁺]_o) to simulate in vivo-like states with asynchronous activity (Figure 8). LFPs recorded at multiple cortical depths and horizontal distances exhibited realistic oscillatory patterns and power spectra, including the experimentally observed distance- and frequency-dependent attenuation (Figure 9).

The S1 model now joins other NetPyNE cortical simulations: generic cortical circuits (Romaro et al., 2021), auditory and motor TC circuits (Sivagnanam et al., 2020; Dura-Bernal et al., 2022a,b), as well as simulations of thalamus (Moreira et al., 2021), dorsal horn of spinal cord (Sekiguchi et al., 2021), Parkinson’s disease (Ranieri et al., 2021) and schizophrenia (Metzner et al., 2020). These large cortical simulations can be extremely computer-intensive, which is a major motivation for NetPyNE’s facilities that allow one to readily simplify the network by swapping in integrate-and-fire or small-compartmental cell models, or by down-scaling to more manageable sizes. CoreNEURON is a state-of-the-art solver optimized for large scale parallel simulations, now included as part of the official NEURON package. The optimization on CPUs and the ability to run across GPUs in CoreNEURON is another key NetPyNE feature enhancing runnability. In the present case, the original S1 model is largely inaccessible, despite the cooperation of its designers, since it requires specialized tools, workflows, and training. Nonetheless, most of the data required to replicate it is available via the NMCP, which we ourselves used to implement the NetPyNE version.

We were able to get substantial speedup (> 2x) for both model setup and run using CoreNEURON despite only using CPUs with no GPU at this time. We note that using CPU cycles/timestep would provide a more direct measure than the total simulation time, which may be affected by other factors such as background processes (Girardi-Schappo et al., 2017). Nonetheless, the speedup obtained is consistent with the 2–7x speedups recently reported when using CoreNEURON on CPUs to simulate large-scale models (Kumbhar et al., 2019; Awile et al., 2022). For example, the NetPyNE-based motor cortex model exhibited a speedup of 3.5x on Google Cloud. When using GPUs, speedups of up to 40x were reported. The differences in firing activity seen with NEURON vs. CoreNEURON are expected due to vectorization of the compute kernels in CoreNEURON and potential differences due to different solvers when using NMODL with sympy. Further differences are to be expected once this is extended to GPUs (Jézéquel et al., 2015; Kumbhar et al., 2019).

We made 2 significant changes in our port to NetPyNE. First, we did not replicate the stochastic K channels that appear in 3.6% of the neurons, making our port somewhat simpler than the original. This channel required writing custom code and made simulations slower, but it could be added to the model in a future iteration. Second, we have not utilized the original cell-to-cell connection mappings that were obtained by BBP from direct microscopic observations of overlap between pre-synaptic axonal fields and post-synaptic dendritic fields (so-called Peter’s principle). In the original BBP S1 model, the use of cell-to-cell connections necessarily limited the simulation to use precisely the original model’s cell morphologies, cell positions and scales. It also required storing and loading large files of connection data. We therefore replaced this connection framework with one based on connection probability based on cell type (including layer), inter-cell distance, and dendritic pattern of post-synaptic locations. Although saving somewhat on space, there is a time-space tradeoff since this requires further calculations on start-up. Despite these limitations, we had excellent agreement with both cell model matching and connection density matching.

Our implementation also incorporates a novel model of thalamic circuitry that recapitulates multiple experimental findings at the single neuron and circuit levels. Thalamic and reticular cell models were adjusted to reproduce the reported resting membrane potential, approximately –60 and –80 mV, respectively (Jahnsen and Llinás, 1984; Destexhe et al., 1996b; Sherman and Guillery, 2009). In the networks, TC cells fired at low frequencies (2–4 Hz), while reticular cells fired at higher rates (6–14 Hz), consistent with values previously reported in the literature (Kim and McCormick, 1998). Thalamic simulations also showed rhythmic rebound bursting when hyperpolarized and regular spiking activity at depolarized potentials (Destexhe and Contreras, 2011). The thalamic network exhibited synchronous activity within and across several populations, as well as synchronous firing with cortical populations, particularly in the in vitro condition. These synchronous patterns likely emerged as a consequence of the implemented intrathalamic and TC connectivity, including the topological organization based on axonal footprints (Lam et al., 2006; Lam and Sherman, 2007, 2011). Taken together, these results make the thalamic circuit a valuable extension to the S1 model, by providing a more realistic input source to the cortical circuit and enabling the study of TC interactions.

Recording the intracellular potential of multiple neurons in vivo requires an elaborate set up and is generally challenging. Extracellular recordings are more accessible and therefore more commonly used in experimental studies. Extracellular potentials are generated by transmembrane currents resulting from neuronal activity. Evidence suggests the main contributor to extracellular signals are synaptic currents (Buzsáki et al., 2012; Reimann et al., 2013). Computational modeling coupled with recordings of field activity in animals can provide insights into the cooperative behavior of neurons and increase our understanding of how these processes contribute to the extracellular signal (Buzsáki et al., 2012; Reimann et al., 2013). The simulated LFPs exhibited a similar range of amplitudes as those recorded experimentally, and reproduced several features of LFPs, including the distance-dependent and frequency-dependent attenuation. This opens the door to future validation of the model by comparing LFPs to those recorded experimentally under different conditions, and to future studies of the biophysical sources of LFPs and the exact contribution of different network populations (Hagen et al., 2018).

To place our model in the context of recent literature, we follow the classification proposed by a recent review of data-driven models structural connectivity at the microcircuit level (Shimoura et al., 2021). Our model can be classified as using conductance-based, morphologically-detailed neurons, with a network size of 38,612 neurons, synaptic plasticity and network spatiality (e.g., distance-based connectivity). Our NetPyNE implementation, together with the original BBP implementation (Markram et al., 2015), constitute the only conductance-based, morphologically-detailed models of S1. These contrast with previous models of S1 (Huang et al., 2022) or of generic sensory cortex (Potjans and Diesmann, 2014) that employ simpler neuron models (leaky integrate and fire point neurons). Models with detailed conductance-based and morphologically-detailed neurons have been developed for other cortical regions, including V1 (Arkhipov et al., 2018; Billeh et al., 2020), M1 (Dura-Bernal et al., 2022b), A1 (Dura-Bernal et al., 2022a), and CA1 (Bezaire et al., 2016; Ecker et al., 2020). Our model is also unique in incorporating thalamic neurons and TC bidirectional topological connectivity. Previous TC circuit models included less biophysically-detailed neuron models and simpler connectivity (Izhikevich and Edelman, 2008), or focused on single cell (Iavarone et al., 2019) or small circuit models (Destexhe et al., 1996a). An impressively detailed model of the thalamoreticular microcircuit has recently been developed, although this is limited to the VPL and RTN somatosensory thalamus regions (Iavarone et al., 2022).

As outlined above, the level of biophysical, morphological and connectivity detail in the model is very high compared to most existing models. Although this makes it harder to simulate and tune, it also enables exploration of a unique set of scientific questions that simpler models cannot address, or at least not with the same level of realism. Here we included two results that require and justify the level of detail of the model. First, we simulated a network state with lower extracellular calcium concentration that more closely resembles the in vivo conditions (Figure 8). Secondly, we calculated realistic LFPs, which critically depend on the sum of transmembrane currents along detailed neuronal morphologies (Figure 9). We also describe the methodology for future model parameter explorations, and provide a basic code set up example to explore the effects of inhibitory GABAergic connections on network dynamics. Examples of parameter explorations in NetPyNE-based biophysically detailed circuit models can be found in our related publications on motor and auditory cortex models (Sivagnanam et al., 2020; Dura-Bernal et al., 2022a,b), including an exploration of the effects of long-range and neuromodulatory inputs.

Consequently, our port of the S1 model provides a quantitative framework that can be used in several ways. First, it can be used to perform in silico experiments to explore sensory processing under the assumption of various coding paradigms or brain disease, including the representation of whisker motion (Bosman et al., 2011; Huang et al., 2022), maximization of sensory dynamic range (Gautam et al., 2015), response to unexpected sensory inputs (Amsalem et al., 2020) schizophrenia (Metzner et al., 2020) and Parkinson’s disease (Ranieri et al., 2021). Second, drug effects can be directly tested in the simulation (Neymotin et al., 2016)—this is an advantage of a multiscale model with scales from molecule to network, which is not available in simpler models that elide these details. Third, the model constitutes a unified multiscale framework for organizing our knowledge of S1 which serves as a dynamical database to which new physiological, transcriptomic, proteomic, and anatomical data can be added. This framework can then be utilized as a community tool for researchers in the field to test hypotheses and guide the design of new experiments.

Data availability statement

The model and data for this study can be found in the following repositories and online platforms: GitHub (https://github.com/suny-downstate-medical-center/S1_Thal_NetPyNE_Frontiers_2022), ModelDB (https://senselab.med.yale.edu/ModelDB/), Open Source Brain (https://www.opensourcebrain.org/), and EBRAINS Model Catalog (https://ebrains.eu/).

Author contributions

FB, JM, WL, and SD-B conceived and designed research and drafted the manuscript. FB, JM, LT, and SD-B implemented and optimized the simulation code. FB and JM prepared the figures. All authors contributed to the manuscript revision, read, and approved the submitted version.

Funding

This work was funded by the following grants: NIH NIBIB U24EB028998, NSF 1904444-1042C, NYS SCIRB DOH01-C32250GG-3450000, and NIH NIDCD R01DC012947. This research was funded in part by the Aligning Science Across Parkinson’s (ASAP-020572) through the Michael J. Fox Foundation for Parkinson’s Research (MJFF).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fninf.2022.884245/full#supplementary-material

Footnotes

1. ^ https://bbp.epfl.ch/nmc-portal
2. ^ https://github.com/BlueBrain/eFEL
3. ^ https://bbp.epfl.ch/nexus/cell-atlas/

References

Amsalem, O., King, J., Reimann, M., Ramaswamy, S., Muller, E., Markram, H., et al. (2020). Dense computer replica of cortical microcircuits unravels cellular underpinnings of auditory surprise response. BioRxiv [Preprint]. doi: 10.1101/2020.05.31.126466

CrossRef Full Text | Google Scholar

Arkhipov, A., Gouwens, N. W., Billeh, Y. N., Gratiy, S., Iyer, R., Wei, Z., et al. (2018). Visual physiology of the layer 4 cortical circuit in silico. PLoS Comput. Biol. 14:e1006535. doi: 10.1371/journal.pcbi.1006535

PubMed Abstract | CrossRef Full Text | Google Scholar

Awile, O., Kumbhar, P., Cornu, N., Dura-Bernal, S., King, J. G., Lupton, O., et al. (2022). Modernizing the NEURON simulator for sustainability, portability, and performance. Front. Neuroinform. 16:884046. doi: 10.3389/fninf.2022.884046

PubMed Abstract | CrossRef Full Text | Google Scholar

Azarfar, A., Calcini, N., Huang, C., Zeldenrust, F., and Celikel, T. (2018). Neural coding: A single neuron’s perspective. Neurosci. Biobehav. Rev. 94, 238–247. doi: 10.1016/j.neubiorev.2018.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Barthas, F., and Kwan, A. C. (2017). “Secondary motor cortex: Where ‘sensory’ meets ‘motor’ in the rodent frontal cortex. Trends Neurosci. 40, 181–193. doi: 10.1016/j.tins.2016.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezaire, M. J., Raikov, I., Burk, K., Vyas, D., and Soltesz, I. (2016). Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent CA1 circuit. Elife 5:e18566. doi: 10.7554/eLife.18566

PubMed Abstract | CrossRef Full Text | Google Scholar

Billeh, Y. N., Cai, B., Gratiy, S. L., Dai, K., Iyer, R., Gouwens, N. W., et al. (2020). Systematic integration of structural and functional data into multi-scale models of mouse primary visual cortex. Neuron 106, 388.e–403.e. doi: 10.1016/j.neuron.2020.01.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonjean, M., Baker, T., Bazhenov, M., Cash, S., Halgren, E., and Sejnowski, T. (2012). Interactions between core and matrix thalamocortical projections in human sleep spindle synchronization. J. Neurosci. 32, 5250–5263. doi: 10.1523/JNEUROSCI.6141-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Borst, J. G. (2010). The low synaptic release probability in vivo. Trends Neurosci. 33, 259–266. doi: 10.1016/j.tins.2010.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosman, L. W., Houweling, A. R., Owens, C. B., Tanke, N., Shevchouk, O. T., Rahmati, N., et al. (2011). Anatomical pathways involved in generating and sensing rhythmic whisker movements. Front. Integr. Neurosci. 5:53. doi: 10.3389/fnint.2011.00053

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., Anastassiou, C. A., and Koch, C. (2012). The origin of extracellular fields and currents — EEG, ECoG, LFP and spikes. Nat. Rev. Neurosci. 13, 407–420. doi: 10.1038/nrn3241

PubMed Abstract | CrossRef Full Text | Google Scholar

Carnevale, N. T., and Hines, M. L. (2006). The neuron book. New York, NY: Cambridge University Press.

Google Scholar

Dai, K., Hernando, J., Billeh, Y. N., Gratiy, S. L., Planas, J., Davison, A. P., et al. (2020). The sonata data format for efficient description of large-scale network models. PLoS Comput. Biol. 16:e1007696. doi: 10.1371/journal.pcbi.1007696

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., and Contreras, D. (2011). “The fine structure of slow-wave sleep oscillations: From single neurons to large networks,” in Sleep and anesthesia: Neural correlates in theory and experiment, ed. A. Hutt (New York, NY: Springer New York), 69–105. doi: 10.1007/978-1-4614-0173-5_4

CrossRef Full Text | Google Scholar

Destexhe, A., Bal, T., McCormick, D. A., and Sejnowski, T. J. (1996a). Ionic Mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070. doi: 10.1152/jn.1996.76.3.2049

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Contreras, D., Steriade, M., Sejnowski, T. J., and Huguenard, J. R. (1996b). In vivo, in vitro, and computational analysis of dendritic calcium currents in thalamic reticular neurons. J. Neurosci. 16, 169–185. doi: 10.1523/JNEUROSCI.16-01-00169.1996

PubMed Abstract | CrossRef Full Text | Google Scholar

Dura-Bernal, S., Neymotin, S. A., Suter, B. A., Dacre, J., Schiemann, J., Duguid, I., et al. (2022b). Multiscale model of primary motor cortex circuits reproduces in vivo cell type-specific dynamics associated with behavior. bioRxiv [Preprint]. doi: 10.1101/2022.02.03.479040

CrossRef Full Text | Google Scholar

Dura-Bernal, S., Griffith, E. Y., Barczak, A., O’Connell, M. N., McGinnis, T., Schroeder, C. E., et al. (2022a). Data-driven multiscale model of macaque auditory thalamocortical circuits reproduces in vivo dynamics. bioRxiv [Preprint]. doi: 10.1101/2022.02.03.479036

CrossRef Full Text | Google Scholar

Dura-Bernal, S., Suter, B. A., Gleeson, P., Cantarelli, M., Quintana, A., Rodriguez, F., et al. (2019). NetPyNE, a tool for data-driven multiscale modeling of brain circuits. Elife 8:e44494. doi: 10.7554/eLife.44494

PubMed Abstract | CrossRef Full Text | Google Scholar

Ecker, A., Romani, A., Sáray, S., Káli, S., Migliore, M., Falck, J., et al. (2020). Data-driven integration of hippocampal ca1 synaptic physiology in silico. Hippocampus 30, 1129–1145. doi: 10.1002/hipo.23220

PubMed Abstract | CrossRef Full Text | Google Scholar

Erö, C., Gewaltig, M. O., Keller, D., and Markram, H. (2018). A cell atlas for the mouse brain. Front. Neuroinform. 12:84. doi: 10.3389/fninf.2018.00084

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuhrmann, G., Segev, I., Markram, H., and Tsodyks, M. (2002). Coding of temporal information by activity-dependent synapses. J. Neurophysiol. 87, 140–148. doi: 10.1152/jn.00258.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Gal, E., London, M., Globerson, A., Ramaswamy, S., Reimann, M. W., Muller, E., et al. (2017). Rich Cell-type-specific network topology in neocortical microcircuitry. Nat. Neurosci. 20, 1004–1013. doi: 10.1038/nn.4576

PubMed Abstract | CrossRef Full Text | Google Scholar

Gautam, S. H., Hoang, T. T., McClanahan, K., Grady, S. K., and Shew, W. L. (2015). Maximizing sensory dynamic range by tuning the cortical state to criticality. PLoS Comput. Biol. 11:e1004576. doi: 10.1371/journal.pcbi.1004576

PubMed Abstract | CrossRef Full Text | Google Scholar

Girardi-Schappo, M., Bortolotto, G. S., Stenzinger, R. V., Gonsalves, J. J., and Tragtenberg, M. H. (2017). Phase diagrams and dynamics of a computationally efficient map-based neuron model. PLoS One 12:e0174621. doi: 10.1371/journal.pone.0174621

PubMed Abstract | CrossRef Full Text | Google Scholar

Gleeson, P., Crook, S., Cannon, R. C., Hines, M. L., Billings, G. O., Farinella, M., et al. (2010). NeuroML: A language for describing data driven models of neurons and networks with a high degree of biological detail. PLoS Comput. Biol. 6:e1000815. doi: 10.1371/journal.pcbi.1000815

PubMed Abstract | CrossRef Full Text | Google Scholar

Gleeson, P., Cantarelli, M., Marin, B., Quintana, A., Earnshaw, M., Sadeh, S., et al. (2019). Open source brain: A collaborative resource for visualizing, analyzing, simulating, and developing standardized models of neurons and circuits. Neuron 103, 395–411.e5. doi: 10.1016/j.neuron.2019.05.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, K., Yamawaki, N., Barrett, J. M., Tapies, M., and Shepherd, G. M. G. (2020). Cortico-Thalamo-cortical circuits of mouse forelimb S1 are organized primarily as recurrent loops. J. Neurosci. 40, 2849–2858. doi: 10.1523/JNEUROSCI.2277-19.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Hagen, E., Næss, S., Ness, T. V., and Einevoll, G. T. (2018). Multimodal Modeling of neural network activity: Computing LFP, ECoG, EEG, and MEG signals with LFPy 2.0. Front. Neuroinform. 12:92. doi: 10.3389/fninf.2018.00092

PubMed Abstract | CrossRef Full Text | Google Scholar

Hennig, M. H. (2013). Theoretical models of synaptic short term plasticity. Front. Comput. Neurosci. 7:154. doi: 10.3389/fncom.2013.00154

PubMed Abstract | CrossRef Full Text | Google Scholar

Hill, S., and Tononi, G. (2005). Modeling sleep and wakefulness in the thalamocortical system. J. Neurophysiol. 93, 1671–1698. doi: 10.1152/jn.00915.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Zeldenrust, F., and Celikel, T. (2022). Cortical representation of touch in silico. Neuroinformatics doi: 10.1007/s12021-022-09576-5 [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Iavarone, E., Simko, J., Shi, Y., Bertschy, M., García-Amado, M., Litvak, P., et al. (2022). Thalamic control of sensory enhancement and sleep spindle properties in a biophysical model of thalamoreticular microcircuitry. bioRxiv [Preprint]. doi: 10.1101/2022.02.28.482273

CrossRef Full Text | Google Scholar

Iavarone, E., Yi, J., Shi, Y., Zandt, B. J., O’Reilly, C., Van Geit, W., et al. (2019). Experimentally-constrained biophysical models of tonic and burst firing modes in thalamocortical neurons. PLoS Comput. Biol. 15:e1006753. doi: 10.1371/journal.pcbi.1006753

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M., and Edelman, G. M. (2008). Large-scale model of mammalian thalamocortical systems. Proc. Natl. Acad. Sci. US.A. 105, 3593–3598. doi: 10.1073/pnas.0712231105

PubMed Abstract | CrossRef Full Text | Google Scholar

Jahnsen, H., and Llinás, R. (1984). Ionic basis for the electro-responsiveness and oscillatory properties of guinea-pig thalamic neurones in vitro. J. Physiol. 349, 227–247. doi: 10.1113/jphysiol.1984.sp015154

PubMed Abstract | CrossRef Full Text | Google Scholar

Jézéquel, F., Lamotte, J.-L., and Saïd, I. (2015). “Estimation of numerical reproducibility on CPU and GPU,” in Proceedings of the 2015 Federated Conference on Computer Science and Information Systems, (Piscataway, NJ: IEEE), 675–680. doi: 10.15439/2015f29

CrossRef Full Text | Google Scholar

Kim, U., and McCormick, D. A. (1998). The functional influence of burst and tonic firing mode on synaptic interactions in the thalamus. J. Neurosci. 18, 9500–9516. doi: 10.1523/JNEUROSCI.18-22-09500.1998

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumbhar, P., Hines, M., Fouriaux, J., Ovcharenko, A., King, J., Delalondre, F., et al. (2019). CoreNEURON : An optimized compute engine for the NEURON Simulator. Front. Neuroinform. 13:63. doi: 10.3389/fninf.2019.00063

PubMed Abstract | CrossRef Full Text | Google Scholar

Lam, Y. W., and Sherman, S. M. (2007). Different topography of the reticulothalmic inputs to first- and higher-order somatosensory thalamic relays revealed using photostimulation. J. Neurophysiol. 98, 2903–2909. doi: 10.1152/jn.00782.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Lam, Y. W., and Sherman, S. M. (2011). Functional organization of the thalamic input to the thalamic reticular nucleus. J. Neurosci. 31, 6791–6799. doi: 10.1523/JNEUROSCI.3073-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lam, Y. W., Nelson, C. S., and Sherman, S. M. (2006). Mapping of the functional interconnections between thalamic reticular neurons using photostimulation. J. Neurophysiol. 96, 2593–2600. doi: 10.1152/jn.00555.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Łęski, S., Lindén, H., Tetzlaff, T., Pettersen, K. H., and Einevoll, G. T. (2013). Frequency Dependence of signal power and spatial reach of the local field potential. PLoS Comput. Biol. 9:e1003137. doi: 10.1371/journal.pcbi.1003137

PubMed Abstract | CrossRef Full Text | Google Scholar

Luczyńska, A., Dziewiatkowski, J., Jagalska-Majewska, H., Kowiański, P., Wójcik, S., Labuda, C., et al. (2003). Qualitative and quantitative analysis of the postnatal development of the ventroposterolateral nucleus of the thalamus in rat and rabbits. Folia Mophol. 62, 75–87.

PubMed Abstract | Google Scholar

Lytton, W. W., Seidenstein, A. H., Dura-Bernal, S., McDougal, R. A., Schürmann, F., and Hines, M. L. (2016). Simulation neurotechnologies for advancing brain research: Parallelizing large networks in NEURON. Neural Comput. 28, 2063–2090. doi: 10.1162/NECO_a_00876

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, P. M. (1991). The barrelettes–architectonic vibrissal representations in the brainstem trigeminal complex of the mouse. I. Normal structural organization. J. Comp. Neurol. 309, 161–199. doi: 10.1002/cne.903090202

PubMed Abstract | CrossRef Full Text | Google Scholar

Markram, H., Muller, E., Ramaswamy, S., Reimann, M. W., Abdellah, M., Sanchez, C. A., et al. (2015). Reconstruction and simulation of neocortical microcircuitry. Cell 163, 456–492. doi: 10.1016/j.cell.2015.09.029

PubMed Abstract | CrossRef Full Text | Google Scholar

McDougal, R. A., Bulanova, A. S., and Lytton, W. W. (2016). Reproducibility in Computational neuroscience models and simulations. IEEE Trans. Biomed. Eng. 63, 2021–2035. doi: 10.1109/TBME.2016.2539602

PubMed Abstract | CrossRef Full Text | Google Scholar

Metzner, C., Mäki-Marttunen, T., Karni, G., McMahon-Cole, H., and Steuber, V. (2020). The effect of alterations of schizophrenia-associated genes on gamma band oscillations. bioRxiv [Preprint]. doi: 10.1101/2020.09.28.316737

CrossRef Full Text | Google Scholar

Meyer, H. S., Wimmer, V. C., Hemberger, M., Bruno, R. M., de Kock, C. P., Frick, A., et al. (2010). Cell type-specific thalamic innervation in a column of rat vibrissal cortex. Cereb. Cortex 20, 2287–2303. doi: 10.1093/cercor/bhq069

PubMed Abstract | CrossRef Full Text | Google Scholar

Migliore, M., Cannia, C., Lytton, W. W., Markram, H., and Hines, M. L. (2006). Parallel network simulations with NEURON. J. Comput. Neurosci. 21, 119–129.

Google Scholar

Moreira, J. V. S., Borges, F. S., Doherty, D., Lytton, W. W., and Dura-Bernal, S. (2021). Topographically detailed computational model of the motor and somatosensory thalamic circuits. Available online at: https://www.abstractsonline.com/pp8/#!/10485/presentation/16321 (accessed February 25, 2022).

Google Scholar

Murray Sherman, S., and Guillery, R. W. (2001). “Chapter II – the nerve cells of the thalamus,” in Exploring the thalamus, eds S. Murray Sherman and R. W. Guillery (San Diego, CA: Academic Press), 19–58. doi: 10.4324/9781315152837-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Murray, J. D., and Anticevic, A. (2017). Toward understanding thalamocortical dysfunction in schizophrenia through computational models of neural circuit dynamics. Schizophr. Res. 180, 70–77. doi: 10.1016/j.schres.2016.10.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Neymotin, S. A., Dura-Bernal, S., Moreno, H., and Lytton, W. W. (2016). Computer modeling for pharmacological treatments for dystonia. Drug Discov. Today Dis. Models 19, 51–57.

Google Scholar

Ohno, S., Kuramoto, E., Furuta, T., Hioki, H., Tanaka, Y. R., Fujiyama, F., et al. (2012). A morphological analysis of thalamocortical axon fibers of rat posterior thalamic nuclei: A single neuron tracing study with viral vectors. Cereb. Cortex 22, 2840–2857. doi: 10.1093/cercor/bhr356

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Reilly, C., Iavarone, E., Yi, J., and Hill, S. L. (2021). Rodent Somatosensory thalamocortical circuitry: Neurons, synapses, and connectivity. Neurosci. Biobehav. Rev. 126, 213–235.

Google Scholar

Parasuram, H., Nair, B., D’Angelo, E., Hines, M., Naldi, G., and Diwakar, S. (2016). Computational modeling of single neuron extracellular electric potentials and network local field potentials using LFPsim. Front. Comput. Neurosci. 10:65. doi: 10.3389/fncom.2016.00065

PubMed Abstract | CrossRef Full Text | Google Scholar

Peña-Rangel, T. M., Lugo-Picos, P. I., Báez-Cordero, A. S., Hidalgo-Balbuena, A. E., Luma, A. Y., Pimentel-Farfan, A. K., et al. (2021). Altered sensory representations in parkinsonian cortical and basal ganglia networks. Neuroscience 466, 10–25.

Google Scholar

Petrof, I., Viaene, A. N., and Sherman, S. M. (2015). Properties of the Primary somatosensory cortex projection to the primary motor cortex in the mouse. J. Neurophysiol. 113, 2400–2407. doi: 10.1152/jn.00949.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Potjans, T. C., and Diesmann, M. (2014). The Cell-type specific cortical microcircuit: Relating structure and activity in a full-scale spiking network model. Cereb. Cortex 24, 785–806. doi: 10.1093/cercor/bhs358

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramaswamy, S., Courcol, J. D., Abdellah, M., Adaszewski, S. R., Antille, N., Arsever, S., et al. (2015). The neocortical microcircuit collaboration portal: A resource for rat somatosensory cortex. Front. Neural Circuits. 9:44. doi: 10.3389/fncir.2015.00044

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranieri, C. M., Montino Pimentel, J., Romano, M. R., Elias, L. A., Romero, R. A. F., Lones, M. A., et al. (2021). A data-driven biophysical computational model of Parkinson’s disease based on marmoset monkeys. IEEE Access 9, 122548–122567.

Google Scholar

Reimann, M. W., Anastassiou, C. A., Perin, R., Hill, S. L., Markram, H., and Koch, C. (2013). A biophysically detailed model of neocortical local field potentials predicts the critical role of active membrane currents. Neuron 79, 375–390. doi: 10.1016/j.neuron.2013.05.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Reimann, M. W., King, J. G., Muller, E. B., Ramaswamy, S., and Markram, H. (2015). An algorithm to predict the connectome of neural microcircuits. Front. Comput. Neurosci. 9:120. doi: 10.3389/fncom.2015.00120

PubMed Abstract | CrossRef Full Text | Google Scholar

Reimann, M. W., Horlemann, A., Ramaswamy, S., Muller, E., and Markram, H. (2017b). Morphological diversity strongly constrains synaptic connectivity and plasticity. Cereb. Cortex 27, 4570–4585. doi: 10.1093/cercor/bhx150

PubMed Abstract | CrossRef Full Text | Google Scholar

Reimann, M. W., Nolte, M., Scolamiero, M., Turner, K., Perin, R., Chindemi, G., et al. (2017a). Cliques of neurons bound into cavities provide a missing link between structure and function. Front. Comput. Neurosci. 11:48. doi: 10.3389/fncom.2017.00048

PubMed Abstract | CrossRef Full Text | Google Scholar

Romaro, C., Najman, F. A., Lytton, W. W., Roque, A. C., and Dura-Bernal, S. (2021). NetPyNE Implementation and rescaling of the potjans-diesmann cortical microcircuit model. Neural Comput. 33, 1993–2032. doi: 10.1162/neco_a_01400

PubMed Abstract | CrossRef Full Text | Google Scholar

Sekiguchi, K., Medlock, L., Dura-Bernal, S., Prescott, S. A., and Lytton, W. W. (2021). Multiscale computer model of the spinal dorsal horn reveals changes in network processing associated with chronic pain. bioRxiv [Preprint]. doi: 10.1101/2021.06.09.447785

CrossRef Full Text | Google Scholar

Shepherd, G. M. G., and Yamawaki, N. (2021). Untangling the cortico-thalamo-cortical loop: Cellular pieces of a knotty circuit puzzle. Nat. Rev. Neurosci. 22, 389–406. doi: 10.1038/s41583-021-00459-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, S. M., and Guillery, R. W. (2009). Exploring the thalamus and its role in cortical function, 2nd Edn. Cambridge, MA: Mit Press.

Google Scholar

Shimoura, R. O., Pena, R. F. O., Lima, V., Kamiji, N. L., Girardi-Schappo, M., and Roque, A. C. (2021). Building a model of the brain: From detailed connectivity maps to network organization. Eur. Phys. J. 230, 2887–2909.

Google Scholar

Sivagnanam, S., Gorman, W., Doherty, D., Neymotin, S. A., Fang, S., Hovhannisyan, H., et al. (2020). Simulating large-scale models of brain neuronal circuits using google cloud platform. PEARC20 (2020) 2020, 505–509. doi: 10.1145/3311790.3399621

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugitani, M., Yano, J., Sugai, T., and Ooyama, H. (1990). Somatotopic organization and columnar structure of vibrissae representation in the rat ventrobasal complex. Exp. Brain Res. 81, 346–352. doi: 10.1007/BF00228125

PubMed Abstract | CrossRef Full Text | Google Scholar

Vázquez, Y., Salinas, E., and Romo, R. (2013). Transformation of the neural code for tactile detection from thalamus to cortex. Proc. Natl. Acad. Sci. U.S.A. 110, E2635–E2644. doi: 10.1073/pnas.1309728110

PubMed Abstract | CrossRef Full Text | Google Scholar