Visual physiology of the layer 4 cortical circuit in silico

Abstract

Despite advances in experimental techniques and accumulation of large datasets concerning the composition and properties of the cortex, quantitative modeling of cortical circuits under in-vivo-like conditions remains challenging. Here we report and publicly release a biophysically detailed circuit model of layer 4 in the mouse primary visual cortex, receiving thalamo-cortical visual inputs. The 45,000-neuron model was subjected to a battery of visual stimuli, and results were compared to published work and new in vivo experiments. Simulations reproduced a variety of observations, including effects of optogenetic perturbations. Critical to the agreement between responses in silico and in vivo were the rules of functional synaptic connectivity between neurons. Interestingly, after extreme simplification the model still performed satisfactorily on many measurements, although quantitative agreement with experiments suffered. These results emphasize the importance of functional rules of cortical wiring and enable a next generation of data-driven models of in vivo neural activity and computations.

About the Authors:

Anton Arkhipov

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0003-1106-8310

Nathan W. Gouwens

Roles Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0001-8429-4090

Yazan N. Billeh

Roles Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0001-5200-4992

Sergey Gratiy

Roles Data curation, Formal analysis, Software, Visualization, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0001-8911-2321

Ramakrishnan Iyer

Roles Data curation, Formal analysis, Software, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Ziqiang Wei

Roles Data curation, Formal analysis, Investigation, Software, Writing – review & editing

Affiliation: Janelia Research Campus, Howard Hughes Medical Institute, Ashburn, Virginia, United States of America

ORCID logo http://orcid.org/0000-0002-8694-1260

Zihao Xu

Roles Data curation, Formal analysis, Investigation, Software, Writing – review & editing

Affiliation: University of California San Diego, La Jolla, CA, United States of America

Reza Abbasi-Asl

Roles Formal analysis, Investigation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0001-7824-4628

Jim Berg

Roles Data curation, Investigation, Resources, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Michael Buice

Roles Formal analysis, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Nicholas Cain

Roles Software, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0002-5848-199X

Nuno da Costa

Roles Data curation, Formal analysis, Investigation, Methodology, Resources, Visualization, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Saskia de Vries

Roles Data curation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0002-3704-3499

Daniel Denman

Roles Data curation, Formal analysis, Investigation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Severine Durand

Roles Data curation, Investigation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

David Feng

Roles Formal analysis, Software, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Tim Jarsky

Roles Data curation, Investigation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Jérôme Lecoq

Roles Data curation, Resources, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Brian Lee

Roles Data curation, Investigation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Lu Li

Roles Data curation, Investigation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Stefan Mihalas

Roles Resources, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Gabriel K. Ocker

Roles Formal analysis, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0001-9627-9576

Shawn R. Olsen

Roles Data curation, Formal analysis, Investigation, Methodology, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

R. Clay Reid

Roles Resources, Supervision, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Gilberto Soler-Llavina

Roles Data curation, Writing – review & editing

Affiliation: Novartis Institutes for BioMedical Research, Cambridge, MA, United States of America

ORCID logo http://orcid.org/0000-0002-1838-401X

Staci A. Sorensen

Roles Data curation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Quanxin Wang

Roles Data curation, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Jack Waters

Roles Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Massimo Scanziani

Roles Methodology, Resources, Supervision, Writing – review & editing

Affiliation: Howard Hughes Medical Institute and Department of Physiology, University of California San Francisco, San Francisco, California, United States of America

Christof Koch

Roles Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing

Affiliation: Allen Institute for Brain Science, Seattle, Washington, United States of America

Introduction

Although our knowledge of the cortex has been improving dramatically thanks to the ongoing revolution in experimental neuroscience methods, the field is still far from an overall understanding of cortical circuits and their specific function. One essential component necessary to address this problem is the development of data-driven quantitative models that integrate experimental information and enable predictive simulations under a wide range of realistic in-vivo-like conditions–following the dictum attributed to Richard Feynman, “What I cannot create, I do not understand” [1]. Whereas detailed data-driven models of cortical tissue have been reported, in particular, [2] and [3], modeling applications at the biophysical level to an in-vivo-like regime have been fewer (although see, e.g., [4] and references therein).

A typical systems neuroscience experiment involves a battery of different stimuli and, ideally, perturbations of the investigated circuit. Reproducing this in simulations of a data-constrained cortical model has proven challenging. To investigate the feasibility of in-vivo-like comprehensive simulations, and to build a platform for further studies, we set out to simulate a set of visual physiology experiments in the mouse primary visual cortex (area V1), with the emphasis on the first step in the cortical processing of visual information–namely, modeling the V1 input layer, the layer 4 (L4). We decided early on to replicate a small set of what we consider to be canonical physiological findings characterizing cells in L4 of V1. Given the thousands or more of published experiments carried out over the years in this region of cortex, our list is small, non-exclusive and may be considered idiosyncratic by some. However, we believe it is critical to start somewhere solid before generalizing indiscriminately. Our list of explored phenomena includes the approximately log-normal distributions of firing rates [5], orientation selectivity [6, 7, 8, 9], oscillatory population dynamics [10, 11, 12, 13, 14, 15], sparsity of responses to natural stimuli [16], amplification of thalamic inputs by recurrent connections [17, 18, 19, 20, 21, 22, 23], preferential connectivity among similarly tuned neurons [24, 25, 26, 27, 28, 29], and a number of others.

The model was constructed in a data-driven fashion from what is known about the L4 circuit organization. Indeed, we think of it as a consistent summary of our collective anatomical and physiological knowledge about this region of the nervous system. Of course, it is not the most compact such summary (e.g. in terms of Kolmogorov complexity) nor was it meant to be that. To the extent that we can reproduce physiology across scales–from post-synaptic potentials to local field potentials–we would argue that we understand the phenomena that we observe experimentally.

Although the model was biophysically and anatomically detailed, we also used simplifications when appropriate, typically choosing computationally inexpensive approximations for biological mechanisms. A crucial component was a set of filters that represented visual information processing from image to the output of the lateral geniculate nucleus (LGN) of the thalamus, which projects to L4. This feature enabled one to use arbitrary movies as visual stimuli. We presented the same or similar sets of stimuli to the model and to mice in experiments and then compared the in silico and in vivo responses.

We asked three major questions: (1) How well does our model reproduce experimentally observed neural responses from the above list? (2) What are the major mechanisms that determine the neuronal activity patterns? And (3), how does the ability to reproduce experimental recordings depend on the level of granularity of the model?

To address (1), we assessed neuronal responses to artificial (e.g., drifting gratings) and naturalistic (e.g., movies) stimuli and selected a number of features of these responses that are generally considered important and interesting in the field. We then benchmarked the model performance on these features against the experimental data. Whereas typically models are developed to explain a specific phenomenon and may aim to reproduce 1–2 observed quantities, the key element in our study was to observe generalization over a wide variety of visual stimuli and response features. We found that our simulations reproduced many of experimental observations (with some exceptions) under a range of different stimuli.

For (2), we performed in silico experiments to investigate how individual neurons process inputs from different sources and how recurrent connections shape the network activity. This approach relied in part on simulated optogenetic experiments paralleling in vivo optogenetic studies. The most striking observation was that tuning properties of neurons were critically affected by the functional connectivity rules.

For question (3), we introduced two much simplified versions of our model, where biophysical neuron models were replaced by point-neuron models with either instantaneous or time-dependent synaptic action, and compared simulations of these simplified networks to the biophysical model. We found that, although the simplified network models qualitatively reproduced the trends observed in the biophysical simulations (as also reported by [30]), the quantitative agreement with experiment suffered. The time-dependent synaptic kinetics in the simplified model allowed for better agreement with the biophysical model and experiment, such as, for example, in terms of producing oscillations in a gamma range.

No circuit in the brain exists in isolation, but including all the brain complexity in the model is impossible at present due to absence of data and inadequate computing resources. Therefore, we attempted to build a model of L4 of V1, which can be handled with available resources and for which a substantial amount of information could be found. The good performance of the model compared to the experiment may indicate a relative compartmentalization of L4 computation (see Discussion); this should not be expected everywhere in the cortex. Our L4 model provides a comprehensive characterization of activity and mechanisms in this cortical circuit and may serve as a stepping-stone for future, more sophisticated studies of all cortical layers. To enable this, we make the software code, the model, and simulation results publicly available (see SI).

Results

Construction and optimization of the model

The network (Fig 1A and 1B) consisted (see Methods) of models of individual neurons [31] from an early version of the Allen Cell Types Database [32], employing compartmental representation of somato-dendritic morphologies (~100–200 compartments per cell) and 10 active conductances at the soma that enabled spiking and spike adaptation. Although recent additions to the Allen Cell Types Database include models of neurons with active conductances in the dendrites as well, those models are very computationally expensive, which was prohibitive for the breadth of our study (see below). In addition, in terms of somatic spike output, the current versions of such models do not exhibit much better performance than the models with active conductances restricted to the soma [32]. Thus, we used these latter, computationally cheaper models. The cells were distributed uniformly in a cylinder ~400 μm in radius, representing the central portion of V1, and 100 μm height with density of 200,000 mm-3 [33]. Five single neuron models represent five “types” of neurons–three major excitatory groups as determined by Cre-lines (Scnn1a, Rorb, Nr5a1, 85% of all cells) and two groups of parvalbumin-positive fast-spiking interneurons (15%, denoted as PV1 and PV2), which form the majority of interneurons in L4 [34]. All 10,000 biophysical cells were exact copies of these five models. The cell models correspond to regular-spiking excitatory cells and fast-spiking interneurons (PV+); whereas non-PV+ interneurons do exist in L4, they are a relative minority [34], and therefore were neglected for simplicity. Furthermore, 35,000 much simpler leaky-integrate-and-fire (LIF) neurons, with only two groups–excitatory and inhibitory–were placed around the biophysically detailed “core” to prevent boundary artefacts (see SI). The complete model accounted for over half of V1 L4 cells (45,000 out of ~70,000). Below, we primarily focus on the properties of the biophysical core circuit.

[Figure omitted. See PDF.]

Fig 1. Construction of the model.

(a) The biophysical and LIF portions of the model on the cortical surface with delineations of cortical areas (top; “VISp” is V1; the other labels starting with “VIS” correspond to higher visual areas) and with individual cells rendered (bottom; only a small subset of cells is shown for clarity). (b) Morphologies and action potential shapes of the five neuron models used to generate the L4 network; numbers of cells of each type are listed. (c-e) Connection probability (c,d) and synaptic weights (e) of excitatory (E) or inhibitory (I) cell targets and sources. The rules incorporate the dependence on the distance between somata (c) or difference of the assigned preferred orientations, Δori (d, e). Rules dependent on Δori were applied to E-to-E connections only, and synaptic weights for all connections were independent of distance. (f) Connection probability computed for cells within 50 μm using actual preferred orientation observed in simulations (cf. [24]). Numbers of connected and total pairs, used to obtain the probability, are shown inside the bars. (g) Three types of LGN filters (ON, OFF, and ON/OFF), superimposed onto an image, providing inputs to a L4 cell. Example filter’s spatial and temporal profiles are at the bottom. (h) The proportion of excitatory thalamocortical synapses (VGLUT2+) in the neuropil of V1 L4, as determined experimentally using EM. Proportions for individual samples of tissue are in gray; mean and s.e.m. for each mouse are in black. Right, an exemplar EM image of a putative VGLUT2+ synapse from the LGN onto an L4 neuron (Sp, spine; Bt, bouton; arrows: postsynaptic densities within the spine).

https://doi.org/10.1371/journal.pcbi.1006535.g001

Three independent instantiations were generated using different random seeds. The connectivity and inputs into these three model instantiations were distinct, but all followed the rules described below. All simulations were performed using the Python 2.7 code (an early version of the BioNet package [35]) employing NEURON 7.4 [36].

Recurrent connections (Fig 1C) were established randomly according to probability decaying with intersomatic distance (e.g., [37]). Recent experiments [25, 26, 27, 29] showed that excitatory neurons in the mouse V1 L2/3 exhibit more likely and stronger connections if neurons prefer similar stimuli (“like-to-like connectivity”). We sought to investigate how such a rule affects neural activity. The rule (Fig 1D) was applied to all excitatory cells by assigning a preferred orientation to each excitatory neuron (which was also used to select LGN inputs, see below) and using the difference between such orientations to compute a probability of connection. A similar like-to-like rule was applied to the amplitude of synaptic weights (Fig 1E). The like-to-like rules were parameterized to correspond approximately to the observations for L2/3 [25] (Fig 1F) and did not apply to pairs containing one or more inhibitory neurons. See Methods (section “Connectivity within the Network” and “Synaptic characteristics”) for further details.

We directly converted visual stimuli to spikes of LGN cells (which provide the input to L4) via space-time separable linear-nonlinear Poisson (LNP) filters (Fig 1G; see, e.g., [38, 39, 40]) of ON, OFF, and ON/OFF types (see Methods), all of which produced relatively transient responses. Such a considerable simplification of the wide variety of response types that the LGN cells exhibit (see, e.g., [41, 42, 17]) was necessary to make the model tractable, especially because rules of connectivity from various functional LGN cell types and L4 in the mouse are largely unknown. Because the retinotopy of the complete model did not cover the entire field of view, we employed 9,000 LGN filters–approximately half of the mouse LGN. Based on the experimental data about receptive fields of L4 neurons [6, 7, 17], we used retinotopy of L4 cells to “pool” inputs from LGN filters with similar retinotopy, and the assigned preferred orientation (see above) to establish the geometry of the ON and OFF subfields (Fig 1G; see Methods). This helped to establish orientation selectivity in the combined LGN inputs to individual L4 cells [17]. Due to the distance dependence of recurrent connections and retinotopy of LGN inputs, the L4 cells that were connected to each other had a higher chance to receive inputs from the same LGN cells, in comparison with unconnected L4 cells (see Methods). This is consistent with recent experimental observations for L4 and L2/3 cells [27].

The numbers of synapses for recurrent connections were chosen based on the literature (see Methods). Using electron microscopy (EM), we found (Fig 1H) that LGN synapses constitute 15–30% of all synapses in the neuropil in L4 (see Methods for details), consistent with observations for mouse S1 and M1 [43]. This corresponds to over 1,000 synapses from the LGN per excitatory L4 neuron (given ~8,000 synapses total per mouse V1 cell [33] on average, we assume relatively small L4 neurons to receive ~6,000 synapses). Thus, the number of LGN-to-L4 synapses is an order of magnitude higher in the mouse than in the cat V1 (for which this number has been long established to be ~100 [44]). These new data determined the numbers of LGN-to-L4 synapses in the model (see Methods for further details). For the majority of connections, multiple synapses per connection were assigned, typically in the 3–7 range (see Methods). Synaptic dynamics was described by a sum of two decaying exponentials. We assumed no short-term plasticity as it is plausible that in vivo synapses operate at steady-state conditions [18]. Long-term plasticity was neglected because simulations covered several seconds at most.

To account for inputs from the rest of the brain and for different brain states, we introduced externally generated waves of background activity sweeping across the L4 model (Fig 2A), inspired by reports of activity propagation over cortical surface in vivo [45]. Poisson spike generators were distributed within the cortical plane and were activated whenever the wave swept through their positions, sending spikes to the nearby L4 cells. This enabled the L4 model (Fig 2B) to operate with an extremely simple analogue of distinct cortical states (see, e.g., [46, 47, 48]), which we call “Bkg. on” and “Bkg. off.” (see Methods). The random generation of the background waves, as well as random generation of background spike trains from the wave profiles and of LGN spike trains from deterministic filter outputs were the sole sources of indeterminacy, leading to differences between individual trials.

[Figure omitted. See PDF.]

Fig 2. Examples of simulated responses to various visual stimuli.

(a) Spontaneous activity in the model was generated by waves of background excitation (“Bkg. on”, yellow arrow denotes the direction of motion of the yellow bar-like region) alternated with intervals of no background excitation (“Bkg. off”). (b) Model activity in a gray screen trial. Examples of membrane potential traces from simulated biophysical cells are in the middle. (c) Spike raster in response to a drifting grating (TF = 4 Hz, at 0 degrees direction). Bottom, example orientation tuning curves for an excitatory and inhibitory cell from simulation. (d) Spikes in response to a 50 ms full-field flash. (e) Spike raster for a single trial of a natural movie (top) and for a temporally scrambled version of the same movie (bottom). The raster in (b) shows all neurons and those in (c-e) for clarity show the 10,000 cells in the biophysical core of the model (inset on top of (c) zooms in on 200 cells). All rasters are examples from one trial; all trials used unique combinations of “Bkg. on” and “Bkg. off” states (shown at the bottom of plots), which overlapped (or not) in different ways with the visual stimuli.

https://doi.org/10.1371/journal.pcbi.1006535.g002

After the steps above, all parameters were fixed except the synaptic weights (conductance amplitudes). The LGN-to-L4 weights were selected to produce experimentally observed excitatory currents from LGN [17]. The weights for recurrent connections were then manually optimized, while constraining post-synaptic potentials and currents (PSPs and PSCs) to the experimentally reported range for mouse cortex [49] (see Methods). The weights were not adjusted individually; instead, we uniformly scaled all weights belonging to a particular connection type, such as, for example, excitatory-to-Rorb connections or inhibitory-to-Nr5a1 connections. Optimization was limited to two training stimuli–a single trial of a drifting grating and one of a gray screen, 500 ms each. The target was to reproduce the mean spontaneous firing rate and the rate in response to the preferred grating (Rmax), while avoiding synchronous epileptic-like activity. Published data from [6] were used for this optimization, before our own experiments were finalized. After optimizing the synaptic weights using this limited training set, we applied the model without any further modification to all the other stimuli, which served as a test set. The exact synaptic weights obtained during the optimization stage may not represent a unique solution, due to degeneracy [50]. Nevertheless, it is reassuring that synaptic properties from our optimized models were consistent with published experimental reports (see S2 Fig), exhibiting, e.g., current amplitudes of ~10–30 pA for excitatory and ~40–60 pA for inhibitory synapses [49] or voltage amplitudes of ~1 mV [51] for L4 excitatory-to-excitatory connections (as measured at the soma).

Benchmarking simulated neuronal activity to the experimental data

