Refractory density model of cortical direction selectivity: Lagged-nonlagged, transient-sustained, and On-Off thalamic neuron-based mechanisms and intracortical amplification

A biophysically detailed description of the mechanisms of the primary vision is still being developed. We have incorporated a simplified, filter-based description of retino-thalamic visual signal processing into the detailed, conductance-based refractory density description of the neuronal population activity of the primary visual cortex. We compared four mechanisms of the direction selectivity (DS), three of them being based on asymmetrical projections of different types of thalamic neurons to the cortex, distinguishing between (i) lagged and nonlagged, (ii) transient and sustained, and (iii) On and Off neurons. The fourth mechanism implies a lack of subcortical bias and is an epiphenomenon of intracortical interactions between orientation columns. The simulations of the cortical response to moving gratings have verified that first three mechanisms provide DS to an extent compared with experimental data and that the biophysical model realistically reproduces characteristics of the visual cortex activity, such as membrane potential, firing rate, and synaptic conductances. The proposed model reveals the difference between the mechanisms of both the intact and the silenced cortex, favoring the second mechanism. In the fourth case, DS is weaker but significant; it completely vanishes in the silenced cortex.DS in the On-Off mechanism derives from the nonlinear interactions within the orientation map. Results of simulations can help to identify a prevailing mechanism of DS in V1. This is a step towards a comprehensive biophysical modeling of the primary visual system in the frameworks of the population rate coding concept.


Introduction
Mathematical models of primary vision, aimed to be comprehensively comparable with experimental data are still under development. In the present paper we implement mechanisms of direction selectivity (DS) of primary visual cortex neurons into a model that has already shown its ability to reflect experimental data on cortical tissue excitability and orientation selectivity.
Visual information processing begins in the retina and subsequently passes through the visual thalamic structures (primarily the lateral geniculate nucleus, or LGN) to the primary visual cortex [1,2,3]. The neurons of the primary visual cortex are selective to various characteristics of a stimulus, such as orientation, the direction of motion, color, etc. [4,5,6]. Neurons preferring a particular orientation (orientational neurons) or direction (directional neurons) of a stimulus are unevenly distributed within the primary visual cortex and are grouped into the so-called functional columns or modules [7,8,9]. The columnar structure constitutes socalled hypercolumns, each including a full set of the orientational or directional columns. The functional maps with orientation hypercolumns have been found in tree shrews, ferrets, cats, and primates. Our modeling is focused to the visual cortex with the orientation hypercolumns. However, we also consider DS mechanisms revealed in rodents, for which experiments with intracellular recordings and optogenetic stimulation have been more developed.
Some important questions regarding the mechanisms underlying DS remain unanswered. One of the questions is whether DS is inherited from the directionally selective subcortical neurons (mainly geniculate cells) or is reconstructed in the cortex de novo. DS of geniculate cells is well-documented for lagomorphs and rodents [10], but data on carnivores and primates point to some directional bias in geniculate cells rather than a prominent DS [11,12]. Recently, Lien and Scanziani [13] used an optogenetic approach to show that even in mice, the cortical DS emerges de novo at the convergence of thalamic synapses on the same cortical neuron, whereas at the geniculate level, DS is negligible.
Most DS models include a time delay between the responses of the spatially shifted neighboring geniculate neurons terminating at the same cortical cell [14], similar to the mechanisms revealed in the retina [15]. These geniculate neurons reveal either an increase or a decrease in their activity in response to the light stimulus (so-called On and Off cells) or differences in the temporal characteristics of their responses. In the second case, LGN neurons are categorized as (i) transient and sustained (T-S), which generate fast-and slow-decaying responses, respectively [13,16,17], or (ii) lagged and nonlagged (L-N), depending on the latency to discharge after spot onset [18,19,20]. Such neurons exist in cats [17,19], monkeys [21,22], and mice [10] and have been proposed to provide the spatiotemporal offset for cortical DS [19,23].
The physiological mechanism of the delay formation was studied in several experimental works [24,25], and it has been demonstrated that the delay between responses to the center and surrounding area of a receptive field is determined by GABA-A (gamma-aminobutyric acid) receptors within the LGN. A complex model of DS has also been proposed based on specific convergent projections of the signals from LGN cells with and without a delay and intracortical interactions [26]. In this model, the temporal difference between the lagged and nonlagged cell responses and the structure of LGN inputs provide DS in V1. An alternative mechanism of DS based on T-S thalamic neurons converging on the same cortical neuron has been recently proposed by Lien and Scanziani [13]. The authors suggested that the T-S neurons generate a fast-and a slow-decaying excitatory postsynaptic current (EPSC), respectively, which combine into a compound EPSC. Thalamic neurons prefer distinct spatial phases, so the decay of the compound EPSC changes with the phase. Time-staggered summation results in large or small F1 modulations of the compound EPSC depending on the direction of the simulated motion. Thus, the DS of the cortical neuron depends on the phase shift between the T-S thalamic neurons.
Alternatively, recently reported experimental data obtained through the use of optogenetics [27] and multielectrode electrophysiological recordings [28] suggest that DS in V1 is determined by a displacement of On and Off subzones of the receptive fields of V1 neurons. However, a mechanism that is based on the On-Off subzone displacement and provides DS for symmetric white-and-black stimuli like gratings is unclear. Summarizing, some alternative mechanisms of DS have been discussed in the literature; however, the dominant mechanism is still unknown.
Besides of feedforward mechanisms, DS is highly dependent on the intracortical interactions, which follows even from the fact of weak correlation between the preferred directions evaluated from the thalamic input and the neuronal spike response [13]. The intracortical interactions provide some important effects of the primary cortex functioning, such as the prolonged responses to brief stimuli [29] or apparent motion [30]. As shown in experiments and models, the balanced intra-layer cortical interactions between inhibitory and excitatory populations play a major role in shaping the dynamic stimulus representations in the early visual cortex [31,32]. The question of interplay between the feedforward DS mechanisms and the intracortical interactions is an open question brought up in the present paper.
In our previous works, we extended our biophysically detailed, conductance-based population model of cortical orientation tuning [33,34]. This approach is alternative to comprehensive network approaches that simulate single neurons explicitly [35,36]. A single population is defined here as a large number of similar uncoupled neurons receiving similar inputs. In these notations, the cortex is considered a layered continuum of coupled populations. Due to the combination of (i) the population level description, which is optimal in the framework of the population coding of information, (ii) a quite accurate population approach, namely, a conductance-based refractory density (CBRD) method [37,38], and (iii) a hypercolumnar structure description of V1, the model was able to reproduce a vast series of experimental observations obtained in V1 in vitro and in vivo. In the present paper, we discuss the reproduction of the abovementioned three DS mechanisms using this model by relying on (i) lagged and nonlagged, (ii) transient and sustained, or (iii) On-and Off-neurons. Moreover, we explain whether DS may occur as an epiphenomenon due to the orientational hypercolumns in the absence of any thalamo-cortical bias. Using these four models, a potential role of intracortical interactions in DS has been identified by mimicking the silencing of the cortex.

