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.

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.

Figure 1

Download asset Open asset

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.

Figure 2 with 3 supplements see all

Download asset Open asset

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

Figure 3 with 2 supplements see all

Download asset Open asset

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

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.

Figure 4 with 1 supplement see all

Download asset Open asset

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: KSS_b = 0.62, KSS_p = 0.63, and KSS_s = 0.54; RS L4: KSS_b = 0.77, KSS_p = 0.60, and KSS_s = 0.63; RS L5 KSS_b = 0.77, KSS_p = 0.77, and KSS_s = 0.78; RS L6: KSS_b = 0.54, KSS_p = 0.45, and KSS_s = 0.47; FS: KSS_b = 0.54, KSS_p = 0.53, and KSS_s = 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 (KSS_b = 0.73, KSS_p = 0.77, and KSS_s = 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).

Figure 5 with 4 supplements see all

Download asset Open asset

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

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

Figure 6 with 9 supplements see all

Download asset Open asset

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 (KSS_b = 0.77, KSS_p = 0.70, and KSS_s = 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.).

Figure 7 with 3 supplements see all

Download asset Open asset

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.

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.

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

1. 1. Rimehaug AR
   2. Stasik AJ
   3. Hagen E
   4. Billeh YN
   5. Siegle JH
   6. Dai K
   7. Olsen SR
   8. Koch C
   9. Einevoll G
   10. Arkhipov A

   (2022) Dryad Digital Repository

   Uncovering circuit mechanisms of current sinks and sources with biophysical simulations of primary visual cortex.

   https://doi.org/10.5061/dryad.k3j9kd5b8

1. 1. Siegle JH
   2. Wakeman W
   3. Jia X
   4. Heller G
   5. Ramirez T
   6. Graddis N
   7. Mei N
   8. Durand S

   (2020) DANDI Archive 000021

   ID 000021. 20191003_AIBS_mouse_ecephys_brain_observatory_1_1.

   https://dandiarchive.org/dandiset/000021

1. 1. Andersen RA
   2. Musallam S
   3. Pesaran B

   (2004) Selecting the signals for a brain-machine interface

   Current Opinion in Neurobiology 14:720–726.

   https://doi.org/10.1016/j.conb.2004.10.005

   • PubMed
   • Google Scholar
2. Conference

   1. Arjovsky M
   2. Chintala S
   3. Bottou L

   (2017)

   Wasserstein generative adversarial networks

   Proceedings of the 34th International Conference on Machine Learnin. pp. 214–223.

   • Google Scholar
3. 1. Arkhipov A
   2. Gouwens NW
   3. Billeh YN
   4. Gratiy S
   5. Iyer R
   6. Wei Z
   7. Xu Z
   8. Abbasi-Asl R
   9. Berg J
   10. Buice M
   11. Cain N
   12. da Costa N
   13. de Vries S
   14. Denman D
   15. Durand S
   16. Feng D
   17. Jarsky T
   18. Lecoq J
   19. Lee B
   20. Li L
   21. Mihalas S
   22. Ocker GK
   23. Olsen SR
   24. Reid RC
   25. Soler-Llavina G
   26. Sorensen SA
   27. Wang Q
   28. Waters J
   29. Scanziani M
   30. Koch C

   (2018) Visual physiology of the layer 4 cortical circuit in silico

   PLOS Computational Biology 14:e1006535.

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

   • PubMed
   • Google Scholar
4. 1. Aru J
   2. Suzuki M
   3. Larkum ME

   (2020) Cellular mechanisms of conscious processing

   Trends in Cognitive Sciences 24:814–825.

   https://doi.org/10.1016/j.tics.2020.07.006

   • PubMed
   • Google Scholar
5. 1. Barthó P
   2. Hirase H
   3. Monconduit L
   4. Zugaro M
   5. Harris KD
   6. Buzsáki G

   (2004) Characterization of neocortical principal cells and interneurons by network interactions and extracellular features

   Journal of Neurophysiology 92:600–608.

   https://doi.org/10.1152/jn.01170.2003

   • PubMed
   • Google Scholar
6. 1. Bastos AM
   2. Vezoli J
   3. Bosman CA
   4. Schoffelen JM
   5. Oostenveld R
   6. Dowdall JR
   7. De Weerd P
   8. Kennedy H
   9. Fries P

   (2015) Visual areas exert feedforward and feedback influences through distinct frequency channels

   Neuron 85:390–401.

   https://doi.org/10.1016/j.neuron.2014.12.018

   • PubMed
   • Google Scholar
7. 1. Bazelot M
   2. Dinocourt C
   3. Cohen I
   4. Miles R

   (2010) Unitary inhibitory field potentials in the CA3 region of rat hippocampus

   The Journal of Physiology 588:2077–2090.

   https://doi.org/10.1113/jphysiol.2009.185918

   • PubMed
   • Google Scholar
8. 1. Belitski A
   2. Gretton A
   3. Magri C
   4. Murayama Y
   5. Montemurro MA
   6. Logothetis NK
   7. Panzeri S

   (2008) Low-frequency local field potentials and spikes in primary visual cortex convey independent visual information

   The Journal of Neuroscience 28:5696–5709.

   https://doi.org/10.1523/JNEUROSCI.0009-08.2008

   • PubMed
   • Google Scholar