How well does the model reproduce experimentally observed neural responses in vivo? The answer will depend on the types of stimuli and features of responses to those stimuli that one chooses for comparison. We reasoned that a successful model should be versatile; in other words, it is not sufficient to be able to model responses to one type of stimuli well, and instead one should strive to reproduce features across many types of stimuli. Therefore, our battery of visual stimuli included gray screen, a variety of drifting gratings, 10 natural images, 3 natural movie clips, 2 types of full-field flashes, and 4 moving-bar stimuli (see example responses in Fig 2C–2E; see also Methods, S3 and S4 Figs, and S2 Table). Altogether, ~3,600 simulations were carried out. Our own experiments used for benchmarking consisted of two series of extracellular electrophysiological recordings in mouse V1 employing multi-electrode silicon probes. In the first [7], responses to drifting gratings with a variety of spatial and temporal frequencies were measured, in both awake and anesthetized mice. In the second (previously unpublished), spontaneous activity (responses during gray screen presentation) as well as responses to natural movies, natural images, full-field flashes, and other stimuli were recorded in awake mice.

We chose several features of visual responses that are generally considered important and are often reported in the field of cortical visual physiology and compared their values between experiment and simulations (Fig 3). It should be noted that, although the measured metrics differed slightly among the three excitatory cell types (Scnn1a, Rorb, and Nr5a1), as well as among the two inhibitory types (PV1 and PV2), we do not ascribe significance to these differences. This is because during model construction and optimization no experimental data was available on distinctions in response features or connectivity between these L4 cell populations, and, in absence of either, we assumed the same optimization targets for all excitatory or inhibitory cells (which, however, allowed for ~1 Hz rate deviation from the target, resulting in different types settling at slightly different activity levels).

[Figure omitted. See PDF.]

Fig 3. Benchmarking the simulation results.

(a) Log-normal like distribution of firing rates. Top, firing rates of all biophysical cells from three models during spontaneous activity and in response to a drifting grating, averaged over 10 trials. Large red dots are firing rates averaged over all cells in 0.1 Hz bins over the spontaneous activity axis. Bottom, examples of firing rate distributions on a log scale for a single trial. Solid lines indicate log-normal fits of the data. (b) Comparison of the spontaneous rates, maximal rates in response to gratings (Rmax), orientation selectivity index (OSI), and direction selectivity index (DSI) between the simulation (by cell type, color) and experimental measurements using extracellular electrophysiology (gray). Here and in all other box plots, the box bottom and top mark the inter-quartile range (IQR), the median is in red, whiskers mark +/-1.5 IQR; mean and standard deviation are shown in black (see Methods for further details). “An.”–experiments in anesthetized mice, “Aw.”–in awake mice. (c) The local field potential (LFP; see Methods) measured at the center of the L4 model, for a drifting grating. The spectra from 10 trials are shown in gray, and the averaged spectrum is in black. (d) The model PSTH in response to a 50 ms flash (average over all biophysical excitatory cells, all models, and all trials, in 2 ms bins). (e) The magnitude and time-to-peak (from flash onset) for the first and second peaks of the response to the 50 ms flash, for both simulation and electrophysiological data. (f) Distributions of lifetime sparsity of simulated and experimental responses of excitatory neurons, computed for three directions of a grating (0, 45, and 90 degrees) and for three movies. See Methods for details on computing all values presented.

https://doi.org/10.1371/journal.pcbi.1006535.g003

The first characteristic we checked is whether neuron firing rates follow a positively-skewed, log-normal-like distribution, which is a ubiquitous hallmark feature of brain activity in vivo [5]. Such distributions were indeed observed in our simulations for spontaneous activity and all types of visually-driven responses (Figs 3A and S5A). Consistent with published reports, our simulated rate distributions spanned 2–3 decades and may widen further in longer simulations (for instance, spontaneous rates were computed from 20 trials of 500 ms each, resulting in the lowest possible rate of 0.1 Hz). Also in agreement with literature [5], individual neurons tended to keep their low or high firing rates across different stimulus conditions (Fig 3A, top). The positive skewness of firing rate distribution (i.e., the third standardized moment, ⟨((f −⟨f⟩)/σf)3⟩, where f is the firing rate and σf is its standard deviation, both computed in the linear (i.e., not log) firing-rate space) was typically between 0 and 4 (S5A Fig). While it is hard to expect an exact match of the experimental skewness distributions from our relatively rough model, it is reassuring that the model reproduced the overall experimental trends of firing rate distributions with positive skewness in the 0–4 range. Note that a normal distribution has zero skewness, whereas log-normal distribution has positive skewness.

The mechanisms underlying log-normal like distributions of activity of individual neurons in vivo are not well understood [5], and may involve both cell-intrinsic and network phenomena. One major mechanism proposed based on point-neuron network simulations [52] involves a transformation of an approximately normal distribution of total inputs to a cortical cell by a rectified and expansive input-output nonlinearity, which results in a heavy-tailed distribution of firing rates. This, combined with the log-normally distributed strength of background inputs in our model (see S1C Fig), likely underlie our observations.

Another aspect of activity that we checked was whether the magnitude of responses was consistent with the experiment; presumably, it is important for network computations that cells of certain types fire at certain rates in response to particular stimuli. For this purpose, we considered the rates of spontaneous activity and maximal rates (Rmax) in response to drifting gratings (Fig 3B; see Methods for definitions). The spontaneous rates were 0.5–1.0 Hz for excitatory and 1.5–2.0 Hz for inhibitory cells, broadly consistent with our experimental measurements and published data [6], and Rmax levels were also similar to experimental ones.

Furthermore, the L4 responses to gratings are known to exhibit orientation and direction selectivity (e.g., [6, 7]). We found that the orientation selectivity index (OSI; see Methods) for excitatory cells was ~0.4–0.5 on average, in excellent agreement with the experiments (Fig 3B). For inhibitory cells, OSI was ~0.05 on average, whereas the experimental value was ~0.3. The relatively low OSI values were due to non-selective thalamic inputs and recurrent connectivity for inhibitory cells (see Methods), which were chosen that way at model construction to conform to an often-expressed view that excitatory cells are tuned and inhibitory fast-spiking cells are not. However, recent experiments [34], as well as our own data (Fig 3B), suggest moderate levels of tuning for these neurons, and therefore future models will need to reflect these observations better. Nevertheless, qualitatively the model already captures the trend of better tuning of excitatory cells compared to inhibitory fast-spiking cells.

Another often-observed phenomenon is that the orientation tuning width of cortical excitatory neurons stays approximately constant with respect to contrast (see, e.g., [6, 8, 38, 53]). This is considered important because simple feedforward models typically produce broadening of tuning curves with contrast, and, thus, contrast invariance may reflect mechanisms inherent to cortical processing. Interestingly, it has been shown that responses in L4 of mouse V1 do not simply maintain tuning width, but actually exhibit sharpening of tuning curves and increase of OSI with contrast for excitatory cells (and broadening of tuning curves/reduction of OSI for inhibitory cells) [9]. We found strong orientation tuning in excitatory cells for contrast as low as 10% in our model (S5C Fig). The changes in tuning curve width (characterized as differences of the half-width at half-height at high and low contrasts, ΔHWHH) were centered relatively narrowly around zero (S5D Fig; 4+/-5 degrees for 80% vs. 10% and –2+/-4 degrees for 80% vs. 30% contrasts on average). Thus, our model of excitatory and inhibitory recurrent connections was adequate for preventing broadening of tuning curves, but not sufficient to produce sharpening of tuning curves at high contrasts [9]. We further characterized OSI changes with contrast and found on average slight increase in OSI of excitatory cells and more substantial decrease in OSI of inhibitory cells (S5E Fig). Whereas this trend is consistent with experimental observations [9], it is not sufficiently strong in our model, most likely due to the insufficiently sharp tuning of inhibitory cells mentioned above. Thus, the model captures the overall tendency of cortical circuits to prevent broadening of tuning curves with contrast, but should be improved in the future to capture the more delicate properties of tuning sharpening and broadening in distinct cell populations.

By contrast to OSI, the model performed poorly on the direction selectivity index (DSI; Fig 3B), which was close to zero in simulations, whereas experimentally it was ~0.4–0.5. This was due to the extreme simplicity of the LGN inputs, which were constructed based on experimental data, but at this point did not include all types of LGN activity observed experimentally [7, 17, 21]. In particular, for simplicity, we did not include the sustained LGN responses (i.e., all LGN filters produced transient responses), whereas the most recent experimental data suggest a critical role for the interplay between sustained and transient LGN inputs in generating direction selectivity in V1 [23]. Proper incorporation of direction selectivity is the subject of ongoing work on the next generation of the model.

Another important aspect of neuronal activity in vivo is the population-level oscillatory rhythms, which are observed in a variety of frequency bands. Whereas oscillations at many of such frequencies are likely caused by non-local interactions in the brain, and therefore cannot be expected to arise in an isolated L4 model, some of them may be generated locally. Indeed, our simulations exhibited oscillations in the 15–50 Hz range (sometimes referred to as “mouse gamma”), with a peak at ~20 Hz (Fig 3C). This is consistent with extracellular electrophysiology data [10], which exhibits a particularly strong peak at similar frequencies in L4.

Global luminance changes present another challenge to models, due to simultaneous engagement of all cells by these stimuli. We studied responses to 50 ms long full-field flashes and observed, in both experiment and simulations, that a sharp and fast peak in activity was evoked by the first white-to-gray transition, followed by a second smaller and wider peak (Figs 2D and 3D). The magnitude of the first peak was about 2–4 times higher in simulation than in the experiment, whereas that of the second peak was approximately the same, and the time course of the response was uniformly ~2-fold faster in simulations (Fig 3D and 3E and S5B Fig). These differences were most likely again caused by the absence of sustained LGN responses in the model. Importantly, however, in both simulations and experiments, the second peak was delayed–instead of 50 ms (the flash duration), it appeared 100 ms after the first (200 ms in the experiment). This reflects a known phenomenon of suppression following luminance change, where the second peak corresponds to release from inhibition (e.g., [54, 55]). Thus, qualitatively the model reproduces well the critical features of cortical responses to global luminance change, which likely affect perception of visual stimuli (such as in, e.g., luminance-induced visual masking).

Finally, artificial stimuli such as gratings are quite different from natural images and movies and it is an important question to ask whether a model reproduces differences in response statistics between artificial and natural stimuli. We did observe in simulations that responses to a natural movie (e.g., Fig 2E) were very distinct from those to a horizontally drifting grating (Fig 2C; both panels use the same ID labels). For a grating, cells with certain IDs respond strongly throughout the simulation (as they prefer the grating’s orientation). By contrast, natural movies evoked episodes of concerted responses in distinct populations of cells (which share similar orientation preference). Notably, such concerted responses were mostly absent when movie frames were shuffled in time (Fig 2E, bottom). To quantify the differences, we computed lifetime sparsity for each cell following Vinje and Gallant [56] (see Methods). Sparsity was higher for movies than for gratings (Fig 3F) in both simulations and experiments, consistent with previous observations (including Ca2+ imaging [16]): 0.77+/-0.16 in simulation vs. 0.71+/-0.22 in experiments for movies and 0.64+/-0.14 vs. 0.65+/-0.25 for gratings. Thus, the model correctly accounts for the distinction in sparsity between these two classes of stimuli.

To investigate how trial-to-trial neural variability in the L4 model responses compares with experimental observations for gratings and natural movies, we computed several metrics (see Methods for details). The coefficients of variation of inter-spike intervals (S6A Fig) were consistent between the experiment and the model, as was the trend that spontaneous activity and drifting gratings produced lower values than natural movies. The model was also consistent with reports of stimulus-induced suppression of variability [57], as the Fano factor (S6B Fig) tended to be lower for stimulus-evoked responses (especially, for gratings) than for spontaneous activity. This effect, however, was weak for excitatory neurons in the model and was not clearly identifiable in our experimental recordings (S6B Fig), where it might have been obscured by the noise due to relatively small number of recorded cells. The signal correlations were widely distributed in both simulation and experiment (S6C Fig), with the median close to 0 for gratings and close to 1 for spontaneous activity. The signal correlations for natural movies tended to be higher than for gratings, although only slightly so in the experiment and more substantially in simulation. The noise correlations (S6D Fig) were more narrowly distributed and on average close to zero in both simulation and experiment.

Mechanisms underlying neural activity and computation in the L4 circuit

We further used the model to shed light on the mechanisms that determine patterns of neural activity and computations performed by the L4 circuit. One important question is to what extent the L4 activity and computations are inherited from the input regions (here, LGN) vs. being shaped by intrinsic, recurrent connectivity (see, e.g., [58]). To investigate this, in addition to regular (“Full”) simulations, we performed a number of simulations where recurrent and background connections were removed, so that neurons received LGN inputs only (“LGN only”). Using in silico voltage clamp (see Methods), we measured (Fig 4A) synaptic excitatory currents at the cell body–from the LGN only, ILGN, and the total, Itot – and found that, for preferred directions of 2 Hz drifting gratings, the fraction of Itot contributed by ILGN was 0.41+/-0.05 for excitatory cells, in good agreement with experiment (0.36+/-0.02 [17]). Also in agreement with experiment [17], the mean ILGN was untuned, since individual LGN filters were mostly untuned [17] (but see [19, 20]), whereas mean Itot and Isub = Itot = ILGN (i.e., the current due to recurrent connections) were well tuned to the grating orientation (Fig 4B), and F1 components (i.e., the amplitude of the mode at the stimulus frequency; see Methods) of all of these currents were tuned. Besides these features that were consistent with experimental recordings [17], we observed one distinction: the F1 component of Isub was substantially smaller than that of ILGN or Itot (Fig 4B) and the temporal dynamics of Isub was in antiphase with that of ILGN (S7A Fig). This turned out to be a space clamp artefact, where a stronger ILGN current increases the membrane voltage at the synapses, resulting in weaker driving force for recurrent excitation and, therefore, antiphase relationship of ILGN and Isub; simulations where LGN current was removed from the recorded cells exhibited Isub with no oscillations at the grating frequency (S7A Fig). Thus, the overall picture was that the ILGN mean was untuned whereas its F1 component was; the recurrent connections added current that was tuned overall, but was not time-modulated (whereas in the experiment [17] it is time-modulated); and the resulting total current was tuned in both mean and F1. The inhibitory currents were mostly untuned at the level of both the mean and F1 (Fig 4B); the F1 did show a preference to orientation, but its magnitude was only a few percent of the mean and the current was not strongly modulated in time (S7A Fig).

[Figure omitted. See PDF.]

Fig 4. Mechanistic characterization of the model.

(a) Cortical amplification of the LGN inputs. The excitatory currents (from the LGN only, as well as total) in biophysical cells were measured using voltage clamp recordings. Top–an example; bottom–distributions of LGN contribution to the total excitatory current across excitatory and inhibitory cells (computed for each cell as the average current over time and over all trials of the preferred orientation). (b) Tuning curves for the mean and F1 component of the total and LGN-only currents, and their difference (“Sub”, i.e., the cortical component), as well as inhibitory current. The data for each cell were normalized to the peak value of the “Total” and shifted so that the preferred direction is at 0 degrees; averages and s.e.m. over all recorded excitatory cells are shown (TF = 2 Hz, contrast 80%). The inhibitory currents were normalized and aligned to their own peak values, since their magnitude is significantly higher than that of excitatory currents. (c) Amplification of excitatory current. Top, the total current vs. the LGN-only current, for an individual Rorb cell (each point is an average over time and over 10 trials). Linear fits (Itot = A ILGN + B) are shown for data aggregated from all grating directions, TFs, and contrasts (black), for one selected direction (yellow), and for a fixed contrast and TF (i.e., representing a sample direction tuning curve; right plot). Bottom, summary of linear fits across all cells analyzed. (d) Tuning curves for mean firing rate in full network simulations (“Full”, red) and in simulations where all connections except the feedforward connections from the LGN were removed (“LGN only”, blue). The data for each cell were normalized to the peak value of the “Full” and shifted so that the preferred direction is at 0 degrees; averages and s.e.m. over all excitatory cells are shown (TF = 2 Hz, contrast 80%). (e) Simulations of responses to a drifting grating, with the LGN activity switched off at 1000 ms. The black curve is the firing rate averaged over all cells, models, and trials; green is the exponential fit. (f) Distribution of the optogenetic modulation index (OMI) by cell type in responses to gratings, for simulations of optogenetic silencing of the Scnn1a population (top). Combined distribution for all biophysical excitatory cells is compared to the experimental result (bottom).

https://doi.org/10.1371/journal.pcbi.1006535.g004

Comparing mean Itot and ILGN of individual cells (Fig 4C) for all grating conditions (i.e., not only the preferred one shown in Fig 4A) and different contrasts, we observed a complex dependence. The relationship of Itot with ILGN was essentially linear along the contrast dimension. The quality of the linear fit,was r2 = 0.76 for excitatory cells for all conditions, and 0.97 for a specific grating direction. This is consistent with the prediction of an earlier, much simpler self-consistent model [58] (see also [59]). For the orientation dimension, the dependence was highly non-linear (r2 = 0.18). The latter was the consequence of the mean of ILGN being untuned to orientation, whereas the tuned current from recurrent connections made the mean of Itot highly tuned (Fig 4B). Due to the same reason, the amplification factor A was approximately 2 across all conditions (since B was very small on average, A can be thought of as the inverse of the LGN contribution described above; across all conditions 1/A = 0.54+/-0.04, Fig 4C), whereas for the preferred orientation it was 1/0.41 = 2.44 (as the LGN contribution at the preferred orientation was 0.41+/-0.05, Fig 4A).

How do these relationships between currents shape the properties of the spiking output? We found that the OSI in L4 at the level of spikes was inherited from the combined LGN inputs to each cell [17], producing weakly selective output, and was strongly increased by recurrent connections (Figs 4D and S7B). Generally, the relationship between the rates in Full and LGN only simulations (S7C Fig) was not close to linear (attempting a linear fitwhere fFull and fLGN are the firing rates in Full and LGN only simulations, and the parameters Af and Bf are marked with the subscript “f” to distinguish them from the parameters of linear fit of currents above, we found r2 = 0.5 for all conditions and 0.65 for one grating direction for excitatory cells). The firing rates at fixed contrast and TF exhibited a much more linear relationship (r2 = 0.88), but only locally–the good linear fit required negative values of parameter Bf (-4+/-3 Hz on average), which is non-physiological given that Bf is the firing rate of the cell under the conditions when the firing rate induced by the LGN-only input is zero. Thus, overall the relationship was non-linear, but close to linear past the firing threshold for LGN-only inputs. Interestingly, the overall effect of recurrent connections on the spiking output in our model was that of suppression, as the Full-network firing rates at non-preferred orientations tended to be smaller than LGN-only rates (Fig 4D and S7C Fig, top right), and only at the preferred orientations Full-network firing reached the same rates as in the LGN-only case. Consequently, Rmax values were approximately the same in the Full and LGN-only cases (S7B Fig).

