We construct a probabilistic decoder, represented by a population of neurons, that reports or estimates the identity of a stimulus from the spiking response of an input population of neurons. We assume an input population with responses that are a function of the stimulus, f(θ_k), plus Gaussian noise, f(θ_k), and the covariance (Q) is equal for all stimulus conditions (Q) such that Q = Q_k = Q_i. Then, the posterior distribution can be written as

where a multivariate Gaussian is

This can be expanded and simplified such that

where

Here, w are the weights over k for each neuron in the decoder population and β is a constant term for each k. Importantly, because we assume Gaussian input, with this formulation, w and β are derived in closed form. More generally, w and β, can be estimated numerically using multinomial logistic regression and this form remains optimal for any input population statistics within the exponential family (e.g., Poisson noise).

To generate input populations (P_IN), we simulated N neurons responding to a stimulus characterized by orientation (θ_k ε [−π/2:π/2K:π/2]). The response of each neuron, r_i, depends on a tuning function, f_i(θ), and an additive noise term, ε_i, describing trial-to-trial variability. Noise is correlated across the population, generated from a multivariate Gaussian distribution with zero mean and covariance C. Orientation tuning functions were defined as:

Here, α is the baseline firing rate, β scales the tuned response, κ scales the tuning bandwidth, and Φ is the orientation preference of each neuron. For homogeneous P_IN all parameters except Φ were fixed: (α, β, κ) = (0, 5, 4). For heterogenous P_IN, we sampled parameters to match measurements from macaque V1 (Ringach et al., 2002) and our ferret V1 data. Tuning bandwidth was generated by converting half-width at 1/√2 height (γ) values from a lognormal distribution (μ = −1, σ = 0.6):

Limited-range correlations were included so neural noise correlation depends on tuning preference difference (Ecker et al., 2011). A correlation matrix, C, was specified by the difference between preferred orientations of neurons and the maximum pairwise correlation, c_o:

where δ is the circular difference and

where I is the identity matrix of size N. We scaled the correlation matrix by the mean firing rate of each neuron to produce Poisson-like noise (Ecker et al., 2011).

Derived weights for a given P_IN were artificially smoothed using the following equation from Park and Pillow (2011):

Here, S⁺ is the pseudoinverse of S, δ(Φ_i−Φ_j) is the circular difference between preferred orientations of neurons, ρ₁ scales the amplitude of smoothing, and ρ₂ scales the functional range of smoothing.

Decoding accuracy was calculated with the mean-squared-error of the maximum a posterior probability (MAP) estimate across t simulated trials of each stimulus (k):

Here, w_k are the weights for a given decoder neuron and θ_k is the true stimulus.

Briefly, female ferrets aged P18–23 (Marshall Farms) were anesthetized with isoflurane (delivered in O2). Atropine was administered and a 1:1 mixture of lidocaine and bupivacaine was administered SQ. Animals were maintained at an internal temperature of 37°C. Under sterile surgical conditions, a small craniotomy (0.8 mm diameter) was made over the visual cortex (7–8 mm lateral and 2–3 mm anterior to lambda). A mixture of diluted AAV1.hSyn.Cre (1:25,000–1:50,000) and AAV1.Syn.FLEX.GCaMP6s (UPenn) was injected (125–202.5 nl) through beveled glass micropipettes (10–15 micron outer diameter) at 600, 400, and 200 microns below the pia. Finally, the craniotomy was filled with sterile agarose (Type IIIa, Sigma-Aldrich) and the incision site was sutured.