9. 1. Billeh YN
   2. Cai B
   3. Gratiy SL
   4. Dai K
   5. Iyer R
   6. Gouwens NW
   7. Abbasi-Asl R
   8. Jia X
   9. Siegle JH
   10. Olsen SR
   11. Koch C
   12. Mihalas S
   13. Arkhipov A

   (2020) Systematic integration of structural and functional data into multi-scale models of mouse primary visual cortex

   Neuron 106:388–403.

   https://doi.org/10.1016/j.neuron.2020.01.040

   • PubMed
   • Google Scholar
10. Book

    1. Buzsáki G

    (2006) Rhythms of the Brain

    Oxford university press.

    https://doi.org/10.1093/acprof:oso/9780195301069.001.0001

    • Google Scholar
11. 1. Buzsáki G
    2. Anastassiou CA
    3. Koch C

    (2012) The origin of extracellular fields and currents--EEG, ECoG, LFP and spikes

    Nature Reviews. Neuroscience 13:407–420.

    https://doi.org/10.1038/nrn3241

    • PubMed
    • Google Scholar
12. 1. Chatzikalymniou AP
    2. Skinner FK

    (2018) Deciphering the Contribution of Oriens-Lacunosum/Moleculare (OLM) Cells to Intrinsic θ Rhythms Using Biophysical Local Field Potential (LFP) Models

    eNeuro 5:ENEURO.0146-18.2018.

    https://doi.org/10.1523/ENEURO.0146-18.2018

    • PubMed
    • Google Scholar
13. 1. Dai K
    2. Gratiy SL
    3. Billeh YN
    4. Xu R
    5. Cai B
    6. Cain N
    7. Rimehaug AE
    8. Stasik AJ
    9. Einevoll GT
    10. Mihalas S
    11. Koch C
    12. Arkhipov A
    13. Marinazzo D

    (2020) Brain Modeling ToolKit: An open source software suite for multiscale modeling of brain circuits

    PLOS Computational Biology 16:e1008386.

    https://doi.org/10.1371/journal.pcbi.1008386

    • PubMed
    • Google Scholar
14. 1. Di S
    2. Baumgartner C
    3. Barth DS

    (1990) Laminar analysis of extracellular field potentials in rat vibrissa/barrel cortex

    Journal of Neurophysiology 63:832–840.

    https://doi.org/10.1152/jn.1990.63.4.832

    • PubMed
    • Google Scholar
15. 1. Drion G
    2. O’Leary T
    3. Marder E

    (2015) Ion channel degeneracy enables robust and tunable neuronal firing rates

    PNAS 112:E5361–E5370.

    https://doi.org/10.1073/pnas.1516400112

    • PubMed
    • Google Scholar
16. 1. Einevoll GT
    2. Pettersen KH
    3. Devor A
    4. Ulbert I
    5. Halgren E
    6. Dale AM

    (2007) Laminar population analysis: estimating firing rates and evoked synaptic activity from multielectrode recordings in rat barrel cortex

    Journal of Neurophysiology 97:2174–2190.

    https://doi.org/10.1152/jn.00845.2006

    • PubMed
    • Google Scholar
17. 1. Einevoll GT
    2. Kayser C
    3. Logothetis NK
    4. Panzeri S

    (2013) Modelling and analysis of local field potentials for studying the function of cortical circuits

    Nature Reviews. Neuroscience 14:770–785.

    https://doi.org/10.1038/nrn3599

    • PubMed
    • Google Scholar
18. 1. Fernández-Ruiz A
    2. Makarov VA
    3. Benito N
    4. Herreras O

    (2012) Schaffer-specific local field potentials reflect discrete excitatory events at gamma frequency that may fire postsynaptic hippocampal CA1 units

    The Journal of Neuroscience 32:5165–5176.

    https://doi.org/10.1523/JNEUROSCI.4499-11.2012

    • PubMed
    • Google Scholar
19. 1. Fernández-Ruiz A
    2. Oliva A
    3. Nagy GA
    4. Maurer AP
    5. Berényi A
    6. Buzsáki G

    (2017) Entorhinal-CA3 Dual-Input Control of Spike Timing in the Hippocampus by Theta-Gamma Coupling

    Neuron 93:1213–1226.

    https://doi.org/10.1016/j.neuron.2017.02.017

    • PubMed
    • Google Scholar
20. 1. Givre SJ
    2. Schroeder CE
    3. Arezzo JC

    (1994) Contribution of extrastriate area V4 to the surface-recorded flash VEP in the awake macaque

    Vision Research 34:415–428.

    https://doi.org/10.1016/0042-6989(94)90156-2

    • PubMed
    • Google Scholar
21. Conference

    1. Głąbska H
    2. Chintaluri HC
    3. Wójcik DK

    (2014) Collection of Simulated Data for Validation of Methods of Analysis of Extracellular Potentials

    Neuroinformatics 2014.

    https://doi.org/10.3389/conf.fninf.2014.18.00035

    • Google Scholar
22. 1. Glickfeld LL
    2. Olsen SR

    (2017) Higher-Order Areas of the Mouse Visual Cortex

    Science 3:251–273.

    https://doi.org/10.1146/annurev-vision-102016-061331

    • Google Scholar
23. 1. Gouwens NW
    2. Berg J
    3. Feng D
    4. Sorensen SA
    5. Zeng H
    6. Hawrylycz MJ
    7. Koch C
    8. Arkhipov A

    (2018) Systematic generation of biophysically detailed models for diverse cortical neuron types

    Nature Communications 9:710.

    https://doi.org/10.1038/s41467-017-02718-3

    • PubMed
    • Google Scholar