Results
We present simulations performed with the model (see Methods), which combines the detailed, conductance-based description of neuronal population activity in area V1 with the simple, filter-based description of the retino-thalamic processing of visual signals. The generalized model retains the advantages of the former model [34], relating simulations to different experimental observations obtained in slices and in vivo in the visual cortex, and allowing the consideration of experiments with moving stimuli. These stimuli, such as moving gratings ( Fig  1A, left) are used to calculate the fields of activity of the LGN neurons having center-surround receptive fields (RFs) (Fig 1A, middle), and the activity of the V1 area having neurons with more sophisticated RFs (Fig 1A, right). The elements are specific to each of DS mechanism.
According to the model, cortical neurons receive inputs from a particular "footprint" in the LGN; this footprint is a domain in the LGN that sends axons to a given V1 neuron [40]. This V1 neuron receives thalamic input only from that footprint domain of the LGN. The orientation and direction selectivity of cortical neurons are determined based on the properties of the LGN neurons and the structure of the footprint. The elongation of the footprint determines the cortical neurons' orientation preference. DS depends on the asymmetry of the V1 neuron projections received from lagged (L) versus nonlagged (N), transient (T) versus sustained (S), or On-versus Off-neurons, for the L-N, T-S, or On-Off mechanisms of DS, respectively. Data on the asymmetry of those connections is lacking in the experimental literature, so for simplicity our simulation assumes that convergent L and N (or T and S) inputs to a cortical cell are spatially segregated, similar to the assumption about adjacent elongated T-and S-cell subfields made in the alternative modeling work [36]. In our simulations, the footprint of a V1 neuron splits into two halves along the axis of elongation, each sending signals from either lagged or nonlagged thalamic cells (Fig 1A, right). A similar footprint is taken for the T-S mechanism (Fig 1B, middle). The elongation of the footprints provides the orientation selectivity while their width determines the optimal wavelength of the stimulating gratings. For the On-Off mechanism, the footprint is shaped like the shifted, oval On-and round Off-subdomains ( Fig  1B, right). The elongation of the On-subdomain provides an orientation selectivity even for the stimuli that evoke only the On-pathway. The magnitude of the shift effects the optimal wavelength of the stimulating gratings. This footprint structure is a prototype of that of typical RFs [28,41].
The sum of the inputs from the thalamic cells depends on the direction of the stimulus in the L-N and T-S mechanisms ( Fig 1B). Together with a threshold-linear-like or sigmoid-like nonlinearity, this sum determines the cortical cell activity, which is significantly different for the preferred and non-preferred directions.
In the On-Off mechanism, the sum is identical for both directions, if the stimulus is as shown in Fig 1B (i.e. the displacement of the white and black subzones in corresponance to the displacement of the On-and Off-subfields of the footprint). This is also the case for the moving gratings of any wavelength. For gratings moving with a unity velocity, this issue can be clarified as follows. We can approximate the thalamic input through the On-subfield of the footprint as |sin(t)| and the Off-input as |sin(t + π)|, and we can introduce a ratio of Offto-On contributions k and a subfield displacement Δφ. The total input when the stimulus is moving in one direction is |sin(t)| + k |sin(t + π + Δφ)|, and for the opposite direction, it is |sin(−t)| + k |sin(−t + π + Δφ)|. The inputs have equal amplitudes for any k and Δφ, which do not provide any preference for any direction. Thus, in the considered cases, DS cannot stem from just the subcortical projections, but it may derive only from nonlinear interactions within the orientation map.
with monkeys [29], a short (50ms) spot of light evokes an activity that lasts more than 200 ms. A similar effect is observed in rodents and is shown to be provided by recurrent cortical excitatory circuits [42]. Our simulation of the response to a round spot produces results very similar to those of the experiment (Fig 2B). This effect is nontrivial because all the processes that provide positive feedback in the dynamical system under consideration are much slower than the response. The response is most sensitive to the recurrent excitatory connection strengths, � g AMPA;E and � g NMDA;E , which reveals the role of intracortical recurrent excitation in maintaining the cortical activity. With all excitatory intracortical conductances reduced by 25% the response duration dramatically decreases. This effect will be a subject of future analysis as it is out of the scope of the present paper. However, here we explore the model with sharply tuned intracortical connections to study the effects of intracortical interactions on DS pre-determined by different thalamo-cortical footprints.

Three mechanisms provide DS
We have examined the mechanism of DS by simulating responses to horizontal gratings of a spatial wavelength 1.2˚moving up and down with the temporal frequency of 2 Hz (Fig 3A). Simulated cortical domain contained 6 orientation hypercolumns located, for simplicity, on a rectangular grid (Fig 3B). In response to the stimulus, the bright spots in the patterns of activity averaged over a large time period correspond with high rates of firing of the V1 excitatory neurons ( Fig 3C). The spots appear in the columns that prefer an orientation similar to that of the stimulus. The patterns are not symmetrical with respect to the central vertical axis, which is due to DS, i.e., different direction preferences for neurons of the neighboring columns with the same orientation preferences belonging to different orientation hypercolumns. The peaks of E-cell activity are located in different hypercolumns depending on the direction of the grating movement. These simulated patterns are similar to those registered with optical imaging, such as those obtained in the cat visual cortex [6] (see their Fig 4A and 4B). The spots have a similar size and smooth shape and shift after a change in the stimulus direction. The comparison of three DS mechanisms in Fig 3 shows that (i) DS is strong in all three cases (Fig 3C, bottom) and characterized with sharp tuning diagrams (Fig 3D), high values of the direction selectivity index (DSIs) for a representative site ( Fig 3D) and considerable values of the DSIs averaged across the entire simulated cortical domain for the firing rate (0.29, 0.46 and 0.28 for the L-N, T-S, and On-Off models, respectively); (ii) DS maps are similar for L-N and T-S mechanisms; (iii) the DS map for the On-Off mechanism is different, though the number of spots is comparable; (iv) in all the cases, the direction response is not well-aligned with the orientation map shown in Fig 3B; (v) the T-S mechanism gives the most stable response on each cycle of the gratings; and (vi) the preferred direction (downward) in the L-N model corresponds with the direction from the lagged to nonlagged subfields of the footprint. The preferred direction in the T-S model corresponds to the direction from the sustained towards transient subfields.

Pinwheel structure and intracortical interactions lead to DS without any subcortical bias
A question arises regarding to what extent DS may be caused by intracortical interactions unbiased with the structure of the footprints. To answer this question, we have considered a model with a footprint of the same elongation but without any discrimination of LGN cells within the footprint. Thus, the thalamo-cortical projections were set to be unbiased towards the stimulus direction. In this model, cortical neurons do not receive any precortical information about the stimulus orientation. At the same time, the pinwheel structure of the thalamic input provides an intracortical bias because a neuron may receive unequal signals from the neighboring neurons with similar and opposite preferences. The model shows that the response patterns are symmetrical ( Fig 4A) according to the orientation map (Fig 3B), as expected; however, DS is still present, as seen from the tuning diagrams and DSIs given for two representative sites in Fig 4C. Though, DS is weaker (averaged over the entire domain DSI for the firing rate is 0.21), transient (Fig 4A, bottom), and the neurons preferring upward or downward directions are located in zones different from those for the three DS models discussed above. These facts along with the asymmetry of the patterns in Fig 3C lead to the conclusion that the direction maps obtained with the explicit DS mechanisms are determined by the precortical bias rather than the intracortical interactions. At the same time, the strength of DS may strongly depend on the intracortical interactions.
To reveal a role of intracortical connections (Fig 2A), the cortex was artificially inhibited. For that purpose we supply the inhibitory neurons with extra depolarizing current 100pA that mimicks the effect of photoactivation of cortical inhibitory interneurons expressing channelrhodopsin-2 [13]. It resulted in increased inhibition, which almost silenced the principle neurons.
In the case without any specific DS mechanism, the artificial inhibition of the cortex results in a significant reduction in the response (Fig 4B). The response is no longer selective to direction (Fig 4B, bottom); however, it is still weakly selective to orientation. Locations of the firing neurons (gray spots in Fig 4B, top) correspond to the horizontal orientation preference regions on the orientation map in Fig 3B. The dramatic decrease in the response in comparison with that for the intact cortex ( Fig 3A) reveals a major role of the intracortical interactions in strengthening the tuning to feature detection.

