Uncovering circuit mechanisms of current sinks and sources with biophysical simulations of primary visual cortex
Abstract
Local field potential (LFP) recordings reflect the dynamics of the current source density (CSD) in brain tissue. The synaptic, cellular, and circuit contributions to current sinks and sources are ill-understood. We investigated these in mouse primary visual cortex using public Neuropixels recordings and a detailed circuit model based on simulating the Hodgkin–Huxley dynamics of >50,000 neurons belonging to 17 cell types. The model simultaneously captured spiking and CSD responses and demonstrated a two-way dissociation: firing rates are altered with minor effects on the CSD pattern by adjusting synaptic weights, and CSD is altered with minor effects on firing rates by adjusting synaptic placement on the dendrites. We describe how thalamocortical inputs and recurrent connections sculpt specific sinks and sources early in the visual response, whereas cortical feedback crucially alters them in later stages. These results establish quantitative links between macroscopic brain measurements (LFP/CSD) and microscopic biophysics-based understanding of neuron dynamics and show that CSD analysis provides powerful constraints for modeling beyond those from considering spikes.
Editor's evaluation
The study demonstrates that utilizing the LFP and/or the CSD in modeling can facilitate model configuration and implementation by revealing discrepancies between models and experiments. The analysis of the biophysical origin of the canonical CSD using the model is an interesting and worthy line of investigation. The dissection of CSD components is detailed and exhaustive. A key novelty of this article is the addition of CSD patterns as another constraint to more accurately infer the model parameters beyond its prior state.
https://doi.org/10.7554/eLife.87169.sa0Introduction
The local field potential (LFP) is the low-frequency component (below a few hundred Hertz) of the extracellular potential recorded in brain tissue that originates from transmembrane currents in the vicinity of the recording electrode (Lindén et al., 2011; Buzsáki et al., 2012; Einevoll et al., 2013; Pesaran et al., 2018; Sinha and Narayanan, 2022). While the high-frequency component of the extracellular potential, the single- or multi-unit activity (MUA), primarily reflects action potentials of one or more nearby neurons, the LFP predominantly stems from currents caused by synaptic inputs (Mitzdorf, 1985; Einevoll et al., 2007) and their associated return currents through the membranes. Thus, cortical LFPs represent aspects of neural activity that are complementary to those reflected in spikes, and as such, they can provide additional information about the underlying circuit dynamics from extracellular recordings.
Applications of LFP are diverse and include investigations of sensory processing (Rall and Shepherd, 1968; Di et al., 1990; Victor et al., 1994; Kandel and Buzsáki, 1997; Mehta et al., 2000a; Mehta et al., 2000b; Henrie and Shapley, 2005; Einevoll et al., 2007; Belitski et al., 2008; Montemurro et al., 2008; Niell and Stryker, 2008; Nauhaus et al., 2009; Bastos et al., 2015; Senzai et al., 2019), motor planning (Scherberger et al., 2005; Roux et al., 2006), navigation (Tort et al., 2008; Makarova et al., 2011; Fernández-Ruiz et al., 2012; Watrous et al., 2013; Fernández-Ruiz et al., 2017), and higher cognitive processes (Pesaran et al., 2002; Womelsdorf et al., 2006; Liu and Newsome, 2006; Kreiman et al., 2006; Liebe et al., 2012). The LFP is also a promising signal for steering neuroprosthetic devices (Mehring et al., 2003; Andersen et al., 2004; Rickert et al., 2005; Markowitz et al., 2011; Stavisky et al., 2015) and for monitoring neural activity in human recordings (Mukamel and Fried, 2012) because the LFP is more easily and stably recorded in chronic settings than spikes. Due to the vast number of neurons and multiple neural processes contributing to the LFP, however, it can be challenging to interpret (Buzsáki et al., 2012; Einevoll et al., 2013; Hagen et al., 2016). While we have extensive phenomenological understanding of the LFP, less is known about how different cell and synapse types and connection patterns contribute to the LFP or how these contributions are sculpted by different information processing streams (e.g., feedforward vs. feedback) and brain states.
One way to improve its interpretability is to calculate the current source density (CSD) from the LFP, which is a more localized measure of activity, and easier to read in terms of the underlying neural processes. The current sinks and sources indicate where positive ions flow into and out of cells, respectively, and are constrained by Kirchoff’s current law (i.e., currents sum to zero over the total membrane area of a neuron). However, the interpretation of current sinks and sources is inherently ambiguous as several processes can be the origin of a current sink or source (Buzsáki, 2006; Pettersen et al., 2006; Einevoll et al., 2007). For example, a current source may reflect an inhibitory synaptic current or an outflowing return current resulting from excitatory synaptic input elsewhere on the neuron. There is no simple way of knowing which it is from an extracellular recording alone (Buzsáki, 2006).
Another approach to uncovering the biophysical origins of current sinks and sources, and by extension the LFP, is to simulate them computationally (Pettersen et al., 2008; Einevoll et al., 2013). Following the classic work by Rall in the 1960s (Rall, 1962), a forward-modeling scheme in which extracellular potentials are calculated from neuron models with detailed morphologies using volume conduction theory under the line source approximation has been established (Holt and Koch, 1999). With this framework, we have achieved a good understanding of the biophysical origins of extracellular action potentials (Koch, 1998; Holt and Koch, 1999; Pettersen and Einevoll, 2008; Hay et al., 2011; Lindén et al., 2010). Expanding on this understanding, models composed of populations of unconnected neurons (e.g., Pettersen et al., 2008; Lindén et al., 2011; Schomburg et al., 2012; Łęski et al., 2013; Sinha and Narayanan, 2015; Hagen et al., 2017; Ness et al., 2018) and recurrent network models (e.g., Traub et al., 2005; Vierling-Claassen et al., 2010; Reimann et al., 2013; Głąbska et al., 2014; Tomsett et al., 2015; Hagen et al., 2016; Hagen et al., 2018; Chatzikalymniou and Skinner, 2018) have been used to study the neural processes underlying LFP.
While interesting insights about CSD and LFP were obtained from these computational approaches, establishing a direct relationship between the biological details of the circuit structure and the electrical signal like LFP remains a major unresolved challenge. One reason is that the amount and quality of data available for modeling the circuit architecture in detail have been limited. This situation improved substantially in recent years, and a broad range of data on the composition, connectivity, and physiology of cortical circuits have been integrated systematically (Billeh et al., 2020) in a biophysically detailed model of mouse primary visual cortex (area V1). In addition, significant improvements were achieved in experimental recordings of the LFP and the simultaneous spiking responses. In particular, the Neuropixels probes (Jun et al., 2017) record LFP and hundreds of units across the cortical depth in multiple areas, with 20 μm spacing between recording channels allowing for an unprecedented level of spatial detail. These developments provide unique opportunities to improve our understanding of circuit mechanisms that determine LFP patterns.
Here, we analyze spikes and LFP from the publicly available Allen Institute’s Visual Coding survey recorded using Neuropixels probes (https://www.brain-map.org; Siegle et al., 2021) and reproduce these using the mouse V1 model developed by Billeh et al., 2020. The model is comprised of >50,000 biophysically detailed neuron models surrounded by an annulus of almost 180,000 generalized leaky-integrate-and-fire units. The neuron models belong to 17 different cell type classes: one inhibitory class (Htr3a) in layer 1, and four classes in each of the other layers (2/3, 4, 5, and 6) where one is excitatory and three are inhibitory (Pvalb, Sst, Htr3a) in each layer. The visual coding dataset consists of simultaneous recordings from six Neuropixels 1.0 probes across a range of cortical and subcortical structures in 58 mice while they are exposed to a range of visual stimuli (about 100,000 units and 2 billion spikes over 2 hr of recording).
In our analysis of this dataset, we identified a canonical CSD pattern that captures the evoked response in mouse V1 to a full-field flash. We then modified the biophysically detailed model of mouse V1 to reproduce both the canonical CSD pattern and laminar population firing rates in V1 simultaneously. We reproduce, in a quantitative manner, the shape and timing of the pattern of current sources and sinks that have been described in considerable detail by experimentalists (e.g., Mitzdorf, 1987; Swadlow et al., 2002; Senzai et al., 2019). This shows that adjustments to synaptic parameters such as weights and placement in addition to a circuit architecture that included feedback are sufficient to reproduce experimental findings on both single-cell measures such as spikes and population-level measures such as CSD. We use this model to explain, in a highly mechanistic manner, the biophysical origins of the various ionic current sinks and sources and their location across the various layers of visual cortex.
In the process of obtaining a model that could reproduce both spikes and CSD, we discovered that the model can be modified by adjusting the synaptic weights to reproduce the experimental firing rates with only minor effects on the simulated CSD, and, conversely, that the simulated CSD can be altered with only minor effects on the firing rates by adjusting synaptic placement. Furthermore, we found that comparing the simulated CSD to the experimental CSD revealed discrepancies between model and data that were not apparent from only comparing the firing rates. Additionally, it was not until feedback from higher cortical visual areas (HVAs) was added to the model that simulations reproduced both the experimentally recorded CSD and firing rates, as opposed to only the firing rates. This bio-realistic modeling approach sheds light on specific components of the V1 circuit that contribute to the generation of the major sinks and sources of the CSD in response to abrupt visual stimulation. Our findings demonstrate that utilizing the LFP and/or the CSD in modeling can aid model configuration and implementation by revealing discrepancies between models and experiments and provide additional constraints on model parameters beyond those offered by the spiking activity. The new model obtained here is freely accessible (https://doi.org/10.5061/dryad.k3j9kd5b8) to the community to facilitate further applications of biologically detailed modeling.
Results
Spikes and LFP were recorded across multiple brain areas, with a focus on six cortical (V1, LM, AL, RL, AM, PM) and two thalamic (LGN, LP) visual areas, using Neuropixels probes in 58 mice (Siegle et al., 2021).
A schematic of the six probes used to perform the recordings in individual mice is shown in Figure 1A, and the spikes and LFP recorded in V1 of an exemplar mouse during presentation of a full-field bright flash stimulus are displayed in Figure 1B, C. The CSD can be estimated from the LFP (averaged over 75 trials) using the delta iCSD method to obtain a more localized measure of inflowing (sinks) and outflowing currents (sources) (Pettersen et al., 2006; Einevoll et al., 2013). The biophysically detailed model of mouse V1 used to simulate the neural activity and the recorded potential in response to the full-field flash stimulus is illustrated in Figure 1E. The model contains 230,924 neurons, of which 51,978 are biophysically detailed multicompartment neurons with somatic Hodgkin–Huxley conductances and passive dendrites, and 178,946 are leaky-integrate-and-fire (LIF) neurons. These neuron models are arranged in a cylinder with a radius 845 μm and a height 860 μm. The multicompartment neurons are placed in the ‘core’ with a radius of 400 μm, while the LIF neurons form an annulus surrounding this core. Cellular models belong to 17 different classes: one excitatory class and three inhibitory (Pvalb, Sst, Htr3a) in each of layers 2/3, 4, 5 and 6, and a single Htr3a inhibitory class in layer 1. The extracellular electric field in the model was recorded on an array of simulated point electrodes (Dai et al., 2020) arranged in a straight line (Figure 1D) and separated by 20 μm, consistent with Neuropixels probes, shown in Figure 1E to scale with the model.

Illustration of experimental data and the biophysical model for mouse primary visual cortex (V1).
(A) Schematic of the experimental setup, with six Neuropixels probes inserted into six cortical (V1, latero-medial [LM], rostro-lateral [RL], antero-lateral [AL], postero-medial [PM], AM) and two thalamic areas (LGN, LP). (B) Top: spikes from many simultaneously recorded neurons in V1 during a single trial. Bottom: spikes from a single neuron recorded across multiple trials. In both cases, the stimulus was a full-field bright flash (onset at time 0, offset at 250 ms). (C) Top: local field potential (LFP) across all layers of V1 in response to the full-field bright flash, averaged over 75 trials in a single animal. Bottom: current source density (CSD) computed from the LFP with the delta iCSD method. (D) Histology displaying trace of the Neuropixels probe across layers in V1, subiculum (SUB) and dentate gyrus (DG). (E) Visualization of the V1 model with the Neuropixels probe in situ. (Image made using VND.)
Uncovering a canonical visually evoked CSD response
We first established a ‘typical’ experimentally recorded CSD pattern to be reproduced with the model. Though there is substantial inter-trial and inter-animal variability in the evoked CSD response, we find that most trials and animals have several salient features in common. In Figure 2A, the trial-averaged evoked CSDs from five individual mice are displayed. In the first four animals (#1–4), we observe an early transient sink arising in layer 4 (L4) ~40 ms after flash onset, followed by a sustained source starting ~60 ms, which covers L4 and parts of layers 2/3 (L2/3) and layer 5 (L5). We also observe a sustained sink covering layers 5 and 6 (L6) emerging around 50 ms, as well as a sustained sink covering layers 1 and 2/3 around 60 ms. An animal that does not fully exhibit what we term the ‘canonical’ pattern is shown in the rightmost plot (#5 in Figure 2A); it has an early L4 sink arising at 40 ms, but this sink is not followed by the sustained sinks and sources from 50 to 60 ms and onward observed in the other animals. The timing and location of sinks and sources are, overall, similar to those described earlier by Givre et al., 1994; Schroeder et al., 1998, Niell and Stryker, 2008, and Senzai et al., 2019.

Variability in experimentally recorded current source density (CSD).
(A) Evoked CSD response to a full-field flash averaged over 75 trials, from five animals in the dataset. (B) The first principal component (PC) computed from the CSD of all n = 44 animals, explaining 50.4% of the variance. (C) Illustration of movement of sinks and sources in the calculation of the Wasserstein distance (WD) between the CSD of two animals in the dataset. The gray lines in the rightmost panels display how the sinks or sources of one animal are moved to match the distribution of sinks or sources of the other animal. (D) Left: WDs from each animal to the PC 1 CSD. Right: pairwise WDs between all 44 animals sorted by their distance to the first PC. (E) CSD from five individual trials in example animal 1. (F) Distribution of pairwise distances between single-trial CSD (red) and pairwise distances between trial-averaged CSD of individual animals (blue). Both are normalized to the maximum pairwise distance between the trial-averaged CSD of individual animals. (G) Pairwise WDs between trials in each of 44 animals (white boxplots), normalized to maximal pairwise WDs between trial-averaged CSD of animals. Gray-colored boxplot shows the distribution of pairwise WDs between trial-averaged CSD of individual animals, and the red stars indicate the n = 5 animals for which the inter-trial variability was greater than the inter-animal variability (assessed with Kolmogorov–Smirnov [KS] tests, p < 0.001 in all cases, see Figure 2—figure supplement 3).
To identify the robust features across animals in this dataset, we performed principal component analysis (PCA) on the trial-averaged evoked CSD from all animals. Out of the 58 animals in the dataset, 5 did not have readable recordings of LFP in V1 during the presentation of the full-field flash stimuli, and the exact probe locations in V1 could not be recovered for 9 other animals due to fading of fluorescent dye or artifacts in the optical projection tomography (OPT) volume (see ‘Materials and methods’). The remaining 44 (out of the 58) animals in the dataset were retained for the CSD analysis. The trial-averaged CSD plots of all these 44 animals are displayed in Figure 2—figure supplement 1. The first principal component (PC 1) (Figure 2B) constitutes a sum of weighted contributions of the CSD patterns from all 44 animals and explains half (50.4%) of the variance. The salient features typically observed in individual animals are also prominent in the PC 1 CSD pattern (Figure 2B), that is, the canonical pattern. In Figure 2—figure supplement 2, the first 10 principal components cumulatively explaining 90% of the variance are plotted.
Quantifying CSD pattern similarity
We use the Wasserstein, or Earth Mover’s, distance (WD) to quantify the differences in CSD patterns (see ‘Materials and methods’), which can then be used to assess how well the simulated CSD matches the CSD typically observed in experiments. The WD reflects the cost of transforming one distribution into another by moving its ‘distribution mass’ around (Rubner et al., 1998; Arjovsky et al., 2017). An often-used analogy refers to the two distributions as two piles of dirt, where the WD tells us the minimal amount of work that must be done to move the mass of one pile around until its distribution matches the other pile (Rubner et al., 1998). In the context of CSD patterns, the WD reflects the cost of transforming the distribution of sinks and sources in one CSD pattern into the distribution of sinks and sources in another pattern, with larger WD indicating greater dissimilarity between CSD patterns. The WDs are computed between the sinks of two CSD patterns and between the sources of two CSD patterns independently, and then summed to form a total WD between the CSD patterns (Figure 2C). The sum of all sinks and the sum of all sources in each CSD pattern are normalized to –1 and +1, respectively, so the WD only reflects differences in patterns, and not differences in the overall amplitude. The WD scales linearly with shifts in space and time.
When computing the WDs between the evoked CSD patterns of individual animals and the canonical pattern, we find that the animals with CSD patterns that, by visual inspection, resemble the canonical pattern (Figure 2A, animals 1–4), are indeed among animals with smaller WD, while the animal with the more distinct CSD pattern (Figure 2A, animal 5) is an outlier (Figure 2D).
The onset of the evoked response is less conspicuous in the single-trial CSD due to pronounced, ongoing sinks and sources, but there is still a visible increase in magnitude from 40 to 50 ms onward (Figure 2E), compatible with the latency of spiking responses to full-field flashes in V1 (Siegle et al., 2021). An oscillation of sinks and sources with a periodicity of ~20 ms, that is, in the gamma range is apparent in the region stretching from L2/3 to the top of L5, which appears to be either partially interrupted or drowned out by more sustained sinks and sources emerging at about 60 ms. At least some of this gamma-range activity derives from the visual flash that covers the entire visual field and that drives retinal neurons and postsynaptic targets in the lateral geniculate nucleus (LGN) in an oscillatory manner (see the pronounced gamma-range oscillation in the LGN firing rate in Figure 3D).

Variability in experimentally recorded spikes.
(A) Trial-averaged laminar population firing rates of regular-spiking (RS) cells, differentiated by layer, and fast-spiking (FS) cells across all layers in response to full-field flash. Black line: average across all animals. Gray shaded area: ±1 standard deviation. (B) Kolmogorov–Smirnov (KS) similarities (see ‘Materials and methods’) between the trial-averaged firing rates of each individual animal and the average firing rate over cells from all animals (black line in A) at baseline (the interval of 250 ms before flash onset), peak evoked response (from 35 to 60 ms after flash onset), and during the sustained period (from 60 to 100 ms). (C) Correlations between trial-averaged firing rates of individual mice and all mice (0–100 ms after flash onset). (D) Baseline-subtracted evoked firing rates for excitatory cells in seven visual areas (average over trials, neurons, and mice). Note the strong, stimulus-triggered gamma-range oscillations in the firing of lateral geniculate nucleus (LGN) neurons (blue). (E) Mean (μ) ± standard deviation (σ) of population firing rates during baseline, peak evoked response, and the sustained period. Averaged across trials, neurons and time windows defined above.
The inter-trial variability is roughly comparable to the inter-animal variability of the trial-averaged responses. By computing the pairwise Wasserstein distances between single trial CSDs within each animal and comparing it to the pairwise WD between the trial-averaged CSD of each animal, we find that inter-trial variability in CSD is significantly lower than the inter-animal variability in trial-averaged CSD (Kolmogorov–Smirnov distance = 0.33; p<0.001) (Figure 2F).
The majority of animals (39 out of 44) have a WD to the first principal component, PC 1, of the CSD that is less than half of the greatest WD between the CSD of individual animals and the PC 1 CSD (Figure 2D); the pairwise WDs between animals are also less than half of the maximum pairwise WD for most animals (921 out of the total 946 pairwise WDs; Figure 2E). This supports the view that most animals exhibit the canonical CSD pattern captured by the PC 1 CSD (Figure 2B). The total inter-trial variability is smaller than the inter-animal variability, both estimated by pairwise WDs (Figure 2F and G), though there are n = 5 animals for which the inter-trial WDs are larger than the inter-animal WDs (Figure 2G, marked by red stars; determined with KS tests on the distribution of pairwise WDs between animals and pairwise WDs between trials in each animal; see Figure 2—figure supplement 3).
Quantifying firing rate variability
For the spike analysis, we distinguish between fast-spiking (FS; putative Pvalb inhibitory) neurons and regular-spiking (RS; putative excitatory and non-Pvalb inhibitory) neurons (see ‘Materials and methods’ and Figure 3—figure supplement 1). All FS-neurons are grouped together into one population across all layers, while the RS-neurons are divided into separate populations for each layer (Figure 3A). The FS-neurons are merged across layers because we set a criterion of at least 10 recorded neurons in any one layer when comparing the population firing rate in individual animals to the average population firing rate in all animals, and only one animal had 10 FS-neurons or more in any layer (Figure 3—figure supplement 2). This criterion was set to have a more reliable estimate of the population firing rates in individual animals.
We use the KS similarity (defined as one minus the KS distance, see ‘Materials and methods’) and correlation to quantify the variability in firing rates. We use the experimental variability as a reference to assess whether the model reproduces firing rates typically observed in experiments. The KS similarity gives the similarity between the distributions of average firing rates across neurons in two populations in selected time windows, with KS similarity = 1 implying identity. As such, KS similarity provides a metric to compare the magnitudes of firing rates in certain time periods. We defined the ‘baseline’ window as the period over 250 ms before the flash onset, the ‘initial peak’ window as 35–60 ms after flash onset, and the ‘sustained’ window as 60–100 ms after flash onset. The KS similarity score during baseline is denoted ‘KSSb,’ during the ‘initial peak’ ‘KSSp,’ and ‘sustained’ ‘KSSs.’ The correlation, on the other hand, is computed between two population firing rates throughout the 100 ms window. The correlation thus gives us a measure of the similarity in the temporal profile of firing rates in this interval, independent of magnitudes. We establish the experimental variability in KS similarities and correlation by computing these metrics between the population firing rates of each individual animal and the average population firing rates of all other animals (averaged over trials for both the individual animals and the average over all other animals) (Figure 3B and C).
The population firing rates for FS neurons are more than twice as high as RS cells during baseline, peak, and sustained. Among the RS populations, the firing rate in L5 is the highest in all periods, followed by L4 and L6, while L2/3 has the lowest firing rates (Figure 3E).
Discrepancy between the original model and experimental observations
We simulated the response to a full-field flash stimulus with the biophysical network model of mouse primary visual cortex as presented in Billeh et al., 2020. As input to the model, we used experimentally recorded LGN spike trains (Figure 4C; see ‘Materials and methods’). A Poisson source, firing at a constant rate of 1 kHz, provides additional synaptic input to all cells, representing the influence from the rest of the brain (‘background’ input). The thalamocortical input consists of spike trains from 17,400 LGN units (Arkhipov et al., 2018; Billeh et al., 2020). The public Neuropixels data contain recordings from 1263 regular-spiking LGN neurons across 32 mice during 75 trials of full-field bright flash presentations, resulting in 94,725 spike trains. To construct the input for each of our 10 simulation trials, we randomly sampled 10 unique subsets of spike trains from this pool until all 17,400 units had been assigned a spike train in each trial.

Local field potential (LFP), current source density (CSD), and spikes from simulations with the original model.
(A) Top: raster plot of all ~50,000 cells in the model’s 400 μm radius ‘core’ region spanning all layers, in a simulation of a single trial with the flash stimulus. Bottom: raster plot and histogram of spikes from 10 trials for an example cell. (B) Top: simulated LFP averaged over 10 trials of flash stimulus. Bottom: CSD calculated from the LFP via the delta iCSD method. (C) Firing rate of experimentally recorded lateral geniculate nucleus (LGN) spike trains used as input to the model. (D) Wasserstein distance between CSD from the original model (blue diamond) and PC 1 CSD from experiments together with the Wasserstein distances from experimental CSD in every animal to PC 1 CSD (boxplot), normalized to maximal distance for animals. (E) Experimentally recorded firing rates (black) and simulated firing rates (blue). (F) Kolmogorov–Smirnov (KS) similarity between firing rates in original model (blue diamond) or individual animals (boxplots) and firing rates in experiments at baseline, peak evoked response, and during the sustained period (defined in Figure 3). (G) Correlation between firing rates of model (blue diamond) or individual animals in experiments (boxplots) and average population firing rates in experiments (0–100 ms). (H) Mean (μ) ± standard deviation (σ) of model firing rates during baseline, peak evoked response, and the sustained period. Averaged across trials, neurons and time windows defined above.
Figure 4A, B displays the resulting spiking pattern across all layers with its associated LFP. The inferred CSD exhibits a strong sink in the L5 and L6 region, matched by a strong source below it, both starting at ~50 ms after flash onset (Figure 4B, bottom). However, the early L4 sink, the later sustained L4 source, and the sustained L2/3 sink typically observed in the experimental CSD (Figure 2A, B) are either absent or too weak compared to the sink and source in L5 and L6. The WD from the simulated CSD to the experimental PC 1 CSD is greater than the WD between the CSD of the farthest outlier animal and the PC 1 CSD (WD = 1.84, normalized to the largest WD between CSD of individual animals and PC 1 CSD). Thus, using experimental variability as a reference, the CSD from this simulation is an outlier (Figure 4C).
The population firing rates of the model, the KS similarities and correlation between the model and the data, are plotted together with the data in Figure 4D–F. The magnitudes of the model firing rates are higher than the experimental firing rates in all populations and time windows (Figure 4H). However, the KS similarities between the model firing rates and the experimental firing rates are still within the minimum to maximum range of the boxplots for the RS L2/3, RS L4, and RS L5 cells in all time windows (Figure 4F), and during baseline for the FS cells. For RS L6 neurons, the KS similarities were among the outliers of the experiments in all time windows, while for FS neurons they were among the outliers during the peak and sustained windows (RS L2/3: KSSb = 0.62, KSSp = 0.63, and KSSs = 0.54; RS L4: KSSb = 0.77, KSSp = 0.60, and KSSs = 0.63; RS L5 KSSb = 0.77, KSSp = 0.77, and KSSs = 0.78; RS L6: KSSb = 0.54, KSSp = 0.45, and KSSs = 0.47; FS: KSSb = 0.54, KSSp = 0.53, and KSSs = 0.49). The temporal profile of the model firing rates is above the minimum of the boxplots for all populations (RS L2/3: r = 0.38***, RS L4: r = 0.62***, RS L5: r = 0.75***, RS L6: r = 0.90***, FS: r = 0.80***; ***p<0.001). For RS L4 and RS L6, the model is in fact outside the experimental distribution, but it is an outlier in a positive sense; the model firing rates for these populations are more similar to the experimental average than the corresponding population firing rates in individual animals.
The original model studied in Figure 4 produced firing rates and orientation and direction tuning consistent with recordings in vivo (Billeh et al., 2020) with some shortcomings, such as relatively slow responses of V1 to the onset of visual stimuli (Arkhipov et al., 2018; Billeh et al., 2020). Here, we see even more inconsistencies reflected clearly in the CSD pattern. This demonstrates the importance of multi-modal characterization of such biologically detailed models. To investigate the properties of the cortical circuit that sculpt the CSD, we manipulated the model and observed how both the CSD and firing rate responses were improved to match the experimental data.
Adjusting the model to fit experimental firing rates
Due to the discrepancy between the magnitudes of the model firing rates and the experimental firing rates, especially with respect to the outliers of the modeled RS L6 and FS neurons, we selectively adjusted the recurrent synaptic weights. We left the synaptic weights between LGN and the V1 model unchanged since they were well constrained by data (Billeh et al., 2020).
We first reduced the synaptic weights from all excitatory populations to the FS PV-neurons by 30% to bring their firing rates closer to the average firing rate in this population in the experiments. This resulted in increased firing rates in all other (RS) populations due to the reduced activity of the inhibitory Pvalb-neurons (Figure 4—figure supplement 1). Therefore, we further applied reductions in the synaptic weights from all excitatory neurons to RS neurons and increases in the synaptic weights from inhibitory neurons to the RS neurons to bring their firing rates closer to the experimental average firing rates. We multiplied the recurrent synaptic weights with factors in the [0.2, 2.5] range until we arrived at a set of weights where none of the model firing rates were among the experimental outliers in any time window (KSSb = 0.73, KSSp = 0.77, and KSSs = 0.70; average across RS populations and the FS population) and temporal profiles (RS L2/3: r = 0.49***, RS L4: r = 0.63***, RS L5: r = 0.71***, RS L6: r = 0.87***, FS: r = 0.86***; *** p<0.001) (Figure 5A–C).

Adjusting the model to fit spikes or current source density (CSD).
(A) Average experimentally (black) and simulated firing rates of experiments in the model with adjusted recurrent synaptic weights (green) and original model (blue). Synaptic adjustments included scaling the weights from all excitatory populations to the PV cells down by 30% to reduce the firing rates in these fast-spiking populations, reducing the synaptic weights from excitatory populations to all others and increasing synaptic weights from all PV cells to all other populations to compensate for the reduced inhibition. (B) Kolmogorov–Smirnov (KS) similarity between firing rates of model versions (markers) or individual animals in experiments (boxplots) and firing rates of experiments at baseline, peak evoked response, and during the sustained (defined in Figure 3). (C) Correlation between simulated firing rates or individual animals (boxplots) and measured firing rates (0–100 ms). (D) Left: PC 1 current source density (CSD) from experiments (see Figure 2). Right: CSD resulting from simulation on model with adjusted recurrent synaptic weights. (E) Wasserstein distance between CSD from model versions and PC 1 CSD from experiments together with Wasserstein distances from CSD in animals to PC 1 CSD (boxplot). (F) Effect of different patterns of placing excitatory synapses onto layer 4 excitatory cells on this population’s contribution to the simulated CSD (left) and to the total simulated CSD (right). These synaptic placement schemes with accompanying inflowing (blue arrows) and outflowing (orange arrows) currents are illustrated in the middle. (G) Effect of synaptic placement on the simulated population firing rate. (H) Contribution of L4 excitatory cells to the simulated CSD in the model where all recurrent connections have been cut (left) and when all active channels have been removed from all cells in the model (right).
The resulting pattern (but not the magnitude) of the CSD, however, was largely unchanged (Figure 5D) compared to the original CSD (Figure 4B). The overall magnitude was reduced, and there were some traces of a sink arising at 40 ms after flash onset, and a L2/3 (and L1) sink after 60 ms, but they were substantially weaker relative to the L5/L6 dipole than they were in the experiments. Furthermore, the large and sustained L4 source after 60 ms was still either absent or too weak to be visible. The WD between the CSD from this version of the model and the experimental PC 1 CSD remained among the outliers of the animals (Figure 5E) (normalized WD = 1.26).
Two-way dissociation between spikes and CSD
Simulations demonstrate that the LFP, and the associated CSD, can be significantly altered by changes to synaptic placement (Einevoll et al., 2007; Pettersen and Einevoll, 2008; Lindén et al., 2010; Lindén et al., 2011; Łęski et al., 2013; Hagen et al., 2017; Ness et al., 2018). As observed in Figure 5A–E and Figure 5—figure supplement 1, adjustments to synaptic weights can modify the population firing rates substantially, yet without substantially changing the pattern of the CSD, that is, the placement and timing of sinks and sources. The inverse can also occur; that is, the CSD pattern can be altered extensively with only minor effects on firing rates (Figure 5F and G, Figure 5—figure supplements 2–4).
In the model’s original network configuration, L4 excitatory neurons received geniculate input from synapses placed within 150 μm from the soma on both basal and apical dendrites, and excitatory, recurrent input from other V1 neurons within 200 μm from the soma on both basal and apical dendrites. We tested the effects of synaptic location by placing all synapses from both LGN and excitatory neurons onto the basal dendrites of L4 excitatory neurons (within the same ranges as in the original configuration). This increased the contribution from the L4 excitatory neurons to the total CSD (Figure 5F, middle row, leftmost plot) by a factor of ~2 and led to a dipole pattern with a single sink at the bottom and a single source at the top, as opposed to having two pairs of sinks and sources like in the case of the original synaptic placement (Figure 5F, top row; leftmost plot). The firing rate of the L4 excitatory cells, however, remained essentially unchanged by this modification (Figure 5G). On the other hand, placing all synapses from LGN and excitatory neurons onto the apical dendrites of L4 excitatory neurons resulted in even greater CSD magnitude from this population (Figure 5F, bottom row; leftmost plot), while the magnitude of its firing rates were reduced (Figure 5G). In this case, the pattern displayed a sink in the middle with a source above and below it. We also find that the somatic Hodgkin–Huxley channels significantly shapes the CSD pattern (Figure 5H), generating a source in the L4 region where the somata are localized.
We quantified the changes in CSD of the full model resulting from changes in recurrent synaptic weights in Figure 5—figure supplement 1 and the changes in firing rates resulting from the changes in synaptic placement onto L4 excitatory cells in Figure 5—figure supplement 2. Our analyses corroborate that the changes in CSD after changing only synaptic weights is minor, and that the changes in firing rate after changing only synaptic placement is minor. More examples of the disparate effects on CSD and spikes from adjusting placement of excitatory synapses on excitatory cells in L2/3 and L5 are displayed in Figure 5—figure supplement 3 and Figure 5—figure supplement 4, respectively.
These results indicate a two-way dissociation that can occur between CSD and firing rates of excitatory neurons. The firing rates can be changed without substantially changing the CSD by modifying the strength of synapses, while the CSD can be changed without substantially changing the firing rates by modifying synaptic location. This suggests that utilizing the CSD in the optimization of the model can provide constraints on the circuit architecture that could not be obtained from spikes alone.
Effects of feedback from higher visual areas to the model
Hartmann et al., 2019 found that feedback from higher visual areas (HVAs) can exert a powerful influence on the magnitude of the evoked LFP response recorded in V1 of macaque monkeys, particularly in the period 80–100 ms after stimulus onset. The sustained L2/3 sink and L4 source we observe in the experimental CSD emerge at 60 ms (Figure 2A and B), which roughly coincides with the peak firing rates in the latero-medial (LM), rostro-lateral (RL), antero-lateral (AL), and postero-medial (PM) cortical areas (Figure 3C). Furthermore, anatomical data indicate that synapses from HVAs terminate on L1 and L2/3 apical dendrites of pyramidal cells (whose cell bodies reside in L2/3 or L5) (Glickfeld and Olsen, 2017; Marques et al., 2018; Hartmann et al., 2019; Keller et al., 2020; Shen et al., 2020). Together, these observations suggest that the sustained L2/3 sink and L4 source might, in part, be induced by feedback from HVAs, where the sink is generated from the input to the apical tufts in L1 and L2/3, and the source may be the return currents of this input.
Of these HVAs, the feedback from LM to V1 is best characterized (Marques et al., 2018; Keller et al., 2020; Shen et al., 2020) and has the highest connection density to V1 (Harris et al., 2019). Based on these considerations, we decided to test the hypothesis that the large sinks and sources in the upper layers were caused, at least in part, by feedback from LM. In addition to the earlier feedforward LGN input and the background input representing the influence of the rest of the brain, we introduced a feedback input constructed from experimentally recorded spike trains in LM. In total, the public Neuropixels dataset has 2075 neurons recorded in LM (simultaneously with the recordings in LGN, V1, and other visual areas) from 42 animals during presentations of the full-field flash stimulus. 1823 of the 2075 neurons were classified as RS, and spike trains from these were used to generate the feedback input to the model (Figure 6A).

Introducing feedback from latero-medial (LM) to V1 in the model.
(A) Firing rate of the experimentally recorded lateral geniculate nucleus (LGN) and LM units used as input to the model. (B) Total current source density (CSD) resulting from simulation with input only from the LM. (C) Left: PC 1 CSD from experiments (see Figure 2). Right: total CSD from simulation with both LGN input and LM input. (D) Wasserstein distance between CSD from model versions and PC 1 CSD from experiments together with Wasserstein distances from CSD in animals to PC 1 CSD (boxplot). (E) Population contributions from populations that receive input from LM. (F) Average population firing rates of experiments (black line) and model versions. (G) Kolmogorov–Smirnov (KS) similarity between simulated firing rates or individual animals (boxplots) and recorded firing rates at baseline, peak evoked response, and the sustained period (defined in Figure 3). (H) Correlation between simulated and experimentally recorded firing rates (0–100 ms).
The synapses from this LM source were placed on the apical dendrites of L2/3 excitatory neurons (within 150 μm from the soma), on the apical tufts (>300 μm from the soma) and the basal dendrites (within 150 μm from the soma) of L5 excitatory neurons, and on the somata and basal dendrites of L2/3, and L5 inhibitory (Pvalb and Sst) neurons (at any distance from the soma). The input onto L2/3 excitatory neurons did generate a sink in L1 and L2/3 and a source below in L4 (Figure 6B–E).
The synaptic weights from LM to the populations targeted by the feedback were initialized at high values (see ‘Materials and methods’), and then adjusted (decreased) by multiplying them with factors in the range [0.05, 0.5] (see ‘Materials and methods’). The weights from the background to the feedback-targeted populations were also multiplied by factors in the range [0.2, 0.5], and the weights of connections from Pvalb neurons to L2/3 excitatory and L5 excitatory neurons were multiplied by factors in the range [0.8, 1.2]. This weight scaling was done until the population firing rates were within the experimental variability. Additionally, the synapses from excitatory populations onto L6 excitatory cells were restricted to be within 150 μm from the soma to reduce the magnitude of the L5/L6 dipole (Figure 6—figure supplement 1; see ‘Materials and methods’). Greater separation between a sink and a source in a dipole moment increases its magnitude, so restricting the range within which synapses can be placed should diminish the L5/L6 dipole’s dominance.
When the model received this feedback input together with the LGN input, the resulting CSD pattern reproduced the main features observed in the experiments (Figure 6C). The WD between the model CSD and the experimental PC 1 CSD was also no longer an outlier (Normalized WD = 0.41; Figure 6D), and the population firing rates remained within the minimum and maximum value of the experimental boxplots for the firing rates in all windows and all populations, both with respect to magnitudes (KSSb = 0.77, KSSp = 0.70, and KSSs = 0.68; average across all populations) and temporal profiles (RS L2/3: r = 0.36***, RS L4: r = 0.64***, RS L5: r = 0.69***, RS L6: r = 0.87***, FS: r = 0.77***, ***p<0.001) (Figure 6F and G). Thus, when average responses to the full-field flash are considered, this final adjusted model exhibits both the CSD and firing rate patterns that are consistent with the experimental observations and are well within animal variability (Figure 6F–H).
Furthermore, we investigated whether the model could reproduce the stereotypical features of the single-unit firing response observed in experiments. To this end, we computed the moments of the distribution of peak firing rates as well as the distribution of latencies to the peak across cells both for individual animals and for the model versions (Figure 6—figure supplements 2 and 3). The original model was an outlier for the first moment of the peak firing rate distributions of RS L2/3, L4, L6, and FS populations, while the intermediate model was within the experimental variability for all populations, and the final model was within the experimental variability for all populations except RS L2/3, where it was just outside the maximum value of the experimental boxplot (Figure 6—figure supplement 2A). All three model versions were outside the experimental variability of the third moment (skewness) of the peak firing rate distributions for the RS L2/3 population, and the original model was just outside for the RS L5 and FS populations as well (Figure 6—figure supplement 2C). Otherwise, all model versions were within the experimental variability for all moments that we considered and for all populations.
With respect to the distribution of latencies to the peak firing rates, all model versions were within the experimental variability for all populations and all moments except the fourth moment (kurtosis) of the FS population, where all model versions were outside the experimental variability. We also computed the relative change in firing rates between neighboring populations (Figure 6—figure supplement 4) and the distributions of greatest curvature of the firing rate across cells (Figure 6—figure supplement 5), and found that all model versions were within the experimental variability on both of these metrics too. Thus, as with the analysis of population firing rates, the original model is an outlier when compared to the experiments on the features of unit firing, while the intermediate and final model versions reproduce most features observed in the experiments.
To check if the model continued to exhibit appropriate orientation and direction tuning after the adjustments made, we ran a simulation with the final model configuration and the same drifting grating stimulus that was utilized in Billeh et al., 2020. We found that the model still displayed firing rates at preferred directions and direction and orientation selectivity indices comparable to those observed experimentally (Figure 6—figure supplement 6).
Identifying the biophysical origins of the canonical CSD
With the canonical CSD (Figure 2B) reproduced, we can use the model to probe the biophysical origins of its sinks and sources. We began by removing all recurrent connections and only feeding the LGN input to the model to find the contribution from the thalamocortical synapses onto excitatory and inhibitory neurons (Figure 7A). The main thalamic contribution to the CSD is from synapses onto excitatory neurons, in line with the expectation that neurons with a spatial separation between synaptic input currents and the return currents dominate the cortical LFP generation (Einevoll et al., 2013). (Neurons without apical dendrites will have largely overlapping synaptic input currents and return currents, resulting in a cancellation of current sinks and sources.).

Biophysical origin of canonical current source density (CSD).
(A) Sinks and sources generated from thalamocortical and (B) feedback synapses. The schematics illustrate which synapses cause the observed sinks and sources. Blue arrows indicate inflowing current (sinks), while orange arrows indicate outflowing current (sources). (C) Total CSD from thalamocortical and feedback synapses (without recurrent connections) with (left) and without (right) active channels in the V1 neurons. (D) Total CSD of model with both thalamocortical and feedback input when inhibitory synapses are removed (cross indicates removed connection). (E) Population contributions to the total CSD in final model with both lateral geniculate nucleus (LGN) and feedback input and recurrent connections. (F) Summary of biophysical origins of the main contributions to the sinks and sources in the canonical CSD in different periods of the first 100 ms after flash onset. More arrows mean more current. Left: before onset of evoked response (0–35 ms). The average inflowing and outflowing current in V1 neurons is zero in this time window. Middle: initial evoked response (35–50 ms). The L4 sink is primarily generated by inflowing current thalamocortical synapses onto L4 excitatory cells. Right: sustained evoked response (50–100 ms). The L5/L6 sink is primarily due to inflowing currents from thalamocortical synapses and recurrent excitatory synapses. Inflowing current at synapses from higher visual areas (HVAs) onto apical tufts of L2/3 and L5 excitatory cells generates, in part, the L2/3 sink, and the resulting return current generates, in part, the L4 source in this time window.
We further observed that the early L4 and the sustained L5/L6 sinks are present in the CSD contributions of excitatory neurons, though the magnitude of the L5/L6 sink is substantially reduced compared to its magnitude when the model is configured with recurrent synapses intact (Figures 4B—6). The sustained L2/3 sink and L4 source, on the other hand, were not visible. This suggests that the early L4 sink and the L5/L6 sink are at least partly generated by thalamocortical synapses. However, the substantially diminished magnitude of the L5/L6 sink indicates that recurrent synapses also contribute significantly to the generation of this sink.
We then removed the LGN input and added the feedback (while keeping the recurrent connections cut), which resulted in a prominent upper layer dipole, with the sink residing in L1 and L2/3, and the source residing in L4 (Figure 7B). Together with their absence when input came from LGN only (Figure 7A), this suggests that the sustained L2/3 sink and the L4 source in the canonical pattern originate at least in part from the feedback synapses onto the apical dendrites of L2/3 and L5 pyramidal cells and the activity this input generates.
To assess the extent to which active channels at the somata contributed to the CSD pattern, we compared the CSD resulting from a simulation with both LGN and feedback input (where the recurrent connections were still cut) when we included or excluded the active channels (NaT, NaP, NaV, h, Kd, Kv2like, Kv3_1, K_T, Im_v2, SK, Ca_HVA, Ca_LVA; only at the soma [see supplementary information in Gouwens et al., 2018 for definitions]) on all neurons in the model. The most prominent discrepancy between the CSD with and without active channels is the magnitude of the L4 source and the L2/3 sink (Figure 7C). In this all-passive setting, the L4 source is significantly attenuated, and the L2/3 sink is either absent or dominated by a source in the same region.
We explored whether the contributions from currents in recurrent connections come primarily from excitatory or inhibitory synapses by removing all connections from inhibitory (Pvalb, Sst, Htr3a) neurons to all other neurons, so that all postsynaptic currents stem from excitatory thalamocortical synapses, excitatory synapses from HVAs, or recurrent excitatory synapses in V1 (Figure 7D and Figure 7—figure supplement 1). Note that inhibitory synaptic currents give rise to sources, while excitatory synaptic currents give rise to a sink. Of course, without inhibition, the network is unbalanced, which limits the conclusions that can be drawn from this simulation. However, the fact that the major sinks and sources are still present is an indication that the currents from excitatory input account for the majority of the sinks and sources observed in the experimental CSD.
The contributions from each population to the total CSD in the final model (Figure 6D) with both LGN and feedback input and intact recurrent connections are displayed in Figure 7D. From this, it is apparent that the L5/L6 dipole is mainly generated by L6 excitatory cells, the L2/3 sink stems from sinks at the apical tufts of the L2/3 and L5 excitatory cells, the L4 sink from both the L4 excitatory and inhibitory cells, while the L4 source is a mix of sources from mainly L2/3, L4, and L5 excitatory cells, as well as the L4 inhibitory cells. (The magnitude of the CSD contribution from L4 inhibitory cells is greater than anticipated. Given their lack of apical dendrites, we would expect their postsynaptic current sinks and sources to largely cancel [Einevoll et al., 2013]. Their contribution can be reduced by scrambling the 3-D orientation of these cells [Figure 7—figure supplement 2]. However, we cannot rule out that L4 inhibitory cells can have a contribution comparable in magnitude to the excitatory cells with the data we have available. We therefore let the L4 inhibitory cells keep their original orientation here.)
We investigated what would happen if we turned the feedback off again at 60 ms for the final model version (Figure 7—figure supplement 3). Most notably, we found that the sustained L4 source was replaced by a sink, and that the sustained L2/3 sink turned more transient as it was significantly reduced in magnitude after about 70 ms. This CSD pattern serves as a prediction that can be tested experimentally, for example, by optogenetic silencing of the HVAs in the sustained periods of the evoked response (e.g., see Keller et al., 2020).
We summarize the main contributions to the canonical CSD in Figure 7F. Before the onset of the evoked response (0–35 ms) there is, on average, no significant net inflow or outflow of current to any neurons. Around 40 ms, an inflow of current from excitatory thalamocortical synapses onto all excitatory neurons and all Pvalb inhibitory neurons appears, with the largest current coming from the synapses targeting basal and apical dendrites of L4 excitatory cells. This is the primary origin of the L4 sink. Following this initial L4 sink, there is a sustained sink in L5/L6 arising at ~50 ms, which originates partly from thalamocortical synapses onto L6 excitatory cells and partly from recurrent synapses from excitatory populations in V1 onto L6 excitatory cells. At ~60 ms, a sustained sink emerges in L1 and L2/3, which partly originates in synapses from HVAs targeting apical tufts of L2/3 and L5 excitatory cells. This feedback results in a stronger return current at the soma and basal dendrites of L2/3 excitatory cells and L5 excitatory cells.
Discussion
In the present study, we analyzed experimentally recorded spikes and LFP during presentation of full-field flashes from a large-scale visual coding dataset derived from mouse visual cortex (Siegle et al., 2021) and simulated the same experimental protocol using a biophysically detailed model of mouse V1 (Billeh et al., 2020). Our analysis of the experimental data focused on the responses in LGN and cortical visual areas V1 and HVAs. We found that the evoked CSD in V1, computed from the LFP, is captured by a canonical pattern of sinks and sources during the first 100 ms after stimulus onset (Figure 2B). This canonical CSD, in response to a flashed, bright field pattern, explains half (50.4 %) of the variance in the trial-averaged CSD responses across animals.
Both the early L4 sink with concurrent sources above and below and the L5/L6 sink with a source below were observed with a similar timing by Senzai et al., 2019. The L4 source and L2/3 sink were also observed in that study but emerge somewhat later than in our data – just after 100ms as opposed to ~60 ms in our canonical pattern. This discrepancy in onset might simply be due to differences in stimuli. In Senzai et al., 2019, the animals were exposed to 100 ms light pulses, while the animals in our data were presented with 250 ms whole-field flashes of a white screen. Nonetheless, the canonical CSD pattern exhibits good overall agreement with the pattern seen in Senzai et al., 2019. Studies in non-human primates where the animals were exposed to flashes also demonstrate good spatio-temporal agreement with the CSD observed here, with sinks and sources occurring not only in the same layers and in the same order, but also at the same time after stimulus presentation (Givre et al., 1994; Schroeder et al., 1998).
We introduced the WD as a method to evaluate the difference between two CSD patterns and used it to quantify the variability in trial-averaged CSD between animals (Figure 2D), the trial-to-trial variability in CSD within animals (Figure 2F and G), and the difference between the model CSD, the trial-averaged CSD of individual animals, and the canonical CSD pattern. This application of the WD to compare CSD patterns comes with certain considerations that are important to note. First, although we compute WD for sinks and sources separately, sinks and sources do not arise independently. Current leaving the extracellular space in one place leads to current entering the extracellular space in another place, so current sinks and sources are inter-dependent. Second, the cost of shifting a sink or a source in space relative to shifting it in time is determined by the relative resolution in space vs. time. This relative cost does not necessarily correspond to the actual cost of changing the underlying physiology such that two distributions of sinks or sources match spatially vs. temporally. Determining the most appropriate relative cost of moving sinks and sources in space vs. time would require more detailed data than currently available and is beyond the scope of this study.
For the firing rate analysis, we utilized KS similarity and correlation to quantify experimental variability and model performance with regard to magnitude and temporal profile, respectively. We also investigated the statistics of unit spike firing by computing the moments of the distributions across cells in each population of peak firing and latencies to peak firing for both individual animals and model versions (Figure 6—figure supplements 2 and 3). Systematic use of quantitative metrics for biophysical modeling at this scale is still relatively uncommon, and our work establishes a set of measures for testing the model on LFP and spiking simultaneously, which can be useful for future studies in the field. Of course, there may well be other metrics that are equally or more suitable, and a systematic investigation into what would be the optimal metrics to apply is an important avenue for future work.
Our aim was to simultaneously reproduce experimentally recorded spikes and the CSD in our simulations. The original model captured spiking responses to gratings well (reproducing, e.g., direction selectivity distributions for different neuronal populations) with variable success when applied to other visual stimuli (Billeh et al., 2020). It was not originally tested on LFP/CSD. We found that, for the full-field flash stimulus, this model did not reproduce the CSD pattern in the upper layers of V1, and the spiking responses for this stimulus also exhibited a number of discrepancies.
After making selective adjustments to the recurrent synaptic weights, the model could reproduce the experimental firing rates (Figure 5A–C), though the discrepancy between the model CSD and the canonical CSD remained (Figure 5D and E), with only minor differences relative to the CSD of the original model (Figure 4B). The fact that the model can capture the experimental firing rates without capturing the experimental CSD and that adjustments to the synaptic weights yielded significant alterations in firing rates with only small changes in the CSD supports the point that LFP/CSD reflects aspects of circuit dynamics that are complementary to those reflected in locally recorded spikes.
Previous simulation studies demonstrated the importance of synaptic placement in shaping the LFP and CSD signature (Einevoll et al., 2007; Pettersen et al., 2008; Lindén et al., 2010; Lindén et al., 2011; Łęski et al., 2013; Hagen et al., 2017; Ness et al., 2018). To uncover the model adjustments that capture firing rates and CSD simultaneously, we explored the effects of changes in the synaptic positioning. In one case, we placed all excitatory synapses onto only basal or apical dendrites of L4 excitatory cells, as opposed to their original placement on both apical and basal dendrites. Moving all excitatory synapses onto basal dendrites resulted in substantial changes in both the pattern and magnitude of the CSD contribution from these L4 excitatory cells, with only minor changes to their firing rates (Figure 5F and G and Figure 5—figure supplement 1). Placing all excitatory synapses on apical dendrites led to somewhat larger changes in firing rates, though still similar to the firing rate of the original model, and to even bigger changes in the CSD magnitude. Performing the same manipulations of synaptic placement on L2/3 or L5 excitatory cells resulted in a similar decoupling of CSD and firing rates (Figure 5—figure supplements 3 and 4). It should, however, be noted that the decoupling is neither perfect nor universal.
This demonstrates a two-way dissociation of the firing rates and the pattern of sinks and sources in the CSD: The firing rates can be substantially altered with small effects on the CSD by adjusting the synaptic weights, and the CSD can be substantially altered with only small effects on the firing rates by adjusting synaptic placement. These results align with findings in, for example, Schroeder et al., 1998 and Leszczyński et al., 2020, where a lack of correlation between MUA and CSD in the upper layers of V1 during flash exposure to non-human primates suggested that these signals can sometimes be decoupled. Our findings imply that the LFP can reveal deficiencies in the model architecture that would not be evident from the firing rates alone, and that, to a certain extent, models can be optimized for firing rates and CSD independently.
Past studies have suggested that LFP can be modulated by attention through feedback during evoked responses (Mehta et al., 2000a; Mehta et al., 2000b). However, their findings indicated that V1 was not significantly affected by this modulation until the period 250 ms or later after stimulus onset. A recent study showed that feedback from higher visual areas can in fact exert a strong influence on the magnitude of LFP already at around 80 ms after stimulus onset (Hartmann et al., 2019). To investigate whether such cortico-cortical influence can contribute to the sinks and sources in the later stages (>50 ms) of the canonical CSD pattern, we added feedback consisting of experimentally recorded spikes from the higher cortical visual area LM (Siegle et al., 2021) impinging on synapses placed onto V1 neurons in our model, using anatomical data (Glickfeld and Olsen, 2017; Marques et al., 2018; Hartmann et al., 2019; Keller et al., 2020; Shen et al., 2020).
We found that the feedback can play a significant role in shaping the sustained sinks and sources (Figure 6B–E). The resulting model CSD reproduced the major sinks and sources identified in the canonical CSD pattern and was no longer among the outliers compared to the experimental variability (Figure 6D). Interestingly, absence of the feedback was not apparent from analysis of the firing rates alone, as the firing rates were already within the experimental variability before adding the feedback, further underscoring the utility of the LFP in illuminating structure–function relations in the circuit. Contributions from other visual cortical areas were not included, even though they too impinge upon neurons in V1 (Harris et al., 2019; Siegle et al., 2021), due to the lack of data characterizing such connections. This awaits future work. Finally, turning off the feedback from LM at 60ms disrupted the CSD pattern in layers 2/3 and 4 during the sustained period (Figure 7—figure supplement 3), which serves as a prediction for what will be observed if HVAs, and particularly LM, are silenced in this period of an evoked flash response. Our findings accord with the view that basal dendrites are the main targets for feedforward input, while the tufts of apical dendrites are the main targets for feedback input, even though basal dendrites can also receive input from long-range feedback connections and apical dendrites receive input from feedforward connections (Aru et al., 2020).
In our analyses, we aligned the CSD plots to depths obtained from histology and provided in Allen Common Coordinate Framework (CCF) coordinates. An alternative approach is to align the CSD to landmarks in the data (Senzai et al., 2019). In Figure 6—figure supplement 7, we explored whether utilizing landmarks rather than histology for alignment would affect our results. We found that this approach did not change our conclusion as the final model CSD was still within the experimental variability while the original and intermediate model CSD were outliers both when using the first principal component and when using the plain average as the target (Figure 6—figure supplements 7B and C and 8). Lastly, using the plain average rather than the first principal component as the target did not significantly affect our results when we aligned to histology (Figure 6—figure supplement 9).
With the major sinks and sources of the canonical CSD pattern reproduced, we explored their biophysical origins. We found that the initial L4 sink originates in the thalamocortical input to L4 excitatory cells, which aligns with suggestions made in Mitzdorf, 1987, Swadlow et al., 2002, and Senzai et al., 2019. The sustained L5/L6 sink comes from postsynaptic currents in L6 excitatory cells triggered by a combination of thalamocortical and recurrent excitatory inputs. The sustained L2/3 sink stems, in part, from input from LM onto the apical tufts of L2/3 and L5 excitatory cells. The sustained L4 source has its origins in a mixture of return currents from L2/3 and L5 excitatory cells resulting from the abovementioned feedback onto the apical dendrites of these cells, as well as contributions from L4 excitatory and inhibitory cells (Figure 7A, B, D and E).
In line with observations made by Reimann et al., 2013, we found that the somatic voltage-dependent membrane currents significantly shape the CSD signature (Figures 5H and 7C). Even so, our findings still emphasize the importance of synaptic inputs in sculpting the CSD, as the addition of synaptic input (Figure 6A–E) and changes to synaptic placement (Figure 5F) substantially altered the CSD pattern.
Recent investigations into the unitary LFP (the LFP generated by a single neuron) from inhibitory and excitatory synapses have suggested that inhibitory inputs exert a greater influence on LFP than excitatory input (Bazelot et al., 2010; Teleńczuk et al., 2017; Telenczuk et al., 2020). While this may be true of unitary effects, the total effect of excitatory input can still be greater if there are significantly more excitatory than inhibitory cells, and, correspondingly, significantly more excitatory synapses. In this V1 model, inhibitory cells make up about 15% and excitatory cells about 85% of the total number of cells, reflecting the cellular composition in mouse V1. Whether the pronounced dominance of excitatory cells is enough to make up for the reduced unitary influence of excitatory cells is an interesting area for further research.
This investigation into the biophysical origins of sinks and sources is limited by the fact that the contributions from recurrent connections are difficult to estimate precisely due to the nonlinear effects of these connections. That is, their contribution cannot simply be found by subtracting the CSD from thalamocortical and feedback synapses with all recurrent connections removed (Figure 7A and B) from the total CSD with the same input and recurrent connections intact (Figure 6C, right). Still, this analysis provides an initial estimate into the biophysical origins of the sinks and sources observed experimentally and demonstrates the insights that can be obtained from modeling of extracellular signals.
There is ample evidence that firing rates and LFP are modulated by the behavioral state of the animal, including measures like the pupil size (considered to be a proxy for arousal level) or running speed (Niell and Stryker, 2010; McGinley et al., 2015; Vinck et al., 2015; Saleem et al., 2017). In this study, the responses averaged over all trials have been the target for the modeling, without regard to any state-dependence of the responses. Our understanding of the state-dependent responses could benefit from the potential to probe the biophysical origins of extracellular signals. Therefore, reproducing these state-dependent responses is an important avenue for future research.
Note that the set of synaptic weights and other parameters that can reproduce the experimental firing rates and CSD is unlikely to be unique. This is a consequence of the degeneracy inherent to biological neural networks, known from both simulation and experimental studies, as many different parameterizations of neuronal networks can perform the same functions (Prinz et al., 2004; Marder and Goaillard, 2006; Drion et al., 2015; O’Leary, 2018). Thus, our network should only be considered an example of a circuit model that can produce firing rates and CSD that match the experimental observations. Obtaining multiple solutions and characterizing their diversity using automatic searches of the parameter space will be an interesting direction for future work. We did not utilize such an approach here because the number of simulations required (typically, many thousands or more for automatic optimization approaches) would currently be prohibitively expensive on a model of this scale and level of complexity: running a 1 s simulation with this model takes ~90 min on 384 CPU-cores (Billeh et al., 2020); a single trial in this study simulates 0.75 s of activity.
The original model used as a starting point here produced firing rates and direction and orientation tuning consistent with recordings during presentations of drifting gratings (Arkhipov et al., 2018; Billeh et al., 2020). In this study, we focused on the analysis and modeling of the response to full-field flashes, but when the final model was tested with the drifting gratings stimulus utilized in Billeh et al., 2020, we found that the present, revised model continued to exhibit orientation and direction tuning even though the adjustments were made with the aim to reproduce the observed CSD and population firing rates during full-field flashes (Figure 6—figure supplement 6). Ideally, the model should reproduce both firing rates and LFP simultaneously not only for flashes or drifting gratings, but for any visual stimulus (out-class generalization). This is a long-term goal and can be called ‘the holy grail’ of visual system modeling.
In this study, we developed a systematic framework to quantify experimental variability in both LFP/CSD and spikes and to evaluate model performance. We identified a canonical CSD pattern observed during presentations of full-field flash stimuli and obtained a bio-realistic model that reproduced both the canonical CSD pattern and spikes simultaneously. This model thus reproduces, in a quantitative manner, the shape and timing of current sinks and sources observed experimentally. We utilized this validated model to explain, mechanistically, the biophysical origins of the various current sinks and sources and their location across the layers of visual cortex. Our models are freely shared and should be useful for future studies disentangling the mechanisms underlying spiking dynamics and electrogenesis in the cortex.
Materials and methods
Experiments
Quality control
Request a detailed protocolOf the 58 mice in the visual coding dataset, 9 were excluded because the exact probe location could not be recovered due to fading of fluorescent dye or artifacts in the OPT volume (Siegle et al., 2021). Another five animals were excluded because they were missing LFP recordings from V1 during presentation of the flash stimulus. Thus, data for 44 animals were retained for the CSD analysis.
For the spike analysis, the same nine animals for which the exact probe location could not be recovered were excluded, and two additional animals were excluded because they did not have any cells recorded in V1, leaving a total of 47 animals for this part of the data analysis.
Neuronal classification
Request a detailed protocolWe distinguished between RS and FS cells by the time from trough to peak of the spike waveforms (Barthó et al., 2004). For cortical cells, the spike duration was bimodally distributed with a dip at ~0.4 ms, while for thalamic cells, it was bimodally distributed with a dip at ~0.3 ms (Figure 3—figure supplement 1). Thus, the cutoff in the classification of cells as RS or FS was set at 0.4 ms for LM and V1, and at 0.3 ms for cells in LGN. Note that several studies have demonstrated that some pyramidal neurons may have spike waveforms short enough to be classified as FS cells (Vigneswaran et al., 2011; Lemon et al., 2021). Thus, some caution is warranted when interpreting the population firing rates.
When comparing the model firing rates to the experimental firing rates, the excitatory and non-Pvalb populations were grouped together in each layer of the model to make up the RS cells in L2/3, L4, L5, and L6, while the Pvalb cells across all layers were grouped together to make up the FS cells of V1. The layer boundaries were taken from the Allen CCF (Oh et al., 2014), allowing for the assignment of each neuron’s position to a specific cortical layer (Siegle et al., 2021).
Model
Request a detailed protocolThe model consists of both biophysically detailed multicompartment neurons and leaky-integrate-and-fire (LIF) point-neurons. In total, there are 51,978 multicompartment neurons with Hodgkin–Huxley conductances at the soma and only passive conductances at the dendrites. These are arranged in a cylinder of radius of 400 μm and height 860 μm (corresponding to the average cortical thickness of V1 taken from the Allen CCF; Billeh et al., 2020; Oh et al., 2014). This cylinder makes up the ‘core’ of the model and is surrounded by an annulus of 178,946 are LIF neurons which has the same height and a thickness of 445 μm. This makes the total number of neurons in the model 230,924 and the radius of the whole cylinder with both biophysically detailed and LIF neurons 845 μm. There are 17 different classes of neuron models. In each layer from 2/3 to 6, there are one excitatory and three inhibitory classes (Pvalb, Sst, Htr3a), while in layer 1 there is a single Htr3a class. The LGN module providing thalamocortical input to the model consists of 17,400 units selectively connected to the excitatory neurons and Pvalb neurons in L2/3 to L6, as well as the non-Pvalb neurons in L1. The background input to all neurons in the model comes from a single Poisson source firing at 1 kHz and represents influence from the rest of the brain. The feedback input to L2/3 and L5 excitatory, Pvalb, and Sst neurons comes from a node representing LM.
Simulation configuration
Request a detailed protocolInstructions on how to run simulations of the model are provided in Billeh et al., 2020. The files and code necessary to run the model versions presented in Figures 4—6 are provided in the directories old_model_fig4 intermediate_model_fig5, and final_model_fig6, respectively, on Dryad (see ‘Data availability’).
Data processing
LFP and CSD
Request a detailed protocolThe LFP in simulations was obtained from the extracellular potential by first downsampling to every other electrode along the probe (resulting in a spatial separation of 40 μm between each recording electrode, equal to the spacing in the public Neuropixels data) and using a low-pass fifth-order Butterworth filter with a cutoff frequency of 500 Hz (utilizing functions scipy.signal.butter and scipy.signal.filtfilt). The same filtering was applied to get the experimental LFP. The CSD was calculated from the experimental and model LFP using the delta iCSD method (Pettersen et al., 2006), where the radius of laterally (orthogonal to the probe axis) constant CSD was assumed to be 400 μm – the radius of the V1 model’s ‘core’ region consisting of biophysically detailed multicompartment neurons. For the experimental CSD, this radius was set to 800 μm, roughly corresponding to the size of mouse V1.
Visual stimulus
The stimulus used to compare the model and the experiments was full-field flashes. In the experiments, the mice were presented with gray screens for 1 s, followed by 250 ms of white screen, and then 750 ms of gray screen over 75 trials. In the simulations, both the stimulus presentation and the pre- and the poststimulus gray screen periods lasted 250 ms, and the number of trials was 10.
Input from LGN)
Request a detailed protocolOriginally, the LGN spike trains used as input to the model were generated with the FilterNet module provided with the model, using 17,400 ‘LGN units’ (Billeh et al., 2020). However, when this input was used for simulations, the onset of the evoked response in V1 was 20–30 ms delayed in comparison with experiments. Therefore, we used experimentally recorded LGN spike trains as input to the model instead. We assigned a recorded spike train to each of the 17,400 LGN units in all trials. In total, the public Neuropixels data contain recordings from 1263 regular-spiking LGN neurons across 32 animals during 75 trials of full-field flash presentations. We divided the total pool of spike trains into 10 subsets, and then randomly sampled spike trains from one subset in each trial until all 17,400 LGN units had been assigned a spike train in all trials.
Input from LM
Request a detailed protocolThe experimentally recorded spike trains in the LM were used to construct the feedback input to V1. In total, the public Neuropixels data contain recordings from 1823 RS LM neurons across 42 animals during presentations of the full-field flash stimulus. Spikes were randomly sampled from the pool of all spike trains to construct a spike train that was used as input to all the cells that were targeted by the feedback in the model. All neurons received the same spike train.
Background input
Request a detailed protocolThe input from the Poisson source firing at 1 kHz was not stimulus dependent. It is a coarse representation of the continuous influence of the rest of the brain on V1.
Dendritic targeting
The rules for placement of synapses were set in Billeh et al., 2020 and were based on reviews of literature on anatomy.
LGN to V1
Request a detailed protocolIn the original model, the synapses from LGN onto excitatory V1 neurons were placed on apical and basal dendrites within 150 μm from the soma, while synapses onto inhibitory V1 neurons were placed on their soma and on their basal dendrites without distance limitations (Billeh et al., 2020). This placement was left unchanged in this study.
V1-V1
Request a detailed protocolThe synapses for recurrent connections were placed according to the following rules in the original model (Billeh et al., 2020):
Excitatory-to-excitatory connections
Request a detailed protocolAll synapses from excitatory V1 neurons onto other excitatory V1 neurons were placed along the dendrites and avoided the soma. In layers 2/3 and 4, the placement of synapses was restricted to be within 200 μm from the somata, while in layers 5 and 6, they could be placed anywhere along the dendrites.
Excitatory-to-inhibitory connections
Request a detailed protocolAll synapses from excitatory V1 neurons onto inhibitory V1 neurons were placed on their somata or dendrites without any distance limitations.
Inhibitory-to-excitatory connections
Request a detailed protocolSynapses from Pvalb neurons onto excitatory V1 neurons were placed on the soma and on the dendrites within 50 μm from the soma. Synapses from Sst neurons were placed only on dendrites and only more than 50 μm from the soma. Synapses from Htr3a neurons were placed on dendrites between 50 and 300 μm from the soma.
Inhibitory-to-inhibitory connections
Request a detailed protocolSynapses from inhibitory neurons to other inhibitory neurons were placed according to the same rules as the inhibitory-to-excitatory connections described above.
These placement rules were kept in this study, except for the synapses from excitatory neurons to excitatory L6 neurons. Here, they were restricted to be within 150 μm of the soma. The purpose of this restriction was to reduce the spatial separation between the current sink and source, and thereby decrease the magnitude of the L6 sink-source dipole.
LM-V1
Request a detailed protocolThe synapses from the node representing LM to V1 were placed on the apical dendrites of L2/3 neurons (within 150 μm from the soma), on the apical tufts (>300 μm from the soma) and the basal dendrites (within 150 μm from the soma) of L5 excitatory cells, and on the somata and basal dendrites of L2/3 and L5 inhibitory cells (at any distance from the soma).
Adjusting synaptic weights
Request a detailed protocolIn the original model, the synaptic weights of thalamocortical connections were based on experimental recordings of synaptic current, while the synaptic weights for recurrent connections were initially set to estimates from literature, then optimized to a drifting gratings stimulus until the model reproduced experimental values of orientation and direction selectivity. In this study, the synaptic weights for thalamocortical connections were left unchanged from the original model. Before the addition of feedback from higher visual areas to the model, the synaptic weights for recurrent connections in V1 were multiplied by factors in the range [0.2, 2.5].
In the original model, the input from the background node represented the influence of the rest of the brain on V1, which included the influence from higher visual areas such as LM. This means that some of the feedback influence from LM on V1 should be present (though coarsely represented) in the input from the background node. When the influence of input from LM to feedback-targeted cells is modeled on its own, the influence of the background node must be updated accordingly. Thus, after the addition of feedback, the synaptic weights from the background node to the populations targeted by feedback (the L2/3 and L5 excitatory, Pvalb, and Sst cells) were multiplied by factors in the range [0.2, 0.5]. The synaptic weights from the node representing LM were initially set equal to the original weights between the background node and the populations targeted by the feedback, but this led to too high firing rates compared to the experimental firing rates in these populations, so they were multiplied by factors in the range [0.2, 0.5]. Finally, the connections from Pvalb neurons in V1 to L2/3 excitatory neurons and L5 excitatory cells were re-scaled in the range [0.8, 1.2] times the weights set prior to the addition of feedback. The ranges reported here were set after experimenting with different ranges to find what would allow the model to reproduce the experimental observations. Only a single value within each range was used in the final model.
Quantification and statistical analysis
Firing rates
Request a detailed protocolThe time-resolved population firing rates (bin size 1 ms, filtered using scipy.ndimage.gaussian_filter with sigma = 2) were computed by averaging the spike count over all cells in a population and over all trials (10 trials in the simulations and 75 trials in the experiments). The distribution of firing rates across cells used in the calculation of the KS similarities was computed by averaging over the time windows baseline, initial peak, and sustained activity (defined in Figure 3) and over all trials.
Kolmogorov–Smirnov similarity
Request a detailed protocolThe KS similarity scores (Billeh et al., 2020) were computed by first calculating the KS distance (using the function scipy.stats.ks_2samp) between two distributions of firing rates across cells, and subtracting this number from 1, such that a KS similarity score of 1 implies identity and a score of 0 implies no overlap between the two distributions. In the comparison of the model to the experimental data, the KS similarity was computed between the distribution of firing rates across cells in each RS and the FS population of the model and the distribution of firing rates across cells from all animals in the corresponding populations. To assess the variability in the experiments, the KS similarity was calculated between the distribution of firing rates across cells in the same RS and FS populations in individual animals, provided there were more than 10 cells recorded in a given population in this animal, and the distribution of firing rates across cells from all other animals.
Correlation
Request a detailed protocolWe computed the similarity in the profile of time-resolved population firing rates with the Pearson correlation coefficient (using the function scipy.stats.pearsonr). The correlation between the model and the experimental firing rates was calculated between model population firing rates and the population firing rates averaged across cells from all animals. The level of experimental variability was assessed by calculating the correlation between population firing rates in each animal and the population firing rates averaged across cells from all other animals.
CSD analysis
Since the number of recording electrodes in V1 are not the same in all animals, we interpolated the CSD of each animal and the CSD from simulations onto dimensions of the same lengths ( points along the depth and points along the time axis for 100 ms time windows) before we quantitatively analyzed the CSD.
PCA
Request a detailed protocolThe trial-averaged CSD of each animal was flattened into a vector of length , and the vectors of all animals were stacked together into a matrix of size 44 × 3000. Then, we performed PCA (using sklearn.decomposition.PCA) on this matrix to obtain the principal components that would constitute sums of weighted contributions of the trial-averaged CSD patterns.
Wasserstein distance (WD)
Request a detailed protocolThe first Wasserstein distance between two distributions and is defined as
where is the cost of moving a unit ‘mass’ from position to following the optimal transport plan in all transport plans (Rubner et al., 1998; Arjovsky et al., 2017).
In the utilization of WD to quantify the similarity between two CSD patterns, the distance between the distribution of sinks in the two patterns and the distance between distribution of sources of the two patterns are calculated separately and summed to form a total WD between the two CSD patterns:
where P1 and P2 refer to the two CSD patterns. The Python Optimal transport library (https://pythonot.github.io/index.html) was used to implement this calculation.
Appendix 1

Effect of adding feedback connections from latero-medial (LM) to L1 inhibitory cells.
(A) Left: current source density (CSD) from simulation of final model version without connections from LM to L1 inhibitory cells. Right: CSD from simulation of final model version with connections from LM to L1 inhibitory cells included. (B) Population firing rates of experiments (black line), final model without connections between LM and L1 inhibitory cells (orange line), final model with connections between LM and L1 inhibitory cells included (purple colored line).
Data availability
The files necessary to run simulations of the different model versions presented in the paper as well as data resulting from simulations of those model versions are publicly available in Dryad: https://doi.org/10.5061/dryad.k3j9kd5b8. The experimental data set utilized is publicly available at: https://portal.brain-map.org/explore/circuits/visual-coding-neuropixels. The code generated for data analysis and producing the figures in this manuscript is publicly available at: https://github.com/atleer/CINPLA_Allen_V1_analysis.git (copy archived at Rimehaug, 2023).
-
Dryad Digital RepositoryUncovering circuit mechanisms of current sinks and sources with biophysical simulations of primary visual cortex.https://doi.org/10.5061/dryad.k3j9kd5b8
-
DANDI Archive 000021ID 000021. 20191003_AIBS_mouse_ecephys_brain_observatory_1_1.
References
-
Selecting the signals for a brain-machine interfaceCurrent Opinion in Neurobiology 14:720–726.https://doi.org/10.1016/j.conb.2004.10.005
-
ConferenceWasserstein generative adversarial networksProceedings of the 34th International Conference on Machine Learnin. pp. 214–223.
-
Visual physiology of the layer 4 cortical circuit in silicoPLOS Computational Biology 14:e1006535.https://doi.org/10.1371/journal.pcbi.1006535
-
Cellular mechanisms of conscious processingTrends in Cognitive Sciences 24:814–825.https://doi.org/10.1016/j.tics.2020.07.006
-
Unitary inhibitory field potentials in the CA3 region of rat hippocampusThe Journal of Physiology 588:2077–2090.https://doi.org/10.1113/jphysiol.2009.185918
-
Low-frequency local field potentials and spikes in primary visual cortex convey independent visual informationThe Journal of Neuroscience 28:5696–5709.https://doi.org/10.1523/JNEUROSCI.0009-08.2008
-
BookRhythms of the BrainOxford university press.https://doi.org/10.1093/acprof:oso/9780195301069.001.0001
-
The origin of extracellular fields and currents--EEG, ECoG, LFP and spikesNature Reviews. Neuroscience 13:407–420.https://doi.org/10.1038/nrn3241
-
Brain Modeling ToolKit: An open source software suite for multiscale modeling of brain circuitsPLOS Computational Biology 16:e1008386.https://doi.org/10.1371/journal.pcbi.1008386
-
Laminar analysis of extracellular field potentials in rat vibrissa/barrel cortexJournal of Neurophysiology 63:832–840.https://doi.org/10.1152/jn.1990.63.4.832
-
Modelling and analysis of local field potentials for studying the function of cortical circuitsNature Reviews. Neuroscience 14:770–785.https://doi.org/10.1038/nrn3599
-
Focal local field potential signature of the single-axon monosynaptic thalamocortical connectionThe Journal of Neuroscience 37:5123–5143.https://doi.org/10.1523/JNEUROSCI.2715-16.2017
-
LFP power spectra in V1 cortex: the graded effect of stimulus contrastJournal of Neurophysiology 94:479–490.https://doi.org/10.1152/jn.00919.2004
-
Electrical interactions via the extracellular potential near cell bodiesJournal of Computational Neuroscience 6:169–184.https://doi.org/10.1023/a:1008832702585
-
Frequency dependence of signal power and spatial reach of the local field potentialPLOS Computational Biology 9:e1003137.https://doi.org/10.1371/journal.pcbi.1003137
-
Intrinsic dendritic filtering gives low-pass power spectra of local field potentialsJournal of Computational Neuroscience 29:423–444.https://doi.org/10.1007/s10827-010-0245-4
-
Local field potential in cortical area MT: stimulus tuning and behavioral correlationsThe Journal of Neuroscience 26:7779–7790.https://doi.org/10.1523/JNEUROSCI.5052-05.2006
-
Parallel readout of pathway-specific inputs to laminated brain structuresFrontiers in Systems Neuroscience 5:77.https://doi.org/10.3389/fnsys.2011.00077
-
Variability, compensation and homeostasis in neuron and network functionNature Reviews. Neuroscience 7:563–574.https://doi.org/10.1038/nrn1949
-
Optimizing the decoding of movement goals from local field potentials in macaque cortexThe Journal of Neuroscience 31:18412–18422.https://doi.org/10.1523/JNEUROSCI.4165-11.2011
-
Inference of hand movements from local field potentials in monkey motor cortexNature Neuroscience 6:1253–1254.https://doi.org/10.1038/nn1158
-
Properties of the evoked potential generators: current source-density analysis of visually evoked potentials in the cat cortexThe International Journal of Neuroscience 33:33–59.https://doi.org/10.3109/00207458708985928
-
Human intracranial recordings and cognitive neuroscienceAnnual Review of Psychology 63:511–537.https://doi.org/10.1146/annurev-psych-120709-145401
-
Stimulus contrast modulates functional connectivity in visual cortexNature Neuroscience 12:70–76.https://doi.org/10.1038/nn.2232
-
h-type membrane current shapes the local field potential from populations of pyramidal neuronsThe Journal of Neuroscience 38:6011–6024.https://doi.org/10.1523/JNEUROSCI.3278-17.2018
-
Highly selective receptive fields in mouse visual cortexThe Journal of Neuroscience 28:7520–7536.https://doi.org/10.1523/JNEUROSCI.0623-08.2008
-
Homeostasis, failure of homeostasis and degenerate ion channel regulationCurrent Opinion in Physiology 2:129–138.https://doi.org/10.1016/j.cophys.2018.01.006
-
Amplitude variability and extracellular low-pass filtering of neuronal spikesBiophysical Journal 94:784–802.https://doi.org/10.1529/biophysj.107.111179
-
Estimation of population firing rates and current source densities from laminar electrode recordingsJournal of Computational Neuroscience 24:291–313.https://doi.org/10.1007/s10827-007-0056-4
-
Similar network activity from disparate circuit parametersNature Neuroscience 7:1345–1352.https://doi.org/10.1038/nn1352
-
Electrophysiology of a dendritic neuron modelBiophysical Journal 2:145–167.https://doi.org/10.1016/s0006-3495(62)86953-7
-
Theoretical reconstruction of field potentials and dendrodendritic synaptic interactions in olfactory bulbJournal of Neurophysiology 31:884–915.https://doi.org/10.1152/jn.1968.31.6.884
-
Encoding of movement direction in different frequency ranges of motor cortical local field potentialsThe Journal of Neuroscience 25:8815–8824.https://doi.org/10.1523/JNEUROSCI.0816-05.2005
-
SoftwareCINPLA Allen V1 analysis, version swh:1:rev:bc32ab40499ba9b7687a1388f461cdf158bba375Software Heritage.
-
The pre-movement component of motor cortical local field potentials reflects the level of expectancyBehavioural Brain Research 169:335–351.https://doi.org/10.1016/j.bbr.2006.02.004
-
ConferenceA metric for distributions with applications to image databasesIEEE 6th International Conference on Computer Vision. pp. 59–66.https://doi.org/10.1109/ICCV.1998.710701
-
The spiking component of oscillatory extracellular potentials in the rat hippocampusThe Journal of Neuroscience 32:11798–11811.https://doi.org/10.1523/JNEUROSCI.0656-12.2012
-
Activation of a cortical column by a thalamocortical impulseThe Journal of Neuroscience 22:7766–7773.https://doi.org/10.1523/JNEUROSCI.22-17-07766.2002
-
A kernel-based method to calculate local field potentials from networks of spiking neuronsJournal of Neuroscience Methods 344:108871.https://doi.org/10.1016/j.jneumeth.2020.108871
-
Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic burstsJournal of Neurophysiology 93:2194–2232.https://doi.org/10.1152/jn.00983.2004
-
Population encoding of spatial frequency, orientation, and color in macaque V1Journal of Neurophysiology 72:2151–2166.https://doi.org/10.1152/jn.1994.72.5.2151
Decision letter
-
Tirin MooreSenior and Reviewing Editor; Howard Hughes Medical Institute, Stanford University, United States
-
Stephanie R JonesReviewer; Brown University, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]
Thank you for submitting the paper "Uncovering circuit mechanisms of current sinks and sources with biophysical simulations of primary visual cortex" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Antonio Fernandez-Ruiz (Reviewer #1); Stephanie R Jones (Reviewer #3).
Reviewer #1 (Recommendations for the authors):
Rimehaug et al., addressed two important questions in this work: the contribution of different circuit components to the generation of network evoked responses and the relationship between single cell firing rates and mesoscopic currents. They used a rigorous quantitative approach combining analysis of experimental data and simulations to address these questions in the mouse visual primary visual cortex. The authors start by characterizing the stereotypical V1 laminar response to full-filed light flashes, which replicates previous studies. They then perform a detailed quantification of within and across animal variability of both mesoscopic and single unit firing rate responses. Next, they attempt to reproduce these responses in a previously published biophysical network model that was shown to reproduce realistic orientation and direction tuning cell response. However, they found that the model did not accurately reproduced responses at the level of local field potentials (LFP) and current source density (CSD), so they proceeded to systematically vary model parameters until obtaining a better fit. One of the main findings of the article is that varying the strength of some synapses in model had a strong effect on firing rates but small on the CSD while varying synapse dendritic placement had strong effects on the CSD but small on firing rates. In the last part of the paper, the author use the fully tuned model to provide mechanistic insights into the biophysical basis of the different component of CSD evoked responses.
The results reported in this manuscript and the methods employed are sound, detailed and rigorous and accurately support the authors interpretations. I found of special relevance the demonstration of the importance of incorporating LFP/CSD features to both understand network dynamics a constrain models. I don't have any major concern with this manuscript but suggest below some points where it can be improved.
– Some of the analysis and simulations can be better motivated. For example, it is not clear, in the section starting in line 244, why it is important to quantify CSD variability and how that relates to the broader goals of the study. The manuscript in quite dense and in occasions is not easy to follow how certain results contribute to the main objective. It would help to state these explicitly at the outset and related each subsequent section to them.
– While the dissection of CSD components is detailed and exhaustive, the analysis of unit firing response fall short. In previous works (like Senzai et al., 2018) a stereotypical timing of single unit firign to this type of stimuli is described in addition to the CSD responses. In the paper there are Can the author's model reproduce this timing of unit activation? In addition to reporting average firing rates, the authors should also describe and interpret the temporal dynamics of unit firing and relate it to their mechanistic description of CSD components in the last section.
– When discussing the biophysical mechanisms of CSD generation a useful concept to introduce is that of dipolar moment. Components with a larger dipolar moment (more separation between sink and source) generate larger LFPs.
– To appeal to a broader audience, it would be useful to related the present findings with both, other mesoscopic network patterns in mouse V1 and in different cortical areas (e.g. sensory evoked γ oscillations, UP-DOWN states, etc.)
The authors should do a better job with citations. Most citations about broad topics like LFP biophysics and modelling or separation and characterization of the generating mechanisms of LFP/CSD components are from the same set of authors, that include those of the present paper. Many other groups have also produced relevant work in these topics, even if in different brain structures (as the hippocampus) or animal models (NHPs) but are not cited here.
Reviewer #2 (Recommendations for the authors):
Using a biophysically detailed model of mouse primary visual cortex, the evoked response to a full-field bright flash is simulated and compared with experimental data.
This paper aims to show how local field potential signal can be helpful in characterizing neural circuitry, and how it provides complementary information to the analysis of neural firing. Spiking patterns alone, they argue, do not uniquely characterize a neural circuit, i.e. there are many models (or variations to the same model) that are capable of reproducing the spiking activity of data, so requiring the model to simultaneously reproduce the CSD would further restrict the model.
The authors try to demonstrate this point by showing that there is a two-way dissociation between CSD pattern and the firing rate pattern. By adjusting synaptic weights they can reproduce the experimental firing rates statistics with minor changes in CSD and by adjusting synaptic placement they can alter the CSD patterns with minor changes in firing rates statistics.
Another support for this claim is that authors can better match both CSD patterns and firing rates by adding input that simulates a feedback from higher visual cortical areas. When feedback is incorporated in the model the performance must improve and it indeed does so based on CSD and firing rate statistics.
In the final section, after having a well-performing model, authors try to identify the pre and postsynaptic populations responsible for different prominent sinks and sources in evoked CSD pattern. They do this by manipulating different synaptic connections in the model: LGN input, cortical top-down feedback, excitatory recurrent connections onto different target neural populations. By comparing the CSD with only one pathway intact and measuring the contribution of that component to the total CSD pattern the authors qualitatively demonstrate contribution of these input-target components to the CSD.
The main approach chosen in the paper – linking detailed experimental data characterizing the evoked response in V1 (firing patterns of a large number of neurons across layers and circuits and LFP/CSD patterns) and a very detailed model of the circuit is a highly ambitious goal. This reviewer is not sharing the authors' optimism that it is achievable at this point beyond illustrations of how different elements of the model could potentially contribute to some elements of the experimental data. The key novelty of this manuscript compared to a previous paper (Billeh et al. 2020) is adding CSD patterns as an additional constraint to more accurately infer the model parameters beyond its prior state. While it is true that adding constraints would undoubtedly be helpful, there is no guarantee that specific hyperparameters used, components of the model simplified or ignored and gross aggregated measures chosen for comparing the patterns in the data and produced by the model would allow identifiability of the model and use it to explain the experimental data. Authors make no justification for many specific choices for a number of unknown high-dimensional parameters e.g. morphology of particular cell types (e.g. L4 pyramids have minimal basal dendrites, Scala et al. 2019, which has a strong effect on Figure 5 results), somato-dendritic distribution of synaptic inputs from different afferents, choice of interneurons (morphology, intrinsic properties, and connectivity) to include and how to wire them in recurrent and feedforward circuit, etc. Given the nonlinear nature of the network dynamics driven by the thalamic input and its nontrivial dependence on very high-dimensional parameter space, I see no way to identify a true model in that space based on very coarse metrics of aggregated neural activity. Care for realism in the model and the use of real data provides a false impression that the outcome of the effort should be the identification of the actual unobserved parameters of the system and the identification of the model. I think that it is a false impression, that might be misleading the readers.
The choice of metrics and dimensionality reduction approaches used to capture the group results on the CSD and firing rate is not suitable for validating the model or might provide poorly interpretable results. The authors note in the discussion, that there might be more suitable metrics and promise to work on it in the future. However, if the choice of metric used in the present paper makes interpretation questionable it puts in question the validity of all the results and conclusions in the manuscript, except for the actual data and model it is based on. Thus, I find it critical that the analysis methods used are chosen carefully and shortcomings discussed below are taken into account.
Inter-animal CSD alignment. It is obvious from the examples and prior work (Senzai et al. 2019) that variation of the location and scale of many CSD dipoles varies greatly. Any further analysis that relies on the CSD has to be thus rescaled and realigned, as Senzai et al. paper did (see their Figure 1). This would allow capturing the true variability due to physiological structure across all the data and not the variability related to the anatomical location of the electrode. I believe that adding this step would very strongly improve the remaining analysis (see more below).
Canonical CSD pattern. The choice of the first principal component of trial-averaged CSDs as the canonical CSD pattern is not at all obvious and has potential ramifications for the remaining analysis that it is based on. The first principal component is capturing the direction of inter-animal variability in the CSD pattern time-depth space affected by the angle of the probe, alignment of the recording sites array to cortical layers, exact retinotopic location, variation of brain states across individual recording sessions etc. This first eigenvector of the covariance matrix does not necessarily capture a common underlying pattern. In fact, mathematically it is not simply a weighted average but depends on the structure of the remaining unexplained variance (ca 50%). The components of the evoked response that strongly covary would contribute more to the first PC, but those that are uncorrelated or partially independent would contribute in an unpredictable way. In light of the focus of the paper on explaining this multi-component structure, it appears that starting from the decomposition of each animal's evoked response, extracting parameters (time, strength) of individual dipole components and characterizing their variability in population would be the most accurate way of characterizing experimental data.
CSD Distance Metrics. The Wasserstein distance used as a measure of distance between two CSD patterns is sensitive to misalignment of the CSD patterns (see above). The WD statistics cannot differentiate the contribution of the single sink and source components as they are aggregated into a single compound measure to which both anatomical variability contributes potentially more than the physiological variability. At the same time, the key elements of the paper are critically linked to this measure and its ability to capture the similarity between the experiment and the model.
KS-similarity. This statistic captures coarse (max of CDF difference) differences of the CDFs and is used to quantify the variability of firing rate magnitude across different populations. As a distance between two distributions, KS distance is not capturing the shape of the distributions and it's not quantifying where they are different. One can see that KS-similarity score is close to 0.75 for all populations in the dataset – the specificity of this measure is low. Importantly, this distribution similarity measure is not sensitive to the nature of change (direction), affected by the emergence of skew or modes in the distributions etc. Critically, the application of the measure to distinct layers doesn't capture joint spatial patterns of the firing rate across layers, which could get it closer to the CSD similarity measure. As the direction of the rate change is not distinguished, very different spatial patterns of firing rate change would be associated with similar KS-similarity values.
Pearson correlation is used to quantify the similarity in the average temporal profile of firing rates of layer-specific populations in experimental and model data. First, In contrast to the KS-similarity measures that strive to capture the distribution of the population firing rates, but ignores temporal dynamics, this measure ignores the variability of neural firing rates and their temporal dynamics. Changes in the shape of the distribution of firing rates will translate into unpredictable changes of the mean firing rate. Variable temporal profiles of neurons would further contribute to the poor interpretability of the mean firing rate evoked response. Given the generally log-normal distribution of the firing rates in many cortical circuits the mean is poorly describing the distribution. Overall, I can't see how mean rate response could be used as an interpretable statistics characterizing complex spatio-temporal dynamics across populations of neurons. A huge spread of confidence interval in all figures showing mean population PSTH is a testament to this point. Second, as the mean firing rate statistics is unlikely normally distributed, given the nonstationary changes of the firing rate, and might not be even unimodally distributed in the 0-100ms time window, using Pearson correlation coefficient is not ideal as it best measures the linear relationship between two normally distributed random variables. It is deviations from a bivariate normal distribution that would influence the value of correlation coefficient introducing a bias. It is hard to interpret what values ranging between 0.3 to 0.6 reflect if anything. Multiple spectral components of the evoked response would contribute in an unaccountable way to the measure further complicating the interpretation.
A small number of simulated trials yields very noisy (high variance) average CSD and PSTH. Given the importance for the interpretation of the results presented in the paper and the stochastic nature of the simulated data, it would be valuable to run longer simulations (more trials) to reach reliable low-variance estimates of the CSD/PSTH.
Analysis of the biophysical origin of the canonical CSD using the model is an interesting and worthy line of investigation. However, from the approach used in this paper that compares CSD due to one pathway to that generated by the full model, it is difficult to make conclusions about the exclusive contribution of that pathway, since this system with many recurrent connections is nonlinear. Figure 7 only shows single examples and relies on very specific, and somewhat arbitrary parameters that are critical for the outcome. Somato-dendritic distribution and density of synapses between a particular input and target neuron, the morphology of that neuron, and the strength of input all make a large impact on the resulting CSD. Most statements are based on a single example (i.e. specific parameters choice) and eye-balling, with no quantification provided. While it could provide some qualitative intuition for a plausible generative model of the CSD of the evoked response in V1, the large number of unknown parameters that were specifically chosen in the model makes it impossible to generalize the results of these single examples to the experimental data.
An appraisal of whether the authors achieved their aims, and whether the results support their conclusions.
Two-way dissociation between neuron's firing rate and CSD
The two-way dissociation between spikes and CSD is not quantitatively evaluated – substantial and minor changes should be demonstrated quantitatively. Claim: It is possible to alter the CSD pattern without substantially changing the firing rates by adjusting synaptic placement. One successful example is provided for this claim: putting all the excitatory inputs to layer four excitatory neurons onto their basal dendrites changes the contribution to CSD without changing the firing rate. Claim: it is possible to alter the firing rates without substantially changing the CSD pattern by adjusting synaptic weights. The example supporting this claim is not quantitatively verified and is not obvious that the change in the CSD pattern is minor. "LFP can reveal deficiencies in the model architecture not evident from firing rates. Models can be optimized for firing rates and CSD independently". Independent optimization is a stronger claim than what can be concluded from the results. The two-way dissociation between spiking and CSD is based on very few simple examples, and the dissociation is not perfect even then. The statistical metrics used to evaluate the change in the rate or CSD patterns are not adequate for such a strong conclusion.
Feedback can shape the sustained sinks and sources in superficial layers
Indeed, the final model simultaneously reproduces the evoked firing rate and CSD patterns in response to bright flash stimuli within animal variability, however, the metrics used cannot be directly interpreted. In essence, this important hypothesis is supported by a very indirect chain of results: experimental high-order firing rate data is fed into the model, which generated CSD/rate patterns, which are compared to the experimental CSD/rate patterns. In doing so the variability of the responses (CSD parameters of the superficial dipole in response to LM firing) provided by single trials and single animals is lost. The direct causal analysis would make this point more convincing.
Utility of CSD in constraining the model. The authors claim that CSD can be utilized to constrain the model. But even in this case, it serves as an after-the-fact confirmation of model performance by incorporating feedback from higher visual areas. Just using preexisting knowledge and experimental data to add feedback from HVA to the model and improving the CSD does not show the utility of CSD in tuning the model parameters. Especially since authors cut recurrent connections when investigating the contribution from external inputs to the CSD pattern it is unclear if feedback dipole cannot be mimicked by stronger excitatory/inhibitory synaptic input from the local circuit or slow delayed currents produced by active conductances. In other words, necessity and sufficiency cannot be demonstrated. Thus the approach chosen here could simply illustrate some ideas that might be further tested with proper analysis and experiments.
Many of the remarks above contained suggestions, few more technical ones are listed here.
An alternative to Wasserstein distance used between CSD patterns, is fitting a mixture model to sinks and sources and comparing the distribution of the parameters (location, time, half-width) of individual components in experimental data to that in the model. ICA decomposition, though now ideal, appears to be a good candidate. The notion of the canonical CSD pattern should rather be replaced by the statistical moments of the distributions of parameters of individual dipoles.
It would be more accurate to capture separately the magnitude of the rate increase and temporal aspects of the rate change (e.g. peak time, curvature of the rate increase curve) for each cell. Such measures, in addition, would replace clumsy aggregate firing rate population statistics computed (currently with KS-similarity) for different slices of time within the evoked response and would reflect spatio-temporal dynamics of the V1 circuit output.
KS-similarity could be replaced by more explicit measures capturing change in the distributions of the peak firing rate across populations. Temporal dynamics in the population can likewise be captured via characteristics of the distribution of the single neuron's peak latency. Statistical moments or direct characteristics of the PDF/CDF can be used for this quantification.
Overall both CSD and neural firing analysis marginalized by layer might not capture the key difference in the V1 circuit response. I believe that a multivariate analysis of the CSD and rate dynamics across layers that emphasizes relative changes of magnitude and latency of neural responses between the layers would provide a more principled measure. While firing rates might vary across experimental animals within each layer, between-layers relative changes in responses might be more consistent. Similar between-layer relative measures can be computed from the model data.
Feedback from higher visual areas. Experimental data provide a possibility for direct analysis of the origin of the superficial dipole: computing spike-CSD coherence (or spike-triggered average CSD) using spikes in higher-order visual areas and CSD in V1 allows for identifying which dipoles are consistent with the causal influence of the higher order visual cortex neural firing. Causal analysis of the firing rates in higher visual areas and CSD sink/source values on a single trial level would potentially allow for more precise conclusions. At any rate, even if correlational, this analysis would allow identifying which components of the CSD are correlated with the top-down inputs firing rate.
After adding feedback, the weights from the background node to feed-back targeted neurons are adjusted, but how is this justified? Please, clarify.
Diminished L4 sink in the final model. "The synaptic weights for thalamocortical connections were left unchanged from the original model" – and maybe that's why the early L4 sink is not prominent in the final model. It is striking that in most simulations the most prominent component in previously published and present evoked CSD profile, the L4 sink, is very weak in the simulation here. As suggested above, it appears that the model should reproduce the balance between various inputs' strengths to different layers manifested in the different CSD sink magnitudes. It would then be important to recapitulate this between-layer relative parameter of the input strength in the model (see comments on that above as well).
It is unclear why 90 minutes for 1 second of simulation presents such a barrier. When compared to experimental work this time seems realistic to provide much more statistically powerfull dataset using simulations. I would suggest increasing the number of trials of simulation.
Recommendations for improving the writing and presentation.
The authors compare the adjusted model to the original model but never mention what the original model was tuned to. "The synapses for recurrent connections were placed according to the following rules in the original model (Billeh et al. 2020)". As a reader, it would be helpful to have an idea of what the rules are based on. Is synaptic placement based on anatomy? Are synaptic weights fitted to evoked response to a particular visual stimulus?
Discussion should include a much more balanced interpretation of the results – it is merely an illustration of some ideas using the model, at the moment, rather than quantitative analysis demonstrating novel results on the composition of the CSD response, importance of the feedback, or recurrency in the actual experimental data. Likewise, the mismatch of the firing rate and CSD is simply illustrated by some examples, rather than revealing a novel result. In fact, the role of the interneurons in this mismatch is likely the key, but it was not discussed. I cannot imagine how a lack of realism in the interneuronal connectivity could be ignored (different PV morphological types, SOM, VIP) if one wants to recapitulate the circuit dynamics in response to the thalamic input. In other systems, for example, hippocampus, the mismatch between the CSD of theta and firing rates of different neuronal groups in entorhino-hippocampal circuits is much better understood and could be compared to here.
Discussion of other inputs not yet accounted for in the model that contribute to the experimental results is missing. What about other thalamic nuclei? Role of cortio-thalamic feedback?
It is clear that further improvements in the model realism will make it prohibitively complex. Authors postpone the future proper multiparametric optimization due to computational complexity. It would thus be useful to discuss the alternative approach – capturing parameters of the low dimensional dynamical system that would recapitulate the observed data or model
Reviewer #3 (Recommendations for the authors):
The authors present analysis of openly available Neuropixel data from the mouse primary visual cortex (V1), and several subcortical areas, provided by the Allen Institute. They seek to explain the detailed neuromechanistic underpinnings of flashed evoked responses in V1 using the previously established large scale V1 model by Billeh et al. 2020, also distributed by the Allen Institute. The authors establish a "canonical CSD pattern" occurring in V1 during the evoked response and work through a sequence of parameter manipulations needed to fit the model to the canonical CSD pattern, while simultaneously accounting for empirically observed firing rate activity in several populations of neurons. The main conclusions from the paper are that (1) feedback from higher order visual areas was necessary to account for both the early-latency (<100 ms) CSD and spiking features in the data, and (2) there is a two way disassociation between firing rates and CSD such that firing rates can be altered with minor effects on the CSD pattern by adjusting synaptic weights, and CSD can be altered with minor effects on firing rates by adjusting synaptic placement on the dendrites. The paper is well written with exciting new findings that are highly impactful to the field. It uses exceptional data and a phenomenally biophysically grounded and detailed large-scale model to make previously unattainable mechanistic interpretations of the data.
There are several directions that could be explored and clarified further that would increase the understanding and impact of the paper. Further, there is a history of earlier literature examining CSD/LFPs from visual evoked responses in VI in non-human primates that comes to similar conclusions on the necessity of "feedback" drive that should be discussed in comparison to the mouse results here.
1) While quantification of the activity of different cell types in different layers is necessarily limited in the empirical data, in the model it is not. The model contains more than 50,000 biophyscially detailed neuron models belonging to 17 different cell type classes: one inhibitory class (Htr3a) in layer 1, and four classes in each of the other layers (2/3, 4, 5, and 6) where one is excitatory and three are inhibitory (Pvalb, Sst, Htr3a) in each layer. In constraining the model to a limited number of features in the data (firing rates of limited cell populations and CSD), the model also makes many more novel and testable predictions about the activity of other cell types and synaptic connections (not fit too) that would be valuable for the readers to see.
For example, while the data is sparse in terms of knowledge of different inhibitory neuron types across populations, the model is not, and can provide testable predictions on the firing activity of the inhibitory neurons in each layer and across three different cell types Pvalb, Sst, Htr3a, along with Htr3a neurons in layer 1. In the final fit model as in Figure 7, it would be useful to detail model predictions about the firing rates in the 4 different interneuron types, and their strength of influence on the RS neurons. These predictions could then be further tested in follow up studies, and perhaps there is a prediction that was not used for fitting but that can be tested in the data set here?
Along this line, what potential mechanisms govern the animals or trials that present atypical CSD or spike rate evoked responses? The methods and results of this study are uniquely positioned to not only account for the circuit mechanisms of a canonical response, but also the cell and circuit mechanisms of a distribution visual flash evoked responses, particularly those categorized as non-canonical. How varied are the non-canonical responses? A supplemental figure showing the CSD from each of the animals, e.g. as shown for other features in Figure S4, would also be helpful to understand the full temporal spatial variability in the CSD. The present study has a distinctly powerful opportunity to make predictions regarding the underlying cell types, circuit elements, connectivity patterns, etc. that control for visual evoked response variability.
Related to this, there are a few important assumptions in the model that should be clarified further:
a. Are the LM to V1 connections known to not target L1 inhibitory neurons? What happens if they do in the simulations?
b. Do the "feedback" inputs in the model also contact the distal dendrites of the L6 RS cells? What about the apical dendrites of the L4 RS cells? What about the L1 cells? Is there any literature supporting the choice to connect the feedback to only a subset of the neurons in the supragranular layers or was it purely because it was necessary to reduce the magnitude of the L5/L6 dipole?
2) There are important details in the comparison made between the model neuron firing rates and empirically recording firing rates that need to be understood more fully. In the empirical recordings, cells are sorted into RS or FS populations, and RS are presumed excitatory with FS presumed inhibitory. The activity of the RS is separated by layer, but the FS is grouped across layers because of the sparsity of recordings. When comparing the model firing rates to the experimental firing rates, the excitatory and "non-Pvalb populations" were grouped together in each layer of the model to make up the RS cells in L2/3, L4, L5, and L6, while the Pvalb cells across all layers were grouped together to make up the FS cells of V1. This makes sense for the experimental data, however, the model has distinct non-Pvalb inhibitory populations including SST and Htr3a. Were these non-Pvalb inhibitory cells grouped with the FS cells or RS cells in the model analysis? Clarification is needed.
Line 557 states that to bring the firing rates of the neurons in the original model closer to the experimental data ["We first reduced the synaptic weights from all excitatory populations to the fast-spiking PV- neurons by 30% to bring their firing rates closer to the average firing rate in this population in the experiments. This resulted in increased firing rates in all other (RS) populations due to the reduced activity of the inhibitory Pvalb-neurons (Figure S6). Therefore, we further applied reductions in the synaptic weights from all excitatory neurons to RS neurons or increases in the synaptic weights from inhibitory neurons to the RS neurons to bring their firing rates closer to the experimental average firing rates." ] The "or" statement here is important, should it be "and"? Was it all inhibitory neurons or only the Pvalb-neurons?
A table with the final parameter values would be immensely helpful in interpreting the microcircuit details contributing to the final model fit in Figure 7.
3) The process of parameter tuning in models, particularly in large-scale models, is challenging. Throughout the results, the authors walk through a series of manipulations to the original model that are necessary to account for the empirical data. While I understand that not all parameters can be estimated, and there is a nice paragraph in the discussion on this starting on Line 981, it would be helpful to have a methods section that describes the overall strategy for parameter manipulations and rationale for the choice of the final parameters. Was it all hand tuning, or was there some estimation method used to fit to data? More specifically, there is a section called "Adjusting synaptic weights' that gives a high level description of a range of multiplication factors over which previous model parameters were adjusted to fit the current data. More details on how this range was determined, and if one value in this range or multiple values were ultimately used in the simulations, and why, is needed.
Details on how the new feedback was implemented in the model is described on Lines 682-690, but I didn't see any description of these inputs in the method section. The reference to Figure S3 seems to be a mistake. I suggest double checking all Figure references. It seems the "Adjusting synaptic weights" section only describes local interactions. Were there feedback inputs in the original Billeh model that were adjusted here? Or are they completely new, and if so, what were the starting weights and how were they chosen?
The modeling results here build from the prior work in Billeh et al. 2020. In this paper, the model is applied to visual flash evoked responses, while in Billeh et al. 2020 the model was tuned to orientation and direction tuning to drifting grating stimuli. It is exciting to see that the updated model reproduces consistent results for firing rates and rate based direction/orientation selectivity, as shown in Supplementary Figure S9. It is stated that for these simulations, the stimulation of the network was as defined in Billeh et al. 2020. One important piece of information to understand is why feedback is needed in the current simulations of flash evoked response, but isn't required for the drifting grading simulations? Would including feedback improve the results? A discussion and/or table summarizing the differences in the two models would be useful to understand what further manipulations may be necessary to accurately account for the LFP/CSD, in addition to firing rates, in Billeh et al. 2020 simulations.
4) One strategy the authors used to examine the importance of various neural mechanisms to the CSD, is to remove various model elements such as synaptic connections or intrinsic ionic currents and see how much the CSD changes. For example, there is a conclusion line 803 that [without inhibition that fact that] "the major sinks and sources are still present is an indication that the currents from excitatory input account for the majority of the sinks and sources observed in the experimental CSD." This could be due to the fact that in the fit model (Figure 7) the inhibitory synapses are very weak (are they?) and hence not playing a major role in the established model configuration. The conclusion that inhibition plays a small role is counterintuitive to the growing literature showing the importance of inhibition to LFP generation (e.g. Telenczuk et al. Scientific Reports 2017), and hence CSD, and further discussion and understanding of the role of the different types of inhibition in this model would be helpful.
A more useful analysis to quantify the impact of specific network features would be a parameter sensitivity analysis that shows the influence of changes in parameters to variance in the resultant CSD; this isn't feasible for all parameters but could perhaps be done for the parameters/network features concluded to account (or not account) for the main feature of the canonical CSD, e.g. as summarized beginning on line 827. Or at least some simulations that show what happens over a range of parameter values for the important parameters.
There is a reference in this paragraph to Figure 7E, which should be Figure 7F. The authors need to double check all Figure references.
5) A similar large scale biophysically detailed cortical column modeling paper (Reimann et al. Neuron 2013) is relevant to discuss further in light of the current results. In particular, that paper showed that "active currents" dominated the LFP, rather than synaptic currents. There is a statement line 890 that "In line with observations made by Reimann et al. (2013), we found that the somatic voltage- dependent membrane currents significantly shape the CSD signature (Figure 5H and 7B)". To my knowledge, Riemann considered the contribution of all active currents, including those in the dendrites, while in this paper they have only explored removing them in the soma. Those in the dendrites may have a greater impact, or maybe not (?). Can the authors show this? Further examination/discussion would be useful.
5) Of particular relevance to these results are earlier laminar recordings to flash evoked responses in VI in non-human primates that suggest the presence of "feedback" and discusses interaction between higher and lower order areas that may contribute to this "feedback": doi: 10.1016/0042-6989(94)90156-2 https://doi.org/10.1093/cercor/8.7.575. A discussion on the conclusions in this paper related to those is needed.
These papers, and many other papers, also show a decoupling between LFP/CSD and firing (MUA). Clarification should be made that such decoupling is not a new concept, but to my knowledge it has yet to be shown how it could be mechanistically implemented using a computational neural model, as the authors have nicely done. Nor has it been shown how simulated LFP/CSD can be altered with only minor effects on the firing rates, as shown here.
Also relevant, is a more recent study of flash evoked responses in humans and non-human primates doi: 10.1126/sciadv.abb0977, which appears to show consistent results to those reported here.
6) There is a paragraph in the methods on line 1644 that begins "This application of the WD to compare CSD patterns comes with certain considerations that are important to note…". I suggest moving this paragraph to the Results section.
https://doi.org/10.7554/eLife.87169.sa1Author response
[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]
Comments to the Authors:
We are sorry to say that, after consultation with the reviewers, we have decided that this work will not be considered further for publication by Life. As you see there are a great number of technical concerns raised by the reviewers (in particular reviewer 2) which we feel, at present, preclude publication. However, we are happy to receive a resubmitted version of a strongly improved manuscript in which it becomes clear that you can convincingly address these technical concerns. In this we recommend writing an action plan of how you would address these concerns in advance.
Reviewer #1 (Recommendations for the authors):
While the dissection of CSD components is detailed and exhaustive, the analysis of unit firing response falls short. In previous works (like Senzai et al., 2018) a stereotypical timing of single unit firing to this type of stimuli is described in addition to the CSD responses. In the paper here can the author's model reproduce this timing of unit activation? In addition to reporting average firing rates, the authors should also describe and interpret the temporal dynamics of unit firing and relate it to their mechanistic description of CSD components in the last section.
The question of how to best characterize spike statistics is an important one, and we thank the reviewer for raising it. We acknowledge that there are relevant aspects of spike statistics that were not captured with the metrics we originally utilized. Reviewer 2 raised similar points, so per the suggestions of reviewers 1 and 2, we have supplemented our existing metrics with the following analyses:
Comparison of the moments of distributions of peak firing rates across cells in the evoked response
Comparison of the moments of distributions of latencies to peak firing rate across cells.
Comparing moments of the relative magnitude of firing rates across laminar populations
Comparing moments of distributions of maximum curvature of firing rate
Summary of the results from (i):
For the first moment (mean) (Figure 6 —figure supplement 2A), the original model was outside the experimental variability for all populations except L5 RS. The final model was just outside of the experimental variability for the first moment of the peak firing rate of L2/3 RS (regular-spiking) cells, but inside for all other populations. The intermediate model was within experimental variability for all populations.
For the second moment (variance) (Figure 6 —figure supplement 2B), all model versions were within experimental variability for all populations.
For the third moment (skewness) (Figure 6 —figure supplement 2C), all three model versions were outside of the experimental variability for L2/3 RS cells, and the original model was outside the experimental variability for L5 RS and FS cells as well. For all other populations, all models were within experimental variability.
For the fourth moment (kurtosis) (Figure 6 —figure supplement 2D), all model versions were within experimental variability for all populations.
Summary of results from (ii):
All model versions were within experimental variability for all populations and all moments except the 4th moment for FS cells, where all model versions were outside the experimental variability (Figure 6 —figure supplement 3)
Summary of results from (iii) and (iv):
For (iii) and (iv), all three model versions were within the experimental variability on all measures (Figure 6 —figure supplement 4 and Figure 6 —figure supplement 5).
Some of the analysis and simulations can be better motivated. For example, it is not clear, in the section starting in line 244, why it is important to quantify CSD variability and how that relates to the broader goals of the study. The manuscript in quite dense and in occasions is not easy to follow how certain results contribute to the main objective. It would help to state these explicitly at the outset and related each subsequent section to them.
The numbers we obtain when quantifying similarity (or distance) between the model and experimental targets cannot necessarily be easily interpreted on their own, so we quantify CSD (and spike) variability in experiments to use as references in our assessment of how similar the signals produced by the model are to those recorded in experiments. Calculating Wasserstein distances from CSD of individual mice to the PC 1 canonical CSD target, for example, gives us an idea of how big of a Wasserstein distance we can expect in modeled CSD that reproduces experimental observations. Calculating pairwise Wasserstein distances between animals provides a reference to use to judge what constitutes a large and a small change in CSD patterns from one model version to the next, which we did in Figure 5 —figure supplement 1. We added a short description of this motivation on lines 247-248 and 351-353.
When discussing the biophysical mechanisms of CSD generation a useful concept to introduce is that of dipolar moment. Components with a larger dipolar moment (more separation between sink and source) generate larger LFPs.
We added the following brief description of this mechanism in the context of adjusting the position of excitatory synapses onto the excitatory L6 cells on lines 702-704:
“Greater separation between a sink and a source in a dipole moment increases its magnitude, so restricting the range within which synapses can be placed should diminish the L5/L6 dipole’s dominance.”
Reviewer #2 (Recommendations for the authors):
1. Spiking patterns alone, they argue, do not uniquely characterize a neural circuit, i.e. there are many models (or variations to the same model) that are capable of reproducing the spiking activity of data, so requiring the model to simultaneously reproduce the CSD would further restrict the model.
We do argue that requiring the model to match CSD in addition to spiking patterns would further constrain the model, that is correct; however, we want to stress that a model that reproduces both spiking activity and CSD is also unlikely to be uniquely characterized, even if it is more restricted than a model that only reproduces spiking activity. Gutenkunst et al., 2007, among others, have demonstrated that many different parametrizations can produce models that fit the data equally well in systems biology (so-called parameter “sloppiness”). This has also been shown for multi-compartment models in neuroscience; e.g. Hay et al., 2011, found that more than 500 parametrizations of their biophysically detailed layer 5 pyramidal neuron could perform the same function equally well.
2. The key novelty of this manuscript compared to a previous paper (Billeh et al. 2020) is adding CSD patterns as an additional constraint to more accurately infer the model parameters beyond its prior state. While it is true that adding constraints would undoubtedly be helpful, there is no guarantee that specific hyperparameters used, components of the model simplified or ignored and gross aggregated measures chosen for comparing the patterns in the data and produced by the model would allow identifiability of the model and use it to explain the experimental data.
We agree that identifiability of a true model, in the sense of a model where all parameters correspond to true experimental values, cannot be guaranteed. Nor can it be expected. Limitations on the data presently available as well as on computational resources necessitates simplifications that restrict the realism of the model. However, this is a feature of all modelling efforts in neuroscience. We would argue that despite these caveats, models, including the one presented here, can still give valuable insight into the mechanistic underpinnings of experimental observations.
3. Authors make no justification for many specific choices for a number of unknown high-dimensional parameters e.g. morphology of particular cell types (e.g. L4 pyramids have minimal basal dendrites, Scala et al. 2019, which has a strong effect on Figure 5 results), somato-dendritic distribution of synaptic inputs from different afferents, choice of interneurons (morphology, intrinsic properties, and connectivity) to include and how to wire them in recurrent and feedforward circuit, etc.
These choices, with the exception of the choice to use L4 pyramids for the results in Figure 5, were not made in this study but in Billeh et al., 2020 and Gouwens et al., 2018, on which this investigation builds. The work reported in those papers was extensive and meticulous, relying as much as possible on systematically obtained experimental data. For example, morphologies of particular cell types were taken from the Allen Cell Types Database, selecting such morphologies from all the cells for which high-quality models were available. The choices that are most relevant to the interpretation of our findings in this paper and their justification are relayed in the Methods section. The choice to use L4 pyramidal cells in Figure 5 was somewhat arbitrary. We have therefore carried out the same simulations and analysis on L2/3 and L5 excitatory cells for comparison and added plots of the results in the supplementary material (Figure 5 —figure supplement 3 and Figure 5 —figure supplement 4).
4. Given the nonlinear nature of the network dynamics driven by the thalamic input and its nontrivial dependence on very high-dimensional parameter space, I see no way to identify a true model in that space based on very coarse metrics of aggregated neural activity. Care for realism in the model and the use of real data provides a false impression that the outcome of the effort should be the identification of the actual unobserved parameters of the system and the identification of the model. I think that it is a false impression, that might be misleading the readers.
In the initial manuscript, we noted how the nonlinear nature of network dynamics limits the precision with which we can estimate the contributions of different pathways to the specific current sinks and sources, and hence, the parameters involved in their generation (see its lines 962-967 or lines 1042-1047 in the one presently submitted). Thus, this is a view we share with the reviewer. We have supplemented our existing metrics with other measures for comparing the neural activity of the model to experiments. These are described in further detail in points 6, 7, and 9 below.
However, even with these considerations, we would argue that our model can still shed light on which parameters are involved in various experimentally observed phenomena, and as more data becomes available in the future, the precision in our estimation of these parameters will improve.
5. The choice of metrics and dimensionality reduction approaches used to capture the group results on the CSD and firing rate is not suitable for validating the model or might provide poorly interpretable results. The authors note in the discussion, that there might be more suitable metrics and promise to work on it in the future. However, if the choice of metric used in the present paper makes interpretation questionable it puts in question the validity of all the results and conclusions in the manuscript, except for the actual data and model it is based on. Thus, I find it critical that the analysis methods used are chosen carefully and shortcomings discussed below are taken into account.
This is addressed in points 6, 7, and 9 below.
6. Inter-animal CSD alignment. It is obvious from the examples and prior work (Senzai et al. 2019) that variation of the location and scale of many CSD dipoles varies greatly. Any further analysis that relies on the CSD has to be thus rescaled and realigned, as Senzai et al. paper did (see their Figure 1). This would allow capturing the true variability due to physiological structure across all the data and not the variability related to the anatomical location of the electrode. I believe that adding this step would very strongly improve the remaining analysis (see more below).
The reviewer correctly points out that the CSD of each animal should be aligned and scaled to be comparable across animals and to ensure that our analysis does not merely reflect variability unrelated to the physiological differences between animals. In our initial manuscript, we did align and scale the CSD of all animals. However, in contrast to the approach in Senzai et al., 2019, we aligned the probes based on the depth of the electrode obtained from histology and provided in Allen Common Coordinate Framework (CCF) coordinates rather than landmarks in the data.
In this revised manuscript, we have implemented an alignment of the CSD based on landmarks in addition to our original alignment based on histology, and we supply the results from a quantitative comparison of the model CSD to the experimental CSD after landmark alignment (Figure 6 —figure supplement 7). Our conclusions regarding which model versions reproduce the experimental CSD patterns and which do not are unchanged by using the landmark alignment approach.
7. Canonical CSD pattern. The choice of the first principal component of trial averaged CSDs as the canonical CSD pattern is not at all obvious and has potential ramifications for the remaining analysis that it is based on. The first principal component is capturing the direction of inter-animal variability in the CSD pattern time-depth space affected by the angle of the probe, alignment of the recording sites array to cortical layers, exact retinotopic location, variation of brain states across individual recording sessions etc. This first eigenvector of the covariance matrix does not necessarily capture a common underlying pattern. In fact, mathematically it is not simply a weighted average but depends on the structure of the remaining unexplained variance (ca 50%). The components of the evoked response that strongly covary would contribute more to the first PC, but those that are uncorrelated or partially independent would contribute in an unpredictable way.
The reviewer is correct that, mathematically, the first principal component is strictly speaking not a weighted average. We have updated the relevant parts of the manuscript such that it reads “weighted contributions” rather than “weighted average”.
The question of how probe alignment could affect this analysis is addressed in point 6 above.
We have added CSD plots for all animals (Figure 2 —figure supplement 1), a plot of the cumulative explained variance of all principal components (Figure 2 —figure supplement 2A), as well CSD plots for the 10 first principal components that cumulatively explain 90 % of the variance (Figure 2 —figure supplement 2B) in the supplementary material. This way, an assessment of whether the first principal component captures a common underlying CSD pattern is available to the reader.
We have also added the analysis where we used a plain average, rather than the first principal component, of trial averaged CSD across animals as the target in our quantitative analysis of model performance. We did this both where the animal CSDs had been aligned to histology, as it was originally, and where they had been aligned to landmarks in the data. In both cases, the model reproduced the experimentally observed CSD pattern and supported our conclusions from the analysis where we used the first principal component of the histology aligned CSD as the target (Figure 6 —figure supplement 7C and Figure 6 —figure supplement 9).
8. The Wasserstein distance used as a measure of distance between two CSD patterns is sensitive to misalignment of the CSD patterns (see above). The WD statistics cannot differentiate the contribution of the single sink and source components as they are aggregated into a single compound measure to which both anatomical variability contributes potentially more than the physiological variability. At the same time, the key elements of the paper are critically linked to this measure and its ability to capture the similarity between the experiment and the model. An alternative to Wasserstein distance used between CSD patterns, is fitting a mixture model to sinks and sources and comparing the distribution of the parameters (location, time, half-width) of individual components in experimental data to that in the model. ICA decomposition, though now ideal, appears to be a good candidate. The notion of the canonical CSD pattern should rather be replaced by the statistical moments of the distributions of parameters of individual dipoles.
It is true that ICA has been successfully applied to distinguish contributions from distinct pathways to recorded LFP in several cases (e. g. Makarov et al., 2010, Makarov et al., 2011, Fernandez-Ruiz et al., 2012, Fernandez-Ruiz et al., 2017). However, in these examples, ICA was applied to LFP recorded in subcortical structures such as hippocampus, and the specific conditions that allow for its successful application in those structures may not hold for other areas of the brain, such as V1. In Makarov et al., 2011, the factors that reduce ICA’s performance are listed and include high synchronicity and a strong imbalance in the relative power of the LFP generated by different biophysical processes. If the various generators of LFP are synchronously activated or their contributions are imbalanced, significant cross-contamination of the LFP attributed to them with ICA may occur. In our data, there is substantial temporal overlap between the biophysical processes during evoked flash response in V1, which would likely produce such cross-contamination. Whether or not there are also strong imbalances in relative power of LFP generated by different biophysical processes is not known for sure and would constitute another potential source of cross-contamination when applying ICA.
Furthermore, the underlying assumption of independence of LFP generators is unlikely to hold true during flash-evoked responses in V1. LGN contributes to all layers in mouse V1, even if L4 is its primary target. This means that there may be significant contributions from LGN to several different sinks in V1. In fact, our own analysis suggests that LGN, together with recurrent connections, is an important contributor to both the L4 and the L5/L6 sink. Thus, the biophysical processes generating the L4 sink and the L5/L6 sink are not independent. Additionally, even if the direct contributions of the LGN input to the L5/L6 sink should be negligible, the recurrent contribution to the L5/L6 sink is dependent on prior LGN input to L4, so the recurrent generation of the L5/L6 sink would in any case not be a pathway that is independent of the LGN input pathway. The same is likely to be true for the input from HVAs such as LM, as they receive input from V1 (Harris et al., 2019, Siegle et al., 2020). Therefore, one cannot expect ICA to successfully separate the different sink/source dipoles into different components.
Note also that Glabska et al., 2014 tested the application of ICA on LFP generated generated by the Traub model (Traub et al., 2005) of somatosensory cortex, and found that no single component could recover the activity of a single population. They could only recover the population activity through combining components, and, even then, only two populations could reliably be recovered, at most three. They also point out that determining which components to combine depended on knowledge of the ground truth population activity in the model, something one does not have access to in the experiments.
We still tested the application of ICA on our experimental and model data. Example results from two animals and the final model are shown in Author response image 1. In accordance with expectations from Makarov et al. 2011, there appears to be significant cross-contamination between the components both for the animals and the model, as similar patterns recur in several components (see e. g. IC 1 and 2 in Author response image 1B), and one component picks up most of the variance (see e. g. IC 2 for the second animal and IC 3 for the final model). It is also not clear which component should be considered to correspond to which biophysical process, nor which components should be compared across animals and between animals and the model. Utilizing all trials instead of trial averaged data or a different number of components did not significantly change the results.