After 3–5 weeks of expression, ferrets were anesthetized with 50 mg/kg ketamine and isoflurane. Atropine and bupivacaine were administered, animals were placed on a feedback-controlled heatingpad to maintain an internal temperature of 37°C, and intubated to be artificially respirated. Isoflurane was delivered throughout the surgical procedure to maintain a surgical plane of anesthesia. An intravenous cannula was placed to deliver fluids. Tidal CO2, external temperature, and internal temperature were continuously monitored. The scalp was retracted and a custom titanium headplate adhered to the skull (Metabond, Parkell). A craniotomy was performed and the dura retracted to reveal the cortex. One piece of custom cover-glass (3 mm diameter, 0.7 mm thickness, Warner Instruments) adhered using optical adhesive (71, Norland Products) to custom insert was placed onto the brain to dampen biological motion during imaging. A 1:1 mixture of tropicamide ophthalmic solution (Akorn) and phenylephrine hydrochloride ophthalmic solution (Akorn) was applied to both eyes to dilate the pupils and retract the nictating membranes. Contact lenses were inserted to protect the eyes. Upon completion of the surgical procedure, isoflurane was gradually reduced and pancuronium (2 mg/kg/h) was delivered IV.

Visual stimuli were generated using Psychopy (Peirce, 2007). The monitor was placed 25 cm from the animal. Receptive field locations for each cell were hand mapped and the spatial frequency optimized (range: 0.04–0.25 cpd). For each soma and dendritic segment, square-wave drifting gratings were presented at 22.5 degree increments (2 s duration, 1 s ISI, 8–10 trials for each field of view).

Two-photon imaging was performed on a Bergamo II microscope (Thorlabs) running Scanimage (Pologruto et al., 2003; Vidrio Technologies) with 940 nm dispersion-compensated excitation provided by an Insight DS+ (Spectraphysics). For spine and axon imaging, power after the objective was limited to <50 mW. Cells were selected for imaging on the basis of their position relative to large blood vessels, responsiveness to visual stimulation, and lack of prolonged calcium transients resulting from over-expression of GCaMP6s. Images were collected at 30 Hz using bidirectional scanning with 512 × 512 pixel resolution or with custom ROIs (frame rate range: 22–50 Hz). Somatic imaging was performed with a resolution of 2.05–10.24 pixels/micron. Dendritic spine imaging was performed with a resolution of 6.10–15.36 pixels/micron.

Imaging data were excluded from analysis if motion along the z-axis was detected. Dendrite images were corrected for in-plane motion via a 2D cross-correlation-based approach in MATLAB or using a piecewise non-rigid motion correction algorithm (Pnevmatikakis and Giovannucci, 2017). ROIs (region of interest) were drawn in ImageJ; dendritic ROIs spanned contiguous dendritic segments and spine ROIs were fit with custom software. Mean pixel values for ROIs were computed over the imaging time series and imported into MATLAB (Sage et al., 2012; Hiner et al., 2017). ΔF/F_o was computed by computing F_o with time-averaged median or percentile filter (10th percentile). For spine signals, we subtracted a scaled version of the dendritic signal to remove back-propagating action potentials as performed previously (Wilson et al., 2016). ΔF/F_o traces were synchronized to stimulus triggers sent from Psychopy and collected by Spike2. Spines were included for analysis if the SNR of the preferred response exceeded two median absolute deviations above the baseline noise (measured during the blank) and were weakly correlated with the dendritic signal (Spearman’s correlation, r < 0.4). Some spine traces contained negative events after subtraction, so correlations were computed ignoring negative values. We then normalized each spine’s responses so that each spine had equal weight. The preferred orientation for each spine was calculated by fitting responses with a Gaussian tuning curve using lsqcurvefit (Matlab). Tuning selectivity was measured as the vector strength index (v) for each neuron’s response:

Here r is the mean responses over the orientations (θ_k) presented for each spine (i). Note, this same index is used to characterize simulated input selectivity.

To compare input tuning (derived synaptic population or measured dendritic spine population) with output tuning (downstream readout or measured somatic tuning) we computed the Pearson Correlation coefficient (Matlab). This correlation was computed on trial-averaged responses across different orientations. For dendritic spines and soma, measured responses across stimulus presentation trials were averaged. For simulated synaptic populations and corresponding downstream readout neurons, we simulated trials by adding noise to each synaptic tuning curve.

A categorical probabilistic decoder reports the probability that a particular stimulus orientation, θ_k, was present given the spiking responses of an input population, R. This can be expressed as a normalized exponential function of the log-likelihood plus the log prior for each θ_k,

where