PLOS COMPUTATIONAL BIOLOGY
oscillations in the case of the preferred direction compared with the non-preferred one was significantly larger for the two initial cycles of the gratings, whereas it was only moderately increased for the later cycles. An even weaker DS was observed in the case without any directional mechanism (Fig 5D, intact cortex). Strong directionality can be observed in two cases-L-N and T-S mechanisms (Fig 5A and 5B, intact cortex)-in which the magnitude of the oscillations is significantly larger for the preferred direction. In all the cases, the responses to the orthogonal, horizontal gratings are smaller than for vertical gratings moving in both preferred and non-preferred directions, revealing strong orientation selectivity (see the tuning diagrams in Figs 3D and 4C).
Following Lien and Scanziani [13], we simulated responses of a V1 neuron to gratings moving in preferred and non-preferred directions under the conditions of an artificially inhibited cortex for each of the three DS mechanisms (Fig 5). The amplitudes of the firing rates are much smaller with the inhibited cortex compared to the intact cortex. The effect of cortical inhibition is different in the three models. A weak direction selectivity could still be observed for the L-N and T-S mechanisms (Fig 5A and 5B, cortex inhibited; Fig 5E), but no difference in the response to the preferred and non-preferred stimuli was evident for the On-Off mechanism ( Fig 5C, cortex inhibited; Fig 5E), such as in the case without a DS mechanism (Fig 5D, cortex inhibited). For the On-Off model, a weak temporally modulated response was observed even when the gratings moved in the orthogonal direction, which was not the case for the other models. The direction tuning diagrams for the inhibited cortex ( Fig 5E) demonstrate weaker tuning than those for the intact cortex (Figs 3D and 4C). The latter ones are more sharp because of intracortical shunting inhibition [43].
Vertical DSI patterns. In the case of the intact cortex, the distribution in the cortical space of the DSI is calculated for the vertical direction, which reveals complex, asymmetrical patterns for the L-N, T-S, and On-Off models ( Mean DSIs. In the case of the intact cortex, mean DSIs for the firing rate are high for the L-N, T-S, and On-Off models (0.29, 0.46, and 0.28, respectively) and are comparable with those observed in animals [44,45]. The inhibition of the cortex results in an expected crucial reduction of the mean firing rate (from 11-13Hz to about 1Hz). DSIs increase in the L-N and T-S models (0.37 and 0.60, respectively), which shows that intracortical tuning has a "contaminating" effect of as a price for amplification. DS vanishes in the On-Off model, as expected and explained at the beginning of the Results section. In the case of the intact cortex, the DSI values based on the firing rates are greater than those for the voltage (Fig 6A, compare the numbers in the bottom two lines), which reveals the threshold nonlinearity as in experiments [13,44]. The effect is opposite for the inhibited cortex where most voltage modulations are subcortically evoked, subthreshold, and thus do not modulate the firing rate.
Direction maps. Complex patterns are also observed in the preferred direction maps for the intact cortex ( Fig 6B). The patterns for the inhibited cortex are more regular, being structured by the thalamic input patterns. (The patterns for the inhibited cortex in the On-Off model and the model without any specific DS mechanism are insignificant because of negligible DSIs.) That the preferred directions of the spike response differ from those of the thalamic input is consistent with experimental data [13]. This difference is revealed to be due to cortex inhibition, which shows the contribution of intracortical interactions.
The model without any specific DS mechanism has a symmetrical vertical DSI pattern in case with an intact cortex ( Fig 6A, left panel in the right pair), in contrast to the other models. The

PLOS COMPUTATIONAL BIOLOGY
DS is considerable (mean DSI 0.21). The direction map ( Fig 6B, left panel in the right pair) is determined mostly by the orientation-tuned thalamo-cortical projections and intracortical interactions. Because of the model's construction, this map is not affected by any specific tuning of input or intracortical connections to the directions of stimuli. That is whythe similarity between this map and those of other models with an intact cortex ( Fig 6B, left panels of each of the left three pairs) reflects a contaminating contribution of the cortex on the direction tuning of neurons, different from the direction tuning caused by the thalamic input. The direction map of the T-S model is the least affected by intracortical interactions.
Variation of parameters. We varied some parameters of the stimulating gratings in the T-S model ( Fig 6C) and compared the simulated results with the control (Fig 6A and 6B, third panels). The comparison show that DS is stronger for faster and finer stimuli (DSIs for the firing rate were 0.29, 0.46, and 0.70 for 1, 2, and 4Hz, respectively; 0.64, 0.46, and 0.35 for 0.9, 1.2, and 2.4 degrees, respectively). This is in contrast to the response amplitudes, which are greater for bigger and slower visual patterns (maximum firing rates were 12, 11, and 8.4 for 1, 2, and 4Hz; 5.4, 11, and 21 for 0.9, 1.2, and 2.4 degrees, respectively). The blockage of NMDA receptors is known to reduce visual responses without significantly changing the degree of DS [46][47][48]. In our simulation, halving the NMDA conductance ( Fig 6C, right panels) consistently increased the mean DSI (from 0.46 to 0.55) with a reduced firing rate (from 11 to 9Hz) and produced avertical DSI pattern much closer to the one for the inhibited cortex (compare with Fig 6A) and the DS map that is an intermediary between the maps for the intact and inhibited cortex.  Fig 7A), T-S (Fig 7B), On-Off ( Fig 7C) and with no specific DS mechanism ( Fig 7D). These simulated signals are similar to experimental data, such as those from [43,44] and [49] (their Fig 5). As expected, the shape of the membrane potential that is traced in response to the preferred stimulus is close to sinusoidal, with maxima and minima corresponding to the oscillations of the gratings and without extra peculiarities. The mean peak amplitude of voltage modulation was about 10mV in the L-N, T-S and On-Off models, which is within the range of the experimental values (mean depolarization 9.5mV and hyperpolarization 3.7mV in [43], voltage modulation about 22 mV for the representative cells in [44] (their Fig 2A), 10mV in [50] and 8mV in [49]). In response to the non-preferred stimulus, this measure is about 2mV in our models, 12mV in [44], 5mV in [50], and 2mV in [49]. Our peak value for the firing rate of excitatory neurons is about 20Hz, which does not exceed the experimental values (80Hz in [44], 40Hz in [43], and 50Hz in [13]). The firing rate bumps are more narrow than those for voltage, similar to those of experiments [43,44]. The input and synaptic conductances are within the range of experimental values. The input conductance of excitatory neurons at the resting state in the model is about 18nS, which is comparable to the mean value of 16nS in the experiment [49]. The total synaptic conductance reaches 30nS in our models. After subtracting the total synaptic conductance in the state of the gray screen stimulation (9nS), the synaptic conductance modulation (21nS) is comparable to the mean experimental value (16nS) from [49].

Synaptic activity underlying DS
As shown in Fig 7A-7C, in all three DS models, the oscillations of the thalamic input to V1 in the case of non-preferred direction are small. The thalamic input oscillations in the case of the preferred direction are similar in all the models. The strongest tuning of the input to the cortex was observed in simulations with the T-S model. The total set of input signals for any neuron consists of NMDA, AMPA, and GABA synaptic conductances that include contributions of both thalamocortical and intracortical connections. NMDA conductance is the largest component in the responses of all the models; however, note that this component is voltagedependent due to magnesium blockage, and thus to exclude an effect of voltage oscillations, the NMDA-signal is plotted with the voltage-dependency factor fixed to 1, as if for zero-magnesium conditions. Therefore, its values are not to be directly compared with the AMPA and GABA conductances. In any case, all three synaptic components react to the change of stimulus direction. Oscillations of glutamatergic components are the most distinct features in the comparison between the responses to preferred and non-preferred stimuli. It is the main source for the voltage and firing rate modulations that produce DS. Comparing all models, the oscillations of the NMDA-component in the case of the preferred stimulus are the largest and most stable for the T-S mechanism. ) and inhibitory GABA (blue) synaptic conductances at the intracortical recurrent connections; the mean membrane potentials of the E-(red) and I-(blue) populations; the E-(red) and I-(blue) population firing rates. The NMDA conductance here does not consider the factor of the voltage-dependent magnesium block. The AMPA conductance was multiplied by 10 to be compared with the GABA and NMDA components. https://doi.org/10.1371/journal.pcbi.1008333.g007