24. 1. Hagen E
    2. Dahmen D
    3. Stavrinou ML
    4. Lindén H
    5. Tetzlaff T
    6. van Albada SJ
    7. Grün S
    8. Diesmann M
    9. Einevoll GT

    (2016) Hybrid scheme for modeling local field potentials from point-neuron networks

    Cereb. Cortex 26:4461–4496.

    https://doi.org/10.1093/cercor/bhw237

    • PubMed
    • Google Scholar
25. 1. Hagen E
    2. Fossum JC
    3. Pettersen KH
    4. Alonso JM
    5. Swadlow HA
    6. Einevoll GT

    (2017) Focal local field potential signature of the single-axon monosynaptic thalamocortical connection

    The Journal of Neuroscience 37:5123–5143.

    https://doi.org/10.1523/JNEUROSCI.2715-16.2017

    • PubMed
    • Google Scholar
26. 1. Hagen E
    2. Næss S
    3. Ness TV
    4. Einevoll GT

    (2018) Multimodal Modeling of Neural Network Activity: Computing LFP, ECoG, EEG, and MEG Signals With LFPy 2.0

    Frontiers in Neuroinformatics 12:92.

    https://doi.org/10.3389/fninf.2018.00092

    • PubMed
    • Google Scholar
27. 1. Harris JA
    2. Mihalas S
    3. Hirokawa KE
    4. Whitesell JD
    5. Choi H
    6. Bernard A
    7. Bohn P
    8. Caldejon S
    9. Casal L
    10. Cho A
    11. Feiner A
    12. Feng D
    13. Gaudreault N
    14. Gerfen CR
    15. Graddis N
    16. Groblewski PA
    17. Henry AM
    18. Ho A
    19. Howard R
    20. Knox JE
    21. Kuan L
    22. Kuang X
    23. Lecoq J
    24. Lesnar P
    25. Li Y
    26. Luviano J
    27. McConoughey S
    28. Mortrud MT
    29. Naeemi M
    30. Ng L
    31. Oh SW
    32. Ouellette B
    33. Shen E
    34. Sorensen SA
    35. Wakeman W
    36. Wang Q
    37. Wang Y
    38. Williford A
    39. Phillips JW
    40. Jones AR
    41. Koch C
    42. Zeng H

    (2019) Hierarchical organization of cortical and thalamic connectivity

    Nature 575:195–202.

    https://doi.org/10.1038/s41586-019-1716-z

    • PubMed
    • Google Scholar
28. Preprint

    1. Hartmann TS
    2. Raja S
    3. Lomber SG
    4. Born RT

    (2019) Cortico-Cortical Feedback from V2 Exerts a Powerful Influence over the Visually Evoked Local Field Potential and Associated Spike Timing in V1

    bioRxiv.

    https://doi.org/10.1101/792010

    • Google Scholar
29. 1. Hay E
    2. Hill S
    3. Schürmann F
    4. Markram H
    5. Segev I

    (2011) Models of neocortical layer 5b pyramidal cells capturing a wide range of dendritic and perisomatic active properties

    PLOS Computational Biology 7:e1002107.

    https://doi.org/10.1371/journal.pcbi.1002107

    • PubMed
    • Google Scholar
30. 1. Henrie JA
    2. Shapley R

    (2005) LFP power spectra in V1 cortex: the graded effect of stimulus contrast

    Journal of Neurophysiology 94:479–490.

    https://doi.org/10.1152/jn.00919.2004

    • PubMed
    • Google Scholar
31. 1. Holt GR
    2. Koch C

    (1999) Electrical interactions via the extracellular potential near cell bodies

    Journal of Computational Neuroscience 6:169–184.

    https://doi.org/10.1023/a:1008832702585

    • PubMed
    • Google Scholar
32. 1. Jun JJ
    2. Steinmetz NA
    3. Siegle JH
    4. Denman DJ
    5. Bauza M
    6. Barbarits B
    7. Lee AK
    8. Anastassiou CA
    9. Andrei A
    10. Aydın Ç
    11. Barbic M
    12. Blanche TJ
    13. Bonin V
    14. Couto J
    15. Dutta B
    16. Gratiy SL
    17. Gutnisky DA
    18. Häusser M
    19. Karsh B
    20. Ledochowitsch P
    21. Lopez CM
    22. Mitelut C
    23. Musa S
    24. Okun M
    25. Pachitariu M
    26. Putzeys J
    27. Rich PD
    28. Rossant C
    29. Sun WL
    30. Svoboda K
    31. Carandini M
    32. Harris KD
    33. Koch C
    34. O’Keefe J
    35. Harris TD

    (2017) Fully integrated silicon probes for high-density recording of neural activity

    Nature 551:232–236.

    https://doi.org/10.1038/nature24636

    • PubMed
    • Google Scholar
33. 1. Kandel A
    2. Buzsáki G

    (1997) Cellular-synaptic generation of sleep spindles, spike-and-wave discharges, and evoked thalamocortical responses in the neocortex of the rat

    The Journal of Neuroscience 17:6783–6797.

    https://doi.org/10.1523/JNEUROSCI.17-17-06783.1997

    • PubMed
    • Google Scholar
34. 1. Keller AJ
    2. Roth MM
    3. Scanziani M

    (2020) Feedback generates a second receptive field in neurons of the visual cortex

    Nature 582:545–549.

    https://doi.org/10.1038/s41586-020-2319-4

    • PubMed
    • Google Scholar