The likelihood, p(R|θ_k), is the probability of the observed responses in an input population given the stimulus k and p(θ_k) is the prior probability of that stimulus class. If p(R|θ_k) is in the exponential family, then L(θ_k) can be written as a weighted sum of the input population response vector plus an offset, which can be estimated numerically via multinomial logistic regression (Ma et al., 2006). For simplicity, we assume the input population has a response that is a function of the stimulus plus Gaussian noise and equal covariance across all stimulus conditions. Although this assumption about the covariance structure deviates from real neural activity, this assumption means the weights and offset can be solved analytically (see “Methods” Section), and as will be shown below, such a simple model makes substantial headway in explaining biological phenomena. Our goal here is to provide a plausible alternative to “somatic selectivity” for the connectivity rules in the cortex. Under the Gaussian assumption, the decoder amounts to:

Here, f(θ_k) is the mean input population response to stimulus orientation, k, K is the total number of orientations, and Q is the covariance matrix. The covariance term captures the influence of each neuron’s response variance (diagonal elements) and the variability shared with other neurons (off-diagonal elements). Intuitively, in the absence of covariability (i.e., off-diagonal elements are zero), the weights are proportional to the signal-to-noise ratio of the neuron (the mean divided by the variance). The term R^Tw_k is the dot product between the population response and weights. The second term, β_k, is an offset for each stimulus. ln p(θ_k) is a constant reflecting the log prior probability of stimulus k. It is worth noting that if the covariance depends on the stimulus, the optimal readout is no longer a linear function of R and is quadratic, which can be interpreted as a complex-cell (Pagan et al., 2016; Jaini and Burge, 2017) and is a potentially fruitful future direction.

In this study, we focus on the weights of this simple Gaussian, equal covariance decoder in order to examine how synaptic tuning from such a simple decoder would arise. Because the optimal weights have an analytic solution (Equation 1), we can see how they depend on the parameters of P_IN. The simplifying assumptions we use to derive the maximum-likelihood weights help build intuitions about what can be expected in biological circuits, and linear weights such as these could be learned by real neural systems (Dayan and Abbott, 2001). A key difference here from prior work is that rather than focus on discrimination (Haefner et al., 2013), we treat orientation estimation as a multiclass identification problem, discretizing θ such that for each possible θ_k, there is a separate weight vector. Thus, in this derivation, the optimal weights depend on the tuning curves themselves, not the derivative.

To understand how synaptic weights depend on input statistics, we derived maximum-likelihood weights for input populations, P_IN, with different tuning and covariance. To generate P_IN we simulated N neurons responding to K oriented stimuli (θ=[−90^o:K/180:+90^o]). We briefly describe the construction of P_IN here (full details are described in the “Methods” Section). Each neuron is defined by a tuning function and noise term, describing trial-by-trial variability, which are summed to generate stimulus-driven responses. We compared two fundamentally different types of input populations that have been used in the literature, homogeneous and heterogenous, as well as the role of correlated variability in shaping readout weights. A homogeneous P_IN consists of shifted copies of a single tuning curve (Figure 2A). Heterogenous P_IN have diverse tuning functions and were generated to match measurements from macaque V1 (Ringach et al., 2002). The heterogenous P_IN consisted of tuning curves closely resembling V1 physiology in terms of the variation in peak firing rate, bandwidth, and baseline firing rate (Figure 2D). Varying amounts of limited-range correlations were included such that the noise correlation between two neurons depends on the difference in their tuning preferences (Ecker et al., 2011; Kohn et al., 2016).

The statistics of P_IN responses, R, will dictate the weight structure for neurons in a decoding population. For a homogeneous P_IN, the weights are smooth across orientation space and exhibit three primary features: a prominent peak about the preferred orientation of the output tuning, slight negative weights for orientations just outside the preferred, and near-zero weights at orthogonal orientations (Figure 2B). With more realistic tuning diversity (heterogenous P_IN), optimal weights are no longer smooth (Figure 2E). While the optimal weights appear to roughly have the same overall shape as for homogeneous P_IN, there is considerable positive and negative weighting across orientation space. Despite substantial changes in optimal weight vectors, the decoder output (i.e., somatic response) tuning was narrow (Figure 2F), similar to the output for the homogeneous case (Figure 2C).