PLOS COMPUTATIONAL BIOLOGY
In contrast to NMDA and AMPA, the GABA-component was rather constant over time (Fig 7). Minor GABA-and larger NMDA-and AMPA-oscillations were in-phase. Both mean voltage and firing rate oscillate during stimulus drift in all the models. In contrast to the T-S model, the oscillations in the L-N and On-Off models were larger for the first and second cycles of the gratings and then stabilized (see Fig 3C and compare Fig 7B to Fig 7A and 7C). A strong two-or three-fold decrement in both mean voltage and firing rate was observed up to the third wave, but thereafter, the oscillations were stable. The oscillations of firing rate are stable on larger time scale (Fig 3C). Therefore, the simulations indicate that the T-S model as the most reliable for DS.
In the case with no specific DS mechanism (Fig 7D), for the particular site of the cortex (S1), shown in Fig 4A, synaptic conductance modulation is similar for both directions of the stimulus according to the tuning diagram ( Fig 4C). The magnitude of the conductance modulations is comparable to those for the L-N and On-Off models, pointing again to the dominant role of intracortical activity evoked not by directional but by orientation-tuned thalamic input. The CBRD-model allows for the reconstruction of a representative neuron's behavior with input variables, such as synaptic conductance, known from the population activity. We simulated a representative neuron of the E-population in the T-S model to compare it with experimental data from [13] (see their Fig 1c and 1e), where the T-S mechanism has been highlighted. The representative E-neuron generates spikes only during the presentation of the preferred direction (Fig 8, top right) whereas, in the case of the non-preferred direction, only subthreshold depolarization was observed (Fig 8, top left). These voltage traces recorded in  1c and 1e). Top to bottom: the membrane potential of a representative excitatory neuron (green); the population firing rate (black); the current as though it was recorded in the voltage-clamp at the holding voltage of -70mV in the cases of an intact (Control, grey) and an inhibited cortex (Thalamic EPSC, black). The signals are from the representative neuron in the cortical site S1 marked by the white circle in Fig 3C. https://doi.org/10.1371/journal.pcbi.1008333.g008 response to moving gratings are consistent with those observed in experiments in vivo when the shape and amplitude of voltage oscillations are compared [13] (see their Fig 1 and present  Fig 8), [49,51]. In the null direction, the voltage and current modulation are significantly larger and spikes are much more numerous. As explained in [13], the temporal modulation of the excitatory postsynaptic current (EPSC) provided by the thalamic input is different for the preferred and non-preferred stimuli, whereas the integrated charge (Q) is maintained (compare the fulfilled areas in Figure 1 from [13] and Fig 8, bottom right and left). Whereas the charges in our intact cortex simulations were different for the preferred, orthogonal, and opposite directions (470, 216, and 259pC, respectively), the charge was the same for the thalamic EPSC evaluated with cortex inhibition (205, 196, and 197pC). The maintenance of the total input charge in the experiment and in our simulations verifies that the transient and sustained neurons are independently active and are not biased toward the stimulus direction. Indeed, the stimulus direction does not affect the footprint of a V1 neuron (Fig 3C top middle) and changes only the timing of the signals from the T-S neurons. This is why the averaged in time thalamic input is independent of the stimulus direction. In contrast, the temporal modulation of the current is determined by the correlation between the firing rates of those T-and S-neurons that contribute to the footprint; it therefore depends upon the stimulus direction.
Thalamic contribution was estimated in [41] as the ratio of the charge due to the thalamic EPSC to the charge due to the EPSC in the intact cortex, Q control /Q inhibited , which was found to be 0.36 for stimulation with gratings of preferred directions. In our simulations, the values of this ratio for responses of the representative neurons located in S1 in the L-N, T-S, and On-Off models were 0.46, 0.43, and 0.47, respectively (Fig 5A-5C), which are comparable with the experimental values. The mean values of the ratio Q control /Q inhibited for all neurons preferring the same stimulus direction in the L-N, T-S, and On-Off models were 0.42 ± 0.05, 0.41 ± 0.07, and 0.50 ± 0.10, respectively. Note that the ratio values for S1 match the mean values.

Identification of a prevailing DS mechanism
The most important question related to DS is likely how a prevailing DS mechanism can be identified. Though responses to drifting gratings have distinct characteristics in different DS models, they are not necessarily robust criteria. For instance, as seen from Fig 7C, the tuning of the thalamic input in V1 is very weak for the On-Off model; however, a moderate directionality was still observed (Figs 3C, 5C and 7C). This is mostly due to the periodicity of the gratings used as a stimulus. As noted at the beginning of the Results section, as weak DS is expected in the On-Off model for stimulation with gratings or for stimuli that excite only the On or Off subzone of the footprint. In contrast, a non-periodic stimulus that excites both subzones is expected to reveal DS in the model. To determine the difference between the On-Off and other DS mechanisms, we propose the application of an experimental protocol with a non-periodic stimulus containing only a moving edge between white and black domains of the screen (i.e. an initially black or white screen changing to white or black, respectively, with a moving edge between the domains of different light intensities [Fig 9]). However, to avoid contamination due to intracortical interactions, the experiment also requires inhibition of the cortex. The simulations showed that an inhibition in the form of the opto-driven excitation of interneurons, as applied in [13], is not sufficient to prevent contamination (Fig 9, gray lines). The principal neurons of the cortex must be silent. For this purpose, a hyperpolarization of neurons must be provided. During experiments, it can be done, presumably by means of halorhodopsin. In simulations, we modelled such hyperpolarization by an additive negative   (Fig 9A and 9C), showing an absence of DS only in the On-Off model. It is in contrast to the response to a bar ( Fig  9B and 9D), which is directionally selective even for the On-Off model with hyperpolarized principle neurons, i.e. with for the thalamic input in this model. This qualitative difference between the responses to an edge and a bar is observed only for the On-Off model. Hence, a similar experimental protocol might help establish the contribution of the On-Off mechanism to DS.

Discussion
In the present modeling study, a previously elaborated model of the neuronal population activity in the neocortex [34], which describes an orientational selectivity, has been adapted for DS as well.

Modeling approach advantages
Responses to visual stimuli, such as gratings, are determined by network activity, to which neuronal populations rather than single neurons substantially contribute. At the same time, this activity is often studied by means of intracellular electrode registrations. Thus, an ideal approach for mathematical modeling should describe the network activity and the activity of recorded representative single neurons, whereas a direct description of all neurons of the network might be considered redundant. This detailed but population-type approach of modeling suits these idealized requirements and contrasts with the other comprehensive models, such as the most prominent ones [35,52] in which the V1 network is directly simulated with a large number of individual neurons. As an advantage, our approach provides an even wider set of experiments captured by the model if those described in our previous paper on orientation tuning and a comparison between in vitro and in vivo experiments are included [34].
Contrary to the detection of orientation based on mainly spatial peculiarities of the receptive fields, the detection of stimulus direction also requires temporal peculiarities of the responses. Thus, a temporal shift between the inputs should be taken into account. A strong DS of the V1 neurons has been observed in the three cases of the integrated thalamic input received from neurons with: (i) lagged and non-lagged, (ii) transient and sustained, and (iii) On and Off receptive fields.