To investigate the L4 circuit mechanisms further, we performed in silico optogenetic silencing of LGN inputs and of a subset of the circuit (see Methods). Despite the significant recurrent amplification, L4 activity shut down rapidly when the LGN spiking was silenced (Fig 4E), in agreement with an analogous experiment [18]. Although the time course of activity decay was faster in the model (2.3 vs. 9+/-2 ms, likely due to the absence of excitatory stimulation from other cortical regions), otherwise the effect of LGN silencing was the same. The widespread and powerful intracortical inhibition appears to be the most plausible driving force behind this effect. We then separately silenced the Scnn1a population of excitatory cells in simulations and conducted analogous experiments in vivo (using Archaeorhodopsin-mediated silencing on randomly selected trials during presentation of TF = 2 Hz drifting gratings, while extracellular multielectrode recordings were performed, see Methods). Results were characterized by converting firing rates with (fp) and without (fcontrol) perturbation to the optogenetic modulation index (OMI),

Simulations qualitatively agreed with experiment (Figs 4F, S8A and S8B) in terms of the OMI distribution over L4 excitatory cells (few inhibitory cells were recorded experimentally), which was approximately bimodal, with one “lobe” concentrated close to -1 (totally silenced neurons) and the other near 0 (weak or no effect). Due to challenges of recording form the thin layer 4, a small number of cells was obtained in the experiment, but they all showed a consistent trend of OMI values belonging to one of these two lobes (S8A Fig). Whereas weak and/or sparse silencing of Scnn1a cells in simulations did not result in a two-lobe distribution, a strong and dense silencing did (S8C Fig), consistent with the experimental perturbation that attempted to silence as many Scnn1a cells as possible. Simulations showed (Fig 4F) that the left lobe in the OMI distribution was due to near-complete silencing of Scnn1a cells. This silencing reduced the amount of excitation in the network and resulted in moderate decrease of firing of inhibitory cells. The net effect is that the firing rates of the other excitatory cell types (Rorb and Nr5a1) remained almost unaffected–yielding the OMI lobe with values close to zero.

We investigated the effect of the logics of connections between L4 excitatory cells on the neural activity and computation. In our regular models, a “like-to-like” (“L”) rule [25, 26, 24, 29] was used for both the connectivity and synaptic weights of excitatory-to-excitatory connections, based on the anticipated orientation tuning of the cells (Fig 1C–1E); because the rule applied to both synaptic connectivity and amplitude distributions, we refer to this set as “LL”. We then studied three alternative sets of models. In one, both the connectivity and weights were randomly (“R”) assigned independently of cell tuning (“RR”). The remaining two sets had the random rule applied to connectivity and like-to-like to weights (“RL”), or vice versa (“LR”). Each set consisted of 3 models. Besides the “R” or “L” rules, everything else was exactly the same between the four sets, including the probabilistic distance-dependent connectivity (Fig 1C).

Synaptic weights for all models were tuned following the same procedure (tuning only involved scaling of weights uniformly across populations, thus not affecting the “L” vs. “R” property, which applies to individual cell pairs), resulting in similar levels of activity for the training stimulus (grating) in all models. We then assessed how functional properties differed. We found that OSIs were extremely reduced in the RR set, with LR being only slightly better tuned, whereas RL came close to the original well-tuned LL set (Fig 5A). Because the mean LGN current was untuned in all models, the critical amplification of orientation selectivity (Fig 4B–4D) came from the lateral L4 connections, which were well tuned in LL and RL sets and poorly tuned in LR and RR sets (Fig 5B). As a result, the total current as well as spiking output was well tuned in LL and RL sets and barely tuned in LR and RR (Fig 5). This suggests a bigger role of synaptic weights than connection probability in shaping network responses [60], which can be easily understood since the “L” rule for weights results in low contributions from non-like-to-like connections, even if such connections are present, thus effectively enforcing the “L” type connectivity.

[Figure omitted. See PDF.]

Fig 5. Like-to-like vs. random connectivity and synaptic weights.

(a) Distribution of OSI for biophysical excitatory cells for the LL, RL, LR, and RR cases. (b) Tuning curves for the mean total and LGN-only currents, and their difference (“Sub”, i.e., the cortical component). The data for each cell were normalized to the peak value of the “Total” and shifted so that the preferred direction is at 0 degrees; averages and s.e.m. over all recorded excitatory cells are shown (TF = 2 Hz, contrast 80%).

https://doi.org/10.1371/journal.pcbi.1006535.g005

Performance of a simplified model

How does the ability to reproduce experimental recordings depend on the level of granularity of the model? To address this question, we built two versions of a radically simplified model, where biophysical neurons and bi-exponential synapses were replaced by LIF neurons and either instantaneous “charge-dump” synapses (using NEURON’s IntFire1 function) or synapses with exponential time dependence of synaptic current (NEURON’s IntFire4 function). These all-LIF models used exactly the same cell-to-cell connections and inputs as the biophysical models and were optimized using the same protocol (see Methods). By contrast to another recent study [30], our coarse-graining procedure did not aim to replicate results of the biophysical simulations, but, rather, aimed to match certain experimental observables using the point-neuron networks in the same way as we did for the biophysical network; nevertheless, both that study and ours agree in that results of the simplified models were quite similar to those of the biophysical models.

The results of the all-LIF simulations were qualitatively similar to the results of biophysical simulations (Fig 6), but specific quantitative distinctions were obvious. The overall appearance of responses (Fig 6A) and mean values of spontaneous rates and Rmax were similar to the biophysical case (Fig 6B). However, there were apparent differences in the way spontaneous rates were distributed (Fig 6B): the bottom three quartiles of cells in all-LIF simulations exhibited no spontaneous firing, whereas in the biophysical case that was true for half or less of the cells; nevertheless, the means across cells were similar. The OSI values were elevated in the IntFire1 case (Fig 6B). Most noticeable, the gamma oscillation was largely absent in the IntFire1 models (Fig 6C). This was likely due to instantaneous synapses in IntFire1, since gamma oscillations [11, 12, 13, 14, 15] are thought to be strongly dependent on the properties of inhibitory perisomatic synapses [11, 14], and especially their time constants. Thus, removal of the appropriate synaptic kinetics from the large and distributed network model like ours may be expected to disrupt gamma oscillations. By contrast, the IntFire4 model employed a simple form of synaptic kinetics (see Methods) and produced oscillations in the range of 10–20 Hz, i.e., close to ~20 Hz gamma oscillation in the biophysical model (Fig 6C; note that the weighted multi-unit activity metric used here is only a proxy for the actual LFP that can be computed for the biophysical model, as shown in Fig 3C, but not for point-neuron models; nevertheless, both metrics show similar trends for the biophysical model, such as the large peak at ~20 Hz).

[Figure omitted. See PDF.]

Fig 6. Comparison of biophysical and all-LIF simulations employing the IntFire1 or IntFire4 models.

(a) An example spike raster in response to a drifting grating in an all-LIF IntFire1 simulation. (b) Spontaneous rate, Rmax, OSI and DSI by cell type. (c) Spectra of multi-unit activity (weighted by 1/r, where r is the distance from the cell to the center of the system). This is used as a proxy to LFP, which cannot be directly computed for the point-neuron models. Note that exact match between this metric and the actual LFP computed for biophysical model (see Fig 3C) is not expected, but they both exhibit similar features, such as a prominent peak at ~20 Hz. (d) Distributions of lifetime sparsity of responses to gratings and movies, averaged over 10 trials. The data are for three directions of a grating (0, 45, and 90 degrees) and for three movies. (e) Distributions of LGN contribution to the total excitatory synaptic inputs across excitatory and inhibitory cells (in the IntFire1 and IntFire4 cases, this is computed for each cell as the average synaptic input over time and over all trials of the preferred orientation; see Methods), for TF = 2 Hz drifting gratings. (f) Distribution of OSI for excitatory cells for the LL, RL, LR, and RR cases.

https://doi.org/10.1371/journal.pcbi.1006535.g006

Further distinctions were observed for sparsity values (Fig 6D). The overall trend of higher sparsity for movies than for gratings was reproduced by both all-LIF models, but they both also exhibited higher sparsity values than the biophysical case and the experiment (Fig 3F). Perhaps the most significant difference with the biophysical model was observed for the magnitude of cortical amplification (Fig 6E) for excitatory cells (interestingly, inhibitory cells did not show substantial differences). The LGN contribution was 0.41+/-0.05 in the biophysical model and 0.36+/-0.02 in experiment [17], whereas it was 0.53+/-0.13 for IntFire1 and even higher, 0.68+/-0.12, for the otherwise more realistic IntFire4. The distributions of amplification values were also quite different in all-LIF models (Fig 6E), exhibiting multiple peaks–apparently for the multiple excitatory cell types–instead of a single peak in the biophysical model. Finally, the functional connectivity had the same consequences in the all-LIF models as in the biophysical case–the orientation tuning was high for models with the LL or RL rules, and low for LR and RR (Fig 6F). On the other hand, the small difference between the LL and RL cases, observed in the biophysical model, was mostly eliminated, and overall the OSIs were higher for the IntFire1 model and, in the RR and LR cases, for IntFire4 as well.

Discussion

The promise of data-driven neuroscience modeling lies in harnessing computing power to establish a platform for discovery that would work hand-in-hand with experiments. For that promise to materialize, models need to be biologically realistic in recapitulating knowledge about brain structure as well as in reproducing in vivo activity. Doing either is difficult, and doing both together is more difficult yet, and is less widely practiced (for some powerful examples, see [2, 3, 4]). We described here a model of the L4 in mouse V1 that was built using a high degree of biological realism and simulated in a framework of an in silico visual physiology experiment. We asked how well the model reproduces activity observed in vivo across a variety of stimuli, which mechanisms underlie the activity and computation in the modeled L4 circuit, and how the levels of simplification used in the model affect its performance.

Our software code, the model, and simulation results are made publicly available (see SI) to enable further efforts in modeling in vivo activity and function.

Despite the small training set of stimuli for which our model was optimized–a gray screen and a single grating presentation for 0.5 s–it generalized well to a large test set of stimuli. The model reproduced major features of in vivo observations with respect to, e.g., the magnitude of responses to gratings, orientation selectivity, prominence of gamma oscillations, long-tailed distributions of firing rates, lifetime sparsity for gratings and movies, trends in variability of neural responses, magnitude of cortical amplification, and effect of optogenetic perturbations of the LGN or the Scnn1a population in L4 (Figs 3 and 4). Such an agreement across stimulus classes and types of observation is remarkable given that, although the model included a significant degree of biological realism, many simplifications were used: neurons possessed active conductances only at the soma, the variety of neuron models was limited to five unique single cell morphologies and sets of membrane conductances, synapses had no short-term plasticity, LGN inputs were simplified, connectivity was established via simple probabilistic rules, most known interneuron cell types [34] were absent, and, finally, and probably most importantly, the influence of most of V1 (L1, L2/3, L5, and L6) and the rest of the brain was reduced to extremely simple background states. This may suggest that many of L4 computations are produced by its local network, and other layers may play a primarily modulatory role–such as gain control exerted by L6 [61]–as far as L4 activity is concerned.

Perhaps even more instructive than successes, a number of deficiencies were observed. In terms of reproducing in vivo activity, the most important issues were the absence of direction selectivity and too fast responses to full-field flashes (Fig 3B and 3E). Both appear to be due to the extreme simplicity of the LGN inputs–in particular, the absence of sustained LGN responses in the model. This provides an immediate direction for improving the model, especially based on the more recent experimental results [23], which is the aim of ongoing work on the next model generation.

In applying the model to study mechanisms underlying the operation of the L4 circuit (Figs 4 and 5), we found the following principles. Tuning in L4 cells arose from combination of convergent LGN inputs, but to reach physiological levels of selectivity, the functional like-to-like connectivity was essential (Figs 4D, 5 and S7B). This suggests that like-to-like rules, experimentally observed in L2/3 [25, 26, 24, 29], are likely to be found in L4 as well, and may play an essential role in determining functional information processing in the cortex. The amplification of excitatory current due to recurrent connections was by 2-3-fold (Fig 4A, 4B and 4C), in quantitative agreement with experiment [17], and was linear [58, 59] along the contrast dimension, but highly non-linear along the orientation dimension (Fig 4C). Despite this strong excitatory amplification, the whole circuit was controlled by even stronger inhibition, which readily shut down the activity in the absence of the LGN input [18] (Fig 4E).

Another insight again followed from a model deficiency: the cortical component of the excitatory current (“Sub”) was orientation-tuned, but poorly modulated in time (Figs 4B and S7A), whereas in the experiment it was modulated at the grating frequency and matched the phase of LGN input [17]. The like-to-like connectivity we used enables orientation tuning of the “Sub” current, but does not take phase into account. Thus, for a given target cell, all source cells supply different phases, explaining why the “Sub” current was not modulated at the grating frequency. The fact that in the experiment it is modulated in phase with LGN current suggests more sophisticated like-to-like connectivity rules that include information about phase [17]. This is consistent with data from L2/3 [25, 24] showing that similarity in orientation preference is a good but not only predictor of connectivity and with theories (e.g., [28]) that lateral connections optimally enhance features such as extended lines.

A further simplification of our network model, replacing all biophysical neurons by LIF units, preserved most general trends in features of neuronal activity, but quantitative agreement with experiment suffered. The levels of orientation selectivity and sparsity of responses were altered (Fig 6B, 6D and 6F). The oscillations at the level of population activity were eliminated when instantaneous synaptic kinetics was used and partially rescued with non-instantaneous synapses (Fig 6C). Although further exploration is necessary to find out how one can match the oscillations spectrum of the biophysical model with a simpler network model, this result supports the importance of synaptic kinetics for generating oscillations (e.g., [11, 14]). Perhaps the most substantial difference with the biophysical model was observed in the magnitude of cortical amplification (Fig 6E), hinting that mechanisms shaping activity and computation in the L4 circuit may not be well represented by these simpler models. Overall, however, it is clear that the simplified models captured well the major aspects of network dynamics.

This suggests that in many cases more complex models may not be necessary to gain important insights about brain networks or to fit the relationship between visual stimuli and neuronal responses (just like in the stationary domain, a single hidden-layer network containing a finite number of neurons can approximate any measurable function under some mild assumptions [62, 63]). However, it is likely that neuronal point-model approximations will be insufficient to capture dendritic nonlinearities, such as synaptic saturation, veto-type inhibition, NMDA or calcium spikes and so on, that may be critical to neuronal selectivity and plasticity.

For the discrepancies that we did observe, many reasons are possible. These include, for example, dendritic filtering and distinct distributions of synapse types over dendritic trees in the biophysical, but not all-LIF models; imperfect match of synaptic kinetics (the IntFire1 model used instantaneous synapses and IntFire4 a mixture of single- and double-exponential synapses, by contrast to double-exponential synapses of the biophysical case, see Methods); and much more sophisticated somatic mechanisms in the biophysical model (mostly Hodgkin-Huxley based [31], as opposed to simple integrate-and-fire mechanisms in both IntFire1 and IntFire4), which enabled spike adaptation, after-spike hyper- and depolarization, and other nonlinear phenomena. Exploring how these different factors contribute to specific differences under systematic, data-driven constrains of our model network is an exciting direction for future studies, which will be enabled by the code and model data that we publicly share (see SI). For example, if IntFire1 and IntFire4 are not sufficient to capture all aspects of dynamics and mechanisms, would the performance improve substantially in the case where all neurons are still represented as points, but are supplied with Hodgkin-Huxley mechanisms? This could shed light, for example, on roles of synapse distributions over dendritic arbors. There are also neuron models that employ much simplified mechanisms that can account for such phenomena as spike adaptation or dendritic computations (e.g., [64, 65, 66]), and such models may furnish a better match with the biophysical simulations and experiments.

Furthermore, although we observed mismatch in certain features between the all-LIF and biophysical models, we cannot claim the all-LIF models are incapable of reproducing these features. A different optimization procedure than our relatively straightforward training on a small dataset (a short trial of grating and another of spontaneous activity, see Methods) could result in a better match. These considerations, and the fact that we did obtain mostly consistent results, support broad applicability of point-neuron models. Similar conclusions were reached in another recent study [30], even though their coarse-graining approach was different (the simpler model was optimized to match results of the complex model, whereas we optimized both the simple and complex models to match experiment) and the point neuron models were more complex than ours (the “Generalized Integrate-and-Fire” model vs. our simple Leaky Integrate-and-Fire model). Thus, while our results exhibiting, e.g., a poor match of cortical amplification, show that it is important to exercise caution in applying simplified models especially when one attempts to gain quantitative understanding of mechanisms, it is likely that point-neuron models will be sufficient to describe major qualitative features of cortical networks.

Ultimately, it is useful to apply a variety of multi-granular modeling techniques to study brain circuits–from abstract multi-layered, machine learning networks to population coding statistical models, point models, and biophysically detailed models reported here [67, 68]. As mentioned above, we here describe the biophysical model as it is the easiest to directly interpret in terms of biological variables that can be queried experimentally and that yields a concise summary of our present-date understanding of the anatomy and the physiology of L4 (and the limits of this knowledge).

Methods

Ethics statement

All experiments, animal treatment and surgical protocols were carried out with authorization from the Institutional Animal care and Use Committee of the Allen Institute in accordance with the Public Health Service (PHS) Policy on Humane Care and Use of Laboratory Animals.

Code availability

The software code used to build models, perform simulations, and perform analysis is included in this published article in its supplementary files and at https://github.com/AllenInstitute/arkhipov2018_layer4.

Model building and simulation details

All simulations were performed with the parallelized code (available in Supplementary Files) written in Python 2.7 and employing NEURON 7.4 [36] as a simulation engine. Simulations were carried out typically on 120 CPU cores (Intel Xeon CPU E5-2620 v2, 2.10GHz; 128 GB RAM per a 24-core node), requiring ~1 hour to simulate 1 s. Visual stimuli for the simulations, as well as experiments (see below), were chosen primarily from the stimulus set in the Allen Brain Observatory [69].