35. Book

    1. Koch C

    (1998) Biophysics of Computation: Information Processing in Single Neurons

    Oxford university press.

    https://doi.org/10.1093/oso/9780195104912.001.0001

    • Google Scholar
36. 1. Kreiman G
    2. Hung CP
    3. Kraskov A
    4. Quiroga RQ
    5. Poggio T
    6. DiCarlo JJ

    (2006) Object selectivity of local field potentials and spikes in the macaque inferior temporal cortex

    Neuron 49:433–445.

    https://doi.org/10.1016/j.neuron.2005.12.019

    • PubMed
    • Google Scholar
37. 1. Lemon RN
    2. Baker SN
    3. Kraskov A

    (2021) Classification of cortical neurons by spike shape and the identification of pyramidal neurons

    Cerebral Cortex 31:5131–5138.

    https://doi.org/10.1093/cercor/bhab147

    • PubMed
    • Google Scholar
38. 1. Łęski S
    2. Lindén H
    3. Tetzlaff T
    4. Pettersen KH
    5. Einevoll GT

    (2013) Frequency dependence of signal power and spatial reach of the local field potential

    PLOS Computational Biology 9:e1003137.

    https://doi.org/10.1371/journal.pcbi.1003137

    • PubMed
    • Google Scholar
39. 1. Leszczyński M
    2. Barczak A
    3. Kajikawa Y
    4. Ulbert I
    5. Falchier AY
    6. Tal I
    7. Haegens S
    8. Melloni L
    9. Knight RT
    10. Schroeder CE

    (2020) Dissociation of broadband high-frequency activity and neuronal firing in the neocortex

    Science Advances 6:eabb0977.

    https://doi.org/10.1126/sciadv.abb0977

    • PubMed
    • Google Scholar
40. 1. Liebe S
    2. Hoerzer GM
    3. Logothetis NK
    4. Rainer G

    (2012) Theta coupling between V4 and prefrontal cortex predicts visual short-term memory performance

    Nature Neuroscience 15:456–462.

    https://doi.org/10.1038/nn.3038

    • PubMed
    • Google Scholar
41. 1. Lindén H
    2. Pettersen KH
    3. Einevoll GT

    (2010) Intrinsic dendritic filtering gives low-pass power spectra of local field potentials

    Journal of Computational Neuroscience 29:423–444.

    https://doi.org/10.1007/s10827-010-0245-4

    • PubMed
    • Google Scholar
42. 1. Lindén H
    2. Tetzlaff T
    3. Potjans TC
    4. Pettersen KH
    5. Grün S
    6. Diesmann M
    7. Einevoll GT

    (2011) Modeling the spatial reach of the LFP

    Neuron 72:859–872.

    https://doi.org/10.1016/j.neuron.2011.11.006

    • PubMed
    • Google Scholar
43. 1. Liu J
    2. Newsome WT

    (2006) Local field potential in cortical area MT: stimulus tuning and behavioral correlations

    The Journal of Neuroscience 26:7779–7790.

    https://doi.org/10.1523/JNEUROSCI.5052-05.2006

    • PubMed
    • Google Scholar
44. 1. Makarova J
    2. Ibarz JM
    3. Makarov VA
    4. Benito N
    5. Herreras O

    (2011) Parallel readout of pathway-specific inputs to laminated brain structures

    Frontiers in Systems Neuroscience 5:77.

    https://doi.org/10.3389/fnsys.2011.00077

    • PubMed
    • Google Scholar
45. 1. Marder E
    2. Goaillard JM

    (2006) Variability, compensation and homeostasis in neuron and network function

    Nature Reviews. Neuroscience 7:563–574.

    https://doi.org/10.1038/nrn1949

    • PubMed
    • Google Scholar
46. 1. Markowitz DA
    2. Wong YT
    3. Gray CM
    4. Pesaran B

    (2011) Optimizing the decoding of movement goals from local field potentials in macaque cortex

    The Journal of Neuroscience 31:18412–18422.

    https://doi.org/10.1523/JNEUROSCI.4165-11.2011

    • PubMed
    • Google Scholar
47. 1. Marques T
    2. Nguyen J
    3. Fioreze G
    4. Petreanu L

    (2018) The functional organization of cortical feedback inputs to primary visual cortex

    Nature Neuroscience 21:757–764.

    https://doi.org/10.1038/s41593-018-0135-z

    • PubMed
    • Google Scholar
48. 1. McGinley MJ
    2. Vinck M
    3. Reimer J
    4. Batista-Brito R
    5. Zagha E
    6. Cadwell CR
    7. Tolias AS
    8. Cardin JA
    9. McCormick DA

    (2015) Waking state: rapid variations modulate neural and behavioral responses

    Neuron 87:1143–1161.

    https://doi.org/10.1016/j.neuron.2015.09.012

    • PubMed
    • Google Scholar
49. 1. Mehring C
    2. Rickert J
    3. Vaadia E
    4. Cardosa de Oliveira S
    5. Aertsen A
    6. Rotter S

    (2003) Inference of hand movements from local field potentials in monkey motor cortex

    Nature Neuroscience 6:1253–1254.

    https://doi.org/10.1038/nn1158

    • PubMed
    • Google Scholar
50. 1. Mehta AD
    2. Ulbert I
    3. Schroeder CE

    (2000a) Intermodal selective attention in monkeys. I: distribution and timing of effects across visual areas

    Cerebral Cortex 10:343–358.

    https://doi.org/10.1093/cercor/10.4.343

    • PubMed
    • Google Scholar