Differences between the three models
Quite similar levels of DS in the cases of L-N and T-S thalamic inputs (see Figs 3C, 5A and 5B) were observed, and a weaker but comparable DS in the case of On-Off inputs was observed. Both the mean voltage and the firing rate oscillate during the stimulus motion. These oscillations were the most stable in the T-S model (Fig 7B). In the On-Off (Fig 7A) and L-N (Fig 7C) models, DS was more evident for the first and second cycles of the gratings, whereas a strong three-fold decrement in both the mean voltage and the firing rate was observed for the third and following cycles. These facts justify favoring the T-S model; however, the transient DS provided by the other models can be also considered functionally sufficient for the secondary visual processing that receives information from the primary visual cortex fragmentarily, separated on the time intervals of eye fixations between saccades. It is interesting that non-periodic stimuli (single pair of black and white stripes, not gratings) result in the most prominent difference between the models (Fig 9). Another stimuli, such as moving dots [53], could also potentially be helpful to distinguish between the models; however, our simulations were limited by computational resources and thus considered only a small area of the cortex, so widespread stimuli could not be processed. Thus, this issue is a matter of future investigations.
The weakest directionality was obtained in the On-Off model. The On-Off division of the visual cells is the first within the visual processing pathway [20,54] and is believed to underlie a DS at the retinal level in lagomorphs and rodents [15,55]. This is not the case for animals with more developed vision, such as carnivores and primates, where the predominant role in DS belongs to higher structures, meaning the visual thalamus and the cortex. Accordingly, the On-Off thalamic mechanisms in these animals (like in our mathematical model) could play a secondary role in the cortical DS. This aspect aligns with the finding that the blockade of an On pathway has no discernible effect on orientation and direction specificities in the cortex of rabbits [56].
For the T-S mechanisms, the question about an origin of DS is related to a separation of the visual system into two main processing streams, so-called X and Y channels, provided by corresponding retinal ganglion cells and thalamic cells. One specific characteristic of the Y and X cells is the duration of their activity after a stimulus onset: fast for the "transient" cells (Y cells) and prolonged for the "sustained" cells (X cells) [16,17,57]. For many years, it has been supposed that only Y cells are responsible for motion processing [58], but D. Mastronarde, A. Humphrey, and A. Saul [19,20] suggested that this is not the case. Instead, the temporal shift between the responses of lagged and non-lagged counterparts irrespective of the X and Y types underlies DS in the visual cortex. The prominent DS in both the L-N and the T-S models is in favor of this suggestion.

Sources of tuning
Tuning to stimulus direction is provided by the cumulative action of the following factors: (i) the projections of LGN to V1, described in Section 2.3 (Fig 1), (ii) the threshold rectification in V1 neurons, and (iii) the recurrent intracortical interactions. Indeed, the LGN input to V1 excitatory neurons is tuned. This tuning is more evident in the T-S model (compare left and right traces in Fig 7A) and is reflected in EPSC recorded in the inhibited cortex in L-N and T-S models (Fig 5A and 5B); it is less evident in the On-Off model (compare left and right traces in Fig 7B) and in EPSC recorded in the inhibited cortex ( Fig 5C). Trivially, it is absent in the model without any DS mechanism (Fig 5D). The second factor results in stronger tuning of the mean membrane voltage of E-neurons than of I-neurons (Fig 7), which is due to the effect of recurrent synaptic inputs. The inhibitory interneuronal activity reflected by the GABAergic conductance on the excitatory neurons contribute to the tuning of voltage with depolarization relative to the spike threshold. The threshold rectification produces the firing rate and thus sharpens voltage tuning. Finally, the firing rate determines the probability for a single neuron of a population to generate spikes, as shown in Fig 8. As a result, the preferred direction is determined by the cumulative action of those factors, not simply by the preferred orientation of the thalamic input, which is consistent with experimental findings [13].

Synaptic conductances underlying cortical tuning
Our simulated conductances are not to be directly compared with experiment-based estimations. In our simulations, the GABA-component during gratings was rather constant over time; minor GABA-and larger NMDA-and AMPA-oscillations were in-phase (Fig 7). On the contrary, known experiment-based estimates of synaptic conductances [44,49,51,59] report that the excitatory and inhibitory conductances oscillate in antiphase. It should be noted that our simulated variables should not be directly compared with those estimates. Rather than three separate (AMPA, NMDA, and GABA) components, the mentioned experimental studies reported estimates of two conductances: excitatory and inhibitory. This method, due to not taking into account the voltage-dependence of NMDA receptor magnesium block, overestimates oscillations of the inhibitory conductance and gives in-phase modulations of the excitatory and inhibitory signals. This limitation of the estimation procedure has been recently revealed [60]. Our current simulations confirm that it is not an antiphase interaction of inhibition and excitation but rather oscillations of the excitatory synaptic drive that underly the temporal modulations of the firing rate response, which contradicts to some earlier models based on the mentioned estimations [61,62] but consistent with more recent work [35] and with the study emphasizing the role of NMDA receptors in cortical temporal tuning [63]. Similar data on the DS of the spike responses and membrane depolarization responses and excitatory and inhibitory interactions underlying DS were obtained using in vivo whole-cell voltage-clamp recording; it was shown that (1) both spike responses and membrane depolarization are sensitive to the stimulus direction, and (2) the bias of the inhibitory responses was much weaker than that of excitation or even not observed [64].

Cortical inhibition and directionality
In L-N and T-S models, the artificial inhibition of the primary visual cortex decreased the directionality level but did not totally disrupt it; this is mostly related to the T-S model. This might mean that cortical interneuronal communication enhances directionality. It is wellknown that cortical neurons with a preference for a particular direction are grouped into directional columns [6,65]. Many studies documented patchy patterns of cortical horizontal interconnections that have approximately the same center-to-center distance, such as between functional columns [66][67][68][69]. In accordance with the wiring economy thesis [70], this connectional pattern, such as cortical columns per se, can be a result of an attempt to minimize the wiring cost for the intracortical connectivity. In this case, a reduction in directionality due to cortical inhibition is expected. Moreover, the major role in strengthening directional tuning belongs to intracortical interactions (Fig 4). This could mean that DS may be caused by the inherent intracortical interactions unbiased by thalamic footprints.
Our modeling helps to estimate the role of pure intracortical interactions, which in the presence of a pinwheel structure of the cortex could lead to a DS that is independent of any thalamocortical input bias to stimulus direction. In fact, we have observed DS in this model, which is much weaker than that in the L-N, T-S, and On-Off models. Naturally, DS disappears after silencing the cortex. The effects of this mechanism cannot be completely excluded for intact cortex functioning. Moreover, these intracortical interactions within the orientation map are critical in allowing DS to emerge in the On-Off mechanism.

Development of the model
Because the detailed simulation presented here is consistent with the experimental data, in particular with those reported in A. Lien and M. Scanziani's paper [13], our CBRD-model is recommended for a prognostic modelling of hypothetic experimental conditions, such as the activation/blockage of a particular thalamic input, changing cortical connections, taking into account particular types of GABAergic neurons, and alterations of neuromodulators and cellular messengers. In a future work, it is also worthwhile to compare model predictions with in vivo data on wild-type and knock-out animals. For example, bursting properties of the cortical pyramidal neurons are determined by the hyperpolarization-activated cationic current [71]. The hyperpolarization-activated cyclic nucleotide-gated channels (HCN) underlie this current; within the CNS, these currents were first described in the retina [72]. In HCN knock-out animals, a special compensatory upregulation of the tonic GABA-A currents was revealed [73,74].
In the retina, it was observed that directionality is dependent on the HCN channels [75]. These conditions can be simulated with the proposed mathematical model.
The synaptic connections in our model were chosen without any optimization tuning. A task-aimed parameter tuning is expected to result in better visual processing performance, as in the model for motion processing in the MST area that uses an unsupervised optimization technique based on multi-cause and principle component analysis [76] or the model that uses Hebbian-like learning to tune the V1 model [35].
The model does not distinguish between different layers of the cortex, which is a reduction of our former model, in which layers 4 and 2/3 were explicitly present [34]. This reduction was made for simplicity. At the same time, we suppose that the transitions from layer to layer provide sharper tuning without a strengthening of recurrent excitatory connections [77]. The effects of multiple layers can be considered in future research.