Network composition.

The network model was constructed from five biophysically detailed models of individual neurons [31] from an early version of the Allen Cell Types Database [32], and two models of leaky-integrate-and-fire (LIF) neurons (one for excitatory type and one for inhibitory type). The biophysical models employed passive conductances in the dendrites and 10 active conductances in the soma, including Na+, K+, and Ca2+ conductances, as described in detail in [31]. The five models represented three types of excitatory L4 neurons (expressing the genes Scnn1a, Rorb, and Nr5a1) and two types of fast-spiking parvalbumin-expressing (PV) interneurons. Cells expressing these markers comprise the majority of neurons in L4 of V1. These models consisted of 264 (Scnn1a), 141 (Rorb), 101 (Nr5a1), 121 (PV1), and 91 (PV2) compartments. The LIF neuron models were implemented using the IntFire1() function of NEURON, which contains two parameters–the time constant and the refractory period–and employs instantaneous “charge-dump” synapses. The time constant of LIF neurons was set to the average time constant from the corresponding biophysical models (from Scnn1a, Rorb, and Nr5a1 for the excitatory LIF neuron and PV1 and PV2 for the inhibitory LIF neuron); the refractory period was 3 ms for both LIF models. All the models are available in the Supplementary Files (see SI).

The 45,000 cells in the network model were distributed as follows: 3700 for Scnn1a, 3300 for Rorb, 1500 for Nr5a1, 800 for PV1, 700 for PV2, 29750 for excitatory LIF, and 5250 for inhibitory LIF. The cells were distributed in a cylinder 100 μm in height; the biophysical cells occupied the inner core with the 400 μm radius, and the LIF cells, the outer shell with the radii from 400 μm to 845 μm.

Three independent models were generated using different random seeds. These three models were used for each connectivity case (LL, LR, RL, and RR, see Main Text). For each model, the recurrent connectivity and synaptic weights (see below) differed across the LL, LR, RL, and RR cases, since different connectivity/weight rules were applied; everything else was identical across these four cases.

Connectivity within the network.

For the purposes of establishing recurrent connections, all neurons were considered to belong to two groups–excitatory (E; Scnn1a, Rorb, Nr5a1, and excitatory LIF) and inhibitory (I; PV1, PV2, and inhibitory LIF). Four types of connections were established: E-to-E, E-to-I, I-to-E, and I-to-I. The probability of connection was chosen as a product of a function dependent on the distance between the somata of the two cells and the function dependent on the difference of the assigned preferred orientation angle (see Fig 1). The latter function was constant, 1.0, for E-to-I, I-to-E, and I-to-I, whereas for E-to-E it differed depending on the model we assumed. The models with like-to-like connection probability (LL and LR, see Main text) employed a linear function, starting at 1 for the zero difference in the preferred orientation angle, and equal to 0.5 at 90 degrees difference; for the models with random connection probability (RL and RR), the function was set to one. The distance dependent function was linear, with the peak at zero distance (0.34 for E-to-E and 0.26 for E-to-I in the case of LL and LR models; 0.255 for both E-to-E and E-to-I in the case of RL and RR models; and 1.0 for I-to-E and I-to-I in all cases). The linear decay was determined by the distance at which the function became zero (300 μm for E-to-E and E-to-I, and 160 μm for I-to-E and I-to-I). For established connections, the number of synapses was selected randomly with uniform probability between 3 and 7 (see, e.g., [2, 70]).

The numbers of synapses from different sources per neuron are not well established for the mouse visual cortex, but are known approximately for cat [44], where, for the L4 excitatory cells (spiny stellate cells projecting to L2/3 or L4, and for pyramidal cells), the number of incoming synapses is ~5,000–7,000. Out of those, e.g., for pyramidal cells, ~1,100 synapses come from L4 excitatory cells, ~2,550 from excitatory cells in other layers, and ~2,000 are excitatory with unassigned source. Assuming the same distribution of sources for unassigned synapses, one can expect that ~1,700 excitatory synapses are from L4 sources. Note that the number of synapses from LGN is ~100–150 in the cat cortex [44, 71]. Converting these numbers to fractions, one finds that ~25–30% of all synapses are from the L4 excitatory sources, 55–60% from other excitatory sources, and 15–20% from inhibitory sources. For the inhibitory cells in L4, the fractions are similar, but the total number of synapses is smaller (~3,250 synapses for basket cells).

We assumed similar fractions for the mouse visual cortex, except the number of synapses from the LGN, for which we experimentally found a much larger number, ~1,000 per cell (see Main Text). Assuming somewhat smaller number of synapses for mouse cells than for cat cells, we posited (Table 1) for the purposes of our model that pyramidal cells in L4 of V1 receive ~6,000 synapses total, with ~1,000 (17%) coming from LGN, ~2,000 (~33%) from L4 excitatory sources, and ~400 (~7%) from L4 PV+ inhibitory sources. The rest was not modeled explicitly (except for the background, see Main Text and below), as it comes from other L4 inhibitory sources (PV- interneurons–presumably, another 7%, i.e., ~400 synapses) and from excitatory sources from other layers and brain regions (~2,200, i.e., 36%), both of which were not represented in the model. For the synapses from the LGN, we did not model explicitly all the LGN cell types (see below), and thus restricted the number of LGN synapses to ~600 out of the expected 1,000 total. For simplicity, we assumed numbers for synapses from LGN and L4 neurons to the PV+ cells to be similar in our model to the numbers used for pyramidal cells. Due to the randomness of connections and of the number of synapses per connection, the actual number of synapses from L4 sources for each cell varied within ~10% from the average (S1A Fig).

[Figure omitted. See PDF.]

Table 1. Assumed average numbers of synapses per an excitatory cell in L4 of mouse V1.

See text for details.

https://doi.org/10.1371/journal.pcbi.1006535.t001

Thus, each cell in our simulations received on average ~3,000 synapses from other excitatory and inhibitory cells in the L4 model and LGN sources. After building our models, including connectivity, we found that, for recurrent connections among the L4 cells, excitatory neurons received connections from 490 +/- 20 and sent connections to 480 +/- 20 neurons on average, whereas inhibitory neurons received connections from 490 +/- 20 and sent connections to 530 +/- 20 neurons on average. In addition, each cell also received a small number of “background” synapses, as described below.

Synapses were distributed on the dendritic trees of the target cells randomly within certain distance constraints, on the soma, basal, or apical processes, according to the literature (e.g., [33, 34, 44, 72, 73, 71, 74, 75, 53]). Namely, all synapses to PV1 and PV2 cells were distributed on the soma and basal dendrites without limitations. For Scnn1a, Rorb, and Nr5a1 cells, excitatory synapses from L4 and background were placed on basal and apical dendrites, 30 to 150 μm from the soma, excitatory synapses from LGN were placed on basal and apical dendrites, 0 to 150 μm from the soma, and inhibitory synapses were placed on the soma and basal and apical dendrites, 0 to 50 μm from the soma.

Model of background inputs.

Whereas the recurrent connections within L4 and the LGN inputs (see below) were modeled explicitly, the connections from the rest of the brain were represented using a simple model of “background” traveling waves (see Main Text). For producing the waves, we distributed 3,000 Poisson spike generators in the x-y space covered by the model (the L4 plane), and drew random connections from these generators to the modeled L4 cells if the distance between a generator and a cell’s soma was within 150 μm, allowing between 18 and 24 such connections per cell. This approximation of a small number of connections (much fewer than what is expected for incoming connections per L4 pyramidal cell from outside L4 or LGN) was chosen for simplicity and for reducing computing expense; in the absence of heterogeneity in the population of sources, providing a much larger number of connections does not appear necessary.

The background sources produced Poisson spike trains from a time dependent firing rate, f(x, y, t), which was controlled by the plane traveling waves (see Fig 2). The waves were pre-generated from a different random seed for each simulation trial. The direction of each wave’s movement in the x-y plane was randomly selected, its width in the direction of movement was kept constant in time, in the perpendicular direction the wave was infinite, and its magnitude was also constant and randomly selected to be between 5 and 15 Hz. In the absence of the wave, f(x, y, t) for each spike generator was set to zero, and it sharply rose to the value given by the wave’s height when the wave moved through the (x, y) location of the generator. The waves were produced so that they rarely overlapped, and thus the firing rate of the spike generators mostly exhibited periods of silence and periods of constant firing output, as waves swept through each generator’s location (see Fig 2).

The utilization of background traveling waves furnished a simple and computationally efficient analogue of different cortical states. To approximate physiological observations with the choice of the wave parameters, we performed a small set of patch-clamp recordings in pyramidal cells (n = 3, targeting L2/3 neurons in V1 for easiness of access) of anesthetized mice. In the neural responses during spontaneous activity (gray screen; S1B Fig), we observed clearly distinct intervals of rest and baseline depolarization, and found that on average the duration of depolarized states was 700+/-300 ms and the interval between the end of one such state and beginning of the next was 1,100+/-600 ms, overall consistent with observations illustrated in the literature (e.g., [47, 46]). Clearly, this is only one measurement under particular conditions, and, generally speaking, the actual characteristics of the cortical states, such as duration of the hyperpolarized state, interval between such states, and also potentially gradations in the degree of depolarization or hyperpolarization, would depend on the type of modulation one considers. These characteristics may exhibit a variety of time scales. Nevertheless, it appears that largely these cortical states may be modulated on the time scale of 1,000 ms or longer, by order of magnitude. Therefore, for simplicity, for the L4 model we generated background waves that lasted from 200 ms to 1,200 ms (i.e., 700 ms on average) and were separated by intervals of 250 ms to 1,750 ms (i.e., 1,000 ms on average)–both parameters being drawn randomly for each wave from uniform distributions. The spatial extent of each wave was set to 2,000 μm, and the propagation speed was set based on that and on the randomly selected duration of the wave.

The remaining parameters were the number of synapses per connection from a spike generator to a L4 cell and the conductance amplitude of these synapses at the synapse site. The amplitudes were set to a single number for cells of one type; the variation in the resulting postsynaptic potentials/currents (PSPs/PSCs) was due to varying synapse number per connection and random placement of synapses on the dendritic tree (see below for a general description of the choice of synaptic weights). First, the weights for the background synapses were set so that the Vm of target cells was elevated by 10–15 mV from rest during the wave condition (which was the amount of elevation we observed between the depolarized and rest states in the patch clamp experiments), in purely feedforward regime in absence of any other connections. The weights were then further adjusted to reproduce the spontaneous firing rates in the context of the full network (see below).

The number of synapses per background connection was the remaining parameter that we used to generate a long-tailed distribution of background input strengths, which we found to be important for producing a skewed (log-normal-like) distribution of firing rates (see Fig 3A), as observed experimentally [5]. Specifically, for each L4 cell, the numbers of synapses received from every background connection were the same, but differed between cells. These numbers were drawn from a long-tailed distribution heavily weighted towards 1 synapse, but permitting up to 16 synapses (S1C Fig). As a result, most L4 cells received a similar amount of background excitation, but small fractions of cells received much higher amounts of background excitation.

Synaptic characteristics.

Bi-exponential synapses (NEURON’s Exp2Syn) were used for biophysical and instantaneous synapses for LIF target cells; the reversal potential was -70 mV for inhibition and 0 mV for excitation. For simplicity, synapses did not employ short-term or other plasticity rules and had a 100% release probability.

Although presence of plasticity or probabilistic release is well documented for cortical synapses, especially in vitro, the role of these phenomena in functional properties in vivo is not well understood. Short-term plasticity may be less prominent in vivo than in vitro (due to saturation to baseline levels) [18], whereas the non-deterministic nature of individual synapses is somewhat alleviated by the fact that multiple such synapses are typically present per connection. In the current study, neglecting plasticity and probabilistic nature of synaptic release lead to significant conceptual simplifications, as well as savings in computing power. The generally good agreement of the simulated neural activity with the in vivo experimental data (see Main Text) suggests that the roles of plasticity and probabilistic release in vivo may be subtle and not straightforward to analyze. Clearly, these roles are important subjects for future investigation.

As a side observation, one should note that, because the modeled synapses were deterministic, the synaptic weights we used are probably underestimated in terms of the conductance values at the synapse location. However, the postsynaptic effects at the soma–peak PSPs and PSCs–were found to be consistent with the experiment (see below).

No complete resource exists yet with a characterization of synaptic weights and kinetics in L4 of mouse V1, and the literature on the subject of synaptic weights in general is large and somewhat disparate in the details of experimental methods used and features recorded and reported (such as post-synaptic potentials, or PSPs, vs. post-synaptic currents, or PSCs, definitions of time constants for PSP or PSC kinetics, etc.). Therefore, choices of synaptic weights and time constants in our models were guided only approximately by the existing literature on the synaptic properties of the cortex in various species (e.g., [25, 49, 51, 76, 77]), and for various cortical areas and layers (e.g., somatosensory [51] and anterior cingulate [49] cortices). Because of the lack of a consistent benchmark, we aimed at restricting the synaptic weights and time constants within approximately an order of magnitude of the values reported across a variety of experimental studies.

The kinetic parameters were fixed in the models. The Exp2Syn’s tau1 and tau2 constants were set to 1 ms and 3 ms for excitatory-to-excitatory synapses, 2.7 ms and 15 ms for inhibitory-to-excitatory synapses, 0.1 ms and 0.5 ms for excitatory-to-inhibitory synapses, and 0.2 ms and 8 ms for inhibitory-to-inhibitory synapses. These constants are related to the temporal characteristics of PSPs and PSCs in complex ways that depend on synapse distributions on the dendritic tree and electrical properties of the cell membrane. Therefore, the specific choices for the values of these constants were based on single-cell modeling runs where we searched for values that provide an approximate match with the literature about kinetics of PSPs and PSCs. Post-hoc analysis of full-network simulations confirmed (S2 Fig) that our choices of tau1 and tau2 resulted in PSP and PSC kinetics that was overall consistent with the literature (e.g., [49, 51]).

The synaptic weights were allowed to be varied globally (such as, e.g., uniformly scaling weights of all excitatory synapses on Scnn1a biophysical cells by a multiplicative constant, etc.) at the optimization stage (see Optimization of synaptic weights below). The optimization stage aimed at reproducing firing rates to a grating, spontaneous activity, and avoid epileptic-like time-locked activity of all cells. Once the weights were scaled so that these aims were approximately reached, they were kept constant throughout all simulations for a given model.

The synaptic weights were given constant values for each connection type (namely, four source types–excitatory from L4, inhibitory from L4, excitatory from LGN, and excitatory from background–to the seven target types–Scnn1a, Rorb, Nr5a1, PV1, PV2, and excitatory and inhibitory LIF), which for the biophysical cells corresponded to a constant value of the peak conductance at the synapse site. For the models where like-to-like synaptic weights were employed (LL and RL, see Main Text), the uniform synaptic weights from above were further multiplied by a factor Fw that depended on the difference between the assigned preferred orientation angle between the two cells:where Δθ is the angle difference (defined within 90°), and σθ = 50° was used in all cases.

Despite weights at the synapse location being constant for most connection types, the actual PSPs or PSCs at the soma were widely distributed, because synapses were randomly distributed on the dendritic trees. The peak somatic PSCs for E-to-E synapses (S1D Fig), for example, were observed in the model to be within 2 orders of magnitude, consistent with experimental findings [5]. Their distribution on the log scale was wide–i.e., similar to a skewed, lognormal distribution (or representable as a sum of a few lognormal distributions) described in the literature [5]. These distributions were very similar between the LL, RL, LR, and RR models, half of which employed constant E-to-E synaptic weights and the other half used Fw-modulated weights.

LGN filters.

To enable simulations with arbitrary movies as visual stimuli, we developed a filter layer, the output of which was used to drive the L4 model, representing inputs from the lateral geniculate nucleus (LGN). The linear-nonlinear Poisson (LNP; see, e.g., [40]) filters accepted movies (x,y,t-arrays) as inputs and produced time series of a firing rate as output. The output firing rates were taken to represent the rates of LGN neurons, which project to L4.

We used 3,000 filters each for three filter types–transient ON, OFF, and ON/OFF [6], i.e., representing 9,000 LGN cells. The transient types correspond to cells that produce an increase in firing rate to transitions from dark to bright (ON), bright to dark (OFF), or both (ON/OFF), such that the firing rate returns relatively quickly back to baseline levels after the change in scene. For simplicity, the sustained LGN cell types [6, 41, 42], for which firing rates may differ substantially from the baseline after the initial transient (depending on the resulting brightness level), were neglected. The major reason for neglecting these types and overall representing the LGN inputs with rather simple filters was that little is known about how different LGN cell types combine their projections onto L4 neurons, and how properties of these projections are correlated with functional properties of the target cells. Including these sustained and other less numerous, but likely important LGN types in the model is an essential objective for future work.

Each LGN filter was represented by one (for ON and OFF filters) or two (for ON/OFF filters) receptive subfields, which were described by the following equations (see Fig 1). The linear response, L(t), of a receptive subfield was given bywhere x and y are the coordinates in the visual space (angles), t is the time, S(x,y,t) is the input signal (a movie), Dt(t) is the temporal kernel, and DS(x,y) is the spatial kernel. For all operations with filters, we used a linear angle approximation for x and y (that is, replacing the tangent of an angle by the angle value in radians, which is approximately correct for small angles). The double-gaussian spatial kernel (center-surround) was used,where we used AS = AC/6,σS = 2σC. The temporal kernel had the form

The input signal was grayscale, represented on a 0 to 255 scale (from black to white), with a time step of 1 ms and the pixel size of 1.25° in x and y (the frames were 192x96 pixels). It was fed directly into ON-type receptive subfields, or converted to 255 − S(x,y,t) and then fed into OFF-type receptive subfields. The linear response L(t) was then combined with the baseline rate R0 and passed through a rectifying nonlinearity, resulting in the firing rate:

For the ON/OFF filter, R(t) from the ON and OFF receptive subfields were summed, and the mean of the corresponding R0 values was subtracted, to produce the final R(t) output.

The resulting firing rate R(t) was converted to spike trains using a Poisson random process, independently for each LGN filter (Fig 7A). The firing rate exhibited a transient at the beginning of each visual stimulus, due to transition from no visual stimulation to a movie. To avoid this artefact, we always prepended 500 ms of full-field gray screen to a movie used as a visual stimulus, and replaced the resulting first 500 ms of R(t) simply by R0.

[Figure omitted. See PDF.]

Fig 7. LGN filters.

