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.

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

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.

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

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.

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

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

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

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.

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.

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

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

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

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

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

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

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.

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.

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.