To explore the importance of decoding weight diversity, we imposed a smoothing penalty on weight vectors (Park and Pillow, 2011). We calculated cross-validated decoder accuracy using the mean-squared error between the maximum a posteriori estimation and true stimulus (see “Methods” Section). Different degrees of smoothing are shown for an example set of weights in Figure 3A. We simulated a range of population sizes (n = 2–2,048) and correlations (c_o= 0, 0.25, 0.50).

Without noise correlations, the accuracy of all decoders increases with population size, with a homogenous P_IN performing best (Figure 3B). In the presence of noise correlations, accuracy saturates for large homogenous P_IN (Figures 3C,D). As previously shown (Ecker et al., 2011), accuracy for heterogeneous populations with limited-range correlations does not saturate (Figures 3C,D). However, this depends on weight diversity. Smoothing the weights for heterogeneous P_IN caused saturation and decreased accuracy (Figures 3C,D), demonstrating that weight amplitude diversity in analytically derived weights distributions are critical for the decoder performance.

In order to compare analytically derived weights with the synaptic inputs onto V1 neurons measured in vivo, we generated excitatory synaptic input populations (P_SYN). Under the assumption that synaptic integration is linear, two synapses of equal weight are the same as one synapse with double that weight. This creates a degeneracy where synapse count and size trade-off. Because current spine imaging techniques typically capture large synapses and there is no relationship between strength and orientation preference (Scholl et al., 2021), we can assume size is fixed and convert the derived weights into a frequency distribution of “synaptic inputs” (Figure 4A). The tuning curve for such a synapse is the tuning curve of the input and thus, the synaptic input population, P_SYN, is the input population resampled with probabilities given by the derived weights. An example P_SYN for a single decoder neuron is shown in Figure 4B (drawn from the heterogeneous P_IN in Figure 2). P_SYN in this example displays some specificity in orientation tuning relative to the somatic output, indicated by a larger proportion of simulated synapses with similar orientation preference as the somatic output (0^o) of the decoder neuron.

We analyzed two-photon calcium recordings from soma and corresponding dendritic spines on individual neurons in ferret V1 during the presentation of oriented drifting gratings (see “Methods” Section). While our model draws from a P_IN matched to measurements from macaque V1, the orientation tuning of layer 2/3 neurons in ferret V1, as measured by two-photon cellular imaging, exhibit a similar range in selectivity (Wilson et al., 2016). Visually responsive and isolated dendritic spines (see “Methods” Section) typically exhibit diverse orientation tuning relative to the somatic output, although some individual cells show greater overall diversity (Figure 5B) than others (Figure 5A). To characterize P_SYN diversity, both for real dendritic spines and simulated inputs, we computed the Pearson correlation coefficient between the tuning curves of individual inputs and the corresponding somatic output (Scholl et al., 2021). For these comparisons, we sampled orientation space in the model spines to match our empirical measurements (22.5 deg increments), and the number of total excitatory inputs recovered for each simulated downstream neuron was set to 100, similar to the average number of visually-responsive spines recorded for each ferret V1 neuron (n = 45, n = 158.9 ± 73.2 spines/cell). Simulations were run 10,000 times, with n = 1,000 for P_IN and c_o= 0.20.