(a) Example responses of a single filter to visual stimuli, as a time-dependent firing rate that is the filter output (blue) and the firing rate computed from generated spike trains, averaged over all trials (green). The panel for images contains responses to 10 images shown in a sequence, 250 ms each. (b) F0 and F1 components of the responses to gratings for two example LGN filters. Tuning curves to orientation, SF, and TF are shown. The data points are averages from generated spike trains over time and over trials. (c) Connecting LGN filters to L4 cells. Geometry in the visual space is illustrated. The left panel shows centers of all filters present in a portion of the visual space around the mapped position of an example excitatory cell from L4. The dashed lines correspond to the “lasso” subfields around one illustrative L4 cell, used to capture input LGN filters of the ON, OFF, and ON/OFF type. The filters that are selected to send inputs to this L4 cells are in deep color; all other filters are dimmed. On the right, the same L4 cell with the filters selected to provide inputs to it are shown. For the filters, the approximate size of their receptive subfields is illustrated (a single subfield for ON or OFF filters and two subfields for ON/OFF filters; the radius of each RF circle is 2σC). (d) Convergence of LGN connectivity onto L4 cells that are not connected to each other (“No con.”), one-way connected (“One-way con.”), and reciprocally connected (“Reciprocal con.”). For each pair of L4 cells, the LGN convergence is defined as the number of LGN filters that connect to both cells divided by the sum of the numbers of LGN filters connected to each of the cells. The data are aggregated from three L4 models.

https://doi.org/10.1371/journal.pcbi.1006535.g007

The centers of the individual filters, (x0,y0), were distributed randomly in the visual space, limited to approximately 130° x 90°. The two receptive subfields of each ON/OFF LGN cell were displaced with respect to each other (in random direction in x, y for each ON/OF cell), so as to enable orientation selectivity. Aside from that, the filters were axially symmetric, resulting in no orientation or direction selectivity. The temporal and spatial frequency (TF and SF) selectivity was thus primarily determined by the constants σC, σS, and k in the spatial and temporal kernels.

A new set of filters has been created for each of the three instantiations of the LL model type. These three instantiations of filter sets were then used for the three models within the LR, RL, and RR types.

The values of parameters for the filters were chosen randomly with a uniform probability within certain range (S1 Table). The ranges were selected to produce approximately the same preferred SF and TF (0.05 cycles per degree, or cpd, and 4 Hz, respectively) across filters, which is approximately in the middle of the typical preferred SF and TF range in real LGN cells [7]. Although real LGN cells exhibit a wide range of preferred SF and TF values, we used the above simplification in order to avoid complexities of how different information processing channels with different SF and TF tuning converge on the L4 cells, as these aspects of feedforward connectivity from LGN to V1 are largely unknown. The specific choice of parameters also allowed for relatively wide tuning over SF and TF for our LGN filters, even though the peak was almost always at SF = 0.05 cpd and TF = 4 Hz. This is generally consistent with experimental observations, where most LGN cells exhibit smooth progression of diminishing responses away from preferred SF and TF, rather than very sharply tuned responses [7].

The above requirements were used to select values of σC and k. The baseline firing rate, R0, was selected to conform approximately to the levels of spontaneous activity exhibited by LGN cells [7] (S1 Table). The remaining free parameter, Ac, was selected so that the maximal observed F0 component (see below) of the responses to gratings was 11–12 Hz (S1 Table), also approximating experimental findings [7]; +/-10% variation was allowed for this parameter.

It should be noted that for this parameterization and for the analysis described below we used the definitions of the F0 and F1 components of responses to drifting gratings following Ref. [7]. Specifically, we used the cycle-averaged firing histograms for each cell (using the cycle period of the drifting grating stimulus); F0 was the mean rate computed from the histogram and F1 was the absolute value of the Fourier component at the frequency of the drifting grating stimulus (i.e., for a sine function, a + b sin(ωt), F0 = a and F1 = |b|).

Response properties of instantiated filters were analyzed and are presented in S1 Table. Example responses of a transient ON and transient ON/OFF filters to gratings, in terms of tuning to orientation, SF, and TF, are also shown in Fig 7B (responses of transient OFF cells are very similar to those of transient ON cells, except being inverted with respect to dark vs. bright). These analyses show that indeed the filters are not selective to the orientation and direction (except for F1 responses for transient ON/OFF cells), and are tuned almost exclusively to SF = 0.05 cycles per degree (cpd) and TF = 4 Hz. However, the SF and TF tuning of individual cells (Fig 7B) is relatively broad, as intended.

Supplying inputs from LGN filters to L4 neurons.

Feedforward connections from groups of LGN filters to the L4 cells were created based on shared retinotopy. Whereas the LGN filters were defined in the visual space, the L4 cells were not, and for connecting them we introduced retinotopy in the L4 model by mapping the x, y coordinates (i.e., the plane of the layer) in the space of L4 cells to the x, y positions in the visual space. The center of the L4 model was assumed to correspond to the center of the visual space (i.e., center of the eye field). The positions were mapped using the conversion factor (from physical to visual space) of 120°/mm in x and 50°/mm in y, following the experimental observation that representation of azimuth and elevation on the surface of mouse V1 differs in length scale approximately by a factor of 2 [78]. With this mapping, each L4 cell was assigned an x,y position in the visual space.

To connect LGN cells to L4 cells, we created separate “lasso” subfields for each of three LGN types (transient ON, OFF, and ON/OFF), independently for each L4 cell (Fig 7C; see also Fig 1 in Main Text). As above, a linear angle approximation was used here. The three subfields were positioned around the visual-space x, y location of the L4 cell. The lasso subfields for inhibitory L4 target cells were bigger, more symmetric, and more overlapping than those for excitatory target cells, since experimental studies show that receptive fields of fast-spiking interneurons tend to be larger than those of excitatory cells in L4 and L2/3, and orientation selectivity of fast-spiking interneurons is also significantly lower (see, e.g., [79]). Specifically, for target L4 cells of PV1, PV2, and inhibitory LIF types, all lasso subfields were centered at the L4 cell’s position and were circular, with the diameter randomly selected from the range [15°, 20°]. For the excitatory target cells, the ON/OFF lasso subfield was centered at the L4 cell’s position and was circular, whereas centers of the ON and OFF subfields were equidistant from the center, separated by an offset distance chosen randomly from the range of [10°, 11°] (thus enabling orientation selectivity and preference for a particular SF). These subfields were ellipsoidal (of equal size for a given target cell), with the ellipse aspect ratio being selected randomly from a range of [2.8, 3.0]. The ellipse minor radius was chosen randomly from the range [3°, 4°], and the major radius was computed by multiplying the minor radius by the aspect ratio. The minor axis was oriented along the line connecting the center of the two subfields. The ON/OFF subfield radius was equal to the minor radius of these ellipses.

For the excitatory L4 target cells, the choice of the direction in which the ON and OFF lasso subfields were positioned (Fig 7C) was determined by the assigned preferred orientation angle that was set for each L4 cell at model construction. The angle between the vector connecting the centers of the ON and OFF subfields and the (1, 0) vector in the visual space was chosen to coincide with this assigned preferred orientation angle. Note that the actual preferred orientation angle was determined by the network activity in the simulations, and generally speaking did differ (more or less) from the assigned one.

Once all the lasso subfields were established for an L4 cell, the corresponding LGN cells were selected based on whether centers of their spatial kernels happened to be inside the subfields or not. Among those inside, a random set of LGN cells was selected, with the total number of each LGN type selected for one L4 target cell being capped at 15 for the inhibitory target cells and at 8 for the excitatory target cells (i.e., because we used 3 types of LGN cells, the maximum number of source LGN cells was 45 for inhibitory targets and 24 for excitatory targets). For the ON/OFF filters, an additional condition was applied during the selection of sources, which required the axis connecting the centers of the ON and OFF subfields of selected filters to be within +/-15° from the assigned preferred orientation angle of the L4 cell. In terms of connections from individual LGN filters to L4 cells, after instantiating models we found broad diversity in the number of target L4 cells per LGN filter, from 0 to several hundred; on average, one LGN filter connected to 100 +/- 170 L4 cells.

Since the LGN filters used did not represent all the LGN cell types (specifically, no sustained types [7]), we assumed that the number of synapses they provided to the L4 cells was smaller than the total expected number of synapses from the LGN (~1,000, see Main Text). We thus set the number of synapses per each LGN filter connected to an L4 cell to be 30, resulting in approximately 600 synapses from the LGN for each L4 excitatory cell.

The geometric constraints above (size and separation of the lasso fields, etc.) were chosen to reflect approximately the characteristics of the receptive fields of the LGN inputs to L4 excitatory cells, based on in vivo measurements of LGN-only currents into L4 cells, resulting from visual stimulation [17], and also taking into account typical sizes of the receptive fields of LGN cells themselves, based on extracellular electrophysiological measurements [6]. Arguably, the constraints we applied were representative of a typical case, but did not reflect considerable diversity of the receptive field sizes and shapes [6, 7]. This was done for simplicity, as we aimed to represent well the responsiveness of the LGN and L4 cells to SFs and TFs that are in the middle of typical values that evoke robust responses in vivo (SF~0.05 cpd and TF~4Hz). Introducing more diversity would require further assumptions about connectivity from LGN to L4 and within L4, in the absence of data. Another note is that the choice of the geometry of “lasso” subfields (high aspect ratios) may appear to contradict observations of the rather symmetrical and highly overlapping subfields for LGN inputs to L4 [17]. However, due to the large size of individual LGN filter receptive fields, the combined receptive subfields for the LGN inputs to L4 cells were not extremely elongated and exhibited a certain degree of overlap (see example on the right in Fig 7C), similar to the experimental observations.

Responses to moving bars and optimization of boundary conditions.

To check whether orientation tuning properties of the model depend on the nature of the stimulus, we used moving bars (S4A Fig) in addition to drifting gratings (see Main Text). Four moving bar stimuli were shown: a vertical bar moving to the right (“Ori 0 degrees”) or a horizontal bar moving up (“Ori 90 degrees”), with a black and a white bar (on a gray background) used for each of the two conditions as separate stimuli. Ten trials were performed for each stimulus. Spikes from each neuron were collected from all trials for a given stimulus and binned into an average time-dependent firing rate (with 50 ms bins). The maximum of the rate over all bins was computed for each neuron to characterize its response to bars.

We found that excitatory neurons had clear preferences to the orientations of the moving bars (S4A Fig), consistent between the white and black bars. Inhibitory neurons exhibited no preference.

We asked whether orientation preference differed between the gratings and bars. To characterize that, we selected excitatory biophysical neurons that preferred 0, 90, 180, or 270 degrees according to responses to gratings (since only horizontal and vertical movements were sampled with bars). The neurons preferring 0 or 180 degrees were combined together as preferring horizontal movement and those preferring 90 or 270 degrees–as preferring vertical movement. We then selected the bar stimulus that evoked the largest response for each of these neurons according to the maximum firing rate metric described above and assigned the neurons to two groups–preferring horizontal or vertical movement based on bar responses. Using those groupings, we computed the difference between the preferred orientation according to gratings and according to bars for each neuron. The difference was zero for every neuron (S4B Fig). Thus, in our model, excitatory neurons exhibited strong orientation preference to moving bars, which was precisely the same as orientation preference for drifting gratings.

By contrast to full-field drifting gratings, a bar is a compact moving stimulus that is expected to affect neural activity unequally across the cortical surface. This property was used at the early stages of constructing the L4 models to test boundary conditions. Models that had biophysical neurons only, employing either a circular or square geometry in the cortical plane, exhibited dramatic differences in the numbers of connections per cell between the center and the periphery (due to the sharp boundary). Using a square geometry with periodic boundary conditions for connectivity, we found that responses to moving bars were smeared across the whole cortical surface. For example, in S4C Fig (top), at t = 500 ms the bar barely reached the first receptive fields of LGN cells projecting to the biophysical L4 cells, and at t = 1000 ms just started moving through the retinotopic locations corresponding to the L4 cells, but these cells started firing soon after t = 500 ms. The firing immediately engulfed the whole model, rather than spreading over the surface following the bar’s movement. To avoid such artifacts, we switched to a model with circular geometry and with added LIF neurons at the periphery. The number of connections per neuron was uniform for the biophysical core in this model (they only dropped off for LIF neurons close to the edge), and responses to bars were more confined. In the example shown in S4C Fig (bottom), this model started responding much later than t = 500 ms, only after the bar moved to the retinotopic location corresponding to biophysical L4 cells, and the responses progressed roughly along the x axis with time. Thus, we subsequently used this geometry for all simulations.

Optimization of synaptic weights.

Once the models have been constructed, including the L4 network, LGN inputs, and background inputs, we fixed all the parameters except the synaptic weights. The weights were then optimized manually in an iterative fashion, where at each iteration a uniform scaling was applied to all weights belonging to a particular connection type. For example, weights of all excitatory-to-Rorb connections would be scaled equally at one iteration, and weights of inhibitory-to-Nr5a1 connections at another. The weights from LGN inputs were adjusted to produce experimentally observed LGN excitatory currents [17] (in the absence of recurrent connections and background). The background weights were scaled so that the firing rates of L4 cell populations in a background-only case were close to (slightly lower than) the target spontaneous firing rates. Finally, the recurrent weights were optimized in the context of a fully connected network with all inputs (the most difficult step). This part of optimization was limited to two training stimuli–a single trial of a drifting grating and one of a gray screen, 500 ms each. The targets for training were the mean spontaneous firing rates and the rates in response to the preferred grating (Rmax) of each cell population; we allowed 0.5 Hz deviation from the target for the spontaneous rate and 1 Hz deviation for Rmax. Another constraint was to avoid synchronous epileptic-like activity. After the optimization was finalized, all the model parameters were frozen, and the model was applied for further simulations of all other stimuli, which were used as a test set.

After models were optimized and simulated, we extracted a subset of synaptic connections (n = 100, each containing 3 to 7 synapses, as explained above) from each model for each connection type (E-to-E, E-to-I, I-to-E, I-to-I). We computed the PSCs and PSPs in voltage clamp and current clamp, respectively, and analyzed their peaks and time course (S2 Fig). For the current clamp measurements, current was not injected, thus allowing for characterization of synaptic properties near cell resting voltage. The voltage clamp measurements employed holding voltages of -70 mV and 0 mV to characterize excitatory and inhibitory synapses, respectively. For the time course characterization, we computed time to peak (time between the presynaptic spike and the peak), rise time (time between 20% and 80% of the peak at the rise stage), time of decay (the time constant from an exponential fit to the time course of decay, from 80% to 20% of the peak value), and the width (width at the half-peak level).

Results of this analysis indicate that the characteristics of all PSPs and PSCs are, on average, very similar between all the models (S2 Fig). The latter observation for the peak values is interesting, as the models were optimized separately based on overall characterization of network dynamics, and in principle one could expect bigger differences in synaptic weights (and, thus, PSP/PSC peaks) for the LL, LR, RL, and RR models. Furthermore, characteristics of PSPs and PSCs in the models are broadly consistent with the values reported in the literature, typically within a factor of 2 or less (e.g., [49, 51]; one should also note differences between the cortical regions and layers, different preparations and experimental conditions, as well as the naturally wide range of observed values). The trends observed in the literature were reproduced as well, such as faster dynamics of postsynaptic events in PV cells than in excitatory cells and a higher amplitude of inhibitory synaptic currents (e.g., [49]).

Stimulation protocol.

All simulations were performed with visual stimuli chosen among drifting gratings, natural movies, natural images, moving bars, and full-field flashes (S3 Fig, the first three are part of the stimulus set in the Allen Brain Observatory [69]). Separate simulations with uniformly gray screen were also carried out to characterize spontaneous activity. The list of all visual stimuli used for all model systems is presented in S2 Table. All stimuli were mapped (linearly subsampled or interpolated) to 192x96 pixels and a time step of 1 ms. Each trial was a separate simulation. Typically, 10 trials were simulated for each stimulus condition, with some exceptions specified below.

Drifting gratings were square-wave alternating white and black stripes (typically at contrast 80%, although simulations with other contrasts were performed as well), drifting in the direction perpendicular to the stripes, at 0, 45, 90, 135, 180, 225, 270, or 315 degrees. For LGN filter characterization, 240 different gratings were used, in combinations of the 8 directions listed above, 6 SF (0.025, 0.05, 0.1, 0.2, 0.4, and 0.8 cpd), and 5 TF (1, 2, 4, 8, and 15 Hz). Due to the high computational expense associated with running simulations of the L4 model, and because the LGN filters were tuned to respond preferentially to SF ~ 0.05 cpd, for the full model simulations we used a subset of these 240 conditions. A single SF = 0.05 cpd and a few TF (typically, 2, 4, and 8 Hz) were used, with all 8 directions (S2 Table). After the initial 500 ms of the gray screen, the gratings were presented for 2500 ms.

The gratings were named from g1 to g240 in the following manner: starting from the minimal TF, SF, and direction angle for g1, the numbering increased first along the TF dimension, then along the SF dimension, and finally along the direction dimension. Thus, for example, g1, g2, g3, g4, and g5 were all at direction 0 degrees, SF = 0.025 cpd, and TF varying from 1 to 15 Hz; or, gratings with a given SF and TF, such as SF = 0.05 cpd and TF = 4 Hz, had numbers that were always separated by 30 –g8, g38, g68, g98, g128, g158, g188, g218, for directions 0, 45, 90, 135, 180, 225, 270, and 315 degrees, respectively. (See Supplementary Files and S2 Table for more details.)

Natural movies were 3 clips from the opening scene of Touch of Evil [80, 69]. After the initial 500 ms of the gray screen, the movies were presented for 4500 ms. The three clips were named "TouchOfEvil_frames_N1_to_N2", where N1 and N2 were the first and last frame based on the sequence from the original movie clip (these frames were 33 ms each, as the movie was encoded at a 30 Hz rate; the clips used for simulations were upsampled to a 1 kHz rate). The N1 and N2 were 1530 and 1680 for the first clip, 3600 and 3750 for the second, and 5550 and 5700 for the third.

For natural images, we used 10 examples taken from the Berkeley Segmentation Dataset [81] and the van Hateren Natural Image Dataset [82] (see also [69]). Each trial consisted of a 3000 ms long simulation, with the first 500 ms being the gray screen, and the rest covered by presentation of the 10 images in a random sequence, for 250 ms each, without interruption. We carried out simulations of 100 such trials, each named "imseq_i", where i was 0 to 99 ("imseq" standing for “image sequence”).