51. 1. Mehta AD
    2. Ulbert I
    3. Schroeder CE

    (2000b) Intermodal selective attention in monkeys. II: physiological mechanisms of modulation

    Cerebral Cortex 10:359–370.

    https://doi.org/10.1093/cercor/10.4.359

    • PubMed
    • Google Scholar
52. 1. Mitzdorf U

    (1985) Current source-density method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomena

    Physiological Reviews 65:37–100.

    https://doi.org/10.1152/physrev.1985.65.1.37

    • PubMed
    • Google Scholar
53. 1. Mitzdorf U

    (1987) Properties of the evoked potential generators: current source-density analysis of visually evoked potentials in the cat cortex

    The International Journal of Neuroscience 33:33–59.

    https://doi.org/10.3109/00207458708985928

    • PubMed
    • Google Scholar
54. 1. Montemurro MA
    2. Rasch MJ
    3. Murayama Y
    4. Logothetis NK
    5. Panzeri S

    (2008) Phase-of-firing coding of natural visual stimuli in primary visual cortex

    Current Biology 18:375–380.

    https://doi.org/10.1016/j.cub.2008.02.023

    • PubMed
    • Google Scholar
55. 1. Mukamel R
    2. Fried I

    (2012) Human intracranial recordings and cognitive neuroscience

    Annual Review of Psychology 63:511–537.

    https://doi.org/10.1146/annurev-psych-120709-145401

    • PubMed
    • Google Scholar
56. 1. Nauhaus I
    2. Busse L
    3. Carandini M
    4. Ringach DL

    (2009) Stimulus contrast modulates functional connectivity in visual cortex

    Nature Neuroscience 12:70–76.

    https://doi.org/10.1038/nn.2232

    • PubMed
    • Google Scholar
57. 1. Ness TV
    2. Remme MWH
    3. Einevoll GT

    (2018) h-type membrane current shapes the local field potential from populations of pyramidal neurons

    The Journal of Neuroscience 38:6011–6024.

    https://doi.org/10.1523/JNEUROSCI.3278-17.2018

    • PubMed
    • Google Scholar
58. 1. Niell CM
    2. Stryker MP

    (2008) Highly selective receptive fields in mouse visual cortex

    The Journal of Neuroscience 28:7520–7536.

    https://doi.org/10.1523/JNEUROSCI.0623-08.2008

    • PubMed
    • Google Scholar
59. 1. Niell CM
    2. Stryker MP

    (2010) Modulation of visual responses by behavioral state in mouse visual cortex

    Neuron 65:472–479.

    https://doi.org/10.1016/j.neuron.2010.01.033

    • PubMed
    • Google Scholar
60. 1. Oh SW
    2. Harris JA
    3. Ng L
    4. Winslow B
    5. Cain N
    6. Mihalas S
    7. Wang Q
    8. Lau C
    9. Kuan L
    10. Henry AM
    11. Mortrud MT
    12. Ouellette B
    13. Nguyen TN
    14. Sorensen SA
    15. Slaughterbeck CR
    16. Wakeman W
    17. Li Y
    18. Feng D
    19. Ho A
    20. Nicholas E
    21. Hirokawa KE
    22. Bohn P
    23. Joines KM
    24. Peng H
    25. Hawrylycz MJ
    26. Phillips JW
    27. Hohmann JG
    28. Wohnoutka P
    29. Gerfen CR
    30. Koch C
    31. Bernard A
    32. Dang C
    33. Jones AR
    34. Zeng H

    (2014) A mesoscale connectome of the mouse brain

    Nature 508:207–214.

    https://doi.org/10.1038/nature13186

    • PubMed
    • Google Scholar
61. 1. O’Leary T

    (2018) Homeostasis, failure of homeostasis and degenerate ion channel regulation

    Current Opinion in Physiology 2:129–138.

    https://doi.org/10.1016/j.cophys.2018.01.006

    • Google Scholar
62. 1. Pesaran B
    2. Pezaris JS
    3. Sahani M
    4. Mitra PP
    5. Andersen RA

    (2002) Temporal structure in neuronal activity during working memory in macaque parietal cortex

    Nature Neuroscience 5:805–811.

    https://doi.org/10.1038/nn890

    • PubMed
    • Google Scholar
63. 1. Pesaran B
    2. Vinck M
    3. Einevoll GT
    4. Sirota A
    5. Fries P
    6. Siegel M
    7. Truccolo W
    8. Schroeder CE
    9. Srinivasan R

    (2018) Investigating large-scale brain dynamics using field potential recordings: analysis and interpretation

    Nature Neuroscience 21:903–919.

    https://doi.org/10.1038/s41593-018-0171-8

    • PubMed
    • Google Scholar
64. 1. Pettersen KH
    2. Devor A
    3. Ulbert I
    4. Dale AM
    5. Einevoll GT

    (2006) Current-source density estimation based on inversion of electrostatic forward solution: effects of finite extent of neuronal activity and conductivity discontinuities

    Journal of Neuroscience Methods 154:116–133.

    https://doi.org/10.1016/j.jneumeth.2005.12.005

    • PubMed
    • Google Scholar
65. 1. Pettersen KH
    2. Einevoll GT

    (2008) Amplitude variability and extracellular low-pass filtering of neuronal spikes

    Biophysical Journal 94:784–802.

    https://doi.org/10.1529/biophysj.107.111179

    • PubMed
    • Google Scholar