Model limitations
The major limitation of this model is the structuring of the visual cortex by functional columns. The population approach is efficient in this case although its usage is not as reasonable in case of salt-and-pepper distribution of orientationally and directionally selective neurons. In the case of the cortex structured by pinwheels, the direction and orientation are implicitly coded by x-y coordinates and do not contribute to the set of independent variables. In contrast, in the salt-and-pepper case, the functional variables would increase the dimensionality of the mathematical problem.
Another model limitation is related to characteristics of thalamic cells as sources of input to cortical cells. The main sources of the geniculate input to the area V1 are layers A and A1 (in carnivores) or magno-and parvocellular layers (in primates). As is known, layers A and A1 are different, and several different features in addition to the ocular dominance of their cells were observed. For example, electrophysiological studies report that On-cells are predominant at the top of the A-layers and Off-cells at the bottom, whereas both types are balanced at the centers of the A-layers. The steepest gradients and maximum differences in the proportions of On-and Off-cells have been obtained in layer A [78]. For the transient and sustained neurons associated with Y and X cells [16,79,80], in layer A, Y cells are concentrated near the interlaminar borders and X cells in the centers of the layers, but these patterns are less obvious for layer A1 [78,81]. Our own data related to the postnatal development of a cat's LGN provides additional evidence of the different dynamics of the maturation of layers A and A1 [82]. These data together allow for assuming there are unequal contributions of the geniculate layers A and A1 to cortical DS. This aspect has not been taken into account in our model because only functional roles and not spatial placements of the thalamic neurons are important in the considered DS mechanisms. Also, the model does not take into account a heterogeneity of either exciting or inhibitory cortical neuron populations. It is well-known that cortical inhibitory networks consist of multiple types of inhibitory neurons [83]; their role in cortical directionality is unclear. To clarify this aspect, the mathematical model can be further developed to describe the behavior of different neuronal types separately.

Main predictions to be verified in experiments
The model predicts that the thalamic input biased toward DS and the intracortical interactions not specifically tuned to DS can together provide an efficient motion detection by cortical neurons in the presence of an orientation map structure. Hypothetically, this statement could be proven or disapproven in experimental settings in which the synaptic plasticity can be switched off specifically for intracortical connections, thus leaving the thalamic projections intact during the development of an animal, knowing that DS appears at a certain age [45]. A lack of strong DS in such treated adult animals would disprove the main model prediction. On the other hand, it is straightforward for the model to account for intracortical connections tuned to the direction preferences of the connected neurons, which is expected to increase DS. Patchy connections primarily between neurons with similar orientation preferences and, presumably, direction preferences are neglected in the model. These connections could drastically change DS, so this aspect should be verified in experiments.
The separation of the thalamic footprints' subfields ( Fig 1B) is another assumption of the model, similar to the assumptions made in alternative modeling studies [35,36]. Disapproval of these assumptions would ruin the mentioned models.
Another model prediction mentioned above concerns the synaptic conductances underlying cortical tuning. This prediction can be tested with experimental reconsideration of synaptic estimations in animals with pinwheels in vivo during stimulations with gratings, as proposed in [60] or with alternative methods. Evidence of considerable oscillations of GABAmediated conductance in antiphase with glutamatergic conductances would demonstrate that the present intracortical tuning of the model is incorrect.
The model emphasizes the importance of strong recurrent excitatory connections and associates this tuning with a prolonged response to short activation. If this association is incorrect, the model should be re-tuned in favor of DS specificity in the connections.
In conclusion, the proposed biophysical model realistically reproduces the activity of the visual cortex, reveals the difference between the DS mechanisms, and helps to elaborate an experimental protocol that would reveal a prevailing mechanism of DS. These findings along with the mentioned previous modeling of orientation selectivity are a step towards a comprehensive biophysical modeling of the primary visual system in the frameworks of the population rate coding concept.

Model of stimulation
Model stimuli are black-and-white sine-wave gratings of various orientations (from 0˚to 180˚) and the spatial frequency 0.83 cycles per degree moving in two opposite directions orthogonal to the orientation at a 2 Hz rate. Stimuli parameters were chosen close to those found optimal in experiments [84].
The gratings are moving relative to the center of the screen (x c , y c ) with some orientation θ 0 , wavelength λ, frequency ν, amplitude of luminance modulation S A and background luminance B: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À x c Þ 2 þ ðy À y c Þ 2 q cos½arctgððx À x c Þ=ðy À y c ÞÞ þ y 0 �=l À n tÞÞ ; ð1Þ A basic set of the parameter values was as follows: θ 0 = 0, λ = 1˚, ν = 2Hz, B = 50, S A = 40.

Model of LGN, and input to V1
We distinguish different types of LGN neurons that innervate V1 neurons, lagged and nonlagged, transient and sustained, and On-and Off-cells, whose activity is characterized by the firing rates L L LGN ðx; y; tÞ, L LGN (x, y, t), L T LGN ðx; y; tÞ, L S LGN ðx; y; tÞ, L On LGN ðx; y; tÞ and L Off LGN ðx; y; tÞ, respectively. We assume that the non-lagged, transient and On-populations are identical, i.e. L T LGN ðx; y; tÞ ¼ L On LGN ðx; y; tÞ ¼ L LGN ðx; y; tÞ. At any given time moment, the firing rate of an LGN neuron is calculated as a convolution of the stimulus with a receptive field (RF). The firing rate L LGN (x, y, t) at any given time moment t is a rectified convolution of the stimulus S (x, y, t) distributed across the retina with a spatial-temporal receptive field (RF) D(x, y, t) [85]. The firing rate is calculated as follows: where the receptive field D(x, y, t)is approximated as the difference of central (excitation for On-neurons) D cen (x, y, t) and surround (inhibition for On-neurons) D sur (x, y, t), each being separable in space and time, and round. The spatial kernel displays a center-surround structure determined by axons from retinal ganglion cells described as a difference of two axisymmetric Gaussian functions [85]: Dðx; y; tÞ ¼ D cen ðx; y; tÞ À D sur ðx; y; tÞ ¼ D cen t ðtÞ À B p s 2 An example of the LGN neuron's RF is shown in Fig 1A, left. The temporal component is determined by alpha-functions, as in [85], with a coefficient controlling the balance between early and late components of LGN cell responses.
Here D cen t ðtÞ and D sur t ðtÞ are the temporal evolution functions D cen t ðtÞ ¼ t=t 2 cen expðÀ t=t cen Þ À t=t 2 late expðÀ t=t late Þ: ð5Þ and D cen t ðtÞ ¼ t=t 2 sur expðÀ t=t sur Þ À t=t 2 late expðÀ t=t late Þ ð6Þ with the time constants τ cen , τ sur and τ late , respectively. The spatial scales of the RF are σ cen and σ sur . The parameters were: τ cen = 10ms, τ sur = 20ms, τ late = 64ms, σ cen = 0.3 deg., σ sur = 1.5 deg. The parameter are from [85,86]. In practice, in order to optimize the memory usage, the convolutions with the time kernels can be substituted by integration of equivalent ordinary differential equations. For example, a variable Qðx; y; tÞ ¼ Z t 0 dt D cen t ðx 0 ; y 0 ; tÞ Sðx 0 ; y 0 ; t À tÞ can be alternatively calculated from the following equation The sustained cells have in 3.5 times slower kinetics [13]. The firing rate L S LGN ðx; y; tÞ is calculated from the Eqs (1)-(6) with the time constant t S cen ¼ 3:5 t cen , t S sur ¼ 3:5 t sur , t S late ¼ 3:5 t late instead of τ cen , τ sur and τ late . The lagged cell activity is L L LGN ðx; y; tÞ ¼ L LGN ðx; y; t À dÞ. For the sake of computational efficiency, the activity L LGN as a function of the delay is calculated approximately as two first terms of the Taylor expansion: L LGN ðx;ỹ; t À dðx; y;x;ỹÞÞ � L LGN ðx;ỹ; tÞ À dðx; y;x;ỹÞ @L LGN ðx;ỹ; tÞ @t For Off-neurons, the firing rate is calculated as L Off LGN ðx; y; tÞ ¼ ½ÀL LGN ðx; y; tÞ� þ . Thalamic input to V1. Orientation and direction selectivity is determined by a footprint of the LGN-to-V1 projections as well as by properties of the LGN neurons, which are different for the considered DS mechanisms. The thalamic input is determined as the firing rate φ Th E (x, y, t), which is a convolution of the LGN neuronal activity of the lagged and non-lagged, or transient and sustained, or On-and Off-cells, A certain form of the footprint function D LGNÀ V1 ðx; y;x;ỹ; iðx; y;x;ỹÞÞ depends on the mechanism of DS.
In the L-N and T-S DS models, the footprint of a direction-selective V1 neuron splits into two halves along the axis of elongation, each sending signals from either lagged or nonlagged (transient/sustained) thalamic cells (Fig 1B). For the On-Off model, the footprint of a V1 cell consists of an elongated On zone and a circular Off zone. The width and elongation of the footprints were 0.3 deg. and 0.8 deg., respectively. The shift was 0.6. deg. The neighboring V1 neurons that belong to different orientational hypercolumns have footprints of similar shapes and prefer the same orientation but opposite directions of the stimulus movement. A kernel expression that determines the DS bias is the equation for the thalamic input into V1 neuron, φ Th E (x, y, t), which is given separately for each of the three mechanisms.
Lagged and nonlagged cell-based mechanism of DS. For the L-N cell-based mechanism of DS (L-N mechanism), two populations of L-N cells were considered. These populations are equally and homogeneously distributed across the LGN (Fig 1A, middle). The L-N cells have round, center-surround receptive fields (Fig 1A, left). In this particular version of the model, only On-cells responding to a bright stimulus in the center of RF and being inhibited in the surrounding area of RF are considered. Lagged cell activity is delayed by 40 ms according to estimations from [25].
In accordance with the footprint of the LGN-to-V1 projections (Fig 1A, middle), the V1 neurons that prefer a certain direction receive an input signal from the nonlagged cells located on one (right) side of the footprint and from the lagged cells located on the other (left) side, whereas V1 neurons preferring the opposite (downward) direction conversely receive an input signal from the nonlagged cells of the second (left) side and from the lagged cells of the first (right) side.
The footprints of V1 neurons are oriented according to the neuronal positions in pinwheels. The pinwheels with clockwise progression of orientation columns are adjacent to the ones with counterclockwise progression. The pinwheel-centers are distributed on the rectangular grid with the pinwheel radius R and indexed by i PW and j PW . The adjacent columns owing to different pinwheels have the same orientation preferences. The coordinates of the pinwheel-center are x PW = (2 i PW − 1)R, y PW = (2 j PW − 1)R. The orientation angle at the point (x, y) of V1 that belongs to the pinwheel (i PW , j PW ) is defined as θ(x, y) = arctan((y − y PW )/(x − x PW )). The progression is determined by the factor ðÀ 1Þ where L LGN is the activity of the LGN neuron (x;ỹ); D LGNÀ V1 ðx; y;x;ỹÞ is the LGN-to-V1 footprint with the width across preferred orientation σ pref and the width across orthogonal orientation σ orth ; dðx; y;x;ỹÞ is the delay that determines contributions of either lagged or nonlagged cells. Because the direction preferences are given by a combination of inputs from the L-N neurons according to the footprint (Fig 1B, left), DS is determined by the delay dðx; y;x;ỹÞ, which is formalized as follows: x 0 ¼ ðx À x cf Þcos y À ðỹ À y cf Þsin y ð13Þ y 0 ¼ ðx À x cf Þsin y þ ðỹ À y cf Þcos y: Here (x cf , y cf ) are the coordinates of the center of the footprint of V1 neuron in LGN, x cf = x cf (x, y), y cf = y cf (x, y).
Transient and sustained cell based mechanism of DS. In the T-S cell based mechanism of DS the input firing rate is where iðx; y;x;ỹÞ is the index of T or S neurons, whose firing rates are calculated with different time constants in the temporal kernel of the convolution. It is defined according to the footprint (Fig 1B, middle): In this case, D LGNÀ V1 ðx; y;x;ỹÞ, x 0 and y 0 are given by Eqs (11), (13) and (14).
On-and Off cell-based mechanism of DS. In the On-and Off cell-based mechanism of DS the input firing rate is x Off with Δ Off as the shift of the Off-subfield, and x 0 and y 0 given by Eqs (13) and (14).