Moving bars were single white or black stripes on gray background, moving with constant velocity in a direction perpendicular to the stripe. Four conditions were assessed–a black or white vertical bar moving from left to right, and a black or white horizontal bar moving from bottom to top. The bars were approximately 4 degrees wide and moved with a velocity of 33.8 degrees per second. These simulations were named "Bbar_v50pixps_hor", "Bbar_v50pixps_vert", "Wbar_v50pixps_hor", and "Wbar_v50pixps_vert".

Two types of full-field flashes were used (named "flash_1" and "flash_2"). For type 1, the screen was gray from 0 to 1000 ms, white from 1000 to 2000 ms, and gray again from 2000 to 3000 ms. For type 2, the screen was gray from 0 to 600 ms, white from 600 to 650 ms, and gray again from 650 to 1500 ms.

Spontaneous activity was measured during 1000 ms long presentation of gray screen. Typically, 20 trials were used. These simulations were named "spont".

All analysis generally disregarded the first 500 ms of each simulation, as an equilibration period (with gray screen presented as a visual stimulus during this time in all cases).

For computational efficiency, a core set of stimuli from those described above was used for all systems (three cases each of LL, LR, RL, and RR), whereas the rest were applied in a targeted manner, primarily for the LL system (S2 Table). Out of those, the richest set of stimuli was applied to model 2 of the LL type (“LL2”).

The LGN filter and background source projections to the L4 cells were generated independently for the 3 models of the LL type and then reused for the 3 models each of the LR, RL, and RR types. The same procedure was applied to recurrent connections for LL and LR cases, as well as for RL and RR. Furthermore, each trial of visual stimulation was accompanied by an independently generated trial of the background activity. These were generated independently for the 3 models of the LL type, and then applied consistently to the 3 models of the LR, RL, and RR types. Thus, for example, the LGN and background feedforward connections, as well as LGN and background spike trains incoming to the L4 cells on a particular trial of a given visual stimulus, were exactly the same between the models LL2, LR2, RL2, and RR2, but different from, e.g., LL1 or RL3.

Modeling optogenetic perturbations.

A number of simulations was performed to mimic the effects of optogenetic manipulations (see below for experimental details). The experimental optogenetic perturbation involved silencing of Scnn1a cells expressing ArchR through application of yellow light. To model such a perturbation, we injected hyperpolarizing current into somata of cells in our simulations (S8 Fig). We selected the cell type to be directly perturbed (typically Scnn1a), and a fraction of the cells in this cell type (all the cells in the cell type that were outside of this fraction did not receive current injections), and applied the current injections to those cells.

Since the light spot used in the experiments was ~2 mm in diameter, in most simulations we applied perturbation over the whole area of our model (which was ~1.7 mm in diameter), including the LIF neurons. In a subset of simulations LIF neurons were not affected (this difference did not qualitatively affect the observed trends). Delivering direct current injections to LIF neurons was not supported in the NEURON implementation we used, and thus we opted for adding inhibitory synapses to LIF neurons that were selected for optogentic perturbations. A single such synapse was added to each silenced LIF neuron, and received a Poisson spike train at 100 Hz frequency. The synaptic weight was calibrated to result in the reduction of the LIF neuron’s firing rate similar to that effected by a -100 pA current injection to the soma of biophysically detailed cells in a benchmark simulation. To represent different current injections, we linearly scaled the synaptic weights for these synapses on the LIF neurons. The LIF neurons were selected for receiving such perturbation based on the fraction of cells represented by the perturbed cell type and the fraction of cells selected for perturbation within the type. For example, Scnn1a constituted 43.5% of all excitatory biophysical neurons; if we were applying optogenetic perturbation to 20% of Scnn1a biophysical cells, then we correspondingly applied perturbation via inhibitory synapses to 8.7% of all LIF excitatory cells.

The modeled optogenetic manipulation was applied at the beginning of the simulation and was maintained throughout the simulated time. For the majority of these simulations, the injection current levels were -30, -50, -100, and -150 pA, whereas the fraction of the Scnn1a population receiving such injections were 20, 50, and 100% (see S2 Table).

Optogenetic silencing of the LGN [18] (Fig 5D) was modeled by eliminating all spikes coming from LGN filters starting at a desired time. For this study, such a time was always set to 1000 ms from the beginning of a simulation.

All-LIF network simulations.

To build an all-LIF version of the L4 model network, we replaced every biophysical cell by a LIF point-neuron model, using either IntFire1() or IntFire4() models of NEURON. The refractory period was kept at 3 ms. The IntFire1() and IntFire4() cell models both describe membrane voltage dynamics with a single-exponential process containing a single parameter, the membrane time constant, which was matched to the membrane constant of the biophysical counterpart cells. All connections (i.e., the binary adjacency matrix) were kept the same as in the biophysical model.

The synapses were instantaneous in the IntFire1() case. In the IntFire4() case, synapses were represented using a single-exponential kinetics for excitatory connections, with the time constant τE, or double-exponential kinetics for inhibitory connections, with time constants τI1 and τI2. Mapping synaptic kinetics for conductances at synapse location on a dendrite in a biophysical cell model to synaptic kinetics at a LIF cell model, with either one- or two-exponential kinetics, is not straightforward. We used a simple approach, in which we averaged the time course of somatic PSPs of each of our biophysical cell model from a large number of synapses (see S2 Fig) and approximated the resulting kinetics using either one- or two-exponential synaptic processes in the corresponding IntFire4() cells. An additional constraint was that in the IntFire4() model τI2 is required to be greater than τE. As a result of this procedure, the values of the IntFire4() synaptic time constants were set to τE = 1 ms, τI1 = 4 ms, τI2 = 17 ms for excitatory cells and τE = 7 ms, τI1 = 1 ms, τI2 = 8 ms for inhibitory cells.

Next, we optimized the synaptic weights in the all-LIF networks. A training procedure included three separate steps, each of which had the training data taken from a single trial with a single visual stimulus. First, we scaled LGN-to-L4 synapses at the level of individual cell types to match the firing rates observed in biophysical simulations for a drifting grating (TF = 4 Hz, SF = 0.05 cpd, orientation 0 degrees) in the fully feedforward regime (and without the background traveling waves). Second, the same was done for the synapses from the background sources to the L4 cells, for a single trial of spontaneous activity (gray screen). These first steps utilized results of biophysical simulations as targets because no experimental data on firing rates in LGN-only or background-only cases are available. Finally, the synapses for recurrent connections were optimized for the fully connected network. This step was done in the same way as for optimizing the biophysical networks, by scaling synaptic weights uniformly within each cell type, the target being the experimental spontaneous activity and Rmax.

Experimental methods

Experimental characterization of LGN-to-L4 synapses.

The synapse counts were performed using disectors on a systematic random sampling scheme. This allows for unbiased sampling and it has been used in other cortical areas and animal models (e.g. [83, 43]). We cannot however identify the cell type of the post synaptic target and therefore we assume that thalamic neurons target cell types according to the available dendrite of each cell. Though exceptions to this rule have been found in other areas or animal models, the contribution of L5 and L6 apical dendrites to the overall dendritic length in L4 is expected to be smaller than the dendrite from L4 neurons and one expects that it will not change substantially the proportion of thalamic synapses onto L4 in ways that would affect the predictions of the model.

We labeled thalamic boutons with immunohistochemistry for vesicular glutamate transporter 2 (VGLUT2), selected sampling sites using rare systematic random sampling and then used the physical disector method introduced by [84] to perform the synaptic counts.

Histology: Wild type mice were perfused with fixative (2% paraformaldehyde and 0.5% glutaraldehyde) and a series of 60 um coronal slices through primary visual cortex (V1) was taken, and alternating slices were processed for electron microscopy (EM) with immunohistochemistry, or processed with Nissl stain for light microscopy (LM).

Slices selected for immunohistochemistry were incubated in anti-VGLUT2 antibody. Nickel intensified diamino benzidin or silver intensified gold particles were used to visualize the labeled elements in the electron microscope (EM). Slices were then osmicated, dehydrated and embedded in Durcupan resin.

We used as a primary antibody anti-VGLUT2 from Millipore, catalog number MAB5504 (lot 2322521). The secondary antibody was biotinylated-goat-anti-mouse IgG, from vector BA-9200 (lot X0623). Anti-VGLUT2 primary antibodies label thalamic terminals in several species including the mouse visual, somatosensory and motor cortices [85, 43] and it had been validated in layer 4 of mouse visual cortex by co-labeling of axons that have been injected in the thalamus [85]. The distribution of labeled terminals in this study matches the one previously observed in the cited works.

Sampling: Slices, grids, sections and sampling sites were chosen using a systematic random sampling (SRS). A random number in the range of the number of slices was generated via computer, and used to select target coronal slices for the EM series. A general region of interest in V1 was then selected based on comparison to the cytoarchitecture of the immediately adjacent Nissl-stained sections. This region of interest, which contained layer 4, was excised from the EM slice, trimmed, and sectioned with an ultramicrotome at 40 nm thickness and serial section ribbons collected onto copper grids. Labeled terminals were then examined and photographed in a Jeol 1200 EX II electron microscope using a digital camera (Gatan).

Low magnification EM images were overlaid and used to co-register the EM sections with the LM (Nissl) images, using vasculature and cell body (soma) landmarks. The Layer 4 region was determined by cyto-architecture, then drawn on the EM image.

Physical disector: A regular sampling grid of 45 x 45 μm was generated within this layer 4 region, and high magnification images were taken at each of these sampling points for disector analysis. Reference and lookup section were separated by one intervening ultrathin section so that the disector ‘z’ dimension was 0.08 μm. Both sections were used as reference and as lookup in order to double the sample [83, 86]. Three hundred and forty one disectors were used from 3 mice. The disectors had a size of 5 × 5 μm and at 3.6 nm per pixel. Synapses and associated structures were classified using conventional criteria [87, 88]. The percentage of thalamic synapses was calculated by dividing the number of synapses formed by labeled boutons by the total number of asymmetric synapses (putative excitatory).

Extracellular electrophysiological recordings in vivo.

All electrophysiological recordings were performed in the left hemisphere of awake adult C57Bl/6 mice (2 to 6 months, males). A few weeks before the recordings, mice were first implanted with a metallic headplate using aseptic conditions and under anesthesia (see details in [7]). The day before or of recording, we performed a craniotomy over V1 (~0.5 x 0.5 mm above the monocular portion of V1, 2.5 mm lateral to lambda) and implanted a reference skull screw; if the day before, then the exposed skull was sealed with Kwik Cast. On the day of recording the exposed cortex and skull were covered with 1% agarose in saline in order to prevent drying and to help maintain mechanical stability. Dexamethasone was given to the animals to avoid brain inflammation (2 mg/kg, sc) and atropine to keep the respiratory tract clear (0.3 mg/kg, ip). Then, awake mice went into the recording setup. The headplate was clamped for stability, while the animal was free to run or remain still on a freely rotating disk. Multi-site silicon electrode arrays (imec Neuropixels) were dipped in DiI allowing post hoc visualization of the electrode path and lowered with a microdrive to a depth of 800–1000 μm. Reference and ground of the electrode were made through screws implanted in the skull. To improve the stability of recorded units, we allowed 20 minutes for the electrodes to settle.

Neurophysiological signals were amplified, band-pass filtered, and acquired continuously at 20 or 30 kHz at 16-bit resolution using an Amplipex system (Amplipex Ltd) or Open Ephys system (open-ephys.org). The spike sorting procedure was described in detail previously [89]. In brief, algorithmic identification and unit assignment of spikes [90, 91] was followed by manual adjustment of the clusters [92, 93]. Only units with clear refractory periods and well-defined cluster boundaries were included in the analyses [90].

Visual stimuli were largely the same as described above for simulations. The number of trials in the experiments was 50 for flashes, 10 for natural movies, and 100 for natural images. Rates of spontaneous activity were calculated by averaging the rates of each cell over 20 non-overlapping intervals, 500 ms long each, that were evenly distributed over 1 minute of continuous gray screen presentation (20 intervals of 500 ms were used for consistency with the 20 trials of gray screen presentation in simulations; see above). Moving bars were not used. Data for drifting gratings was taken from a previously published dataset [7].

After the recording, mice were perfused and brains fixed, sectioned, and stained with DAPI (see details of histology and imaging in [7]) for reconstruction of the electrode path to assign a layer to each recording site.

The laminar location of neurons was identified using histology and confirmed, when possible, by investigating the CSD depth profile of responses to full-field flashes. Putative excitatory neurons belonging to L4 according to this annotation were used for analysis. Since L4 is thin, and the fraction of inhibitory neurons is only ~15%, we could not obtain a sizeable number of such neurons for L4 in our recordings. Because of that, putative fast-spiking inhibitory neurons from all layers were combined for analysis.

Whole-cell patch-clamp in vivo recordings.

Data of in vivo intracellular membrane potential (Vm) was obtained from genetically defined neuron populations in cortical layer 2/3 of mouse V1 via two-photon targeted whole-cell recordings. Adult (2–4 months old, n = 2) transgenic mice expressing a red fluorescent cytosol marker tdTomato (tdT) under the control of the Cux2-CreERT2 promoter were used. Recordings were performed under Isoflurane (1–1.5% in O2) anesthesia with ~7 MΩ glass micropipettes filled with K-Gluconate based internal solution (containing, in mM: potassium gluconate 125, NaCl 10, HEPES 20, Mg-ATP 3, Na-GTP 0.4, in ddH2O; 290 mOsm; pH 7.3; Alexa 488,50 μg/ml). The exposed V1 region was covered with Artificial Cerebrospinal Fluid (ACSF) containing (in mM: NaCl 126, KCl 2.5, NaH2PO4 1.25, MgCl2 1, NaHCO3 26, glucose 10, CaCl2 2, in ddH2O; 290 mOsm; pH 7.3). tdT expressing neurons within 100 μm– 300 μm underneath the pia surface were targeted and patched, as described previously [94, 95, 96].

Two-photon imaging was performed using a Sutter Moveable Objective Microscope (Sutter, CA USA) coupled with a tunable Ti:Sapphire femtosecond laser (Chameleon Ultra II, Coherent USA), controlled by the open-source software ScanImage 3.8 (Janelia Research Campus/ Vidrio Technologies), to visualize L2/3 neurons in V1. Neurons were imaged with a 40x water-immersion objective (LUMPLFLN 40XW, Olympus USA) at the excitation wavelength of 920 nm. Electrophysiology signals were amplified with a Multiclamp 700B, digitized with a Digidata 1440B at 20kHz, using pClamp software (Molecular Devices) and stored on a PC (Dell, USA). Quality of the pipette was checked both electrically and optically. Pipette would be discarded if disproportional increase of pipette resistance (Rp) occurred, which indicates occlusion of a pipette tip. Optically checked ejection of dye under two-photon imaging was also used to indicate whether the pipette was free of occlusion. If both the electrical and optical measures agreed the tip is clogged, the penetration would be terminated and the pipette was retracted and discarded.

To conduct in vivo whole-cell recordings, a clean pipette was manually advanced towards the target neuron with a low positive pressure (~20 mbar above ambient) under two-photon imaging with Rp constantly monitored. Pipette capacitance was compensated before patching the cell. Standard procedures were used to rapidly form a gigaseal [94, 95, 96]. After gigaseal formation, a brief negative pressure was used to rupture the membrane inside the tip to achieve the whole-cell configuration. Vm was recorded under current clamp mode with series resistance (Rseries) appropriately compensated. Spontaneous Vm activity was recorded for at least 1–2 mins. The recordings were not corrected for the liquid junction potential.

Optogenetic manipulations and extracellular electrophysiological recordings in vivo.

We performed in vivo multichannel silicon probe recordings from the primary visual cortex of mice as previously described [61]. Briefly, mice were anesthetized with a combination of 5 mg/kg chlorprothixene (intraperitoneal) and 1.2 g/kg urethane (intraperitoneal); during surgery mice were supplemented with 0.5–1.0% isoflurane. A head frame with an ~2 mm central opening was mounted over V1 and the skull was thinned using a dental drill. PBS was then applied to the thinned skull and sharpened forceps were then used to remove a small portion of skull wide enough to permit insertion of a NeuroNexus 32-channel linear probe (A1x32-Edge-5mm-20-177). The probe was inserted at a depth of 800–1,000 μm. PBS was then applied to keep the craniotomy moist. Following a recovery period following probe insertion of at least 20 minutes, visual stimuli were presented using a gamma-corrected, Dell 52 3 32.5 cm LCD monitor (60 Hz refresh rate, mean luminance 50 cd/m2, 25 cm from contralateral eye). We used the Psychophysics Toolbox [97] to generate and present full-field sinusoidal drifting gratings (temporal frequency, 2 Hz; Spatial frequency, 0.04 cycles per degree; 2 orientations, 0 and 90 degrees; contrasts 0, 10, 18, 32, 64, and 100% contrast). Drifting gratings were presented for 1500 ms and optogenetic stimulation trials were interleaved with non-optogenetic trials. Each trial was separated with a 3–6 s inter-trial interval using a gray screen. To photostimulate archaerhodopsin (optogenetic trials) we used a 1 mm LED fiber with ~20 mW power output at the fiber tip, which was placed 0.5 mm from skull above recording location (590 nm, 1 mm diameter, Doric Lenses). The optogenetic light delivery lasted 1.9 s and started 200 ms prior to the visual stimulus and ended 200 ms following the stimulus. Black foil (Thor Labs) was shaped into a shield around LED and headframe in order to prevent the LED light from reaching the eyes. Recordings were amplified 1,000X and band-pass filtered between 0.3 Hz and 5 kHz using an AM systems 3500. Acquisition was done at 20 kHz with a NIDAQ PCIe-6239 board using custom MATLAB software (MathWorks).

Data were analyzed with custom-written software using MATLAB. Single units were isolated using software provided by D.N. Hill, S.B. Mehta, and D. Kleinfeld [98]. Signals were high-pass filtered at 500 Hz and waveforms were extracted from four adjacent electrode sites. Spikes were defined as events exceeding 4–5 standard deviations of the noise. Waveforms were clustered using a k-means algorithm and further aligned using a graphical user interface. Fisher linear discriminant analysis and refractory period violations were used to assess unit isolation quality. Units were assigned a depth based on the channel in which they showed the strongest signal. In this work, only the data corresponding to drifting gratings with 100% contrast was used.