66. 1. Pettersen KH
    2. Hagen E
    3. Einevoll GT

    (2008) Estimation of population firing rates and current source densities from laminar electrode recordings

    Journal of Computational Neuroscience 24:291–313.

    https://doi.org/10.1007/s10827-007-0056-4

    • PubMed
    • Google Scholar
67. 1. Prinz AA
    2. Bucher D
    3. Marder E

    (2004) Similar network activity from disparate circuit parameters

    Nature Neuroscience 7:1345–1352.

    https://doi.org/10.1038/nn1352

    • PubMed
    • Google Scholar
68. 1. Rall W

    (1962) Electrophysiology of a dendritic neuron model

    Biophysical Journal 2:145–167.

    https://doi.org/10.1016/s0006-3495(62)86953-7

    • PubMed
    • Google Scholar
69. 1. Rall W
    2. Shepherd GM

    (1968) Theoretical reconstruction of field potentials and dendrodendritic synaptic interactions in olfactory bulb

    Journal of Neurophysiology 31:884–915.

    https://doi.org/10.1152/jn.1968.31.6.884

    • PubMed
    • Google Scholar
70. 1. Reimann MW
    2. Anastassiou CA
    3. Perin R
    4. Hill SL
    5. Markram H
    6. Koch C

    (2013) A biophysically detailed model of neocortical local field potentials predicts the critical role of active membrane currents

    Neuron 79:375–390.

    https://doi.org/10.1016/j.neuron.2013.05.023

    • PubMed
    • Google Scholar
71. 1. Rickert J
    2. Oliveira S de
    3. Vaadia E
    4. Aertsen A
    5. Rotter S
    6. Mehring C

    (2005) Encoding of movement direction in different frequency ranges of motor cortical local field potentials

    The Journal of Neuroscience 25:8815–8824.

    https://doi.org/10.1523/JNEUROSCI.0816-05.2005

    • PubMed
    • Google Scholar
72. Software

    1. Rimehaug AE

    (2023) CINPLA Allen V1 analysis, version swh:1:rev:bc32ab40499ba9b7687a1388f461cdf158bba375

    Software Heritage.

    https://archive.softwareheritage.org/swh:1:dir:2ace4a17f8fdf874a86821913488320804ffcc8a;origin=https://github.com/atleer/CINPLA_Allen_V1_analysis;visit=swh:1:snp:f1b3d0f86648ca18ab6280411c9cc704e015f61f;anchor=swh:1:rev:bc32ab40499ba9b7687a1388f461cdf158bba375
73. 1. Roux S
    2. Mackay WA
    3. Riehle A

    (2006) The pre-movement component of motor cortical local field potentials reflects the level of expectancy

    Behavioural Brain Research 169:335–351.

    https://doi.org/10.1016/j.bbr.2006.02.004

    • PubMed
    • Google Scholar
74. Conference

    1. Rubner Y
    2. Tomasi C
    3. Guibas LJ

    (1998) A metric for distributions with applications to image databases

    IEEE 6th International Conference on Computer Vision. pp. 59–66.

    https://doi.org/10.1109/ICCV.1998.710701

    • Google Scholar
75. 1. Saleem AB
    2. Lien AD
    3. Krumin M
    4. Haider B
    5. Rosón MR
    6. Ayaz A
    7. Reinhold K
    8. Busse L
    9. Carandini M
    10. Harris KD
    11. Carandini M

    (2017) Subcortical source and modulation of the narrowband gamma oscillation in mouse visual cortex

    Neuron 93:315–322.

    https://doi.org/10.1016/j.neuron.2016.12.028

    • PubMed
    • Google Scholar
76. 1. Scherberger H
    2. Jarvis MR
    3. Andersen RA

    (2005) Cortical local field potential encodes movement intentions in the posterior parietal cortex

    Neuron 46:347–354.

    https://doi.org/10.1016/j.neuron.2005.03.004

    • PubMed
    • Google Scholar
77. 1. Schomburg EW
    2. Anastassiou CA
    3. Buzsáki G
    4. Koch C

    (2012) The spiking component of oscillatory extracellular potentials in the rat hippocampus

    The Journal of Neuroscience 32:11798–11811.

    https://doi.org/10.1523/JNEUROSCI.0656-12.2012

    • PubMed
    • Google Scholar
78. 1. Schroeder CE
    2. Mehta AD
    3. Givre SJ

    (1998) A spatiotemporal profile of visual system activation revealed by current source density analysis in the awake macaque

    Cerebral Cortex 8:575–592.

    https://doi.org/10.1093/cercor/8.7.575

    • PubMed
    • Google Scholar
79. 1. Senzai Y
    2. Fernandez-Ruiz A
    3. Buzsáki G

    (2019) Layer-specific physiological features and interlaminar interactions in the primary visual cortex of the mouse

    Neuron 101:500–513.

    https://doi.org/10.1016/j.neuron.2018.12.009

    • PubMed
    • Google Scholar
80. Preprint

    1. Shen S
    2. Jiang X
    3. Scala F
    4. Fu J
    5. Fahey P
    6. Kobak D
    7. Tan Z
    8. Reimer J
    9. Sinz F
    10. Tolias AS

    (2020) Distinctorganization of two cortico-cortical feedback pathways

    bioRxiv.

    https://doi.org/10.1101/2020.02.27.968792

    • Google Scholar
