Abstract
The spatial organization and dynamic interactions between excitatory and inhibitory synaptic inputs that define the receptive field (RF) of simple cells in the cat primary visual cortex (V1) still raise the following paradoxical issues: (1) stimulation of simple cells in V1 with drifting gratings supports a wiring schema of spatially segregated sets of excitatory and inhibitory inputs activated in an opponent way by stimulus contrast polarity and (2) in contrast, intracellular studies using flashed bars suggest that although ON and OFF excitatory inputs are indeed segregated, inhibitory inputs span the entire RF regardless of input contrast polarity. Here, we propose a biologically detailed computational model of simple cells embedded in a V1-like network that resolves this seeming contradiction. We varied parametrically the RF-correlation-based bias for excitatory and inhibitory synapses and found that a moderate bias of excitatory neurons to synapse onto other neurons with correlated receptive fields and a weaker bias of inhibitory neurons to synapse onto other neurons with anticorrelated receptive fields can explain the conductance input, the postsynaptic membrane potential, and the spike train dynamics under both stimulation paradigms. This computational study shows that the same structural model can reproduce the functional diversity of visual processing observed during different visual contexts.
SIGNIFICANCE STATEMENT Identifying generic connectivity motives in cortical circuitry encoding for specific functions is crucial for understanding the computations implemented in the cortex. Indirect evidence points to correlation-based biases in the connectivity pattern in V1 of higher mammals, whereby excitatory and inhibitory neurons preferentially synapse onto neurons respectively with correlated and anticorrelated receptive fields. A recent intracellular study questions this push–pull hypothesis, failing to find spatial anticorrelation patterns between excitation and inhibition across the receptive field. We present here a spiking model of V1 that integrates relevant anatomic and physiological constraints and shows that a more versatile motif of correlation-based connectivity with selectively tuned excitation and broadened inhibition is sufficient to account for the diversity of functional descriptions obtained for different classes of stimuli.
Introduction
Hubel and Wiesel (1962) hypothesized that orientation selectivity emerges because of specific alignment of feedforward excitatory ON/OFF inputs from the lateral geniculate nucleus (LGN) onto contrast-matched receptive field (RF) subfields of layer 4 (L4) neurons in the primary visual cortex (V1). Simultaneous recordings of LGN and V1 connected pairs provided support for this wiring rule (Toyama et al., 1981; Alonso et al., 2001; Sedigh-Sarvestani, 2017). The principles of intracortical connectivity are less clear. In cat layer 4 (L4) simple cells, intracellular recordings revealed the characteristic push–pull behavior, whereby presentation of a sign-matched stimulus in the RFs subfield evokes depolarization (push), whereas stimulus of opposite polarity evokes hyperpolarization (pull) of the membrane potential (Ferster, 1986; Hirsch et al., 1998). These observations lead to the hypothesis of an anticorrelated arrangement of excitation versus inhibition across the RF space, whereby stimulation of a RF subfield with a sign-matched stimulus evokes predominantly excitation, whereas opposite stimulus polarity evokes predominantly inhibition (Troyer et al., 1998; Ferster and Miller, 2000; Hirsch, 2003).
A way for such an anticorrelated interaction between excitation and inhibition across the RF to be implemented is for the connectivity in L4 cells to adhere to a push–pull schema; excitatory and inhibitory neurons preferentially connect other neurons with, respectively, correlated and anticorrelated RFs (Troyer et al., 1998; Anderson et al., 2000; Ferster and Miller, 2000; Hirsch, 2003). Such push–pull connectivity has been examined in numerous computational models, either on its own (Troyer et al., 1998; Ferster and Miller, 2000; Kremkow et al., 2016b) or in combination with untuned inhibition (Lauritzen and Miller, 2003; Teich and Qian, 2006).
Both extracellular and intracellular studies support a correlation-based rule for excitatory connectivity showing that co-oriented pairs of V1 cells have a higher probability of connection (Michalski et al., 1983; Monier et al., 2003; Ko et al., 2011; Denman and Contreras, 2014; Lee and Bonin et al., 2016; Wilson et al., 2016), and that excitatory neurons preferentially connect other neurons with correlated (in phase) RFs (Cossell et al., 2015).
The anticorrelation rule for inhibition would predict maximal inhibitory synaptic inputs from inhibitory cells; the RF of these cells is cotuned in orientation but in antiphase. Such an antiphase arrangement of excitation and inhibition is supported by intracellular studies showing that drifting grating elicits excitatory and inhibitory inputs that are in antiphase (Ferster, 1988; Anderson et al., 2000; Monier et al., 2003, 2008; Tan et al., 2011). Furthermore, inhibitory inputs on excitatory cells in the ferret layer 2/3 are biased toward opposite direction of movement (Wilson et al., 2018).
However, a recent study revealed that when the RFs of simple cells are mapped with optimally oriented flashed bars, inhibition is evoked broadly and is colocalized across the RF regardless of stimulus polarity (Taylor et al., 2018), questioning the idea of a strict antiphase arrangement of excitation and inhibition across the space of RF. Models using push–pull connectivity (whether alone or in combination with untuned inhibition) do not account for these findings as they predict the inhibition to elicit spatially offset peaks across the RF field for bars of opposite polarity. Models relying solely on phase-unspecific inhibition (Ben-Yishai et al., 1995; Douglas et al., 1995; Carandini and Ringach, 1997) cannot explain the antiphase relationship because of sinusoidal grating stimuli. Thus, we currently lack a mechanistic explanation consistent with the key experimental findings on the stimulus-evoked interaction between excitation and inhibition in V1.
To resolve this issue, we have integrated our experimental findings into a unified large-scale model of cat V1 granular layer. Our model reproduces both the antiphase behavior of excitation versus inhibition for drifting gratings and the in-phase excitation, but broad inhibition, for flashed bars. Exploration of the connectivity parametric space of excitatory versus inhibitory tuning selectivity shows that a wide range of parametrizations induce neural dynamics comparable to those in cat simple cells. This suggests that elementary functional properties such as orientation tuning, contrast invariance, and the push–pull organization can arise robustly across a variety of connectivity schemes. However, the underlying conductance dynamics change significantly across this parameter space, demonstrating that very diverse motifs of excitatory/inhibitory interactions can underlie the same functional properties (Baudot et al., 2013; Fournier et al., 2011).
Materials and Methods
The basic architecture of the model used in this study is derived from a recent data-driven model of V1 (Antolík et al., 2019), based on a consensus of the experimental electrophysiological and anatomic literature in cat area 17. However, here, we have restricted simulations to the thalamocortical input recipient simple cells originating from layer 4 (the original model also comprised complex cells in supragranular layers 2/3). Below, we present a summarized description of the model used here and a full list of modifications from the original. We refer the reader to the parent model for further details (Antolík et al., 2019). The model has been implemented in the Mozaik framework (Antolík and Davison, 2013), whereas the NEST simulator (version 2.1.1; Gewaltig and Diesmann, 2007) was used as the back end for all simulations described in this article.
V1 model
The cortical model corresponds to a 1.2 × 1.2 mm patch of L4 of cat primary visual cortex centered at 5 degrees of visual field eccentricity (see Fig. 2). It contains 3600 cortical neurons and ∼4 million synapses which are driven by spikes representing LGN input. We simulated a population of excitatory neurons (corresponding to spiny stellate neurons in Layer 4) and one population of inhibitory neurons (representing all subtypes of inhibitory interneurons) in a 4:1 ratio (Gabbott and Somogyi, 1986; Beaulieu et al., 1992; Markram et al., 2004; Beaulieu and Colonnier, 1985).
All neurons were modeled as single-compartment integrate and fire units. Specifically, we used the adaptive exponential (AdExp) integrate-and-fire model, which is computationally efficient, offers a broad range of firing dynamics, and offers more realistic membrane potential time courses than simpler integrate and fire schemes (Brette and Gerstner, 2005; Naud et al., 2008; Destexhe, 2009). We set the membrane resistance (Rm) of all cortical neurons to 250 MΩ (Monier et al., 2008). We set the membrane time constant of excitatory neurons to 20 ms and of inhibitory neurons to 10 ms, close to values observed experimentally in cat V1 in vivo (Monier et al., 2008). The refractory period during which the membrane potential is held at −72 mV for all simulations was set to 2 ms and 0.5 ms for excitatory and inhibitory neurons, respectively. Overall, these neural parameter differences between excitatory and inhibitory neurons reflect the experimentally observed greater excitability and higher maximum sustained firing rates of inhibitory neurons (McCormick et al., 1985; Contreras and Palmer, 2003). The excitatory and inhibitory reversal potentials Eexc and Einh were set to 0 mV and −80 mV, respectively, in accordance with values observed experimentally (Monier et al., 2008). The threshold slope factor and the adaptation time constant of the AdExp model were set to 2.0 mV and 88 ms, respectively, for all neurons (Naud et al., 2008). We did not consider the subthreshold and spike-triggered adaptation in this study, and thus the corresponding parameters of the AdExp model were set to zero.
Retino-thalamic pathway model
We did not explicitly model the retinal circuitry but instead used the retinal spike pattern, predicted by the center-surround model of RFs, to feed LGN neurons (Fig. 2). The centers of both ON and OFF LGN neuron RFs were uniformly randomly distributed in the visual space, with a density of 100 neurons per square degree. Each LGN neuron had a spatiotemporal receptive field, with a difference of Gaussians spatial profile and a biphasic temporal profile defined by a difference of γ functions. The exact spatial and temporal parameters have been adopted from Allen and Freeman (2006).
To obtain the spiking output of a given LGN neuron, the visual stimulus was sampled into 7 ms frames and convolved with the spatiotemporal receptive field of the cell. In addition, saturation of the LGN responses with respect to local contrast and luminance was modeled (Papaioannou and White, 1972; Bonin et al., 2005). For simplicity, the local luminance (ll) was calculated as the mean of luminance values, and local contrast (lc) as the SD of the luminance values sampled within the RF of the given neuron. The response of the linear receptive field was separated into a direct current (DC; mean luminance) component (rl) and a contrast component (rc). The saturation of the two components is modeled with two Naka-Rushton functions as follows:
The resulting luminance and contrast temporal traces were then summed and injected into integrate-and-fire neurons as a current, inducing stimulus-dependent spiking responses. In addition to the stimulus-dependent drive, neurons were also injected with white noise current. The magnitude and variance of this current was such that neurons fire at an average rate of 10 Hz in the no-stimulus condition (Troyer et al., 1998). This artificially elicited spontaneous discharge was calibrated to reproduce experimentally observed spontaneous rates.
Thalamo-cortical model pathway
All cortical neurons in the model receive connections from the model LGN (Cruikshank et al. (2007). For each neuron, the spatial pattern of thalamocortical connectivity was determined by a Gabor distribution, inducing the elementary RF properties in cortical neurons (Jones and Palmer, 1987a; Tusa et al., 1978; Troyer et al., 1998; Fig. 2) as follows:
For individual neurons the orientation θ, phase ψ, size σ, spatial frequency λ, and aspect ratio γ of the Gabor distribution were selected as follows. To induce functional organization in the model, we used an existing model of stimulus-dependent orientation map development (Antolík and Bednar, 2011) that uses Hebbian learning to compute a stabilized link map that conditions an orientation map. Such a precomputed orientation map, corresponding to the 1.2 × 1.2 mm of a simulated cortical area, was overlaid onto the modeled cortical surface, thereby assigning each neuron an orientation preference of θ. For the sake of simplicity, the phase ψ of the Gabor distribution was assigned randomly, but see the recent findings on the topological organization of the ON/OFF pathways in V1 (Jin et al., 2011; Wang et al., 2015; Lee and Huang et al., 2016; Kremkow et al., 2016a). The remaining parameters were set to constant values, matching the average of measurements in cat V1 RFs located in the parafoveal area (Jones and Palmer, 1987b). Specifically, the size σ was set to 0.15 degrees of visual field, the spatial frequency λ to 0.8 cycles per degree, and the aspect ratio γ to 0.57 (Pei et al., 1994). For any given neuron, the thalamocortical connections were generated by overlaying its Gabor template over the model's sheets of LGN ON and OFF neurons, then drawing the connections randomly (with replacement) from the resulting probability distribution over the LGN neurons. Each cortical neuron received between 60 and 180 (drawn uniformly within these bounds) thalamocortical synapses (da Costa and Martin, 2011).
Corticocortical connectivity
We modeled 1000 corticocortical synaptic inputs per modeled excitatory cell. Inhibitory neurons received 20% fewer synapses than excitatory neurons to account for their smaller size, but otherwise synapses were formed proportionally to the two cell type densities. The synapses were drawn probabilistically with replacement. The geometry of the corticocortical connectivity was determined based on two main principles, the connection probability falls off with increasing cortical distance between neurons (Budd and Kisvárday, 2001; Buzás et al., 2006; Stepanyants et al., 2009; Fig. 2), and connections are formed preferentially between neurons with similar functional properties (Ko et al., 2011). The two principles were each expressed as a connection-probability density function, then multiplied and renormalized to obtain the final connection probability profiles from which the actual corticocortical synapses were drawn. The following sections describe how the two probability density profiles of connectivity were obtained.
Spatial extent of local intracortical connectivity
The exact parameters of the spatial extent of the model local connectivity were established based on a reanalysis of data from the cat published in Stepanyants (Stepanyants et al., 2008). The details of the analysis can be found in the publication of the parent model (Antolík et al., 2019), but briefly, the probability of potential connectivity data from Stepanyants et al. (2008) was averaged to obtain a 1D profile expressing probability of connectivity as a function of distance between two cortical neurons of a specified type (excitatory or inhibitory).
To obtain a parametric representation of the distance connectivity profiles we fitted them with a zero mean hyperbolic distribution as follows:
The resulting values of the parameters of the fitted hyperbolic distributions for all combinations of presynaptic and postsynaptic neuron types, which were used to generate the local connectivity distance-dependent profiles in the model, can be found in Table 1.
The parameters of hyperbolic profiles of potential connectivity derived from connectivity data in Stepanyants et al. (2008) used to generate the distance-dependent profile of local connectivity in the model
Functionally specific connectivity
Recent studies of local connectivity in mice have shown that this iso-iso preference might be a general connectome feature independent of the species-specific presence of orientation preference maps, because despite the salt-and-pepper V1 organization, local connections in rodents also showed a weak bias toward connecting neurons with similar receptive field properties (Cossell et al., 2015) or a similar orientation (Denman and Contreras, 2014). Finally, the dominant antiphase relationship between excitatory and inhibitory conductances in cat V1 simple cells (Anderson et al., 2000; Monier et al., 2008; Tan et al., 2011) has traditionally been interpreted as further evidence for functionally specific inputs (Anderson et al., 2000; Kremkow et al., 2016b). Overall, these results point to a moderate tendency of excitatory neurons toward connecting nearby neurons of similar receptive field properties, although this bias increases somewhat for more distant postsynaptic neurons. In this study we restrict our investigation to the local connectivity.
Because of the lack of clarity and specificity of experimental data on the functional connectivity in V1, we use previously hypothesized schemes of functional connectivity and adapted them to be compatible with the above experimental findings. Among cortical neurons, we assume push–pull connectivity (Troyer et al., 1998; Fig. 2) of varying strength of the functional biases. For each pair of cortical neurons, the correlation c between their RFs was calculated. The connectivity likelihood for a given pair of neurons is given by the following:
Synapses
Synaptic inputs were modeled as transient conductance changes, with exponential decay with time-constant τe = 1.1 ms for excitatory synapses and constant τi = 1.9 ms for inhibitory synapses. Relatively little data exists on the exact strength of synapses between different neural types and layers. The synaptic weights were thus selected to achieve an overall balance between excitation and inhibition that supports reasonable levels of both spontaneous and evoked activity while being compatible with the limited physiological findings. Specifically, we have set each excitatory synapses onto an inhibitory neuron to 2.4 nS, whereas all remaining synapses in the model were set to 1.4 nS. But note that the connection generation algorithm allows multiple synapses to form between pairs of neurons. We have also modeled synaptic depression for thalamocortical and excitatory corticocortical synapses (Abbott et al., 1997) using the model from Markram (Markram et al., 1998). We did not model short-term plasticity for inhibitory synapses as it is not well studied, but see the review in Kripkee and Froemke (2017). For all excitatory synapses we assume parameters corresponding to moderate depression (U = 0.75, τrec = 150, τpsc = 3.0, and τfac = 0; Markram et al., 1998).
Delays
The delays in the feedforward thalamocortical pathway were drawn from a uniform distribution bounded between 1.4 and 2.4 ms. For all intracortical connectivity, a distance-dependent delay with a propagation constant of 0.3 m/s (Bringuier et al., 1999) was used, corresponding to the slow propagation of action potentials along the intra-areal (lateral) unmyelinated processes. Furthermore, we have included a constant additive factor in all synaptic delays, specifically 1.4 ms to excitatory to excitatory synapses, 0.5 ms to excitatory to inhibitory synapses, 1.4 ms to inhibitory to excitatory synapses, and 1.4 ms to inhibitory to inhibitory synapses, in line with experimental observations (Ohana et al., 2012).
Visual stimuli
We used two types of visual stimuli, flashed stationary bars and drifting gratings. The bars were vertically oriented, matching the preferred orientation of the 30 excitatory cells, the responses of which are quantified in this manuscript, and either bright or dark contrast (100%) at varying locations across the visual field. Sinusoidally varying contrast gratings (spatial frequency, 0.8 cycles per degree; temporal frequency, 2 Hz; 100% contrast) were presented at various orientations, including optimal (0°) and orthogonal (90°) for the recorded cells.
Data analysis
Analysis of the model output was performed using customized code in MATLAB. Responses [Vm, excitatory conductance (gE), and inhibitory conductance (gI)] from multiple repeated trials were averaged for each neuron. We measured peak responses (Vm depolarization from rest or increase in synaptic conductance) as the average within a 50 ms window across all bar positions for each cell, allowing for some variability in the peak response times among bar positions. For visualization purposes, these responses were fit to a Gaussian distribution for each RF subregion. Visual space coordinates were rescaled and translated so that 0 corresponded to the peak of the Vm response to bright and 1 to the peak Vm response to dark, and the distance between the centers of adjacent Vm-derived RF subregions defined 1 receptive field unit.
In the experimental studies that form the basis of the present model (Monier et al., 2003; Taylor et al., 2018); excitatory and inhibitory synaptic conductances were estimated using current clamps. A full replay of such current-clamp experiments would be extremely resource demanding in our large-scale detailed model and in the context of the extensive parameter search that we performed. But rather than estimating indirectly the underlying simulated synaptic conductances from voltage traces recorded under a current clamp, the full network simulation allows us to record them directly. To compare our virtual recordings with those estimated experimentally, we simulated current-clamp responses to drifting sinusoidal gratings and estimated the underlying conductances from the voltage traces. On the whole, the estimated synaptic conductances (Fig. 3A) are found similar to the underlying simulated model synaptic conductances (Fig. 3B) but with somewhat smoother variations. The level of variability of the estimated conductance traces in the simulations matches that of experimental data in cat V1 (Fig. 3C, showing data from Baudot et al. 2013). The key measure of our model (the temporal correlation between the excitatory and inhibitory conductance traces) remains the same, whether based on virtual recordings or experimental current-clamp measurements (Fig. 4D). We thus conclude, that the use of simulated conductance traces in this computational study is fully justified.
To quantify the selectivity to stimulus orientation, we calculated the half-width at half-height (HWHH) as in (Alitto and Usrey, 2004), using a Gaussian fit to the orientation tuning data as follows:
Results
Simple cell characteristics and classification criteria
To evaluate the physiological fidelity of our model, we first established a set of criteria for which it has to match experimental data. Simple cells in V1 L4 have many well-characterized features in response to visual stimuli, and although many of these have been demonstrated in a previous version of this model (Antolík et al., 2019), here we focused on responses to drifting gratings and to flashed bars (Fig. 1A). When stimulated with drifting gratings, simple cells respond with maximal Vm depolarization and spike output to a preferred orientation, and this orientation selectivity is invariant to contrast (Fig. 1B, top). The Vm and firing rate are temporally modulated in phase with the grating (Fig. 1A,B, bottom). Optimally oriented bright bars flashed in ON subregions evoke depolarization and an increase in firing rate, whereas dark bars in the same ON subregion evoke hyperpolarization and an absence of spikes (Fig. 1C). This characteristic behavior in simple cell RFs, observed at the Vm and spike levels (Fig. 1C, bottom), are classically interpreted as the landmarks of push–pull connectivity. The measure of a negative (anticorrelated) spatial correlation of the Vm response to light versus dark bar stimuli is used here to characterize the simpleness of V1 receptive fields.
Schematics of V1 response characteristics to bars and gratings. A, Left, Example of a 2D RF of a simple cell. Middle, Types of stimuli used in our simulations, grating (top), bars (bottom). Right, Vm and conductance responses. B, Schematized spiking, Vm, and conductance response profiles to drifting gratings. Left, The FR and Vm responses to a grating at the preferred orientation of the cell are modulated in time. Right, Excitatory (gE, red) and inhibitory (gI, blue) conductances have been shown to be modulated in temporal antiphase in response to this stimulus. C, Spiking and Vm responses to flashed bright or dark bars at the preferred orientation of the cell. Top, The 1D RF calculated from the response of the cell to bright minus the response to dark reveals the three RF subunits (ON, OFF, ON), also shown in 2D (left inset, with overlaid bar positions in dotted lines). Middle, Spiking responses to bright (gray thick line) or dark (black thin line) bars at different positions across the RF of the cell in visual space. Bottom, Vm responses to flashed bars have a similar pattern and are spatially anticorrelated. D, Conductance responses over space to flashed bars, as shown by recent experimental results (Taylor et al., 2018). Although the gE responses over space are anticorrelated, the gI responses are broad and highly spatially correlated. E, Conductance responses over space to flashed bars, as predicted based on the anticorrelated conductance responses to gratings. In this scheme, gE and gI responses to bright and dark bars are spatially anticorrelated.
When it comes to synaptic conductances under these stimulation protocols, the literature has produced seemingly conflicting results. In response to drifting gratings at the optimal orientation, the current-clamp observation of spatial anticorrelation between the temporal patterns of dominant excitation and inhibition driven by opponent contrast have led some authors to infer that the push–pull organization observed for Vm and spikes extends into the conductance domain and implies that gE and gI are separated in space, or spatially anticorrelated (Fig. 1E). This spatial segregation predicts that stimulation of an RF subregion with signed matched flashed bar (bright bar on ON subregion or dark bar on OFF subregion) will evoke mostly gE, whereas the opposite contrast bar will evoke mostly gI. However, experimental evidence has suggested that when measured with flashed bars, the spatial footprint of gI is broad across the entire RF, with colocalized peaks for the two stimulus polarities (Taylor et al., 2018; Fig. 1F; Borg-Graham et al. 1998, their Fig. 2), which is in agreement with voltage-clamp measurements of spatial overlap of excitation and inhibition in simple cell RFs (Borg-Graham et al., 1998; Monier et al., 2008). This suggests that although gE responses to bright versus dark bars across the RF are anticorrelated over space, gI responses to bright versus dark bars have a positive correlation. We used these observations to guide our exploration of connectivity parameter space in the model, described next.
Model construction
We modified a model of the thalamocortical LGN-L4 V1 circuit (Antolík et al., 2019), the details of which are described above in Materials and Methods and schematized in Figure 2. Briefly, the model contains center-surround LGN cells that are connected to a network of excitatory and inhibitory neurons in L4, all of which have simple RFs. Feedforward connections from the LGN follow RF correlation rules (Hubel and Wiesel, 1962; Alonso et al., 2001; Sedigh-Sarvestani et al., 2017; Fig. 2D), whereas local connections between excitatory and inhibitory cells were determined based on a combination of distance between the cell pair (Fig. 2A) and the RF correlations of the cell pair (Fig. 2C,E). The connectivity shown in this first iteration of the model, which best reproduces the experimental data, can be described as a push–pull-like organization, that is, connections from excitatory cells to other neurons were more likely if the RFs were highly correlated, and much less likely if they were anticorrelated (Fig. 2C,E, red line), whereas connections from inhibitory cells were only mildly biased toward anticorrelated RFs (Fig. 2C,E, blue line). Later, we further explore a subset of this intracortical connectivity parameter space.
Schematic of LGN-V1 model architecture. A, Cortical lateral connectivity. Functional specificity of local connections is not shown, and connection ranges are not to scale. B, Extent of modeled visual field and example of RFs of one ON- and one OFF-center LGN relay neuron. The model is retinotopically organized, and the extent of the modeled visual field is larger than the corresponding visuotopic area of modeled cortex to prevent clipping of LGN RFs. C, Local connectivity cortex follows a biased push–pull organization; excitatory connections are biased toward correlated RFs, whereas inhibitory connections are mildly biased toward anticorrelated RFs (E). D, Afferent RFs of cortical neurons are formed by sampling synapses from a probability distribution defined by a Gabor function overlaid on the ON and OFF LGN sheets, where positive parts of the Gabor function are overlaid on ON and negative on OFF sheets. The ON regions of RFs are indicated by the white color, and OFF regions by black. E, Connection bias factor as a function of RF correlation between presynaptic and postsynaptic cells. Excitatory cells (red) are much more likely to connect to other cells with similar RFs, and inhibitory cells (blue) are slightly more likely to connect to other cells with anticorrelated RFs (Monier et al., 2003, their Fig. 6). Later, we explore a parameter space of excitatory and inhibitory connection specificity by varying these values of σEXC and σINH.
We characterized the behavior of this network in response to two sets of visual stimuli—flashed stationary bright and dark bars at preferred orientation of the neurons at different locations across the RF, and full-field drifting gratings at multiple orientations. The spiking output, Vm, gE, and gI, was recorded from 30 excitatory simple cells estimated to be in L4 in response to these stimuli.
Broad inhibition in response to flashed bars
We first tested the model's response to flashed bars of bright or dark contrast presented across the width of simple cell RFs, at their preferred orientation. Vm and spike responses to a single bright or dark bar presented within the ON subregion of an example cell in the network (Fig. 4A, left) show the expected push–pull relationship; a bright bar evoked depolarization and spikes, whereas a dark bar evoked hyperpolarization and spiking suppression. The underlying synaptic conductances also matched experimental data, with both stimuli evoking large magnitudes of gI increase but only the bright bar evoking a significant increase in gE (Fig. 4A, right).
Comparison of true and estimated mean synaptic excitatory and inhibitory conductance in model neurons. A–C, The mean over 10 trial of excitatory (red) and inhibitory (blue) conductance trace in response to 2 s presentation of sinusoidal grating of optimal orientation, drifting at 2 Hz. A, Conductance estimated through a simulated current-clamp method in three representative model simple cells. B, True conductance in the same three model neurons as in A, C, An example of excitatory and inhibitory conductance estimated through a current clamp in a simple cell recorded in cat V1, adopted from Baudot et al. (2013). D, Pearson correlation between estimated excitatory and inhibitory conductance for 30 model neurons (abscissa) plotted against the same measure derived from the true recorded conductance of the same neurons (ordinate).
Responses of model neurons to flashed bars. A, Vm, center, four example single trials overlaid (top rows) and average ±SD (bottom row) and conductance (right, average ±SD) responses of an example simple cell from the model to a bright bar flashed in its ON subregion (top row) and a dark bar flashed in the same position (middle row). The cell fired in response to the bright bar in the ON subregion, and this response was underlain by a combination of gE and gI (circles indicate maximum response). A dark bar in the ON subregion evoked hyperpolarization of the Vm, underlain by minimal gE and a strong gI response. B, Average Vm, gE, and gI responses to flashed bright and dark bars across space for the same example cell. Circles represent raw values at each location, and lines show a Gaussian best fit for illustration only. Spatial correlation values between bright and dark responses (right) were calculated from the data points. C, Histograms of spatial correlations from n = 30 cells from the simulation (black) and from n = 17 cells from a previously published dataset (pink; Taylor et al., 2018). Vm and gE responses to bright versus dark bars are negatively correlated, whereas gI responses are positively correlated. Arrowheads indicate medians of each distribution.
We plotted the peak evoked values of Vm, gE, and gI to each bar position (Fig. 4A, peak values indicated by circles) over space in Figure 4B. The Vm responses to bright versus dark bars over space were anticorrelated (Fig. 4B, top; Vm spatial correlation = −0.95), reflecting the RF of the cell and following the push–pull characterization of simple cells. The gE responses also had a negative correlation (Fig. 4B, middle; gE spatial correlation = −0.23). The gI responses, however, produced a broad footprint in response to both bright and dark stimuli, resulting in a positive correlation over space (Fig. 4B, bottom; gI spatial correlation = +0.93). The correlations were performed on the peak evoked values in response to bright and dark for each position, not on the Gaussian fits to the data, which are shown in the figure for clarity.
The example cell was a representative sample of the population of 30 excitatory cells analyzed in the simulation. The distributions of spatial correlation values for Vm, gE, and gI are plotted in Figure 4C for both the cells in this simulation (black) and for the experimental dataset (purple) from Taylor et al. (2018). The spatial correlation values for Vm in the simulation are significantly more negative than those from the data (Fig. 4C, top; simulation median Vm spatial correlation = −0.93, p = 1.7e-6, Wilcoxon signed rank test; experiment median Vm spatial correlation = −0.49, p = 1.4e-3, Wilcoxon signed rank test). The stronger anticorrelation in the simulation arises from strong hyperpolarizing responses to bars flashed in the RF subregion of opposite polarity, which reflects the more depolarized Vrest of the cells in the simulation relative to the experimental neurons. In any case, the anticorrelated Vm responses to bright versus dark bars show that the Vm of neurons behave in a push–pull manner typical of simple cells in L4.
The peak gE responses to bright and dark bars were spatially anticorrelated, whereas the peak gI responses were strongly spatially correlated (Fig. 4C; simulation median gE spatial correlation = −0.35, p = 4.7e-6, Wilcoxon signed rank test; simulation median gI spatial correlation = +0.81, p = 1.6e-6, Wilcoxon signed rank test). The simulation matches the data very well for the spatial correlation values of gE and gI (experiment median gE spatial correlation = −0.43, p = 0.017, Wilcoxon signed rank test; experiment median gI spatial correlation = +0.71, p = 2.9e-4, Wilcoxon signed rank test; calculated from data in Taylor et al., 2018).
In summary, despite only a weak bias toward antiphase connectivity between inhibitory and excitatory neurons, this network produced spatially restricted excitation and broad inhibition in response to flashed bars across the RF, consistent with experimental findings (Taylor et al., 2018).
Antiphase inhibition in response to drifting gratings
Having replicated the experimental findings of spatially broad inhibition evoked by flashed bars across the width of simple cell RFs, we asked whether this network would also produce antiphase (temporally anticorrelated) modulation of excitatory and inhibitory conductances in response to drifting sinusoidal gratings. A drifting grating at the preferred orientation produced modulated spike output and Vm (Fig. 5A), whereby spikes occur at the peaks of the modulated Vm trace when the RF subregions are aligned with the matching luminance regions of the drifting sinusoidal grating stimulus. The conductance responses to the preferred orientation grating were noisy but did appear to be modulated in temporal antiphase (Fig. 5B, top). By smoothing the conductance traces, we revealed a clear anticorrelated relationship between gE and gI (Fig. 5B, bottom), with a correlation value of −0.74. The population of cells analyzed in the simulation showed a strongly negative distribution of temporal correlations, with a median at −0.61 (Fig. 5C; n = 30 cells, p = 5.3e-4, Wilcoxon signed rank test).
Model neurons responses to drifting gratings. A, Vm and spiking responses over time from an example cell (same cell as in Figure 3) to drifting gratings at its preferred orientation (top, single trials; middle, average without spikes ±SD; bottom, spike rasters from 10 trials). B, The conductance responses from the same cell over time to the same grating. Top, Average gE (red) and gI (blue) from 10 trials. Bottom, Smoothed conductances (see above, Materials and Methods) clearly illustrate a negative temporal correlation (correlation coefficient = −0.74) between gE and gI. C, Histogram of temporal correlations from n = 30 cells from the simulation. As in the example cell, the majority of cells show temporal anticorrelation between gE and gI to the drifting grating. Arrowhead indicates the median correlation from this distribution.
In summary, a strong bias of excitatory connectivity and a weak bias of inhibitory connectivity among anticorrelated neurons in L4 produced a strong antiphase relationship between excitation and inhibition, consistent with experimental findings (Ferster, 1988; Anderson et al., 2000; Monier et al., 2003; Tan et al., 2011).
Exploring parameter space of intracortical connectivity
The simulation shown in Figures 4 and 5 had specific values determining the RF-correlation-based likelihood of connectivity between any pair of neurons in L4 (Fig. 2E; σE = 1, σI = 2). As shown above, this simulation reproduced the temporal and spatial conductance dynamics shown in previous experiments from different labs. We next assessed the importance of these intracortical connectivity rules in determining the behavior of the network by exploring the parameter space of σE and σI. We independently varied σE and σI from 0.6 to 4.0, with lower values representing less selective connections and vice versa (Fig. 6A). Throughout Figure 6, 2D plots illustrate inhibitory connectivity in the abscissa and excitatory connectivity in the ordinates. The sigmoidal scale functions of probability of connection as a function of RF correlation for the four extreme points in parameter space are illustrated on the right (Fig. 6A, a–d). We note that this parameter space explores the classical push–pull connectivity scheme, which at the extreme would be characterized by the lower left portion (Fig. 6A, c). We explored from this extreme to almost no bias, independently for excitatory and inhibitory connections.
Spiking, Vm, and conductance characterizations across a parameter space of intracortical connectivity schemes. A, Illustration of the parameter space explored in this figure. The sigmoidal function parameters σE (for excitatory connections) and σI (for inhibitory connections) were varied independently between 0.6 (highly selective) and 4 (virtually unselective) values. Throughout this figure, selectivity of excitatory connections grows to the left (of the x-axis) and that of inhibitory connections grows to the bottom (of the y-axis) of each matrix. Right, RF correlation-based connectivity curves for the four extreme combinations (corners of parameter space, a–d; Figure 2). B, Variation of spiking and Vm characteristics across explored parameter space. Left, The median HWHH (degrees) of the spiking responses to high-contrast drifting gratings at varying orientations. Center, The median change in HWHH between the tuning curves in response to gratings at high versus low contrast. Right, The median spatial correlation of Vm responses to bright versus dark bars (as in Fig. 4) at the preferred orientation of the cells. Aside from the lower/lower-left region of parameter space (corresponding with highly selective inhibitory connectivity), these metrics did not vary substantially across the connectivity parameter space. Black squares indicate parametrizations where orientation tuning could not be well fit with Gaussian curve. C, Variation of conductance characteristics across explored parameter space. Left, The median temporal correlation of gE and gI in response to high-contrast drifting gratings at the preferred orientation of the cells as in Figure 5C. Center, The median spatial correlation of gE responses to bright versus dark bars flashed in different locations across each cell's RF as in Figure 4. Right, The median spatial correlation of gI responses to bright versus dark bars flashed in different locations across each cell's RF as in Figure 4. D, Illustration of the region of parameter space that yields simulations with values of spiking, Vm, and conductance response characteristics that are consistent with experimental findings. This region of parameter space corresponds with relatively selective excitatory connections (0.6 < σE < 2.0) and much less selective inhibitory connections (σI > 1.2). Note that the simulation illustrated in Figures 4 and 5 is shown within this parameter space with an asterisk (*).
We first explored the role of connection selectivity on the spiking and Vm characteristics of model L4 neurons in response to drifting gratings and flashed bars (Fig. 6B). We quantified the median orientation tuning width measured as HWHH of the spiking response of the neurons to drifting gratings (Fig. 6B, left), the change in HWHH between low- and high-contrast gratings (Fig. 6B, center), and the spatial correlation of the Vm responses to bright and dark bars at the preferred orientation of the neurons (Fig. 6B, right).
Overall, we found only a small effect, with the exception of a section of parameter space where excitatory and especially inhibitory connections were highly selective. This was reflected in a less physiologically realistic behavior; for example, orientation tuning curves of most neurons could not be well fit by a Gaussian function (Fig. 6B, black squares). Immediately surrounding this area in parameter space, the neurons exhibited large contrast-dependent changes in orientation tuning (Fig. 6B, center, yellow squares), indicating that highly selective inhibitory connections preclude contrast-invariant orientation tuning. This portion of parameter space also featured neurons with less anticorrelated Vm in response to flashed bars (Fig. 6B, right). Closer inspection of population dynamics of models in this subdomain revealed unstable, often oscillatory, behavior explaining the breakdown of elementary functional features reported in Figure 6B. In summary, we conclude that with the exception of highly selective connection schemes, a relatively wide range of connectivity parameters are sufficient to reproduce experimentally determined spiking and Vm characteristics of simple cells in L4 of V1.
We next examined the differential effect of the selectivity of the excitatory and inhibitory connectivity profiles on the relationship between excitatory and inhibitory input conductances. In contrast to spiking and Vm responses (Fig. 6B), the temporal and spatial correlations of synaptic conductances showed strong variations (Fig. 6C). However, gE and gI responses to drifting gratings remained robustly anticorrelated even when the specificity of inhibitory connectivity was decreased, approaching zero or positive values for extremely broad tuning values (Figs. 6C, left, 7A,B). This striking observation occurred despite the decreasing magnitude of the phase-dependent modulation of the inhibitory synaptic input (Fig. 7A,B). Thus, even weakly biased inhibitory connectivity can give rise to experimentally observed levels of temporal push–pull between excitatory and inhibitory conductances in response to drifting gratings. We also observed a breakdown of the anticorrelation between the gE and gI when selectivity of the inhibitory connections was very high (Fig. 6C, bottom edge). But it is important to point out that in this region of the parameter space, the model dynamics were unstable, and thus the correlation measure of the push–pull was dominated by large stimulus-independent runaway events characterized by concomitant rise and fall of excitation and inhibition throughout the cortical network.
Push–pull measures for different model parametrizations. A, B, C, D, Both refer to parametrizations marked with a corresponding letter in Figure 5D. Each shows the average gE (top left) and gI (middle left) responses to flashed bars across space for an example cell; gE and gI responses to drifting grating stimulus (bottom left) for the same example cell; histograms of spatial correlations between gE reponses to dark and bright bars (top right), between gI responses to dark and bright bars (top middle) and between gE and gI responses to drifting sinusoidal grating stimulus. All plotting conventions match corresponding plots in Figures 4 and 5.
Increasing the selectivity of inhibitory connections thus decreased the spatial correlations between gI responses to bright versus dark bars (Fig. 6C, right) and at the same time decreased the temporal anticorrelation between excitatory and inhibitory conductances in response to gratings (Fig. 6C, left). Crucially, the temporal anticorrelation of gE and gI to gratings broke down only at very broad selectivity of inhibitory connectivity (σI greater than ∼3), whereas the inhibition evoked by bars became broad and overlapping already at moderate inhibitory connection selectivity values (σI greater than ∼1.2). Thus, one can find regions of the parameter space where the model behaves in line with experimental evidence for both drifting grating and flashed bars (Fig. 6D).
The following question remains: How can there exist model parametrizations where selectivity of inhibitory connections is sufficiently broad to induce overlapping inhibition across space in response to flashed bars but sufficiently specific to induce the temporally anticorrelated excitation and inhibition observed in response to drifting gratings? We hypothesize that the explanation lies in the difference of temporal dynamics of the input drive respectively fed by drifting grating versus flashed bar stimulus. A flashed bar is an abrupt step stimulus that evokes strong initial response (onset) followed by rapid decline to a lower sustained level (adaptation). This justifies why, as most experimenters do, Taylor et al. (2018) focused their analysis only on the early 50 ms of the response unaffected by the adaptation. The early response is, however, dominated by the stimulus onset dynamics, characterized by concomitant rise of excitation and inhibition (Borg-Graham et al., 1998, their Fig. 2). Such a short flashed stimulation protocol does not allow enough time for the cortical network to self-stabilize into a dynamic regime dictated by the biases in the functional circuitry unless the biases are very strong, as revealed in our parameter search (Fig. 6). This nonstationarity hypothesis is supported by the experimental evidence showing sharpening of orientation tuning of spike response over time from stimulus onset (Ringach et al., 1997; Schummers et al., 2007). This is in contrast to a steady-state response to drifting grating stimuli that engages the cortical network in a steady slowly changing fashion, which reduces adaption over period of seconds. This delayed and long integration time allows the network dynamics to settle in a dynamic regime dictated by connectivity biases, making even weak biases in the cortical connectivity robustly functionally affecting both excitatory and inhibitory input conductances. Incidentally, the onset dynamics with concomitant increase in excitation and inhibition are also present for the drifting grating stimuli, both in our model (Fig. 5B; Antolík et al., 2019, their Figs. 5, 10) or in cat (Monier et al., 2008, their Fig. 5). But because of the short duration of the onset dynamics, this has little influence on the measure of correlation between excitation and inhibition in response to the long stimulation periods covering several cycles of drifting gratings.
The parameter search also revealed that the selectivity of excitatory and inhibitory connections did not act independently, as is evidenced by a moderate slope in the correlation levels between gI responses to bright versus dark bars (Fig. 6C, right). This is illustrated in the two-model parametrizations shown in Figure 7, C and D, and marked by corresponding letters in Figure 6D. Both parametrizations correspond to models with the same inhibitory but different excitatory selectivity. Of the two models, the model with more selective excitation showed lower correlations between gI to bright versus dark bars (Fig. 7C), compared with the model with broader excitatory selectivity (Fig. 7D). This can be explained by the fact that the excitatory selectivity parameter dictates also the selectivity of connections from excitatory to inhibitory neurons, which, however, influences the strength of the phase-dependent modulation of the output of inhibitory neurons.
An important property of the model that emerged from the parameter search was the higher sensitivity to changes in inhibition than excitation. Figure 6B shows that the transition of the model dynamics to unstable happened mostly along the bottom edge of the plots, indicating that smaller change in selectivity of inhibitory than excitatory connectivity can push the model into the unstable state. Similarly, the temporal correlation between gE and gI (Fig. 6C, left panel) changes more rapidly along the inhibitory selectivity axis. Such greater sensitivity to changes in inhibition are likely because of the thalamocortical excitatory connections imposing strong bias on excitation in the cortex. On the other hand, there are no thalamocortical inhibitory connections. This means that the same relative change to the excitatory cortical circuitry has a less relative impact on the overall excitation in the cortex than an equal change to inhibitory cortical circuitry has on the overall inhibitory cortical dynamics. This hypothesis is consistent with the finding that the correlation (across RF space) between synaptic conductance evoked by ON and OFF bar stimuli exhibits greater magnitude of change when excitation selectivity varies (Fig. 6C, middle) than when inhibition selectivity varies (Fig. 6C, right).
Together, the parameter search showed that experimentally observed simple cell characteristics of spiking, Vm, and synaptic conductances in response to both drifting gratings and flashed bars arise when excitatory connections are relatively selective (presynaptic and postsynaptic neurons have correlated RFs) and when inhibitory connections are only mildly selective (presynaptic and postsynaptic neurons tend to have slightly anticorrelated RFs). This area of the parameter space is illustrated in Figure 6D.
Response to flashed bright or dark squares
The Hirsch et al. (1998) experiments have provided grounding support to the push–pull concept by showing that the suppression (pull) evoked by stimulus flashed in the RF subfield of mismatched polarity is because of am increase of inhibition, rather than just a withdrawal of excitation. As expected, our model reproduces the characteristic push–pull response to flashed squares (Hirsch et al., 1998), just as it does to flashed bar stimulus (Taylor et al., 2018). Our simulations show Vm depolarizations for flashed point-like stimuli, the contrast (light/dark) of which matches the polarity of the RF subregion (Fig. 8A,B), whereas opponent hyperpolarization is induced by a stimulus of the opposite contrast in the same RF location (Fig. 8C,D).
Response of model neurons to flashed square stimuli under different current injection levels. The four columns of the figure show data from four representative excitatory cortical model neurons. A, C, Receptive field maps (light for ON subfield, dark for OFF subfield) with the position of the contrast-matched square stimulus (A) and contrast-mismatched square (C) stimulus overlaid on the stimulation grid. Grid spacing was 0.3°. B, D, Top, Traces indicate conductance changes (gE, red; gI, blue). Bottom, Traces indicate membrane potential (black line). B, response to contrast-matched stimulus. D, response to contrast-mismatched stimulus. E–G, Voltage response to contrast-mismatched stimulus during + 0.1 nA (E), −0.1 nA (F), and −0.2 nA (G) current injection. H, Total synaptic conductance change from baseline calculated following the same procedure as in Hirsch et al. (1998). Stimulus presentation is indicated by black bar under each voltage trace.
Whereas evoked hyperpolarization grows with the intensity of intracellularly injected positive current (Fig. 8E), the amplitude diminishes (Fig. 8F) and eventually reverts into depolarization (Fig. 8G) for negative currents when reaching the apparent reversal potential of dominant inhibition, which is in line with Hirsch et al. (1998). When measuring the membrane resistance during the visually evoked hyperpolarization in response to sign-mismatched flashed square stimulus, we observe roughly a doubling of the membrane conductance (Fig. 8H), also in agreement with Hirsch et al. (1998). Finally, although excitatory conductance input is present only in response to flashed point stimulus with contrast matching the polarity of the RF subfield, inhibition is present regardless of the RF subfield polarity throughout the whole extent of the RF, just as predicted by the flashed bar experiments in Taylor et al. (2018). This shows that our model predictions account for the experimental observations in both stimulation paradigms, reported independently by Hirsch et al. (1998) and Taylor et al. (2018).
Relationship between the spike and input orientation tuning of conductances
To elucidate the generation of orientation tuning in the model we investigated the relationship between the orientation tuning of the spiking response of the neuron and of its excitatory and inhibitory input conductances. First, let us point out that in our model, the preset push–pull bias imposed in intracortical connectivity (Fig. 2A) implies that both excitatory and inhibitory synaptic inputs originate predominantly from neurons of similar orientation preference. Thus, we would expect orientation tuning of the spike response of a neuron and both excitatory and inhibitory inputs to coincide. Although it is generally admitted that neurons in V1 receive inputs preferentially from excitatory and inhibitory neurons that are cotuned in the orientation domain (Michalski et al., 1983; Smith and Kohn, 2008; Ko et al., 2011; Denman and Contreras, 2014; Lee and Bonin et al., 2016), more direct voltage-clamp measurements reveal a diversity of the relative tuning preferences of the neuron's excitatory and inhibitory inputs (Monier et al., 2003, 2008), much larger than that assumed in the spike-based literature or inferred by current-clamp experiments (Anderson et al., 2000).
To validate our model against these experimental findings, we applied the analysis from Monier et al. (2003) to the responses to drifting sinusoidal gratings. Figure 9 shows the results from our model (Fig. 9A–C) as well as replotted analogous data from Monier et al. (2003; Fig. 9D,E). The response magnitude of the Vm of a simple cell to a drifting grating can be measured by the DC component of the subthreshold signals (F0 component) and by the magnitude of the modulation at the frequency of the drift (F1 component).
Tuning properties of the excitatory and inhibitory input conductances in excitatory simple cells. A, B, Scatter plot of preferred orientations of gE (abscissa) and gI (ordinates) input conductances in response to the drifting sinusoidal grating stimuli. The preferred orientation of the synaptic input is plotted relative to the orientation preference of the spiking response for each model cell. The iso condition corresponds to the bottom left corner of each scatter plot, where gE, gI and spike output preferences are coaligned. The F0 component was used to compute the orientation preference (A), whereas we used the F1 component of the input conductance (B). C, The HWHH of the excitatory versus inhibitory input conductance tuning curves. D, As in A and B but taken from the in vivo data from simple V1 cells in cat (Monier et al., 2003; Fig. 5, bottom left). Because the test stimulus was a moving bar, the separation into F1 and F0 components was not possible. We have excluded the plot cells shown in Monier et al. (2003), which had low selectivity (Monier et al., 2003, their Fig. 5) because for such cells, the preferred direction preference cannot be determined reliably. E, As in C but data from cat (Monier et al., 2003; Fig. 5, top left). To make the Monier et al. (2003) data directly comparable to the present study, we have collapsed in D and E opposite directions corresponding to the same orientation within the (iso, cross) range. All model data correspond to the condition marked with * in Figure 6D.
As Figure 9A shows, the tuning of the F0 component of both excitatory and inhibitory inputs was similar to that of the spiking response. Paradoxically, a significant number of neurons exhibited different tuning of input conductance and spiking response. For such neurons, excitation and inhibition showed similar tuning. The same behavior has been identified in a subgroup of neurons reported in Monier et al. (2003; Fig. 9D, points along the identity axis). On the other hand, if we examine the F1 component of the conductance signals, we see a different pattern. There are few neurons for which the orientation preference of the excitatory conductance differs from that of the spiking response, but in contrast, the preference of orientation tuning of the inhibitory inputs varies widely and is independent of the excitatory preference. Again, such behavior has been noticed in a second subgroup of neurons reported in Monier et al. (2003; Fig. 9D, points along y-axis). Finally, we find that the mean tuning width of the excitatory and inhibitory inputs in model neurons is remarkably similar (Fig. 9C) and in very good agreement with Monier et al. (2003; Fig. 9E), although we do observe lower overall variability of the tuning width in comparison to the experimental data. Such lower variability could reflect the greater regularity that our model construction process induces in comparison to the biological V1 circuitry.
We should note that there are, however, some differences between the present computational study and Monier et al. (2003). First, Monier et al. used moving bars rather than drifting sinusoidal gratings, which precludes the separate analysis of the F1 and F0 components in their study. Second, their study was not restricted to layer 4 simple cells, and we cannot exclude the existence of systematic tuning differences between the different layers and functional cell types. Regardless, the presented analysis on subthreshold orientation tuning genesis shows a close fit between our simulations and these experimental data, which strengthens the validity of the working hypothesis of our model.
Discussion
We have presented a biologically constrained model of the thalamocortical visual circuit that replicates two sets of seemingly contradictory experimental findings. The model employs correlation-based local connectivity among simple V1 cells in which excitatory projections are strongly biased toward cells with spatially correlated RFs, whereas inhibitory projections are weakly biased toward cells with spatially anticorrelated RFs. This connectivity regime generates broad inhibition evoked by flashed bars across the RF (Taylor et al., 2018) as well as antiphase modulation of excitation and inhibition in response to drifting gratings (Anderson et al., 2000; Monier et al., 2008; Tan et al., 2011), resolving the seeming contradiction in the experimental data.
Generally, these findings support the view that V1 simple cell connectivity in cat V1 does not adhere to a strict push–pull schema. It is often assumed that antiphase modulation of gE and gI in response to drifting gratings requires an underlying connectivity rule in which the inhibitory inputs to a simple cell must arise from an inhibitory cell cotuned in orientation but antiphase in the spatial dimension (Troyer et al., 1998; Anderson et al., 2000). This idea was recently challenged by our experimental findings showing that inhibition could be triggered from flashed bar stimuli across simple cell RF, independent of phase, suggesting a largely indiscriminate scheme for inhibitory connectivity not biased by RF correlations (Monier et al., 2008; Taylor et al., 2018). The simulations presented here support a compromise between the two strict connectivity schemes; the parameter search (Fig. 6) showed that neither a strictly push–pull nor a fully indiscriminate inhibitory connectivity scheme is compatible with the experimental findings.
The tendency of excitatory neurons synapsing onto neurons with correlated functional properties has been well established (Michalski et al., 1983; Monier et al., 2003; Smith and Kohn, 2008; Ko et al., 2011; Denman and Contreras, 2014; Gerard-Mercier et al., 2016; Lee and Bonin et al., 2016), supporting the correlated part of the push–pull-like scheme. However, evidence on connectivity rules of inhibitory neurons, especially on the specific implication of the push–pull-like scheme that the inhibitory neurons should preferentially target neurons with anticorrelated functional properties was lacking. Importantly, a series of recent studies from the Fitzpatrick lab gave support to this notion, showing that inhibitory neurons in layer 2/3 receive selective inputs from potential presynaptic pools (Wilson et al., 2017) and preferentially target other inhibitory neurons with an opposite direction preference (Wilson et al., 2018). Here, we focus on the relationship between the phase of RFs of pairs of neurons and the likelihood of a connection between them, whereas Wilson et al. (2018) investigated the direction preference, not phase. But if we extend a spatial RF kernel into the temporal domain, two neurons preferring the same orientation but an opposite direction preference will be perfectly anticorrelated. The Wilson et al. (2018) study thus offers important indirect evidence supporting the existence of the anticorrelated inhibition part of the push–pull-like scheme as hypothesized here (see also example of directionally anticorrelated inhibition in Monier et al. (2003, their Fig. 3, cell 11).
In the present study, we assume the putative correlation-based canonical circuit only for simple cells receiving direct thalamic input. A wider repertoire of gE/gI configurations was observed electrophysiologically in cells with simple-like spike-based behavior but an unidentified layer of origin (Monier et al., 2003, 2008; Baudot et al., 2013). Interestingly, our results show that despite the regularities imposed by our correlation-based circuit scheme, the model also generates gE/gI configurations that are not straightforwardly predicted from the correlation-based connectivity scheme (Fig. 9), and, crucially, specific relationships between the relative gE/gI tunings emerge that are in a remarkable accordance with the experimental data (Fig. 9).
A question remains on the source of this wider repertory of configurations if they cannot be linked directly to the underlying connectivity scheme in our model. We hypothesize that the probabilistic nature of connection generation used in the model is likely to induce small random biases in individual neurons, which are further amplified by the recurrent cortical dynamics. Another likely contribution, already proposed by Monier et al. (2003), is local neighborhood heterogeneity in the orientation maps, whereby a neuron, depending on location of the neuron in the map (close to pinwheel singularity vs in the center of an iso-orientation domain), could recruit patterns of excitatory and inhibitory inputs of varying orientation preference. Because of local averaging, the direct impact of these factors is likely small, but it is further amplified by the recurrent interactions of the cortical network. In this study, the size of the model, spanning approximately a single hyper-column, and the need to record only neurons in the central area of the model to avoid border effects (see above, Materials and Methods), did not allow us to collect sufficiently diverse data to test this hypothesis.
The majority of the model parameters were set based on existing experimental data. But some, where relevant data or quantification are missing, could not be fully constrained. Among these are the functional biases of the local connections between excitatory and inhibitory cells, which, as extensively discussed above, constitute key elements for the results of the model. However, for these two parameters, we have explored extensively the parameter space, determining the range where observed data could be reproduced (Fig. 6). It is unfortunately computationally infeasible to perform a systematic exploration of all remaining, not fully constrained, model parameters. But it is important to point out that all model parameters are further indirectly constrained by the battery of functional tests that the model has been already validated against in the parent manuscript (Antolík et al., 2019). However, our empirical experience with the model identified several parameters, to which the model dynamics are more sensitive, such as the relative contribution of the feedforward versus intracortical excitatory drive to cells, or generally parameters that influence the overall balance of excitation versus inhibition.
As the discovery of poorly orientation-tuned inhibitory neurons in the granular layer of V1 (three reconstructed smooth cells in Hirsch et al. 2003), untuned inhibition has been welcome in computational models, the use of which enabled simulation of contrast invariant orientation tuning and low-pass temporal frequency tuning in simple cells (Lauritzen and Miller, 2003). But later studies (Cardin et al., 2007; Nowak et al., 2008) that investigated this issue in 55 fast-spiking cells, failed to identify inhibitory neurons that lack orientation tuning. Cardin et al. (2007) found that when measured as HWHH, fast-spiking neurons were on average only moderately more broadly tuned than excitatory ones. They did not observe any untuned fast-spiking cells, and only one with an orientation tuning >100 degrees. Consistently, Nowak et al. (2008) report only two (of 19) inhibitory neurons with relative unselective response amplitude above 50%, and none above 70%. They also show that the orientation tuning width of fast-spiking cells forms a continuum across the population. The present model does produce some poorly tuned inhibitory cells (Fig. 6D, in model condition marked with *, 16.6% of inhibitory cells exhibit HWHH of orientation tuning of spike response >45° and only 2.5% of excitatory cells do so), which contribute to the moderately selective inhibition that has been identified here as a prerequisite for explaining the experimental data. Note, however, when pooling all layer recordings, voltage-clamp conductance decomposition and nonlinear spike triggered analysis of intracellular responses reveal a similar orientation tuning width in the presumed excitatory and inhibitory sources of simple cells (Monier et al., 2003, their Fig. 6B; Fournier et al., 2014, their Fig. 8B), which is also reproduced in our model (Fig. 9C). Thus, although our model is consistent with the statistical evidence of some poorly tuned inhibitory neuron in V1, a separate class of untuned inhibitory cells was not added. The important finding is that even without this built-in tuning constraint, the model did generate contrast-invariant orientation tuning as experimentally reported and quantified in Figure 6B.
Other attempts have been made to constrain neural modeling of V1 with a wealth of biological data. At the structural level, large-scale integrative approaches have been gaining traction recently. In mice, several studies explored extensive morphologically detailed circuit reconstructions of the sensory cortex (Markram et al., 2015; Hawrylycz et al., 2016; Arkhipov et al., 2018). In higher mammals, a number of less anatomically detailed studies investigate visually evoked subthreshold network dynamics in biologically realistic spiking networks (Tao et al., 2006; Wielaard and Sajda, 2006; Rangan et al., 2009; Chariker et al., 2016). But rarely did models try to reconcile structure and function assessed from the microscopic to the mesoscopic scales in a way to confront predictions with the stimulus dependency of visual responses reported in V1 in vivo.
Although we strive to gradually introduce the full breath of known anatomic and functional properties of the early visual system into a single model (Antolík et al., 2019), necessarily this has to be a gradual, incremental process. It is thus important to identify the key limitations of the present model iteration, pertinent to the subject of this study, to guide future work. The presented model focuses on simple cells receiving direct thalamic input. Although it is well established that the majority of simple cells are confined to layer 4 and layer 6 in cat V1 (Gilbert, 1977; Hirsch et al., 1998; Martinez et al., 2005), it should be acknowledged here that the key experimental studies constraining this computational work (Anderson et al., 2000; Monier et al., 2003, 2008; Tan et al., 2011; Taylor et al., 2018) did not provide the exact layer origin of all the recorded simple cells. However, as the exact laminar location is of lesser importance here than the simple push–pull behavior of the spiking RF, it remains fair, in our view, to accept that the present model best applies to the organization of excitation and inhibition in simple cells in a virtual thalamocortical recipient layer (which may collapse layer 4 and layer 6 in cat V1). This model has still other pending structural unknowns. The topological organization of the ON/OFF thalamic pathway (Jin et al., 2011; Lee and Huang et al., 2016; Wang et al., 2015; Kremkow et al., 2016a) should be incorporated into the model in future because it can have interesting implications for the organization of excitation and inhibition in relation to the RF of the neuron and explain, for instance, the OFF bias of inhibition shown by Taylor et al. (2018). Finally, a more refined biophysical dissection should be engineered, concerning the subcomponents of the excitatory and inhibitory synaptic composite drive because of the various receptor types (NMDA, AMPA, GABAA, GABAB), which may be recruited differentially by the different visual input regimes.
In summary, this study shows that the explicit multiscale integration of the structural and biophysical properties of the underlying neural substrate can be a powerful approach for accounting for the diversity of electrophysiological measurements of V1 physiology. Our model demonstrates that in complex dynamical systems such as the visual cortex, using restricted results from a limited number of experimental conditions (and visual contexts) is insufficient to make far-reaching, generalized conclusions on the synaptic nature of the underlying connectivity. Only systematic, comprehensive computational models of the primary visual cortex that incorporate a broad range of established constraints from different experimental designs can lead to reconciliation of the often seemingly contradictory diversity of experimental findings and, in turn, to an accurate characterization of the system under study.
Footnotes
This work was supported by National Institutes of Health Grant R01EY027205 (D.C.), National Science Foundation Grant Graduate Research Fellowship Program DGE-1321851 (M.M.T.), National Eye Institute Vision Training Grant T32 EY007035-38 (M.M.T., D.C.), Centre National de la Recherche Scientifique (A.D., Y.F.), European Community Human Brain Project (H2020-720270 and H2020-785907; A.D.), Agence National de la Recherche (ANR-Horizontal-V1 and ANR-17-CE37-0006; Y.F.), Charles University Improvement of Internationalization in the Field of Research and Development project (CZ.02.2.69/0.0/0.0/17_050/0008466; J.A.), Charles University Primus Research Program 20/MED/006, France-United States Fellowship from the Chateaubriand Foundation (M.M.T), Agence National de la Recherche 353 PARADOX, and the ICODE excellence network. M.M.T was hosted by the European Institute of Theoretical Neuroscience in Paris. We thank Madineh Sedigh-Sarvestani and Larry Palmer for suggestions and discussion.
The authors declare no competing financial interests.
- Correspondence should be addressed to Jan Antolik at antolik{at}ksvi.mff.cuni.cz