ICA decomposition of experimental CSD and simulated CSD.
(A) Two examples of applying ICA on trial averaged CSD of animals. Four resulting independent components are displayed from left to right, and the trial averaged CSD of the example animal is shown in the rightmost plot. (B) Four independent components resulting from applying ICA to trial averaged CSD from the simulation on the final model version.
Lastly, we want to note that we are aware that in Senzai, Fernandez-Ruiz, and Buzsáki, 2019, ICA was successfully applied to decompose LFP recorded in V1. However, in that study, ICA was applied after filtering in the γ band (30-100 Hz), as its purpose was to identify layer-specific γ oscillations. Filtering in the γ band first is not appropriate in our study, since we are equally interested in signals outside the γ range. Thus, even though ICA was applicable to V1 data in Senzai, Fernandez-Ruiz, and Buzsáki, 2019, that does not imply a suitable decomposition of LFP/CSD in our study.
The Wasserstein distance metric is not without limitations, which we note in the manuscript. However, it does not require any human characterization of which sinks and sources should be paired and compared across animals and between animals and the model. We wanted to employ a metric that could compare whole CSD patterns; that is, the relative strength, shape, timing and position of sinks and sources without such, often subjective, pairing. Thus, a compound measure like the Wasserstein distance was in our view the most appropriate.
9a. KS-similarity. This statistic captures coarse (max of CDF difference) differences of the CDFs and is used to quantify the variability of firing rate magnitude across different populations. As a distance between two distributions, KS distance is not capturing the shape of the distributions and it's not quantifying where they are different. One can see that KS-similarity score is close to 0.75 for all populations in the dataset – the specificity of this measure is low. Importantly, this distribution similarity measure is not sensitive to the nature of change (direction), affected by the emergence of skew or modes in the distributions etc. Critically, the application of the measure to distinct layers doesn't capture joint spatial patterns of the firing rate across layers, which could get it closer to the CSD similarity measure. As the direction of the rate change is not distinguished, very different spatial patterns of firing rate change would be associated with similar KS-similarity values.
b. Pearson correlation is used to quantify the similarity in the average temporal profile of firing rates of layer-specific populations in experimental and model data. First, In contrast to the KS-similarity measures that strive to capture the distribution of the population firing rates, but ignores temporal dynamics, this measure ignores the variability of neural firing rates and their temporal dynamics. Changes in the shape of the distribution of firing rates will translate into unpredictable changes of the mean firing rate. Variable temporal profiles of neurons would further contribute to the poor interpretability of the mean firing rate evoked response. Given the generally log-normal distribution of the firing rates in many cortical circuits the mean is poorly describing the distribution. Overall, I can't see how mean rate response could be used as an interpretable statistics characterizing complex spatio-temporal dynamics across populations of neurons. A huge spread of confidence interval in all figures showing mean population PSTH is a testament to this point. Second, as the mean firing rate statistics is unlikely normally distributed, given the nonstationary changes of the firing rate, and might not be even unimodally distributed in the 0-100ms time window, using Pearson correlation coefficient is not ideal as it best measures the linear relationship between two normally distributed random variables. It is deviations from a bivariate normal distribution that would influence the value of correlation coefficient introducing a bias. It is hard to interpret what values ranging between 0.3 to 0.6 reflect if anything. Multiple spectral components of the evoked response would contribute in an unaccountable way to the measure further complicating the interpretation.
c. KS-similarity could be replaced by more explicit measures capturing change in the distributions of the peak firing rate across populations. Temporal dynamics in the population can likewise be captured via characteristics of the distribution of the single neuron's peak latency. Statistical moments or direct characteristics of the PDF/CDF can be used for this quantification.
d. It would be more accurate to capture separately the magnitude of the rate increase and temporal aspects of the rate change (e.g. peak time, curvature of the rate increase curve) for each cell. Such measures, in addition, would replace clumsy aggregate firing rate population statistics computed (currently with KS-similarity) for different slices of time within the evoked response and would reflect spatio-temporal dynamics of the V1 circuit output.
e. Overall both CSD and neural firing analysis marginalized by layer might not capture the key difference in the V1 circuit response. I believe that a multivariate analysis of the CSD and rate dynamics across layers that emphasizes relative changes of magnitude and latency of neural responses between the layers would provide a more principled measure. While firing rates might vary across experimental animals within each layer, between-layers relative changes in responses might be more consistent. Similar between-layer relative measures can be computed from the model data.
The question of how to best characterize spike statistics is an important one, and we thank the reviewer for raising it. We acknowledge that there are relevant aspects of spike statistics that were not captured with the metrics we originally utilized. Reviewer 1 raised similar points, so per the suggestions of reviewers 1 and 2, we have supplemented our existing metrics with the following analyses:
Comparison of the moments of distributions of peak firing rates across cells in the evoked response
Comparison of the moments of distributions of latencies to peak firing rate across cells.
Comparing moments of the relative magnitude of firing rates across laminar populations
Comparing moments of distributions of maximum curvature of firing rate
Summary of the results from (i):
For the first moment (mean) (Figure 6 —figure supplement 2A), the original model was outside the experimental variability for all populations except L5 RS. The final model was just outside of the experimental variability for the first moment of the peak firing rate of L2/3 RS (regular-spiking) cells, but inside for all other populations. The intermediate model was within experimental variability for all populations.
For the second moment (variance) (Figure 6 —figure supplement 2B), all model versions were within experimental variability for all populations.
For the third moment (skewness) (Figure 6 —figure supplement 2C), all three model versions were outside of the experimental variability for L2/3 RS cells, and the original model was just outside the experimental variability for L5 RS and FS cells as well. For all other populations, all models were within experimental variability.
For the fourth moment (kurtosis) (Figure 6 —figure supplement 2D), all model versions were within experimental variability for all populations.
Summary of results from (ii):
All model versions were within experimental variability for all populations and all moments except the 4th moment for FS cells, where all model versions were outside the experimental variability (Figure 6 —figure supplement 3)
Summary of results from (iii) and (iv):
For (iii) and (iv), all three model versions were within the experimental variability on all measures (Figure 6 —figure supplement 4 and Figure 6 —figure supplement 5).
10. A small number of simulated trials yields very noisy (high variance) average CSD and PSTH. Given the importance for the interpretation of the results presented in the paper and the stochastic nature of the simulated data, it would be valuable to run longer simulations (more trials) to reach reliable low-variance estimates of the CSD/PSTH.
Since we use experimentally recorded spike trains as input to the model, the number of trials we can run in a simulation is limited by the number of experimentally recorded spike trains we have in the structures providing input to primary visual cortex. If we increase the number of trials, we have fewer spike trains to use to create unique inputs in each trial. And the fewer spike trains we use to create input for each trial, the more variable the input in each trial will be. Thus, increasing the number of trials will not necessarily reduce the variability in the simulations. We have experimented with different trial numbers and found that our current number gives a reasonable trade-off between the variability due to limited spike trains and variability due to a limited number of trials.
11. Analysis of the biophysical origin of the canonical CSD using the model is an interesting and worthy line of investigation. However, from the approach used in this paper that compares CSD due to one pathway to that generated by the full model, it is difficult to make conclusions about the exclusive contribution of that pathway, since this system with many recurrent connections is nonlinear. Figure 7 only shows single examples and relies on very specific, and somewhat arbitrary parameters that are critical for the outcome. Somato-dendritic distribution and density of synapses between a particular input and target neuron, the morphology of that neuron, and the strength of input all make a large impact on the resulting CSD. Most statements are based on a single example (i.e. specific parameters choice) and eye-balling, with no quantification provided. While it could provide some qualitative intuition for a plausible generative model of the CSD of the evoked response in V1, the large number of unknown parameters that were specifically chosen in the model makes it impossible to generalize the results of these single examples to the experimental data.
We agree that it’s difficult to make conclusions about the exclusive contributions of specific pathways due to the nonlinear effects of recurrent connections (see lines 1042-1047). Indeed, this difficulty was among the reasons why we limited our analysis to one that could provide a qualitative intuition for a plausible generative model of the evoked CSD response in V1 in the first place. Attempting to give precise and generalizable numbers for the contributions of specific pathways lies beyond what is possible with the model at this stage. The abovementioned limitation, as well as other simplifications of the model (e. g., that only the feedback from LM was modeled, while V1 also receives feedback from other higher visual areas), and limitations in the data presently available prohibits exact enumeration of the exclusive contribution from specific pathways to specific sinks and sources. Intuition for a plausible generative model of CSD is, in our view, a more appropriate representation of what can be interpreted from our results at this time.
12. The two-way dissociation between spikes and CSD is not quantitatively evaluated – substantial and minor changes should be demonstrated quantitatively. Claim: It is possible to alter the CSD pattern without substantially changing the firing rates by adjusting synaptic placement. One successful example is provided for this claim: putting all the excitatory inputs to layer four excitatory neurons onto their basal dendrites changes the contribution to CSD without changing the firing rate. Claim: it is possible to alter the firing rates without substantially changing the CSD pattern by adjusting synaptic weights. The example supporting this claim is not quantitatively verified and is not obvious that the change in the CSD pattern is minor. "LFP can reveal deficiencies in the model architecture not evident from firing rates. Models can be optimized for firing rates and CSD independently". Independent optimization is a stronger claim than what can be concluded from the results. The two-way dissociation between spiking and CSD is based on very few simple examples, and the dissociation is not perfect even then. The statistical metrics used to evaluate the change in the rate or CSD patterns are not adequate for such a strong conclusion.
We have now added quantitative evaluations of both the change in spiking with changes in synaptic placement alone and changes in CSD with changes in synaptic weights (Figure 5 – —figure supplement 2 and Figure 5 —figure supplement 1, respectively), and these evaluations support our original conclusions. The added quantification strengthens our analysis, and we thank the reviewer for making this suggestion. Furthermore, we have provided more examples of the effects of manipulating synaptic placement in Figure 5 —figure supplement 3 and Figure 5 —figure supplement 4. We want to note that our manuscript did not state that our results demonstrated a perfect dissociation and that our comment that models can be optimized for firing rates and CSD independently was preceded by the qualifier “to a certain extent”. We do not hold the view that CSD and spikes can be completely dissociated and be utilized entirely independently in optimization, as this standpoint can neither be supported by our own data and analysis nor previous work on this topic.
13. Feedback can shape the sustained sinks and sources in superficial layers
Indeed, the final model simultaneously reproduces the evoked firing rate and CSD patterns in response to bright flash stimuli within animal variability, however, the metrics used cannot be directly interpreted. In essence, this important hypothesis is supported by a very indirect chain of results: experimental high-order firing rate data is fed into the model, which generated CSD/rate patterns, which are compared to the experimental CSD/rate patterns. In doing so the variability of the responses (CSD parameters of the superficial dipole in response to LM firing) provided by single trials and single animals is lost. The direct causal analysis would make this point more convincing.
This is an interesting idea, and we will consider it for future extensions of our work.
14.Utility of CSD in constraining the model. The authors claim that CSD can be utilized to constrain the model. But even in this case, it serves as an after-the-fact confirmation of model performance by incorporating feedback from higher visual areas. Just using preexisting knowledge and experimental data to add feedback from HVA to the model and improving the CSD does not show the utility of CSD in tuning the model parameters. Especially since authors cut recurrent connections when investigating the contribution from external inputs to the CSD pattern it is unclear if feedback dipole cannot be mimicked by stronger excitatory/inhibitory synaptic input from the local circuit or slow delayed currents produced by active conductances. In other words, necessity and sufficiency cannot be demonstrated. Thus the approach chosen here could simply illustrate some ideas that might be further tested with proper analysis and experiments.
We cannot, at this point, definitively rule out that the feedback from HVAs does not contribute to the late dipole in upper layers and that the origin of this dipole is predominantly local (stemming from activity driven by recurrent connections within V1). However, a few observations make this less plausible:
The late dipole is roughly 20 ms delayed relative to the onset of the initial response (which occurs about 40 ms after stimulus presentation). That’s a substantially longer delay than most local circuit mechanisms can explain. For example, the conduction delay of monosynaptic connections between granular and supragranular processes would not be longer than a few ms, and we can observe transient sinks in supragranular layers already at about 45 ms for many animals (Figure 2 —figure supplement 1), which we believe is more likely to reflect the local input from L4 to L2/3.
The sustained L4 source could in principle originate in input from inhibitory populations. However, the firing rate of FS cells is decreasing in the time window after 60 ms (Figure 3A), which does not support the possibility of strong inhibition in this time window.
All attempts to mimic the late dipole with changes in the excitatory/inhibitory synaptic input were unsuccessful, whereas providing feedback inputs from HVA readily helped to produce the late dipole resembling the one seen in experiments.
Of course, this does not mean that local circuit explanations of the late dipole are completely ruled out, but they do make them less parsimonious than the explanation by feedback from HVAs.
15. After adding feedback, the weights from the background node to feed-back targeted neurons are adjusted, but how is this justified? Please, clarify.
In the original model, the input from the background node represented the influence of the rest of the brain on V1, which included the influence from higher visual areas such as LM. This means that some of the feedback influence from LM on V1 should be present (though coarsely represented) in the input from the background node. When the influence of input from LM to feedback targeted cells is modelled on its own, the influence of the background node must be updated accordingly. This justification has been incorporated in the manuscript on lines 1703-1708.
16. Diminished L4 sink in the final model. "The synaptic weights for thalamocortical connections were left unchanged from the original model" – and maybe that's why the early L4 sink is not prominent in the final model. It is striking that in most simulations the most prominent component in previously published and present evoked CSD profile, the L4 sink, is very weak in the simulation here. As suggested above, it appears that the model should reproduce the balance between various inputs' strengths to different layers manifested in the different CSD sink magnitudes. It would then be important to recapitulate this between-layer relative parameter of the input strength in the model (see comments on that above as well).
It is unclear why 90 minutes for 1 second of simulation presents such a barrier. When compared to experimental work this time seems realistic to provide much more statistically powerfull dataset using simulations. I would suggest increasing the number of trials of simulation.
The reviewer is right that the L4 sink appears too weak relative to the other sinks and sources when compared to what’s typically observed in our own and previously published experimental data. However, we would like to point out that there are several examples in the experiments we analyzed where the L4 sink is weak compared to the other sinks and sources (see session ids 816200189, 778240327, 799864342, and 794812542 in Figure 2 —figure supplement 1), so the model isn’t necessarily producing a CSD pattern that’s not biologically plausible, even if it is less stereotypical. Nonetheless, we intend to explore ways to make the L4 sink stronger relative to other sinks and sources in future work.
Regarding the relative strengths of different CSD sink magnitudes: this is already a part of the Wasserstein distance metric. For example, if the relative strengths of the sinks and sources in infragranular layers are too strong compared to sinks and sources in the granular and supragranular layers (as was the case in the original model), this is reflected in a greater Wasserstein distance (Figure 4B and 4D). The Wasserstein distance also allows us to explore how the sinks of one pattern are moved to match the sinks of another pattern (illustrated in Figure 2C) without any human characterization of what constitutes a L4 sink, a L2/3 sink, or a L5/L6 sink. This means that we already have a metric that allows us to assess the relative strengths of sinks in each layer.
The question of the number of trials in a simulation is addressed in point 10 above.
17. The authors compare the adjusted model to the original model but never mention what the original model was tuned to. "The synapses for recurrent connections were placed according to the following rules in the original model (Billeh et al. 2020)". As a reader, it would be helpful to have an idea of what the rules are based on. Is synaptic placement based on anatomy? Are synaptic weights fitted to evoked response to a particular visual stimulus?
A description of this has been added to the “Dendritic targeting” (on lines 1650-1651) and “Adjusting synaptic weights” (on lines 1695-1699) sections in Methods.
On lines 1650-1651:
“The rules for placement of synapses in Billeh et al. 2020 were based on reviews of literature on anatomy.”
On lines 1695-1699:
“In the original model, the synaptic weights of thalamocortical connections were based on experimental recordings of synaptic current, while the synaptic weights for recurrent connections were initially set to estimates from literature, then optimized to a drifting gratings stimulus until the model reproduced experimental values of orientation and direction selectivity.”
18a. Discussion should include a much more balanced interpretation of the results – it is merely an illustration of some ideas using the model, at the moment, rather than quantitative analysis demonstrating novel results on the composition of the CSD response, importance of the feedback, or recurrency in the actual experimental data.
b. Likewise, the mismatch of the firing rate and CSD is simply illustrated by some examples, rather than revealing a novel result. In fact, the role of the interneurons in this mismatch is likely the key, but it was not discussed. I cannot imagine how a lack of realism in the interneuronal connectivity could be ignored (different PV morphological types, SOM, VIP) if one wants to recapitulate the circuit dynamics in response to the thalamic input. In other systems, for example, hippocampus, the mismatch between the CSD of theta and firing rates of different neuronal groups in entorhino-hippocampal circuits is much better understood and could be compared to here.
a) With the quantitative analysis added to the manuscript (see points 6, 7, 9, and 12 above), we consider the current interpretation of our results to be fair.
b) It is not clear why lack of realism is mentioned here. Our model includes distinct and experimentally constrained connectivity patterns for PV, SOM, and VIP neurons, parameterized individually for each layer based on experimental data. This realistic interneuronal connectivity was implemented in the model in Billeh et al., 2020, which we are using here as a starting point, as was different PV, SOM, and VIP morphological types.
19. Discussion of other inputs not yet accounted for in the model that contribute to the experimental results is missing. What about other thalamic nuclei? Role of cortio-thalamic feedback?
The effects of cortico-thalamic feedback should, at least partly, be present in the LGN input to the model, as we are using experimentally recorded spike trains as input. The role of other thalamic nuclei is beyond the scope of our study but would be an interesting avenue for future work.
20. It is clear that further improvements in the model realism will make it prohibitively complex. Authors postpone the future proper multiparametric optimization due to computational complexity. It would thus be useful to discuss the alternative approach – capturing parameters of the low dimensional dynamical system that would recapitulate the observed data or model.
We and other members of our research groups are working on approaches that can produce low dimensional dynamical systems that recapitulate the observed data or model, and the findings of these projects will be presented in future papers.
Reviewer #3 (Recommendations for the authors):
1. While quantification of the activity of different cell types in different layers is necessarily limited in the empirical data, in the model it is not. The model contains more than 50,000 biophyscially detailed neuron models belonging to 17 different cell type classes: one inhibitory class (Htr3a) in layer 1, and four classes in each of the other layers (2/3, 4, 5, and 6) where one is excitatory and three are inhibitory (Pvalb, Sst, Htr3a) in each layer. In constraining the model to a limited number of features in the data (firing rates of limited cell populations and CSD), the model also makes many more novel and testable predictions about the activity of other cell types and synaptic connections (not fit too) that would be valuable for the readers to see.
For example, while the data is sparse in terms of knowledge of different inhibitory neuron types across populations, the model is not, and can provide testable predictions on the firing activity of the inhibitory neurons in each layer and across three different cell types Pvalb, Sst, Htr3a, along with Htr3a neurons in layer 1. In the final fit model as in Figure 7, it would be useful to detail model predictions about the firing rates in the 4 different interneuron types, and their strength of influence on the RS neurons. These predictions could then be further tested in follow up studies, and perhaps there is a prediction that was not used for fitting but that can be tested in the data set here?
Along this line, what potential mechanisms govern the animals or trials that present atypical CSD or spike rate evoked responses? The methods and results of this study are uniquely positioned to not only account for the circuit mechanisms of a canonical response, but also the cell and circuit mechanisms of a distribution visual flash evoked responses, particularly those categorized as non-canonical. How varied are the non-canonical responses? A supplemental figure showing the CSD from each of the animals, e.g. as shown for other features in Figure S4 (in the original submission, Figure 2 —figure supplement 3 in the current submission), would also be helpful to understand the full temporal spatial variability in the CSD. The present study has a distinctly powerful opportunity to make predictions regarding the underlying cell types, circuit elements, connectivity patterns, etc. that control for visual evoked response variability.
This is an interesting suggestion, and one we can pursue in future explorations of the model. Done thoroughly, this could create enough material to warrant another article on its own.
We have included a supplemental figure displaying the CSD from each animal (Figure 2 —figure supplement 1) to give the reader a better idea of the variability in CSD responses.
2. Related to this, there are a few important assumptions in the model that should be clarified further:
a. Are the LM to V1 connections known to not target L1 inhibitory neurons? What happens if they do in the simulations?
b. Do the "feedback" inputs in the model also contact the distal dendrites of the L6 RS cells? What about the apical dendrites of the L4 RS cells? What about the L1 cells? Is there any literature supporting the choice to connect the feedback to only a subset of the neurons in the supragranular layers or was it purely because it was necessary to reduce the magnitude of the L5/L6 dipole?
a) The literature does not suggest that LM to L1 connections are known not to target L1 inhibitory neurons. However, to our knowledge, only apical dendrites of neurons whose somata reside in L2/3 and L5 are mentioned as targets of connections that terminate in L1. Still, we have run simulations that included connections between LM and L1 inhibitory neurons. We have added plots that display the results from these simulations in Appendix 1 – Figure 1.
b): No, the feedback does not target L4, L6, or L1 cells. We connected the feedback only to L2/3 and L5 cells because those were the only cells that were mentioned as targets in the literature. We cannot rule out that LM targets cells in the other layers as well, but since the literature does not indicate these layers as important targets for LM afferents, we chose to leave them out for now.
3. There are important details in the comparison made between the model neuron firing rates and empirically recording firing rates that need to be understood more fully. In the empirical recordings, cells are sorted into RS or FS populations, and RS are presumed excitatory with FS presumed inhibitory. The activity of the RS is separated by layer, but the FS is grouped across layers because of the sparsity of recordings. When comparing the model firing rates to the experimental firing rates, the excitatory and "non-Pvalb populations" were grouped together in each layer of the model to make up the RS cells in L2/3, L4, L5, and L6, while the Pvalb cells across all layers were grouped together to make up the FS cells of V1. This makes sense for the experimental data, however, the model has distinct non-Pvalb inhibitory populations including SST and Htr3a. Were these non-Pvalb inhibitory cells grouped with the FS cells or RS cells in the model analysis? Clarification is needed.
The non-PV inhibitory cells were grouped together with excitatory cells to make up the RS populations in the model analysis, which corresponds to the current understanding of what the experimentally observed RS and FS populations are: FS are mostly PV inhibitory cells, whereas RS are the excitatory AND non-PV inhibitory cells. This is explained on lines 340-342 and 1572-1575 but we acknowledge that this could have been made clearer.
4. Line 557 states that to bring the firing rates of the neurons in the original model closer to the experimental data ["We first reduced the synaptic weights from all excitatory populations to the fast-spiking PV- neurons by 30% to bring their firing rates closer to the average firing rate in this population in the experiments. This resulted in increased firing rates in all other (RS) populations due to the reduced activity of the inhibitory Pvalb-neurons (Figure S6 in the original submission, figure 4 —figure supplement 1 in the current submission). Therefore, we further applied reductions in the synaptic weights from all excitatory neurons to RS neurons or increases in the synaptic weights from inhibitory neurons to the RS neurons to bring their firing rates closer to the experimental average firing rates."]
a. The "or" statement here is important, should it be "and"?
b. Was it all inhibitory neurons or only the Pvalb-neurons?
a): It should be “and”. The relevant part of the manuscript has been corrected.
b): It was only the Pvalb-neurons.
5. The process of parameter tuning in models, particularly in large-scale models, is challenging. Throughout the results, the authors walk through a series of manipulations to the original model that are necessary to account for the empirical data. While I understand that not all parameters can be estimated, and there is a nice paragraph in the discussion on this starting on Line 981, it would be helpful to have a methods section that describes the overall strategy for parameter manipulations and rationale for the choice of the final parameters.
a. Was it all hand tuning, or was there some estimation method used to fit to data?
b. More specifically, there is a section called "Adjusting synaptic weights' that gives a high level description of a range of multiplication factors over which previous model parameters were adjusted to fit the current data. More details on how this range was determined, and if one value in this range or multiple values were ultimately used in the simulations, and why, is needed.
a) It was all hand tuning. We postpone automatic estimation to future projects because the computational cost of running the number of simulations required for such methods is currently prohibitively high with this model.
b) The ranges were set after experimenting with different ranges to find the ones that would allow the model to reproduce the experimental observations. Only a single value within each range was used in the final model. We have added this information to the relevant part in the Methods section of the manuscript.
6a. Details on how the new feedback was implemented in the model is described on Lines 682-690, but I didn't see any description of these inputs in the method section.
b. It seems the "Adjusting synaptic weights" section only describes local interactions. Were there feedback inputs in the original Billeh model that were adjusted here? Or are they completely new, and if so, what were the starting weights and how were they chosen?
c. The modeling results here build from the prior work in Billeh et al. 2020. In this paper, the model is applied to visual flash evoked responses, while in Billeh et al. 2020 the model was tuned to orientation and direction tuning to drifting grating stimuli. It is exciting to see that the updated model reproduces consistent results for firing rates and rate based direction/orientation selectivity, as shown in Supplementary Figure S9 (in the original submission, Figure 6 —figure supplement 6 in the current submission). It is stated that for these simulations, the stimulation of the network was as defined in Billeh et al. 2020. One important piece of information to understand is why feedback is needed in the current simulations of flash evoked response, but isn't required for the drifting grading simulations? Would including feedback improve the results? A discussion and/or table summarizing the differences in the two models would be useful to understand what further manipulations may be necessary to accurately account for the LFP/CSD, in addition to firing rates, in Billeh et al. 2020 simulations.
a) A description of the feedback input can be found under “Input from lateromedial area (LM)” under “Visual stimulus” in the Method section.
b) There was no feedback input to the original Billeh et al. 2020 model. In the revised version of the model that we report in this manuscript, the starting weights were set equal to the weights from the background node to the feedback targeted population. This was done because the input from the background node was a representation of the influence from the rest of the brain before the feedback input from LM had been added, which means that in the original model background represented the influence from HVAs (in addition to influences from all other areas of the brain), and could therefore be seen as a simple, first approximation of the weights from LM.
c) Feedback was included in the drifting gratings simulation (Figure S9 in the original submission, Figure 6 —figure supplement 6 in the current submission) too. Though the simplicity of the connectivity and input from LM in this model meant that we did not expect the direction and orientation selectivity of the model improve, we included it to make this simulation more similar to the final model simulation of the evoked flash response, so that we could test if any of the previous features of the model had been “broken” by the inclusion of feedback.
7. One strategy the authors used to examine the importance of various neural mechanisms to the CSD, is to remove various model elements such as synaptic connections or intrinsic ionic currents and see how much the CSD changes. For example, there is a conclusion line 803 that [without inhibition that fact that] "the major sinks and sources are still present is an indication that the currents from excitatory input account for the majority of the sinks and sources observed in the experimental CSD."
a. This could be due to the fact that in the fit model (Figure 7) the inhibitory synapses are very weak (are they?) and hence not playing a major role in the established model configuration.
b. The conclusion that inhibition plays a small role is counterintuitive to the growing literature showing the importance of inhibition to LFP generation (e.g. Telenczuk et al. Scientific Reports 2017), and hence CSD, and further discussion and understanding of the role of the different types of inhibition in this model would be helpful.
a) We do not have reason to believe that the inhibitory synapses are weak relative to excitatory synapses in the model, as the model reproduces the experimentally observed firing rates of both fast-spiking PV-cells and regular-spiking cells. If the inhibitory synapses were too weak, one would expect the firing rate in excitatory (RS) populations to be significantly higher than the experimental firing rates in these populations. It is possible that both the excitation and inhibition are weak, which could produce firing rates compatible with experimental observation, but that should in any case not result in a disproportionate contribution to the CSD from the excitatory populations, since both excitatory and inhibitory synapses are weak in this scenario.
b) This is an important point. The following discussion of our findings in light of this literature has been added on lines 1031-1040:
“Recent investigations into the unitary LFP (the LFP generated by a single neuron) from inhibitory and excitatory synapses have suggested that inhibitory inputs exert a greater influence on LFP than excitatory input (Bazelot et al., 2010; Teleńczuk et al., 2017; B. Teleńczuk, M. Teleńczuk, and Destexhe, 2020). While this may be true of unitary effects, the total effect of excitatory input can still be greater if there are significantly more excitatory than inhibitory cells, and, correspondingly, significantly more excitatory synapses. In this V1 model, inhibitory cells make up about 15 % and excitatory cells about 85 % of the total number of cells, reflecting the cellular composition in mouse V1. Whether the pronounced dominance of excitatory cells is enough to make up for the reduced unitary influence of excitatory cells is an interesting area for further research.”
8. A more useful analysis to quantify the impact of specific network features would be a parameter sensitivity analysis that shows the influence of changes in parameters to variance in the resultant CSD; this isn't feasible for all parameters but could perhaps be done for the parameters/network features concluded to account (or not account) for the main feature of the canonical CSD, e.g. as summarized beginning on line 827. Or at least some simulations that show what happens over a range of parameter values for the important parameters.
This is a good suggestion, but the analysis in this paper is already extensive in scope. We therefore postpone such an analysis to a future paper.
9. A similar large scale biophysically detailed cortical column modeling paper (Reimann et al. Neuron 2013) is relevant to discuss further in light of the current results. In particular, that paper showed that "active currents" dominated the LFP, rather than synaptic currents. There is a statement line 890 that "In line with observations made by Reimann et al. (2013), we found that the somatic voltage- dependent membrane currents significantly shape the CSD signature (Figure 5H and 7B)". To my knowledge, Riemann considered the contribution of all active currents, including those in the dendrites, while in this paper they have only explored removing them in the soma. Those in the dendrites may have a greater impact, or maybe not (?). Can the authors show this? Further examination/discussion would be useful.
The number and level of detail of biophysically detailed neurons in this model unfortunately makes simulations with active conductances in the dendrites prohibitively expensive at this time.
10. Of particular relevance to these results are earlier laminar recordings to flash evoked responses in VI in non-human primates that suggest the presence of "feedback" and discusses interaction between higher and lower order areas that may contribute to this "feedback": doi: 10.1016/0042-6989(94)90156-2 https://doi.org/10.1093/cercor/8.7.575. A discussion on the conclusions in this paper related to those is needed.
These papers, and many other papers, also show a decoupling between LFP/CSD and firing (MUA). Clarification should be made that such decoupling is not a new concept, but to my knowledge it has yet to be shown how it could be mechanistically implemented using a computational neural model, as the authors have nicely done. Nor has it been shown how simulated LFP/CSD can be altered with only minor effects on the firing rates, as shown here.
Also relevant, is a more recent study of flash evoked responses in humans and non-human primates doi: 10.1126/sciadv.abb0977, which appears to show consistent results to those reported here.
Discussions of these papers have been added on lines 897-901 and 966-970. We see no explicit mention of feedback in these papers, but we do discuss whether CSD patterns in our data align with the patterns described in those papers as well as their observations on decoupling between LFP/CSD and MUA. Other studies by the same group that focus on the effect of feedback in evoked V1 responses are discussed on lines 974-979.
On lines 897-901:
“Studies in non-human primates where the animals were exposed to flashes also demonstrate good spatio-temporal agreement with the CSD observed here, with sinks and sources occurring not only in the same layers and in the same order, but also at the same time after stimulus presentation (Givre, Schroeder, and Arezzo, 1994; Schroeder, Mehta, and Givre, 1998).”
On lines 966-970:
“These results align with findings in e.g. Schroeder, Mehta, and Givre, 1998 and Leszczyński et al. 2020, where a lack of correlation between MUA and CSD in the upper layers of V1 during flash exposure to non-human primates suggested that these signals can sometimes be decoupled.”
On lines 974-979:
“Past studies have suggested that LFP can be modulated by attention through feedback during evoked responses (Mehta, Ulbert, and Schroeder, 2000a; Mehta, Ulbert, and Schroeder, 2000b). However, their findings indicated that V1 was not significantly affected by this modulation until the period 250 ms or later after stimulus onset. A recent study showed that feedback from higher visual areas can in fact exert a strong influence on the magnitude of LFP already at around 80 ms after stimulus onset (Hartmann et al., 2019).”
There were many interesting suggestions for further analyses in our study, and though we address many of them, as described above, several would, in our view, generate enough material to warrant another paper, which appears excessive. Our manuscript is already extensive and ambitious in scope and level of detail, and, we believe, has provided insights that no other study so far has been able to offer. We believe the revisions described above have substantially improved this work and address the most important concerns of the reviewers.
https://doi.org/10.7554/eLife.87169.sa2Article and author information
Author details
Funding
Simula School of Research
- Atle E Rimehaug
European Union Horizon 2020 Research and Innovation program (785907)
- Espen Hagen
European Union Horizon 2020 Research and Innovation program (945539)
- Espen Hagen
Research Council of Norway (COBRA - project number 250128)
- Alexander J Stasik
IKTPLUSS-IKT and Digital Innovation (300504)
- Alexander J Stasik
National Institute of Neurological Disorders and Stroke (R01NS122742)
- Kael Dai
- Josh H Siegle
- Shawn R Olsen
- Christof Koch
- Anton Arkhipov
- Yazan N Billeh
National Institute of Biomedical Imaging and Bioengineering (R01EB029813)
- Kael Dai
- Josh H Siegle
- Shawn R Olsen
- Christof Koch
- Anton Arkhipov
- Yazan N Billeh
Allen Institute
- Kael Dai
- Josh H Siegle
- Shawn R Olsen
- Christof Koch
- Anton Arkhipov
- Yazan N Billeh
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
Research reported in this publication was supported the Simula School of Research and Innovation, CINPLA, the European Union Horizon 2020 Research and Innovation Program under Grant Agreement No. 785907 and No. 945539 (Human Brain Project [HBP] SGA2 and SGA3), the National Institute of Neurological Disorders and Stroke of the National Institutes of Health under Award Number R01NS122742, and by the National Institute of Biomedical Imaging and Bioengineering of the National Institutes of Health under Award Number R01EB029813, as well as by the Allen Institute. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. We acknowledge the use of Fenix Infrastructure resources, which are partially funded from the European Union’s Horizon 2020 research and innovation program through the ICEI project under the grant agreement No. 800858. We thank the Allen Institute founder, Paul G Allen, for his vision, encouragement, and support.
Senior and Reviewing Editor
- Tirin Moore, Howard Hughes Medical Institute, Stanford University, United States
Reviewer
- Stephanie R Jones, Brown University, United States
Version history
- Preprint posted: February 25, 2022 (view preprint)
- Received: February 23, 2023
- Accepted: July 10, 2023
- Accepted Manuscript published: July 24, 2023 (version 1)
- Version of Record published: August 1, 2023 (version 2)
- Version of Record updated: August 11, 2023 (version 3)
Copyright
© 2023, Rimehaug et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 872
- Page views
-
- 154
- Downloads
-
- 3
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Computational and Systems Biology
- Immunology and Inflammation
T cells are required to clear infection, and T cell motion plays a role in how quickly a T cell finds its target, from initial naive T cell activation by a dendritic cell to interaction with target cells in infected tissue. To better understand how different tissue environments affect T cell motility, we compared multiple features of T cell motion including speed, persistence, turning angle, directionality, and confinement of T cells moving in multiple murine tissues using microscopy. We quantitatively analyzed naive T cell motility within the lymph node and compared motility parameters with activated CD8 T cells moving within the villi of small intestine and lung under different activation conditions. Our motility analysis found that while the speeds and the overall displacement of T cells vary within all tissues analyzed, T cells in all tissues tended to persist at the same speed. Interestingly, we found that T cells in the lung show a marked population of T cells turning at close to 180o, while T cells in lymph nodes and villi do not exhibit this “reversing” movement. T cells in the lung also showed significantly decreased meandering ratios and increased confinement compared to T cells in lymph nodes and villi. These differences in motility patterns led to a decrease in the total volume scanned by T cells in lung compared to T cells in lymph node and villi. These results suggest that the tissue environment in which T cells move can impact the type of motility and ultimately, the efficiency of T cell search for target cells within specialized tissues such as the lung.
-
- Computational and Systems Biology
Biomedical single-cell atlases describe disease at the cellular level. However, analysis of this data commonly focuses on cell-type centric pairwise cross-condition comparisons, disregarding the multicellular nature of disease processes. Here we propose multicellular factor analysis for the unsupervised analysis of samples from cross-condition single-cell atlases and the identification of multicellular programs associated with disease. Our strategy, which repurposes group factor analysis as implemented in multi-omics factor analysis, incorporates the variation of patient samples across cell-types or other tissue-centric features, such as cell compositions or spatial relationships, and enables the joint analysis of multiple patient cohorts, facilitating the integration of atlases. We applied our framework to a collection of acute and chronic human heart failure atlases and described multicellular processes of cardiac remodeling, independent to cellular compositions and their local organization, that were conserved in independent spatial and bulk transcriptomics datasets. In sum, our framework serves as an exploratory tool for unsupervised analysis of cross-condition single-cell atlases and allows for the integration of the measurements of patient cohorts across distinct data modalities.