Across all simulated inputs, input-output tuning correlation was higher for homogeneous P_SYN compared to heterogeneous P_SYN (median r_hom = 0.60, median r_hom = 0.18, n = 900,000; Figure 5C). Tuning correlation between all imaged dendritic spines and soma was low (median r_cell = 0.31, n = 7,151 spines from 45 cells), more closely resembling our model with a heterogenous P_IN. As somatic orientation selectivity (i.e., tuning bandwidth) varies for single cells in ferret V1 (Goris et al., 2015; Wilson et al., 2016), we next examined the average input-output tuning correlation across individual cells (Figure 5D). Here, the homogeneous model exhibited greater specificity then the heterogeneous model (median r_hom = 0.52; median r_het = 0.18, n = 90,000). For ferret V1 cells, we observed similar spine-soma correlation as the heterogeneous simulation (median r_cell = 0.20, n = 45). Ferret V1 cells were not statistically different from neural decoders with a heterogeneous P_SYN (p = 0.19, Mann-Whitney test), while neural decoders with a homogeneous P_SYN were significantly more correlated with the inputs (p < 0.0001, Mann-Whitney test). A small percentage of imaged cells (11.1%, n = 5/45) had synaptic populations whose mean tuning correlation was within the 95% confidence interval of the homogeneous model distribution. This observation may be influenced by contamination from back-propagating action potentials, the limited sensory space probed, or indicate a true functional specialization of a subpopulation of excitatory neurons. Additionally, some cells had negative average correlations with their spines, which never occurred in the models—potentially indicating linearities between the spines and soma. It is also important to emphasize that both synaptic populations and the heterogeneous model exhibit a positive bias in tuning correlations, illustrating that while inputs are functionally diverse, they are, on average, more similarly tuned to the cell/decoder output.

Given the differences between ferret V1 neurons, we quantified the degree to which synaptic populations on each neuron matched tuning correlation distributions from models of homogenous and heterogenous P_IN, by calculating the Kullback-Leibler divergence (D_KL, bin size = 0.05, see “Methods” Section). Across our population, imaged neurons more closely resembled simulations with heterogenous, compared to homogenous, P_IN (93.3%, n = 42/45; Figure 5E) and D_KL from a heterogenous model was consistently larger (p < 0.0001, sign rank Wilcoxon test). This trend held for a range of histogram bin sizes (0.001–0.20). Importantly, the models are not fit to data. They are derived entirely from the statistics of the input population, so this correspondence between the heterogenous model and the data results from no free parameters.

In addition to the similarity in input-output tuning correlation, we observed several trends predicted by the heterogeneous model that were evident in synaptic populations imaged in vivo. Simulated excitatory inputs correlated with the decoder output were not more selective for orientation (see “Methods” Section; bootstrapped PCA slope = 0.001 ± 0.004 s.e., n = 10,000 simulation runs). For two-photon data, a minuscule, but significant, trend was evident (bootstrapped PCA slope = 0.03 ± 0.1 s.e., n = 7,151). So while selective inputs are proposed to provide more information about encoded stimulus variables (Seriès et al., 2004; Shamir and Sompolinsky, 2006; Zavitz and Price, 2019) and unselective (or poorly selective) inputs could convey information through their covariance with selective neurons (Zylberberg, 2017), our model and experimental data suggest co-tuned and orthogonally-tuned inputs exhibit a wide range of tuning selectivity. Response variability (i.e., standard deviation) across trials for simulated excitatory inputs was significantly smaller for “null” orientations (± 90 deg) than at the “preferred” (median = 0.30 and IQR = 0.14, median = 0.38 and IQR = 0.22, respectively; p < 0.001, Wilcoxon rank-sum test) and this trend was also observed in our two-photon data (null: median = 0.13 and IQR = 0.15; preferred: median = 0.23 and IQR = 0.31, respectively; p < 0.001, Wilcoxon rank-sum test). As both modeled and imaged neurons had “null”-tuned excitatory inputs that exhibited less response variability, these inputs may carry useful information about when the preferred stimulus is not present.

Taken together, our decoding framework with realistic (i.e., heterogenous orientation tuning), noisy input populations suggests the collection of orthogonally-tuned excitatory inputs in cortical neurons in vivo is not unexpected. Instead, the synaptic architecture of layer 2/3 neurons in the ferret visual cortex is likely optimized for the readout of upstream populations tuned to orientation. As functional heterogeneity and correlated noise are ubiquitous in sensory populations (Jia et al., 2010, 2011; Chen et al., 2011, 2013; Wertz et al., 2015; Wilson et al., 2016, 2018; Iacaruso et al., 2017; Scholl et al., 2017; Kerlin et al., 2019; Ju et al., 2020), the fact that simple linear model utilizes orthogonally-tuned excitatory inputs suggests this is a general principle across mammal species.