Model of V1
The model of V1 cortical area activity is an extension of a previous model describing orientation processing [34]. Neurons in V1 show more complex RFs than LGN cells. The V1 neurons vigorously respond to gratings moving throughout their RFs. It was assumed that most of neurons in V1 are directionally selective [87,88] despite the fact that some neurons are tuned only to orientation [4]. V1 is modeled as a 2-dimensional continuum of neuronal populations. In simulations, the modeled cortical area of the cortex is 1 mm x 1.5 mm and includes six orientation hypercolumns. (Spatial coordinates as arguments of variables are omitted here.) Each point of the cortical continuum contains two neuronal populations, excitatory (E) and inhibitory (I), connected by AMPA (α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid), NMDA (N-methyl-d-aspartate), and GABA-A mediated synapses providing recurrent intracortical interactions and AMPA and NMDA for the geniculate input. The strengths of the external connections correspond to the pinwheel architecture, and thus neurons receive inputs in accordance with their orientation and direction preferences. The profile of the intracortical connections is isotropic, i.e., the maximum conductances depend on the distance between the pre-and postsynaptic populations. According to the definition of a neuronal population [89], neurons of one population located at a given point both receive a common input that comes from presynaptic populations and an individual noise that takes into account any differences in the intrinsic and input parameters of neurons as well as a spontaneous activity of ionic channels. The membrane potentials and ionic channel states of these neurons are dispersed due to the noise, and thus they are distributed in a space of neuronal refractority. Neurons that fire contribute to the population firing rate, which is the output measure of the population activity. The mathematical description of each population is based on the probability density approach [89,90], namely, the CBRD approach [37,38,91,92], where the neurons within each population are distributed according to their refractority states, which are characterized with the time elapsed since the last spike, t � . Single population dynamics are governed by the equations for neuronal density, the mean over noise realizations voltage, and the gating variables. The model for each E-and I-neuronal populations takes into account two neuronal compartments and a set of voltagegated ionic currents, including the adaptation currents.
According to the structure of interneuronal connections, the neuronal population firing rate determines the presynaptic firing rate, which in turn controls the dynamics of synaptic conductances. The presynaptic firing rate predetermined by the excitatory population firing rate determines the dynamics of AMPA and NMDA synaptic conductances. The inhibitory population controls GABA conductance. The synaptic conductances are the input signals received by the postsynaptic neuronal populations. The membrane voltage distribution across t � determines the output firing rate and so on.
Model neurons have two compartments: somatic and dendritic. The inhibitory synapses are assumed to be located at the soma, whereas the excitatory synapses are at the dendrites. Because in experiments the synaptic currents are usually measured at the soma, only somatic conductances are calculated in the model. The two-compartment model from [34,93] implicitly solves a reverse voltage-clamp problem, thus estimating the dendritic synaptic current. Each of the somatic AMPA, GABA, and NMDA conductances (the AMPA conductance mediated by the external thalamic terminals is marked by the index with a dash, i.e. g AMPA 0 ,E ) is calculated with a second order ordinary differential equation [34]. An input function to the equation is the presynaptic firing rate. Equations connecting the somatic firing rate with the presynaptic firing rate are as given in [34], whereas the presynaptic thalamic input has been modified as described. The full CBRD model for interacting adaptive regular spiking pyramidal cells and fast spiking interneurons is given below.
CBRD-approach. The conductance-based refractory density approach [37,38,91,94] considers a population of an infinite number of Hodgkin-Huxley-like neurons receiving both a common input and an individual for each neuron noise. In any arbitrary case of transient or steady-state stimulation the firing rate of such population can be quite precisely and computationally effectively calculated by solving a system of equations in partial derivatives, 1-d transport equations. The equations govern an evolution of neuronal states distributed in a phase space of the time elapsed since last spikes, t � . They contain the Hodgkin-Huxley equations for the membrane voltage and gating variables (except sodium ones), parameterized by t � , as well as the equation for the neuronal density in t � -space, ρ p (t, t � ), where the index of a population p substitutes for E or I. The output characteristic of the population's activity is the firing rate ν p (t), which is equal to ρ p in the state of a spike, t � = 0, i.e. ν p (t) � ρ p (t, 0).
Basic neurons have 2-compartments with the somatic and dendritic voltages U p (t, t � ) and U p d ðt;t � Þ. In comparison with one-compartment model, the extra parameters is the ratio of dendritic to somatic conductances γ and the dendritic length. The inhibitory synapses are located at soma, contributing into the somatic synaptic current I soma , whereas the excitatory synapses are at dendrites, determining the dendritic synaptic current I dendr . Due to the construction of the 2-compartment model [93], both type synaptic conductances are imposed to be somatic, in spite of the localization, in order to be compared with experimental whole-cell somatic registrations. Approximations of voltage-gated ionic currents I p voltageÀ gated ðU p ; t; t � Þ differ for excitatory and inhibitory neurons. Parameterized by t � , the governing equations are as follows: @r p @t þ @r p @t � ¼ À r p HðU p ;g p tot Þ; ð24Þ Voltage-dependent channels of excitatory neurons. For E-neurons, the approximations for the components of the total (except sodium), voltage-gated current I E voltageÀ gated ðU E ; t; t � Þ ¼ À I DR À I A À I M À I AHP are mainly based on the CA1 pyramidal cell model from [95], where instead of full description of calcium dynamics and calcium-dependent potassium currents a cumulative after-spike hyperpolarization (AHP) current, that provides an effect of slow spike timing adaptation [96]. The set of ionic currents includes the voltage-dependent potassium currents I DR and I A responsible for spike repolarization, the slow potassium current I M that contributes to spike frequency adaptation and the potassium current I AHP , implicitly dependent on calcium dynamics and contributing to spike frequency adaptation. Approximating formulas for the currents I Na , I DR , I A and I M are taken from [95]; the approximation for I AHP is given in [96].
The voltage-dependent potassium current I DR : @x @t þ @x @t � ¼ @y @t þ @y @t � ¼ t x ¼ 1=ða þ bÞ þ 0:8 ms; x 1 ¼ a=ða þ bÞ; a ¼ 0:17 exp ððU E þ 5Þ � 0:090Þ ms À 1 ; b ¼ 0:17 exp ðÀ ðU E þ 5Þ � 0:022Þ ms À 1 ; t y ¼ 300 ms; The voltage-dependent potassium current I A : @x @t þ @x @t � ¼ @y @t þ @y @t � ¼ The voltage-dependent potassium current I M : @x @t þ @x @t � ¼ @y @t þ @y @t � ¼ t x ¼ 1=ða þ bÞ þ 8 ms; t y ¼ 1000 ms; boundary conditions for eqs.(A2-A14) at t � = Δt AP which are as follows: The reset values for the fast gating variables in Eqs (45) and (46) were obtained with the basic single neuron model. With a rather arbitrary input providing a spike, these values were measured at the moment of a voltage maximum at the spike. The reset level for each slow conductance in the CBRD model was calculated as a sum of its value at a peak of spike-release distribution in the t � -space and an increment at spike: where t �� is such that rðt; t � Þ Hðt; t � Þ: The increment values for the slow gating variables in Eqs (47) and (48) were also measured at a single spike of the single neuron model.
Parameters of basic neurons are as follows: � g DR ¼ 0:76 mS=cm 2 ; � g A ¼ 4:36 mS=cm 2 ; � g M ¼ 0:76 mS=cm 2 ; � g AHP ¼ 0:6 mS=cm 2 ; Here S is the membrane area. The dependence of V th (t � ) is taken from a full single neuron model [37], allowing to take into account the effect of sodium channel inactivation on the threshold dynamics [98]. σ ν is the noise amplitude meaning the dispersion of individual neuron's voltage fluctuations in a stationary state. Its scaling with g syn approximately reflects the fact of the synaptic noise increase with the increase of mean synaptic drive [99]. Stochastic input to E2-neurons I noise was modeled as Ornstein-Uhlenbeck process with the time correlation 10 ms and the dispersion 20pA for the ID-regime simulation and 40pA for the IID-regime simulation.
The equations for the input synaptic conductances are given below, as well as the values of the reversal potentials. When calculating the dynamics of a neural population, the integration of Eqs (25), (26) and (28)- (39) determines the evolution of the distribution of voltage U E across t � . Then, the effect of crossing the threshold and the diffusion due to noise are taken into account by H-function, Eq (27), substituted into the equation for neuronal density, Eq (24). The integral Eq (40) results in the output firing rate ν E (t).
Lognormal distribution of synaptic weights within each population. The CBRDapproach is generalized to the case of lognormal distribution of synaptic weights within each p-population [100]. In this case, instead of equal total synaptic current, neurons receive lognormally distributed current. For the current scaled by its mean across the distribution, η, the distribution is The membrane potential of neurons parameterized with η, U p Z , can be found as U p Z ðt; t � Þ ¼ ðU p ðt; t � Þ À U p free ðt � ÞÞ x þ U p free ðt � Þ; ð50Þ where U p free ðt � Þ is the unperturbed potential defined for zero synaptic input. The density of neurons parameterized by η and distributed in the phase space t � is denoted as r p Z ðt;t � Þ. Calculation of r p Z ðt;t � Þ requires solving of a continuum of Eq (1) (or Eq (22)) for r p Z instead of ρ p with HðU p Z ;dU p Z =dtÞ. The output firing rate is defined as In numerical simulations, we set the parameter of the lognormal distribution σ LN = 0.75 and discretized the η-space by 15 intervals.