81. 1. Siegle JH
    2. Jia X
    3. Durand S
    4. Gale S
    5. Bennett C
    6. Graddis N
    7. Heller G
    8. Ramirez TK
    9. Choi H
    10. Luviano JA
    11. Groblewski PA
    12. Ahmed R
    13. Arkhipov A
    14. Bernard A
    15. Billeh YN
    16. Brown D
    17. Buice MA
    18. Cain N
    19. Caldejon S
    20. Casal L
    21. Cho A
    22. Chvilicek M
    23. Cox TC
    24. Dai K
    25. Denman DJ
    26. de Vries SEJ
    27. Dietzman R
    28. Esposito L
    29. Farrell C
    30. Feng D
    31. Galbraith J
    32. Garrett M
    33. Gelfand EC
    34. Hancock N
    35. Harris JA
    36. Howard R
    37. Hu B
    38. Hytnen R
    39. Iyer R
    40. Jessett E
    41. Johnson K
    42. Kato I
    43. Kiggins J
    44. Lambert S
    45. Lecoq J
    46. Ledochowitsch P
    47. Lee JH
    48. Leon A
    49. Li Y
    50. Liang E
    51. Long F
    52. Mace K
    53. Melchior J
    54. Millman D
    55. Mollenkopf T
    56. Nayan C
    57. Ng L
    58. Ngo K
    59. Nguyen T
    60. Nicovich PR
    61. North K
    62. Ocker GK
    63. Ollerenshaw D
    64. Oliver M
    65. Pachitariu M
    66. Perkins J
    67. Reding M
    68. Reid D
    69. Robertson M
    70. Ronellenfitch K
    71. Seid S
    72. Slaughterbeck C
    73. Stoecklin M
    74. Sullivan D
    75. Sutton B
    76. Swapp J
    77. Thompson C
    78. Turner K
    79. Wakeman W
    80. Whitesell JD
    81. Williams D
    82. Williford A
    83. Young R
    84. Zeng H
    85. Naylor S
    86. Phillips JW
    87. Reid RC
    88. Mihalas S
    89. Olsen SR
    90. Koch C

    (2021) Survey of spiking in the mouse visual system reveals functional hierarchy

    Nature 592:86–92.

    https://doi.org/10.1038/s41586-020-03171-x

    • PubMed
    • Google Scholar
82. 1. Sinha M
    2. Narayanan R

    (2015) HCN channels enhance spike phase coherence and regulate the phase of spikes and LFPs in the theta-frequency range

    PNAS 112:E2207–E2216.

    https://doi.org/10.1073/pnas.1419017112

    • Google Scholar
83. 1. Sinha M
    2. Narayanan R

    (2022) Active dendrites and local field potentials: biophysical mechanisms and computational explorations

    Neuroscience 489:111–142.

    https://doi.org/10.1016/j.neuroscience.2021.08.035

    • PubMed
    • Google Scholar
84. 1. Stavisky SD
    2. Kao JC
    3. Nuyujukian P
    4. Ryu SI
    5. Shenoy KV

    (2015) A high performing brain-machine interface driven by low-frequency local field potentials alone and together with spikes

    Journal of Neural Engineering 12:036009.

    https://doi.org/10.1088/1741-2560/12/3/036009

    • PubMed
    • Google Scholar
85. 1. Swadlow HA
    2. Gusev AG
    3. Bezdudnaya T

    (2002) Activation of a cortical column by a thalamocortical impulse

    The Journal of Neuroscience 22:7766–7773.

    https://doi.org/10.1523/JNEUROSCI.22-17-07766.2002

    • PubMed
    • Google Scholar
86. 1. Teleńczuk B
    2. Dehghani N
    3. Le Van Quyen M
    4. Cash SS
    5. Halgren E
    6. Hatsopoulos NG
    7. Destexhe A

    (2017) Local field potentials primarily reflect inhibitory neuron activity in human and monkey cortex

    Scientific Reports 7:40211.

    https://doi.org/10.1038/srep40211

    • PubMed
    • Google Scholar
87. 1. Telenczuk B
    2. Telenczuk M
    3. Destexhe A

    (2020) A kernel-based method to calculate local field potentials from networks of spiking neurons

    Journal of Neuroscience Methods 344:108871.

    https://doi.org/10.1016/j.jneumeth.2020.108871

    • PubMed
    • Google Scholar
88. 1. Tomsett RJ
    2. Ainsworth M
    3. Thiele A
    4. Sanayei M
    5. Chen X
    6. Gieselmann MA
    7. Whittington MA
    8. Cunningham MO
    9. Kaiser M

    (2015) Virtual Electrode Recording Tool for EXtracellular potentials (VERTEX): comparing multi-electrode recordings from simulated and biological mammalian cortical tissue

    Brain Structure & Function 220:2333–2353.

    https://doi.org/10.1007/s00429-014-0793-x

    • PubMed
    • Google Scholar
89. 1. Tort ABL
    2. Kramer MA
    3. Thorn C
    4. Gibson DJ
    5. Kubota Y
    6. Graybiel AM
    7. Kopell NJ

    (2008) Dynamic cross-frequency couplings of local field potential oscillations in rat striatum and hippocampus during performance of a T-maze task

    PNAS 105:20517–20522.

    https://doi.org/10.1073/pnas.0810524105

    • PubMed
    • Google Scholar