The following mouse lines were used for in vivo extracellular recordings: Scnn1a-Tg3-Cre (stock number: 009613; B6;C3-Tg(Scnn1a-cre)3Aibs/J) and Ai35(RCL-Arch/GFP) (stock number: 012735; B6;129S-Gt(ROSA)26Sortm35.1(CAG-aop3/GFP)Hze/J).

Data analysis

Data analysis was performed on all simulations carried out for a given stimulus, as described in S2 Table. That typically involved at least two independent models (usually three models) and multiple trials. Data from these multiple models and trials were typically combined together to produce summary plots.

Summary figures that present box plots of various analyzed quantities adhere to the following standard. The box bottom and top represent the lower boundary of the second quartile of the data and the top boundary of the third quartile, respectively, thus marking the inter-quartile range (IQR); the median is shown in red. Whiskers mark +/-1.5 IQR from the bottom and top of the plot, as long as there are data points in that range. The rest of the data (outliers) are shown as separate dots or other symbols. In cases when the plot range does not include the whole span of the data (which was sometimes necessary due to the log-normal-like distribution of the data points over several orders of magnitude), the outliers outside of the plot range are indicated by an arrow. The numbers next to such arrows (“N1/N2”) indicate the number of such outliers (N1) and the total number of data points (N2). The mean and s.e.m. of the data are shown on these plots as circles with error bars (thick black lines).

Computation of the features of firing rate responses.

Spontaneous activity was characterized by analyzing the responses of cells from 20 gray screen trials, each 1,000 ms long. For each cell, the spikes from the second 500 ms of every trial were used to compute the firing rate, which was then averaged over trials.

The maximal response to gratings (Rmax), the orientation selectivity index (OSI) and the direction selectivity index (DSI) were computed using responses to drifting gratings (disregarding the first 500 ms of the gray screen). First, the firing rate for each cell was obtained, using spikes from the 2,500 ms of grating presentation, and averaging over all trials. This rate was computed separately for every grating condition. The Rmax metric for each cell was the maximal firing rate among all conditions.

Following [7], we used the firing rate directly to compute the OSI and DSI (i.e., without subtracting the spontaneous rate). The analysis was performed for each cell independently. For each orientation, an average response across all SF and TF was computed, and the orientation that corresponded to the greatest averaged response was assigned as the preferred orientation of the cell. The SF and TF that evoked the strongest response at the preferred orientation were then considered as preferred frequencies. OSI and DSI were computed using the responses at all orientations for these preferred SF and TF. The OSI was computed as one minus the circular variance,where k is the index for summation over all directions θk, and fk is the firing rate response for each direction. The DSI was computed aswhere fpref is the firing rate response at the preferred direction (θpref) and fnull is the response at the opposite direction (θpref + π).

Responses to flashes were characterized by averaging firing rates of multiple cells in multiple trials, in 2 ms bins. The number of trials was 10 for simulations and 50 for experiments. For simulations, the averaging was done separately for each of the two models for which 50 ms flash stimuli were simulated, and for each of the three biophysical excitatory cell populations, resulting in 6 firing rate vs. time traces. For experiments, the averaging was done separately for each mouse over all L4 excitatory cells, resulting in 8 traces. The peak magnitude and time to peak from the flash onset were then computed on each of those traces for the 1st and 2nd peaks.

Spectra of the local field potential (LFP) and of spiking activity.

The extracellular potential Φ was calculated along the central axis of the model (perpendicular to the L4 plane), at locations distributed with increments of 10 μm along the axis. The extracellular potential at the j-th location (i.e., electrode site) for a particular cell was computed as Φj = ∑kRjkIk, where Ik is the membrane current through the k-th neuronal compartment and Rjk is the transfer resistance between the k-th neuronal compartment and the j-th electrode site. The transfer resistances were computed using the line-source approximation [99] assuming homogeneous and isotropic extracellular conductivity of 0.3 S/m. The contributions of the individual cells to the recordings are then summed up to find the total recorded potential at each recording electrode site. The membrane currents were obtained with NEURON’s cvode.use_fast_imem() function which returns the transmembrane currents without needing to utilize the computationally expensive extracellular mechanism. The computation of the extracellular potential was implemented by modifying NEURON’s advance() function—called on each time step—to include the function call to our Python module computing the extracellular potential.

The local field potential (LFP) was obtained by low-pass filtering the simulated extracellular potential with the 6th-order Butterworth filter, with the cutoff frequency of 200 Hz. The power spectra of the LFP (computed employing the fast Fourier transform) were then used for analysis of oscillations in population activity. The LFP recorded at the site closest to the center of the model was utilized for this analysis.

Power spectra of spiking activity (S7 Fig) were computed as follows. A position in space was selected for recording; in all cases in this work we chose the center of the model for this purpose. The signal was then computed as the inverse-distance weighted multiunit activity, i.e., the time-binned firing rate accumulated from all neurons with the weight 1/ri, where ri is the distance from the neuron’s soma to the recording position. The time bin of 1 ms was used. The power spectrum was then computed using fast Fourier transform.

Both LFP and spiking activity spectra were averaged over all trials of a given stimulus.

Computation of sparsity.

Sparsity of responses was computed following the definition in [56]. The formula used waswhere S is the sparsity, fi is the response within time bin i (we used the firing rates averaged for each cell within each time bin over all trials for a given visual stimulus), and n is the number of bins. This defines the lifetime sparsity for each cell, for a given visual stimulus. We used time bins of 33.3 ms to match the frame rate (30 Hz) of natural movie stimuli.

For purposes of illustrating sparsity in our simulations, we computed lifetime sparsity for the three natural movies and three drifting grating stimuli (SF = 0.05 cpd, TF = 4 Hz, orientation of 0, 45, and 90 degrees). Since gratings were presented for 2.5 seconds vs. 4.5 seconds for the movies, only the first 2.5 seconds of movie responses were used for this analysis (to enable equal comparison).

Measurement of the excitatory currents in the model.

To measure the currents flowing through the somata of biophysical cells, we used NEURON’s function SEClamp(), which provides a simple voltage clamp. Out of all 10,000 biophysical cells in a model, every 200th cell was voltage-clamped at the soma at -70 mV, to measure the excitatory current [17]. Thus, the currents were recorded for 50 cells in such a simulation (e.g., cell ID 2, 202, 402, 602, 802, …, 9802). The currents were saved from the SEClamp() objects, at each simulation time step.

For the measured current, the F0 (mean) and F1 (amplitude at the frequency of the stimulus) components were obtained according to Ref. [17]. The current was averaged over trials and then cycle-averaged. F0 was computed as the mean (over the cycle time) of the cycle-averaged current. F0 was then subtracted from the cycle-averaged current, and the result was fit with a sinusoid function at the stimulus frequency, b sin(ωt + φ), where ω = 2πν (ν being the stimulus frequency). The F1 component value was taken as F1 = 2|b|.

For all-LIF models, we measured the contribution of LGN synaptic inputs to the total excitatory synaptic inputs as a proxy to the excitatory synaptic currents measured for the biophysical model. Specifically, for each spike arriving to a synapse, that synapse contributed its weight to the accumulating sum for the target cell; the sum was computed for synapses from LGN only or for all synapses and was normalized by time and number of trials.

Computation of the coefficient of variation of inter-spike intervals, Fano factor, and noise and signal correlations.

To study the variability of spike counts, we report coefficient of variation, Fano factor, signal and noise correlation of spike counts for both simulation and experiment. Each of these factors were computed separately for gratings, natural movies and spontaneous activity.

The coefficient of variation was defined as the squared ratio of standard deviation to the mean of the inter-spike intervals. The data points from all trials were included.

To compute the Fano factor, we found the variance and mean of the spike counts in the 50 ms sliding window moving in 10ms steps. The Fano factor was defined as the slope of the linear regression line relating the variances to the means of spike counts in these sliding windows in all trials [57].

The signal correlation was computed as the Pearson correlation coefficient between the trial-averaged spike counts for each pair of neurons. For gratings, we computed the correlation over spike counts in 8 different orientations. For natural movies and spontaneous activity, the correlation was computed for binned spike counts in non-overlapping windows, each 333 ms long (corresponding to 10 frames of a natural movie). The noise correlation was computed as the Pearson correlation coefficient between single-trial spike counts for each pair of neurons, and then averaged over stimuli (8 orientations for gratings and non-overlapping 333 ms windows for natural movies and spontaneous activity).

Supporting information

S1 Fig. Model construction.

(a) Distribution of the number of synapses from L4 sources onto the five cell types in the network model. Three curves are shown for each type, corresponding to three independent model instantiations. (b) An example trace of the membrane voltage Vm from patch clamp recordings from V1 L2/3 cells in vivo (under anesthesia) during spontaneous activity. Blue portion of the trace corresponds to the data points removed from the analysis of the state characteristics (clipped spikes). The gray stripes mark the identified rest states. (c) Number of synapses formed on an L4 cell by a connection from a single background source. The distribution of synapse numbers over all background connections is shown, combined from the three independent models. (d) Distribution of peak somatic EPSC values for E-to-E recurrent connections in all simulated models (for biophysical neurons), on the log scale. Data for three models per system type (LL, LR, RL, or RR) are shown (see Main Text).

https://doi.org/10.1371/journal.pcbi.1006535.s001

(TIF)

S2 Fig. Features of PSPs and PSCs for the recurrent connections.

See Online Methods for details of box plots. The features are voltage or current peak (“peak”), time from spike to peak (“t_to_peak”), rise time (“t_rise”), decay time (“t_decay”) and the PSP or PSC width (“t_width”). (a) Somatic PSP features. (b) Somatic PSC features. For each feature and each model type (i.e., LL, LR, RL, or RR), the sample sizes are n = 900 for “E-to-E” and “I-to-E” and n = 600 for “E-to-I” and “I-to-I”.

https://doi.org/10.1371/journal.pcbi.1006535.s002

(TIF)

S3 Fig. Visual stimuli.

Examples of visual stimuli used for simulations and experiments are shown, such as (a) drifting gratings, (b) natural movies, (c) static natural images, (d) moving white or black bars, and (e) full-field flashes.

https://doi.org/10.1371/journal.pcbi.1006535.s003

(TIF)

S4 Fig. Responses to moving bars.

(a) Responses of each biophysical neuron in one model to black and white bars; either a vertical bar was moving in a horizontal direction (“Ori 0 degrees”) or a horizontal bar was moving in a vertical direction (“Ori 90 degrees”). The responses shown were obtained from time-dependent firing rates (in 50 ms bins) averaged over all trials of a given stimulus; the maximum over all bins is computed for each neuron. The neuron IDs for each type are arranged according to the neurons’ assumed direction preference for gratings (see Online Methods), from 0 degrees for the first ID of a type to 360 degrees for the last (hence the pseudo-periodicity apparent in the plots). The types are Scnn1a (IDs 0 to 3699), Rorb (3700 to 6999), Nr5a1 (7000 to 8499), PV1 (8500 to 9299), and PV2 (9300 to 9999). (b) The difference ΔOri between the preferred orientations of a neuron according to responses to gratings and to bars, averaged over all excitatory neurons that prefer 0, 90, 180, or 270 degrees for gratings. The averages and standard deviations are exactly zero for all three models tested. (c) Spike rasters (left) for biophysical neurons from pilot simulations of responses to a horizontally moving white bar, using different model layouts illustrated on the right. For each spike, the position of the neuron along the x dimension (which coincides with the direction of the moving bar) is plotted versus spike time. Top, a model without LIF neurons, with biophysical neurons confined to a rectangular area, and using periodic boundary conditions for connectivity. Bottom, a model with biophysical neurons confined to a cylinder, with LIF neurons distributed at the periphery (no periodic boundary conditions)–that is, the model layout chosen for all simulations reported in the Main Text. The approximate extent of the receptive fields (RFs) of LGN cells that feed into the biophysical portion of the model are marked by white dashed lines. Note that in these preliminary test simulations, the parameters of the moving bar (its width and speed) were somewhat different from those chosen later for production simulations.

https://doi.org/10.1371/journal.pcbi.1006535.s004

(TIF)

S5 Fig. Additional characteristics of visual responses.

(a) Distributions of skewness of firing rates. Left, simulation; right, electrophysiological experimental recordings. (b) The PSTHs from experimental electrophysiological recordings in response to a 50 ms flash (average over all L4 excitatory cells or all inhibitory cells recorded, and all trials, in 2 ms bins). (c) Example tuning curves of a single Scnn1a or PV1 cell to drifting gratings at contrasts C = 80% and C = 10%. For the Scnn1a cell, responses normalized to the peak of the tuning curve are also shown (middle). The data are averages over 10 trials. Error bars: standard deviation. Dashed lines: spontaneous rate (it is close to zero for the example Scnn1a cell shown). (d) Summary of responses to the gratings at different contrasts (C = 30% or 10% vs. C = 80%). The distributions of differences ΔHWHH = HWHH(C = 80%)—HWHH(C = 30%) (top) and ΔHWHH = HWHH(C = 80%)—HWHH(C = 10%) (bottom), are shown for all excitatory cells, with the average +/- standard deviation indicated. (e) Same as (d) for the differences of OSI for excitatory (red) and inhibitory (blue) neurons.

https://doi.org/10.1371/journal.pcbi.1006535.s005

(TIF)

S6 Fig.

Comparison of variability and correlations between simulations (left) and experiment (right). Results of the analysis are shown for gratings (magenta), natural movies (green) and spontaneous activity (beige). (a) Coefficient of variation of inter-spike intervals. (b) Fano factor. (c) Signal correlations. (d) Noise correlations.

https://doi.org/10.1371/journal.pcbi.1006535.s006

(TIF)

S7 Fig. Features of the responses of the L4 circuit.

(a) Temporal dynamics of current amplification (for a TF = 2 Hz grating). The total and LGN-only currents are shown, along with their difference (“Sub”, i.e., the cortical component); the top plot shows the inhibitory current. Data for individual cells for the preferred direction are averaged over 10 trials; the traces are then shifted in time (periodically with respect to the trial start and end, i.e., 500 ms and 3000 ms) so that the phases of the LGN component are aligned across all cells. The data are then averaged over all biophysical excitatory cells. Additional tests were performed, where LGN and background inputs to the voltage-clamped cells were eliminated (gray). In this case, the measured current is from the L4 circuit only. This demonstrates that in simulations the oscillation of the “Sub” component (black) is due to a space clamp artifact. (b) Comparison of the maximal responses to drifting gratings (Rmax) and OSI between the full network simulations (“Full”) and the purely feedforward simulations with the inputs to L4 coming only from the LGN filters (“LGN only”). Note that the OSIs of PV1 and PV2 cells are nominally high for the feedforward case because the firing rates of these cells in that situation are very close to zero. An occasional rare spike results in a “strong” response in comparison with zero responses for most trials of most orientations, which leads to elevated OSI values. (c) Comparison of the firing rates in the full network vs. firing rates in a purely feedforward model receiving only inputs from the LGN. Top, an individual Rorb cell (each point is an average over time and over 10 trials). Linear fits are shown for data aggregated from all grating directions, TFs, and contrasts (black), for one selected direction (yellow), and for a fixed contrast and TF (i.e., representing a sample direction tuning curve; right plot). Bottom, summary of linear fits across all cells analyzed.

https://doi.org/10.1371/journal.pcbi.1006535.s007

(TIF)

S8 Fig. Optogenetic suppression of the Scnn1a cells in the L4 circuit.

(a) Summary of experimental manipulations. The OMI index data are shown for each cell in all 6 experiments (colors, m1 to m6) in the top plot. The bottom plot shows the same data normalized to the total number of cells from all 6 experiments and the combined normalized distribution (gray). Note that in the bottom plot, the curves for m1, m3, m5, and m6 fully or partially overlap. (b) Spike rasters from two simulations–with and without optogenetic perturbation of the Scnn1a population. In both cases, the same drifting grating is presented (TF = 4 Hz), with the stimulation illustrated at the bottom. The hyperpolarizing current at the level of -100 pA, representing the optogenetic perturbation, was injected in the somata of 100% of the Scnn1a cells. (c) OMI computed across all recorded layer 4 excitatory neurons in experiments and simulations. The experimental curve (gray; it is the same in all plots and same as in (a)) is used as a benchmark. The simulation conditions differ between the plots (in the amount of current injected and the proportion of cells that receive the injected current), and thus the simulation curves (black) are different.

https://doi.org/10.1371/journal.pcbi.1006535.s008

(TIF)

S1 Table. Parameters and properties of LGN filters.

The ranges of parameters from which values are randomly selected for each instantiated filter are shown. Properties of the filter responses to visual stimuli are indicated for mean and standard deviation based on spike train responses from 10 trials for gratings and 20 trials for all other stimuli. Orientation selectivity index (OSI) and direction selectivity index (DSI) are computed separately for the F0 and F1 components of the responses to gratings.

https://doi.org/10.1371/journal.pcbi.1006535.s009

(DOCX)

S2 Table. Simulations of the L4 models.

The simulated systems (of the LL, RL, LR, and RR types) are indicated together with the path name in the file repository. The five types of visual stimuli, as well as spontaneous activity (“Spont.”) are listed in columns. If simulations of a particular type were performed for a given system, they are indicated together with the corresponding file name suffix and the number of trials.

https://doi.org/10.1371/journal.pcbi.1006535.s010

(DOCX)

SI 1. A zip archive containing the python/NEURON code for running simulations, along with a simple example of a model and a simulation using that model (see instructions in the README.txt file within the archive).

https://doi.org/10.1371/journal.pcbi.1006535.s011

(ZIP)

SI 2. A zip archive containing software code for building L4 network models (see instructions in the README.txt file within the archive).

https://doi.org/10.1371/journal.pcbi.1006535.s012

(ZIP)

SI 3. A zip archive containing software code for distributing LGN filter models in the visual space and producing LGN spike trains in response to movies, as well as for instantiating background inputs (see instructions in the README.txt file within the archive).

https://doi.org/10.1371/journal.pcbi.1006535.s013

(ZIP)

SI 4. A zip archive containing software code for generating visual stimulus movies–gratings, moving bars, and full field flashes, as well as ready natural movies and natural images (see instructions in the README.txt file within the archive).

https://doi.org/10.1371/journal.pcbi.1006535.s014

(ZIP)

SI 5. A zip archive containing software code for analysis of simulation results (see instructions in the README.txt file within the archive).