Modeling of V1 silencing and recordings from single neurons
To reveal the contribution of LGN-to-V1 projections into DS, Lien and Scanziani [13,41] optogenetically inhibited the visual cortex of mice with interneurons expressing channelrhodopsin. Monochrome light led to the activation of interneurons and the subsequent inhibition of their target neurons, including the excitatory neurons in layer 4. In this way, the contribution of cortico-cortical interconnections was minimized, thus revealing the contribution of LGN-to-V1 connections. In this model, we substituted the effect of light by the effect of the depolarizing current 100 pA applied to the interneurons.
To register the behavior of a single representative neuron, we simulated the neuron with the basic Hodgkin-Huxley-like model that was used for the construction of the CBRD model. Because in the cortex each single neuron receives synaptic inputs determined by the activity of neuronal populations, the representative neuron was simulated with synaptic conductances obtained from the equations of the population activity. In addition, the neuron can be controlled by a patch-clamp electrode that injects a current in the current-clamp mode or holds the membrane potential in the voltage-clamp mode. For instance, experimental registrations by Lien and Scanziani [13,41] were done with the voltage clamped at the level of the GABAergic channel reversal potential, thus measuring EPSC. In simulations, it is equal to -77mV. Because the representative neuron is simply used to observe the population dynamics, it does not affect the network. Locations of the representative neurons are marked with white dots in Figs 3C, 4A and 4B. For the sake of statistical measurements of the ratio of the charges supplied with the EPSC for intact versus inhibited cortex Q control /Q inhibited (Fig 8), we simulated patchclamp recordings in the voltage-clamp mode in all nodes of the spatial grid. Specifically, for each of the DS models we: (i) calculated EPSCs in all nodes of the spatial grid, (ii) calculated Q control and Q inhibited by integrating over the interval from 0 to 1600ms the EPSCs with the subtracted background level before the stimulation (-135 and -125pA for the intact and inhibited cortex, respectively), and (iii) calculated the mean value and the dispersion of Q control /Q inhibited across all nodes with the same preferred orientation as in the site shown in Fig 3C.

Data analysis
For each point of the computational domain the field of the excitatory population firing rate ν E (t, x, y) and the membrane potential � U E ðt; x; yÞ were averaged over the time interval from 600 to 1600ms since the stimulus onset. Based on modulations of these signals, the field of the preferred directions was found. At each point, the direction selectivity index (DSI) was calculated as in [13]: DSI = (Pref − Nonpref)/Pref, where Pref and Nonpref are the amplitudes of the modulated component of the response (firing rate or voltage) to gratings of the preferred and opposite directions, The mean DSIs were calculated by averaging over the entire computational domain.
The charge Q was calculated as an integral of EPSC(t) − EPSC 0 over the time period from 0 to 1600ms, where EPSC 0 is the steady-state response to the initial gray screen. With the basic parameters it was equal to -135 and -125pA for the intact and inhibited cortex, respectively.