90. 1. Traub RD
    2. Contreras D
    3. Cunningham MO
    4. Murray H
    5. LeBeau FEN
    6. Roopun A
    7. Bibbig A
    8. Wilent WB
    9. Higley MJ
    10. Whittington MA

    (2005) Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts

    Journal of Neurophysiology 93:2194–2232.

    https://doi.org/10.1152/jn.00983.2004

    • PubMed
    • Google Scholar
91. 1. Victor JD
    2. Purpura K
    3. Katz E
    4. Mao B

    (1994) Population encoding of spatial frequency, orientation, and color in macaque V1

    Journal of Neurophysiology 72:2151–2166.

    https://doi.org/10.1152/jn.1994.72.5.2151

    • PubMed
    • Google Scholar
92. 1. Vierling-Claassen D
    2. Cardin JA
    3. Moore CI
    4. Jones SR

    (2010) Computational modeling of distinct neocortical oscillations driven by cell-type selective optogenetic drive: separable resonant circuits controlled by low-threshold spiking and fast-spiking interneurons

    Frontiers in Human Neuroscience 4:198.

    https://doi.org/10.3389/fnhum.2010.00198

    • PubMed
    • Google Scholar
93. 1. Vigneswaran G
    2. Kraskov A
    3. Lemon RN

    (2011) Large identified pyramidal cells in macaque motor and premotor cortex exhibit “thin spikes”: implications for cell type classification

    The Journal of Neuroscience 31:14235–14242.

    https://doi.org/10.1523/JNEUROSCI.3142-11.2011

    • PubMed
    • Google Scholar
94. 1. Vinck M
    2. Batista-Brito R
    3. Knoblich U
    4. Cardin JA

    (2015) Arousal and locomotion make distinct contributions to cortical activity patterns and visual encoding

    Neuron 86:740–754.

    https://doi.org/10.1016/j.neuron.2015.03.028

    • PubMed
    • Google Scholar
95. 1. Watrous AJ
    2. Lee DJ
    3. Izadi A
    4. Gurkoff GG
    5. Shahlaie K
    6. Ekstrom AD

    (2013) A comparative study of human and rat hippocampal low-frequency oscillations during spatial navigation

    Hippocampus 23:656–661.

    https://doi.org/10.1002/hipo.22124

    • PubMed
    • Google Scholar
96. 1. Womelsdorf T
    2. Fries P
    3. Mitra PP
    4. Desimone R

    (2006) Gamma-band synchronization in visual cortex predicts speed of change detection

    Nature 439:733–736.

    https://doi.org/10.1038/nature04258

    • PubMed
    • Google Scholar

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.

[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:

1. Comparison of the moments of distributions of peak firing rates across cells in the evoked response

2. Comparison of the moments of distributions of latencies to peak firing rate across cells.

3. Comparing moments of the relative magnitude of firing rates across laminar populations

4. 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 4^th 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.

Author response image 1

Download asset Open asset

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:

1. Comparison of the moments of distributions of peak firing rates across cells in the evoked response

2. Comparison of the moments of distributions of latencies to peak firing rate across cells.

3. Comparing moments of the relative magnitude of firing rates across laminar populations

4. 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 4^th 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:

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

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

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

Author details

1. Atle E Rimehaug

   Department of Informatics, University of Oslo, Oslo, Norway

   Contribution

   Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

   For correspondence

   atleeri@ifi.uio.no

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-8312-9875

2. Alexander J Stasik

   Department of Physics, University of Oslo, Oslo, Norway

   Contribution

   Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-1646-2472

3. Espen Hagen

   1. Department of Physics, University of Oslo, Oslo, Norway
   2. Department of Data Science, Norwegian University of Life Sciences, Ås, Norway

   Contribution

   Conceptualization, Resources, Software, Formal analysis, Investigation, Methodology, Writing – review and editing

   Competing interests

   No competing interests declared

4. Yazan N Billeh

   MindScope Program, Allen Institute, Seattle, United States

   Contribution

   Conceptualization, Software, Methodology

   Competing interests

   No competing interests declared

5. Josh H Siegle

   MindScope Program, Allen Institute, Seattle, United States

   Contribution

   Data curation, Investigation, Visualization, Writing – review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-7736-4844

6. Kael Dai

   MindScope Program, Allen Institute, Seattle, United States

   Contribution

   Resources, Software, Methodology

   Competing interests

   No competing interests declared

7. Shawn R Olsen

   MindScope Program, Allen Institute, Seattle, United States

   Contribution

   Resources, Supervision

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-9568-7057

8. Christof Koch

   MindScope Program, Allen Institute, Seattle, United States

   Contribution

   Conceptualization, Resources, Supervision, Funding acquisition, Methodology, Project administration, Writing – review and editing

   Competing interests

   No competing interests declared

9. Gaute T Einevoll

   1. Department of Physics, University of Oslo, Oslo, Norway
   2. Department of Physics, Norwegian University of Life Sciences, Ås, Norway

   Contribution

   Conceptualization, Resources, Supervision, Funding acquisition, Methodology, Project administration, Writing – review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-5425-5012

10. Anton Arkhipov

    MindScope Program, Allen Institute, Seattle, United States

    Contribution

    Conceptualization, Resources, Supervision, Funding acquisition, Methodology, Project administration, Writing – review and editing

    For correspondence

    antona@alleninstitute.org

    Competing interests

    No competing interests declared

    Additional information

    Lead contact

    "This ORCID iD identifies the author of this article:" 0000-0003-1106-8310

• 872

  Page views

• 154

  Downloads

• 3

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.