https://doi.org/10.1371/journal.pcbi.1006535.s015

(ZIP)

Acknowledgments

We are grateful to Gabe J. Murphy for helpful discussions. We thank the Allen Institute founders, Paul G. Allen and Jody Allen, for their vision, encouragement and support.

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

References

1. California Institute of Technology. Caltech Library Archives [Internet]. Available from: http://archives-dc.library.caltech.edu/islandora/object/ct1:483

2. Markram H. et al. Reconstruction and Simulation of Neocortical Microcircuitry. Cell 163, 456–492 (2015). pmid:26451489

3. Traub R.D. et al. Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts. J. Neurophysiol. 93, 2194–2232 (2005). pmid:15525801

4. Zhu W, Shelley M, and Shapley R. A neuronal network model of primary visual cortex explains spatial frequency selectivity. J. Comput. Neurosci. 26, 271–287 (2009). pmid:18668360

5. Buzsaki G., and Mizuseki K. The log-dynamic brain: how skewed distributions affect network operations. Nat Rev Neurosci. 15, 264–278 (2014). pmid:24569488

6. Niell CM, Stryker MP. Highly selective receptive fields in mouse visual cortex. J Neurosci 28, 7520–7536 (2008). pmid:18650330

7. Durand S. et al. A comparison of visual response properties in the lateral geniculate nucleus and primary visual cortex of awake and anesthetized mice. J Neurosci 36, 12144–12156 (2016). pmid:27903724

8. Li Y. et al. Broadening of Inhibitory Tuning Underlies Contrast-Dependent Sharpening of Orientation Selectivity in Mouse Visual Cortex. J. Neurosci. 32, 16466–16477 (2012). pmid:23152629

9. Finn I.M., Priebe N.J., and Ferster D. The emergence of contrast-invariant orientation tuning in simple cells of cat visual cortex. Neuron 54, 137–152 (2007). pmid:17408583

10. Welle CG and Contreras D. Sensory-driven and spontaneous gamma oscillations engage distinct cortical circuitry. J. Neurophysiol. 115, 1821–1835 (2016). pmid:26719085

11. Brunel N., and Wang X.-J. What Determines the Frequency of Fast Network Oscillations with Irregular Neural Discharges? J. Neurophysiol. 90, 415–430 (2003). pmid:12611969

12. Borgers C., and Kopell N. Synchronization in Networks of Excitatory and Inhibitory Neurons with Sparse, Random Connectivity. Neural Computation 15, 509–538 (2003). pmid:12620157

13. Tiesinga P, and Sejnowski TJ. Cortical enlightenment: are attentional gamma oscillations driven by ING or PING? Neuron 63, 727–732 (2009). pmid:19778503

14. Buzsaki G., and Wang X.-J. Mechanisms of gamma oscillations. Annu Rev Neurosci. 35, 203–225 (2012). pmid:22443509

15. Bos H., Diesmann M., and Helias M. Identifying Anatomical Origins of Coexisting Oscillations in the Cortical Microcircuit. PLoS Comput Biol. 12, e1005132 (2016). pmid:27736873

16. Froudarakis E. et al. Population code in mouse V1 facilitates readout of natural scenes through increased sparseness. Nat. Neurosci. 17, 851–857 (2014). pmid:24747577

17. Lien AD, Scanziani M. Tuned thalamic excitation is amplified by visual cortical circuits. Nat Neurosci 16, 1315–1323 (2013). pmid:23933748

18. Reinhold K., Lien A.D., and Scanziani M. Distinct recurrent versus afferent dynamics in cortical visual processing. Nature Neurosci. 18, 1789–1797 (2015). pmid:26502263

19. Kondo S, and Ohki K. Laminar differences in the orientation selectivity of geniculate afferents in mouse primary visual cortex. Nat Neurosci. 19, 316–319 (2016). pmid:26691830

20. Sun W., et al. Thalamus provides layer 4 of primary visual cortex with orientation- and direction-tuned inputs. Nature Neurosci. 19, 308–315, (2016). pmid:26691829

21. Li Y. et al. Linear Transformation of Thalamocortical input by Intracortical Excitation. Nat Neurosci. 16, 1324–1330 (2013). pmid:23933750

22. Morgenstern N.A., Bourg J., and Petreanu L. Multilaminar networks of cortical neurons integrate common inputs from sensory thalamus. Nature Neurosci. 19, 1034–1040 (2016). pmid:27376765

23. Lien A.D., and Scanziani M. Cortical direction selectivity emerges at convergence of thalamic synapses. Nature 558, 80–86 (2018). pmid:29795349

24. Ko H. et al. Functional specificity of local synaptic connections in neocortical networks. Nature 473, 87–91 (2011). pmid:21478872

25. Cossell L. et al. Functional organization of excitatory synaptic strength in primary visual cortex. Nature 518, 399–403 (2015). pmid:25652823

26. Lee W.C. et al. Anatomy and function of an excitatory network in the visual cortex. Nature 532, 370–374 (2016). pmid:27018655

27. Yoshimura Y., Dantzker J.L.M., and Callaway E.M. Excitatory cortical neurons form fine-scale functional networks. Nature 433, 868–873 (2005). pmid:15729343

28. Iyer R. and Mihalas S. Cortical circuits implement optimal context integration. bioRxiv, 158360. (2017).

29. Wertz A., Trenholm S., Yonehara K., Hillier D., Raics Z., Leinweber M., Szalay G., Ghanem A., Keller G., Rózsa B., Conzelmann K.-K., and Roska B. Single-cell–initiated monosynaptic tracing reveals layer-specific cortical network modules. Science 349, 70–74 (2015). pmid:26138975

30. Rössert C. et al. Automated point-neuron simplification of data-driven microcircuit models. https://arxiv.org/abs/1604.00087

31. Gouwens N. W. et al. Systematic generation of biophysically detailed models for diverse cortical neuron types. Nature Communications 9, 710 (2018). pmid:29459718

32. Allen Institute for Brain Science. Allen Cell Types Database [Internet]. Available from: http://celltypes.brain-map.org.

33. Braitenberg V. and Schuz A. Cortex: Statistics and Geometry of Neuronal Connectivity. Berlin: Springer (1998).

34. Tremblay R., Lee S., and Rudy B. GABAergic Interneurons in the Neocortex: From Cellular Properties to Circuits, Neuron 91, 260–292 (2016). pmid:27477017

35. Gratiy S.L. et al. BioNet: a Python interface to NEURON for modeling large-scale networks. PLoS ONE 13: e0201630 (2018). Package available from: https://alleninstitute.github.io/bmtk/bionet.html. pmid:30071069

36. Hines ML, Carnevale NT. The NEURON simulation environment. Neural Comput. 9, 1179–1209 (1997). pmid:9248061

37. Perin R, Berger TK, Markram H. A synaptic organizing principle for cortical neuronal groups. Proc Natl Acad Sci USA 108, 5419–5424 (2011). pmid:21383177

38. Troyer T.W., Krukowski A., Priebe N.J., Miller K.D. Contrast-invariant orientation tuning in cat visual cortex: Feedforward tuning and correlation-based intracortical connectivity. J. Neurosci. 18, 5908–5927 (1998). pmid:9671678

39. Krukowski A.E. and Miller K.D. Thalamocortical NMDA Conductances and Intracortical Inhibition Can Explain Cortical Temporal Tuning. Nature Neurosci., 4, 424–430 (2001). pmid:11276234

40. Simoncelli E. P., Paninski L., Pillow J. & Swartz O. Characterization of Neural Responses with Stochastic Stimuli. The Cognitive Neurosciences 3rd edn., 327–338 (Ed. Gazzaniga M.) MIT press (2004).

41. Grubb M.S. and Thompson I.D. Quantitative characterization of visual response properties in the mouse dorsal lateral geniculate nucleus. J. Neurophysiol. 90, 3594–3607 (2003). pmid:12944530

42. Piscopo D.M. et al. Diverse Visual Features Encoded in Mouse Lateral Geniculate Nucleus. J. Neurosci. 33, 4642–4656 (2013). pmid:23486939

43. Bopp R., et al. An Ultrastructural Study of the Thalamic Input to Layer 4 of Primary Motor and Primary Somatosensory Cortex in the Mouse. J. Neurosci. 37, 2435–2448 (2017). pmid:28137974

44. Binzegger T., Douglas R.J., and Martin K.A. A quantitative map of the circuit of cat primary visual cortex. J. Neurosci. 29, 8441–53 (2004).

45. Sato TK, Nauhaus I, Carandini M. Traveling waves in visual cortex. Neuron 75, 218–229 (2012). pmid:22841308

46. Tan AY, et al. Sensory stimulation shifts visual cortex from synchronous to asynchronous states. Nature 509, 226–229 (2014). pmid:24695217

47. McGinley MJ, et al. Waking State: Rapid Variations Modulate Neural and Behavioral Responses. Neuron 87, 1143–1161 (2015). pmid:26402600

48. Harris K.D., and Thiele A. Cortical state and attention. Nature Reviews Neurosci. 12, 509–523 (2011).

49. Wu L-J, Li X, Chen T, Ren M, Zhuo M, Characterization of intracortical synaptic connections in the mouse anterior cingulate cortex using dual patch clamp recording. Mol Brain 2, 32 (2009). pmid:19828050

50. Prinz A. A., Bucher D., and Marder E. Similar network activity from disparate circuit parameters. Nature Neurosci. 7, 1345–1352 (2004). pmid:15558066

51. Lefort S., Tomm C., Floyd Sarria J.-C., and Petersen C.C.H. The Excitatory Neuronal Network of the C2 Barrel Column in Mouse Primary Somatosensory Cortex. Neuron 61, 301–316 (2009). pmid:19186171

52. Roxin A. et al. On the Distribution of Firing Rates in Networks of Cortical Neurons. J. Neurosci. 31, 16217–16226 (2011). pmid:22072673

53. Banitt Y., Martin K.A., and Segev I. A biologically realistic model of contrast invariant orientation tuning by thalamocortical synaptic depression. J. Neurosci. 27, 10230–10239 (2007). pmid:17881529

54. Tucker T.R. and Fitzpatrick D. Luminance-Evoked Inhibition in Primary Visual Cortex: A Transient Veto of Simultaneous and Ongoing Response. J. Neurosci. 26, 13537–13547 (2006). pmid:17192437

55. Macknik S. L. and Martinez-Conde S. The spatial and temporal effects of lateral inhibitory networks and their relevance to the visibility of spatiotemporal edges. Neurocomputing 58–60, 775–782 (2004).

56. Vinje WE and Gallant JL. Natural stimulation of the nonclassical receptive field increases information transmission efficiency in V1. J. Neurosci. 22, 2904–2915 (2002). pmid:11923455

57. Churchland M.M. et al. Stimulus onset quenches neural variability: a widespread cortical phenomenon. Nat. Neurosci. 13, 369–378 (2010). pmid:20173745

58. Douglas R.J., et al. Recurrent excitation in neocortical circuits. Science 269, 981–985 (1995). pmid:7638624

59. van Vreeswijk C, and Sompolinsky H. Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science 274, 1724–1726 (1996). pmid:8939866

60. Schaub M.T. et al. Emergence of Slow Switching Assemblies in Structured Neuronal Networks. PLOS Comp. Bio. 11, e1004196 (2015).

61. Olsen S.R., Bortone D.S., Adesnik H., and Scanziani M. Gain control by layer six in cortical circuits of vision. Nature 483, 47–52 (2012). pmid:22367547

62. Cybenko G. Approximations by superpositions of sigmoidal functions. Mathematics of Control, Signals, and Systems 2, 303–314 (1989).

63. Hornik K., Stinchcombe M., and White H. Multilayer feedforward networks are universal approximators. Neural Networks 2, 359–366 (1989).

64. Teeter C. et al. Generalized leaky integrate-and-fire models classify multiple neuron types. Nature Communications 9, 709 (2018). pmid:29459723

65. Shai A. S. et al. Physiology of Layer 5 Pyramidal Neurons in Mouse Primary Visual Cortex: Coincidence Detection through Bursting. PLoS Comput Biol 11, e1004090 (2015). pmid:25768881

66. Poirazi P., Brannon T., and Mel B.W. Pyramidal Neuron as Two-Layer Neural Network. Neuron 37, 989–999 (2003). pmid:12670427

67. Hawrylycz M. et al. Inferring cortical function in the mouse visual system through large-scale systems neuroscience. PNAS 113, 7337–7344 (2016). pmid:27382147

68. Cain N. et al. The Computational Properties of a Simplified Cortical Column Model. PLOS Comp. Bio. 12, e1005045 (2016).

69. Allen Institute for Brain Science. Allen Brain Observatory–Visual Coding [Internet]. Available from: http://observatory.brain-map.org/visualcoding/.

70. Silver R.A. et al. High-Probability Uniquantal Transmission at Excitatory Synapses in Barrel Cortex. Science 302, 1981–1984 (2003). pmid:14671309

71. da Costa N. and Martin K.A.C. How Thalamus Connects to Spiny Stellate Cells in the Cat’s Visual Cortex. J. Neurosci. 31, 2925–2937 (2011). pmid:21414914

72. Richardson R.J., Blundon J.A., Bayazitov I.T., and Zakharenko S.S. Connectivity Patterns Revealed by Mapping of Active Inputs on Dendrites of Thalamorecipient Neurons in the Auditory Cortex. J. Neurosci. 29, 6406–6417 (2009). pmid:19458212

73. Schoonover C.E. Comparative Strength and Dendritic Organization of Thalamocortical and Corticocortical Synapses onto Excitatory Layer 4 Neurons. J. Neurosci. 34, 6746–6758 (2014). pmid:24828630

74. Douglas R.J., and Martin K.A.C. Inhibition in cortical circuits. Current Biology 19, R398–R402 (2009). pmid:19467204

75. Spruston N. Pyramidal neurons: dendritic structure and synaptic integration. Nat. Rev. Neurosci. 9, 206–221 (2008). pmid:18270515

76. Pfeffer C.K., et al. Inhibition of inhibition in visual cortex: the logic of connections between molecularly distinct interneurons. Nature Neurosci. 16, 1068–1076 (2013). pmid:23817549

77. Beierlein M., Gibson J.R., and Connors B.W. Two Dynamically Distinct Inhibitory Networks in Layer 4 of the Neocortex. J. Neurophys. 90, 2987–3000 (2003).

78. Niell C.M., Bonin V., and Andermann M.L. Functional Organization of Circuits in Rodent Primary Visual Cortex. In The New Visual Neurosciences, Werner J.S.and Chalupa L.M., Eds. MIT Press (2013).

79. Liu B. et al. Visual Receptive Field Structure of Cortical Inhibitory Neurons Revealed by Two-Photon Imaging Guided Recording. J. Neurosci. 29, 10520–10532 (2009). pmid:19710305

80. Welles O. Touch of Evil (United States: Universal—International). (1958).

81. Martin D., Fowlkes C., Tal D., and Malik J. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. Proc. Eighth IEEE Int. Conf. Comput. Vision. ICCV 2001 2, 416–423 (2001).

82. van Hateren J.H., and van der Schaaf A. Independent component filters of natural images compared with simple cells in primary visual cortex. Proc. Biol. Sci. 265, 359–366 (1998). pmid:9523437

83. da Costa N. M., Hepp K., & Martin K. A. C. A systematic random sampling scheme optimized to detect the proportion of rare synapses in the neuropil. Journal of Neuroscience Methods, 180, 77–81. (2009). pmid:19427532

84. Sterio D. The unbiased estimation of number and sizes of arbitrary particles using the disector. J Microsc, 134 (Pt 2), 127–136 (1984). pmid:6737468

85. Coleman J.E., et al. Rapid Structural Remodeling of Thalamocortical Synapses Parallels Experience-Dependent Functional Plasticity in Mouse Primary Visual Cortex. J. Neurosci. 30, 9670–9682 (2010). pmid:20660250

86. Howard, C., and Reed, M. Unbiased stereology (2nd ed.). New York. (2005).

87. Gray E. G. Electron microscopy of synaptic contacts on dendrite spines of the cerebral cortex. Nature, 183, 1592–1593 (1959). pmid:13666826

88. Colonnier M. Synaptic patterns on different cell types in the different laminae of the cat visual cortex. An electron microscope study. Brain Research, 9, 268–287 (1968). pmid:4175993

89. Mizuseki K, Sirota A, Pastalkova E, Buzsaki G. Theta oscillations provide temporal windows for local circuit computation in the entorhinal-hippocampal loop. Neuron 64, 267–280 (2009). pmid:19874793

90. Harris KD, Henze DA, Csicsvari J, Hirase H, Buzsaki G. Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements. J Neurophysiol 84, 401–414 (2000). pmid:10899214

91. Pachitariu M, Steinmetz NA, Kadir S, Carandini M and Harris KD. Kilosort: realtime spike-sorting for extracellular electrophysiology with hundreds of channels. bioRxiv (2016).

92. Hazan L, Zugaro M, Buzsaki G. Klusters, NeuroScope, NDManager: a free software suite for neurophysiological data processing and visualization. J Neurosci Methods 155, 207–216 (2006). pmid:16580733

93. Rossant C. et al. Spike sorting for large, dense electrode arrays. Nat. Neurosci. 19, 634–641 (2016). pmid:26974951

94. Kitamura K, Judkewitz B, Kano M, Denk W, Häusser M. Targeted patch-clamp recordings and single-cell electroporation of unlabeled neurons in vivo. Nat Methods 5, 61–67 (2008). pmid:18157136

95. Judkewitz B, Rizzi M, Kitamura K, Häusser M. Targeted single-cell electroporation of mammalian neurons in vivo. Nat Protoc. 4, 862–869 (2009). pmid:19444243

96. Long B, Li L, Knoblich U, Zeng H, Peng H. 3D Image-Guided Automatic Pipette Positioning for Single Cell Experiments in vivo. Sci Rep. 5, 18426 (2015). pmid:26689553

97. Brainard D.H. The Psychophysics Toolbox. Spat Vis. 10, 433–436 (1997). pmid:9176952

98. Fee M.S., Mitra P.P., Kleinfeld D. Automatic sorting of multiple unit neuronal signals in the presence of anisotropic and non-Gaussian variability. J. Neurosci. Meth., 69, 175–188 (1996).

99. Holt G.R., and Koch C. Electrical interactions via the extracellular potential near cell bodies. J Comput Neurosci. 6, 169–184 (1999). pmid:10333161