ABSTRACT Biophysically detailed description of the mechanisms of orientation and direction selectivity of visual cortex neurons is still in development. We have incorporated a simplified, filter-based description of retino-thalamic visual signal processing including a mechanism of direction selectivity into a conductance-based model of neuronal population activity of the primary visual cortex. Our simulations of the cortical response to drifting gratings have verified that the mechanism based on asymmetrical projections of lagged and non-lagged thalamic neurons to the cortex provides the direction selectivity to the extent consistent with experimental evidences, and that the biophysical model realistically reproduces such characteristics of the visual cortex activity as membrane potential, firing rate, synaptic conductances etc. INTRODUCTION Primary visual cortex neurons are selective to various characteristics of the stimulus: orientation, direction of motion, color, etc. [14]. Up to now, some important questions about the underlying mechanisms of these properties remain unanswered. Mathematical modeling is a perspective method of the V1 studies. Previously we have offered a biophysically detailed model reproducing the activity of V1 and its individual neurons in response to the stationary gratings (the so-called orientation tuning) [8, 11]. Here, we extend this model to the case of moving stimuli taking into account the direction selectivity (DS) of V1 neurons. Direction-selective neurons exist in almost all animals that have a visual apparatus. As discussed, for example, in [3], there are three principal models of DS, namely: 1) Hassenstein-Reichardt; 2) Barlow-Levick, and 3) fully opponent Hassenstein-Reichardt model. Most models of DS include a time delay between the spatially shifted inputs in a cortical cell [1-2]. The physiological mechanism of this delay formation is disclosed, firstly, in [4] and then, in more details, in [24], where with the help of intracellular _in vivo_ registration it was demonstrated that the lateral geniculate nucleus (LGN) neurons fall into two classes: lagged and non-lagged cells; and a delay of the lagged neurons is determined by the effects of the inhibitory-excitatory synaptic complexes formed on synaptic axonal terminals of retinal ganglion cells in LGN. [17] and [23] proposed a complex schematic model of DS on a basis of the specific convergent projections of the signals from lagged and non-lagged LGN cells, as well as on the intracortical interactions. So, the temporal and spatial differences of LGN inputs appear to provide DS in V1. We have taken into account these mechanisms so as to reproduce DS in our biophysically detailed model of V1. For modeling V1 we used the conductance-based refractory density approach (CBRD) (Chizhov and Graham 2007), which allows us to verify the role of not only network structure or interactions but also cellular properties, such as ionic currents, in V1 functioning. METHODS The structure of the model providing the direction selectivity in V1 on account of geniculate timing Figure 1 shows schematically the proposed model implementing the lagged/non-lagged distinction of LGN inputs to a V1 cell as the mechanism of DS. We consider the following three spatially two-dimensional levels: 1) Screen showing visual stimuli, see the left panel on Figure 1. We used a visual stimulus in the form of the drifting sinusoidal luminance gratings, which is characterized by the spatial and temporal frequencies and by the orientation. 2) LGN, consisted of neurons with the center-surround receptive field (RF), see the central panel on Figure 1. RF of a neuron is a specific location in the visual field, the stimulus in which excites this neuron the most strongly. Neurons with the center-surround RF structure respond strongly to a stimulus in the center of RF and are inhibited by a stimulus in the surround of RF. The center-surround structure, which is shown schematically on the left panel in Figure 1, is typically described by a difference of Gaussians, see [12, eq. 2.46]. Its spatial and temporal terms are usually assumed to be independent. We describe the temporal component by a function as in [12, eq. 2.47], but adding a coefficient controlling the balance between early and late components of LGN cell response. The firing rate of an LGN neuron at any given time is expressed as a convolution of RF with the stimulus, eq. 2.1 in [12], rectified then by a non-linear function which zeroed the firing rate if the convolution is negative. The model of LGN cells is described in detail in [25]. We also distinguish between lagged (with time delay 30 μs) and non-lagged (without delay) cells in LGN. 3) V1, consisted of neurons with more sophisticated RF, see the right panel on Figure 1. We consider V1 neurons being into the orientation hypercolumns, at that the neurons of neighboring hypercolumns preferring opposite directions of movement. The details of the model which provides DS of V1 neurons are as follows. V1 neurons receive inputs from LGN cells which are located within the footprint described by asymmetric Gaussian. The neurons of two neighboring hypercolumns preferring the same orientation but opposite directions have the same footprint. The difference of preferred directions is provided by the following property of input. The footprint is split along more elongated axis into two halves. V1 neurons preferring one direction get input signal from non-lagged cells located in one (right) side of the footprint and from lagged cells located in the other (left) side, whereas V1 neurons preferring the opposite (downward) direction conversely get input signal from non-lagged of the second (left) side and from lagged of the first (right) side. Neurons in V1 have more complex RF structures than LGN cells and respond stronger to moving (or steady) gratings within their RFs then to spots. We assume all neurons in V1 to be tuned to direction, in spite of the fact that some neurons are tuned only to orientation [15]. Biophysically detailed mathematical model of V1 The model of V1 is based on the CBRD model describing orientation processing [11]. V1 is modeled as a continuum in the two-dimensional space along the cortical surface. Each point of the continuum contains 2 populations of neurons, excitatory and inhibitory. We model AMPA, NMDA and GABA-mediated synaptic currents in neurons in order to describe recurrent interactions, whereas for the LGN input we consider, certainly, only AMPA and NMDA-mediated currents. The strengths of the external connections are distributed according to the pinwheel architecture, thus neurons receive inputs according to their orientation and direction preferences. The strengths of the intracortical connections, the maximum conductances, are distributed according to the locations of pre- and postsynaptic populations. The most strong connections are local and isotropic, thus depending only on the distance between the pre- and postsynaptic populations and decaying with the distance. Patchy long-distance connections are also taken into account, but they are much weaker. Within one population located at a given point the neurons receive the common inputs from presynaptic populations and an individual noise. The membrane potentials and ionic channel states of these neurons are dispersed due to the noise, being rather distributed in a space of neuron refractority characterized by a time passed since the last spike of each neuron. Neurons with the voltages crossing a spike threshold contribute into the population firing rate, which is the output characteristic of the population activity. According to the structure of connections described above, the firing rate determines the presynaptic firing rate, which, in turn, controls the dynamics of synaptic conductances. The presynaptic firing rate pre-determined by the excitatory population firing rate determines the dynamics of AMPA and NMDA synaptic conductance, the inhibitory one controls GABA conductance. The synaptic conductances are the inputs for the postsynaptic neuronal populations. Subjected to the electrotonic effect provided by the dendritic compartment, the membrane voltage distribution across the refractory time space is controlled by the synaptic conductances. Again, it determines the output firing rate and so on. The external presynaptic firing rate is distributed according to retinotopic projection. The orientation preferences are set by the pinwheel architecture. 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 iPW and jPW. The adjacent columns owing to different pinwheels have the same orientation preferences. The coordinates of the pinwheel-center are xPW = (2 iPW − 1)R,  yPW = (2 jPW − 1)R. The orientation angle for the point (x, y) of V1 which belongs to the pinwheel (iPW, jPW) is defined as θ(x, y)=arctan((y − yPW)/(x − xPW)). The progression is determined by the factor ( − 1)iPW + jPW. Finally, the input firing rate is \varphi _{Th\; E} (x,y,t)=\iint dd\, D^{LGN-V1} (x,y,,)\, L_{LGN} (,,t-\delta(x,y,,)) , where $D^{LGN-V1} (x,y,,)={\pi \; \sigma _{pref} \sigma _{orth} } \exp \left(- }{\sigma _{pref}^{2} } - }{\sigma _{orth}^{2} } \right)$, where LLGN is the activity of the non-lagged LGN neurons at ($,$), (xcf, ycf) are the coordinates of the centre of the footprint of V1 neuron in LGN, xcf = xcf(x, y), ycf = ycf(x, y); δ is the delay of activity of the lagged neurons. The direction preferences are given by a combination of inputs from lagged and non-lagged LGN neurons: L_{LGN} &&(,,t-\delta(x,y,,)) \approx \nonumber \\ &&\left[L_{LGN} (,,t)-\delta(x,y,,) (,,t)}{\partial t} \right]_{+} \delta(x,y,,) = 30~ms,\quad \quad (-1)^{i_{PW} +j_{PW} } x'>0; \quad 0, \quad , x'=(-x_{cf} )\cos \theta -(-y_{cf} )\sin \theta , y'=(-x_{cf} )\sin \theta +(-y_{cf} )\cos \theta . The mathematical description of each population is based on the probability density approach, namely CBRD [5-7], where the neurons within each population are distributed according to their phase variable, the time elapsed since their last spikes, t*. The CBRD for adaptive regular spiking pyramidal cells and fast spiking interneurons is given in [11]. The model of an excitatory neuron takes into account a set of voltage-gated ionic currents. It is based on the CA1 pyramidal cell model from [13], where instead of full description of calcium dynamics and calcium-dependent potassium currents the cumulative afterhyperpolarization (AHP) current has been introduced according to the work [16], which provides the effect of slow spike adaptation. The inhibitory synapses are assumed to be located at soma, whereas the excitatory synapses are at dendrites. As the synaptic conductances are usually measured at soma, the two-compartment model from [10, 11] implicitly solves the reverse voltage-clamp problem, estimating the dendritic synaptic current. To describe the synaptic AMPA, GABA and NMDA dynamics the second order ordinary differential equations for the non-dimensional synaptic conductances were used, see [9]. The input function to these equations are the presynaptic firing rate. The calculations of presynaptic firing rates for horizontal connections are presented in [11], whereas the presynaptic thalamic input has been modified as described above. The basic parameters are given in [11]. The specific parameters of our simulation as follows. Synaptic parameters: $_{ampa\; Th\; E} =0.4$ mS/cm², $_{ampa\; E\; E} =0.4$ mS/cm², $_{ampa\; E\; I} =0.4$ mS/cm², $_{gaba\; I\; E} =0.8$ mS/cm², $_{gaba\; I\; I} =0.2$ mS/cm², $_{nmda\; Th\; E} =0.8$ mS/cm², $_{nmda\; E\; E} =1.6$ mS/cm², $_{nmda\; E\; I} =1.6$ mS/cm², τrAMPA = 1.7ms, τdAMPA = 8.3 ms, τrGABA = 0.5 ms, τdGABA = 20 ms, τrNMDA = 6.7 ms, τdNMDA = 100 ms, Mg = 2 mM, VAMPA = VNMDA = 0, VGABA = −77 mV. Region geometry: cortex region is 1 mm x 1.5 mm, containing 2x3 regularly distributed pinwheels. Numerical parameters: time step 0.1 ms, spatial grid 23x35, the t*-space was limited to 100 ms and discretized by 100 points.

INTRODUCTION Primary visual cortex neurons are selective to various characteristics of the stimulus: orientation, direction of motion, color, etc. [14]. Up to now, some important questions about the underlying mechanisms of these properties remain unanswered. Mathematical modeling is a perspective method of the V1 studies. In our previous works it was developed a biophysically detailed V1 model that reflects orientation tuning [11], [8]. Here we develop this model, supplying it with direction selectivity. Direction-selective (DS) neurons exist in almost all animals that have a visual apparatus. Principal DS models (Hassenstein-Reichardt model, Barlow-Levick model and fully opponent Hassenstein-Reichardt model) are well known and discussed in [3]. Most of the DS models include a time delay between the spatially shifted inputs in cortical cell [1-2]. There is an important question about the source of the time delay. The physiological mechanism of this delay formation is disclosed in experimental work [4] and in more details in [24], where with the help of intracellular _in vivo_ registration it was demonstrated that LGN neurons fall into two classes of lagged and non-lagged cells, and the delay is determined by the effects of the inhibitory-excitatory synaptic complexes formed on synaptic axonal terminals of retinal ganglion cells in LGN. In the papers [17], [23], it was proposed a complex schematic model of DS, on a basis of a specific convergent projections of the signals from LGN cells with and without a delay and intracortical interactions. The temporal difference and the structure of LGN inputs provides DS in V1. We have taken into account these mechanisms to reproduce DS in our biophysically detailed visual cortex model. METHODS Lagged and non-lagged cell based mechanism of DS A simple scheme of the proposed model is presented in Figure 1. The model reproduces the direction selectivity in V1, in response to moving gratings. The model contains three two-dimensional levels: the screen with moving gratings (left panel); LGN (central panel) consisted of neurons with the center-surround receptive field (RF); V1 (right panel) consisted of neurons with more sophisticated RF. The center-surround RF structure is typically described as a difference of Gaussian functions. An example of the LGN neuron’s RF is shown in Fig. 1 (on the left). There are two types of LGN cells – lagged (with time delay 30 μs) and non-lagged (without delay). At the level of V1, neurons are organized into orientation hypercolumns, and the model was constructed so that the neurons of neighboring hypercolumns would prefer opposite directions of movement. The details of the model which provides DS in V1 are as follows. V1 neurons receive inputs from LGN cells which are located within the footprint described by asymmetric Gaussian. The neurons of two neighboring hypercolumns preferring the same orientation but opposite directions have the same footprint. The difference of preferred directions is provided by the following property of input. The footprint is splitted along more elongated axis into two halves. V1 neurons preferring one direction get input signal from non-lagged cells located in one (right) side of the footprint and from lagged cells located in the other (left) side, whereas V1 neurons preferring the opposite (downward) direction conversely get input signal from non-lagged of the second (left) side and from lagged of the first (right) side. Stimulus Presented stimulus is a moving gradient characterised by the orientation, spatial frequency and temporal frequency. LGN Neurons in LGN display a center-surround receptive field (RF) structure, responding most strongly to visual input in the center of a specific location of the visual field and being inhibited by visual input in the surround. This center-surround RF structure is typically described as a difference of Gaussians (see [12] eq. 2.46). The spatial and temporal aspects of RF are usually assumed to be independent. We describe the temporal component by a function as in [12], eq. 2.47, but add a coefficient controlling the balance between early and late components of LGN cell response. The firing rate of an LGN neuron at any given time is expressed as a convolution of RF with the stimulus, eq. 2.1 in [12], rectified then by a non-linear function which zeroed the firing rate if the convolution is negative. The model of LGN cells is described in detail in [25]. V1 Neurons in V1 show more complex RFs than LGN cells and respond vigorously to moving (or steady) oriented gratings within their RFs. We assume all neurons in V1 to be tuned to direction, in spite of the fact that some neurons are tuned only to orientation [15]. The model of V1 is based on the CBRD model describing orientation processing [11]. V1 is modeled as a continuum in 2-dimensional space along the cortical surface. Each point of the continuum contains 2 populations of one-type’s neurons, excitatory and inhibitory. The populations are interconnected by synapses of AMPA, NMDA and GABA types and have AMPA- and NMDA-synapses formed by external axons. The strengths of the external connections are distributed according to a pinwheel architecture, thus neurons receive inputs according to their orientation and direction preferences. The strengths of the intracortical connections, the maximum conductances, are distributed according to the locations of pre- and postsynaptic populations. The most strong connections are local and isotropic, thus depending only on the distance between the pre- and postsynaptic populations and decaying with the distance. Patchy long-distance connections are also taken into account, but they are much weaker. Within one population located at a given point the neurons receive the common inputs from presynaptic populations and an individual noise. The membrane potentials and ionic channel states of these neurons are dispersed due to the noise, being rather distributed in a space of neuron refractority characterized by a time passed since the last spike of each neuron. Neurons with the voltages crossing a spike threshold contribute into the population firing rate, which is the output characteristic of the population activity. According to the structure of connections described above, the firing rate determines the presynaptic firing rate, which, in turn, controls the dynamics of synaptic conductances. The presynaptic firing rate pre-determined by the excitatory population firing rate determines the dynamics of AMPA and NMDA synaptic conductance, the inhibitory one controls GABA conductance. The synaptic conductances are the inputs for the postsynaptic neuronal populations. Subjected to the electrotonic effect provided by the dendritic compartment, the membrane voltage distribution across the refractory time space is controlled by the synaptic conductances. Again, it determines the output firing rate and so on. The external presynaptic firing rate is distributed according to retinotopic projection. The orientation preferences are set by the pinwheel architecture. 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 iPW and jPW. The adjacent columns owing to different pinwheels have the same orientation preferences. The coordinates of the pinwheel-center are xPW = (2 iPW − 1)R,  yPW = (2 jPW − 1)R. The orientation angle for the point (x, y) of V1 which belongs to the pinwheel (iPW, jPW) is defined as θ(x, y)=arctan((y − yPW)/(x − xPW)). The progression is determined by the factor ( − 1)iPW + jPW. Finally, the input firing rate is \varphi _{Th\; E} (x,y,t)=\iint dd\, D^{LGN-V1} (x,y,,)\, L_{LGN} (,,t-\delta(x,y,,)) , where $D^{LGN-V1} (x,y,,)={\pi \; \sigma _{pref} \sigma _{orth} } \exp \left(- }{\sigma _{pref}^{2} } - }{\sigma _{orth}^{2} } \right)$, where LLGN is the activity of the non-lagged LGN neurons at ($,$), (xcf, ycf) are the coordinates of the centre of the footprint of V1 neuron in LGN, xcf = xcf(x, y), ycf = ycf(x, y); δ is the delay of activity of the lagged neurons. The direction preferences are given by a combination of inputs from lagged and non-lagged LGN neurons: L_{LGN} &&(,,t-\delta(x,y,,)) \approx \nonumber \\ &&\left[L_{LGN} (,,t)-\delta(x,y,,) (,,t)}{\partial t} \right]_{+} \delta(x,y,,)=30~ms\quad \quad (-1)^{i_{PW} +j_{PW} } x'>0; \quad 0, \quad , x'=(-x_{cf} )\cos \theta -(-y_{cf} )\sin \theta , y'=(-x_{cf} )\sin \theta +(-y_{cf} )\cos \theta . The mathematical description of each population is based on the probability density approach, namely, the conductance-based refractory density (CBRD) approach [5-7], where the neurons within each population are distributed according to their phase variable, the time elapsed since their last spikes, $\ts$. The CBRD for adaptive regular spiking pyramidal cells and fast spiking interneurons is given in [11]. The model of an excitatory neuron takes into account a set of voltage-gated ionic currents. It is based on the CA1 pyramidal cell model from [13], where instead of full description of calcium dynamics and calcium-dependent potassium currents the cumulative after-hyperpolarization (AHP) current has been introduced according to the work [16], which provides the effect of slow spike adaptation. The inhibitory synapses are assumed to be located at soma, whereas the excitatory synapses are at dendrites. As the synaptic conductances are usually measured at soma, the two-compartment model from [10, 11] implicitly solves the reverse voltage-clamp problem, estimating the dendritic synaptic current. To describe the synaptic AMPA, GABA and NMDA dynamics the second order ordinary differential equations for the non-dimensional synaptic conductances were used, see [9]. The input function to these equations are the presynaptic firing rate. The calculations of presynaptic firing rates for horizontal connections are presented in [11], whereas the presynaptic thalamic input has been modified as described above. The basic parameters are given in [11]. The specific parameters of our simulation as as follows. Synaptic parameters: $_{ampa\; Th\; E} =0.4$ mS/cm, $_{ampa\; E\; E} =0.4$ mS/cm, $_{ampa\; E\; I} =0.4$ mS/cm, $_{gaba\; I\; E} =0.8$ mS/cm, $_{gaba\; I\; I} =0.2$ mS/cm, $_{nmda\; Th\; E} =1.6$ mS/cm, $_{nmda\; E\; E} =1.6$ mS/cm, $_{nmda\; E\; I} =1.6$ mS/cm, τrAMPA = 1.7ms, τdAMPA = 8.3 ms, τrGABA = 0.5 ms, τdGABA = 20 ms, τrNMDA = 6.7 ms, τdNMDA = 100 ms, Mg = 2 mM, VAMPA = VNMDA = 0, VGABA = −77 mV. Region geometry: cortex region is 1 mm x 1.5 mm, containing 2x3 regularly distributed pinwheels. Numerical parameters: time step 0.1ms, spatial grid 24x36, the $\ts$-space was limited to 100 ms and discretized by 100 points. RESULTS The constructed model proposes the detailed, conductance-based description of neuronal population activity in V1 and simple, filter-based description of retino-thalamic visual signal processing. The generalized model keeps the advantages of the former model [11] to relate different experimental observations obtained in slices and in vivo in visual cortex and allows us to consider experiments with moving stimuli. Here we verify by simulations (Fig.2) that the mechanism based on asymmetrical projections of lagged and non-lagged thalamic neurons to the cortex does allow the biophysical model to reproduce the direction selectivity in an extent consistent with experimental evidences, e.g. [19]. Also, we consider effects of surround suppression, stimulus locking, apparent motion and divisive normalization. Directional selectivity In our numerical experiments we presented moving gratings and reproduced the visual cortex activity. The modeled area of the cortex was as large as 1mm x 1.5mm and included 6 orientation hypercolumns. Spatio-temporal distribution of neuronal firing activity reveals the tuning of neurons to the direction of grating movement. As seen in Fig.2, the spike response of a neuron to the preferred and non-preferred directed stimuli were different (compare bottom plots in Fig.2C and D). This excitatory neuron located in the orientation column that prefers horizontal orientation. Its preferred direction is upward. A nearby neuron of the neighbor, lower orientation hyper-column prefers opposite direction (traces are not shown). Such tuning is provided by the cumulative effects of the described in Section 2.4 projection of LGN to V1 (Fig.1), the threshold rectification in V1 neurons and the recurrent intracortical inter-actions. Indeed, the light intensity was the same for the both stimuli (top frames in Fig.2A,B and top plots in Fig.2C,D). According to the model construction, the signal in LGN is not tuned (not shown) but its projection to V1 excitatory neurons is tuned (compare second to bottom plots in Fig.2C and D). The mean membrane voltage of the excitatory neurons is even more tuned (compare third to bottom plots in Fig.2C and D), which is due to the effect of recurrent synaptic inputs (forth to bottom plots in Fig.2C and D). The inhibitory interneuronal activity is reflected in the GABAergic conductance at the excitatory neurons. The threshold rectification producing the firing rate (fifth to bottom plots in Fig.2C and D) sharpens the voltage. The firing rate determines the probability for a single neuron of a population to generate spikes. With known synaptic conductances, a behavior of a single representative neuron is recon-structed and shown in bottom plots of Fig.2C and D. Stimulus locking The effect of long response to short stimulus called here as stimulus locking is the most challenging. It observed in experiments in cats and monkey . The response to 50 ms-spot of light lasts as long as 300 ms. At the same time, it is not trivial to find such long characteristic time scale in the mechanisms involved in generation of the response, except the adaptation which should lead rather to damping than prolongation of the response. However, our model is able to reproduce the effect ([Fig_SL], bottom trace). We have found that the most critical factor for this effect is the recurrent excitation. The effect was observed for the $_{ampa\; E\; I}$ from the range of 0.4 ÷ 0.5 mS/cm². The mechanism might be explained with the vector plot ([Fig_VectorPlot]). For the particular parameters chosen, the tuning of the cortical excitability is characterized with a basin of slow change of states in the phase plot of νE, νI (green bounded region). A stimulus abruptly shifts the state of the system from low rates (left bottom point) to high rates (right top point). Then the system gradually travels from this point down to the initial state (orange dotted line). Because the trajectory goes through the low change rate basin, the excitatory and inhibitory firing rates change slowly in the interval between 100 and 300 ms, as seen in Fig.[Fig_SL]. Apparent motion The stimulus locking effect allows the cortex to keep excitation for a while and thus link the stimulus with the following stimuli. In particular, this property is expressed in apparent motion, which can be observed, for example, if compare the activity in response to a moving spot (Fig.[Fig_MovSpot]) with that to a jumping spot (Fig.[Fig_3Spots]). The responses are similar in the sense that the focus of activity gradually travels from one cite to another without jumps. Surround suppression The surround suppression is revealed in experiments with different size stimuli. The response of the cortex non-monotonically depends on stimulus size [Rubin et al. 2015], such that the amplitude of the response to a larger stimulus might occur smaller than that to smaller stimulus. It is reproduced in our model (Fig.[Fig_SS]). Divisive normalization Sharpening by saccades As observed in , saccades increase the precision of the cortical response to stimulus. We reproduce the effect in simulations, as seen in comparison of Fig.[Fig_Saccades] with Fig.[Fig_noSaccades]. High amplitude bursts of the excitatory firing rate occur after rapid change of the “eye” position (Fig.[Fig_XYsaccades]), if the new position corresponds to bright stripe of the gratings. DISCUSSION We have introduced the direction selectivity in the model by implying the existence of a direction selectivity map, however such map is not observed in all animals. It is present in rabbits and ferrets. In monkeys, the directional map is absent in V1 but present in V2 [18]. With or without maps in V1, DS is a common feature of visual cortex neurons for various animals. It is quite possible that the mechanism of DS in animals without the directional map is similar to the mechanism in animals with such a map, like the one implemented in our model. The main mechanism of DS, according to our model, is the same as in the principle Hassenstein-Reichardt model based on the delayed excitation. Namely, the time delay between lagged and non-lagged neurons of LGN provides DS in V1 neurons. At the same time, Barlow-Levick mechanism also contributes to the direction tuning in V1 neurons because of the intracortical inhibition that is taken into account in our model. Concluding, the proposed model might be recommended as a tool to systematize various experimental observations focused to study mechanisms of the primary vision. ACKNOWLEDGEMENTS This work was supported by Grant No. 15-04-06234a from the Russian Foundation for Basic Research to AVC. REFERENCES 1. Adelson, E.H., Bergen, J.R.: Spatiotemporal energy models for the perception of motion. J. Opt. Soc. Am. A. 2, 284–299 (1985) 2. Barlow, H.B., Levick, W.R.: The mechanism of directionally selective units in the rabbit’s retina. J. Physiol. (Lond.). 178, 477–504 (1965) 3. Borst A., Helmstaedter M.: Common circuit design in fly and mammalian motion vision. Nat Neurosci. 18(8), 1067-1076 (2015) 4. Cai D., DeAngelis G.C., Freeman R.D.: Spatiotemporal receptive field organization in the lateral geniculate nucleus of cats and kittens. J. Neurophysiol. 78(2), 1045-1061 (1997) 5. Chizhov A.V., Graham L.J., Turbin A.A.: Simulation of neural population dynamics with a refractory density approach and a conductance-based threshold neuron model. Neuro-computing. 70(1), 252-262 (2006) 6. Chizhov A.V., Graham L.J.: Efficient evaluation of neuron populations receiving colored-noise current based on a refractory density method. Phys. Rev. E 77, 011910. (2008) 7. Chizhov A.V., Graham L.J.: Population model of hippocampal pyramidal neurons, linking a refractory density approach to conductance-based neurons. Phys. Rev. E 75, 011924 (2007) 8. Chizhov A.V., Smirnova E.Y., Graham L.J.: Mapping between V1 models of orientation selectivity: From a distributed multi-population conductance-based refractory density model to a firing-rate ring model. BMC Neuroscience. 10, 181 (2009) 9. Chizhov A.V.: A model for evoked activity of hippocampal neuronal population (in Rus-sian). Biofizika. 47(6), 1086-1094 (2002) 10. Chizhov A.V.: A two-compartment model for the dependence of a postsynaptic potential on a postsynaptic current, measured by the patch-clamp method (in Russian). Biofizika. 49(5), 877-880 (2004) 11. Chizhov A.V.: Conductance-based refractory density model of primary visual cortex. J. Comput. Neurosci. 36, 297–319 (2014) 12. Dayan P., Abbott L.F.: Theoretical Neuroscience: Computational And Mathematical Mod-eling of Neural Systems. The MIT Press (2001) 13. Graham L.J.: Interpretations of Data and Mechanisms for Hippocampal Pyramidal Cell Models. In: Ulinski, P.S., Jones, E.G., Peters, A. (eds.) Cerebral Cortex, vol.13, pp. 19-138. Plenum Press, New York (1999) 14. Hubel, D. H., Wiesel, T. N.: Receptive fields of single neurones in the cat’s striate cortex. J. Physiol. 148, 574-591 (1959) 15. Hubel, D. H., Wiesel, T. N.: Receptive Fields, Binocular Interaction and Functional Archi-tecture in the Cat’s Visual Cortex. J. Physiol. 160, 106-154 (1962) 16. Kopell N., Ermentrout G.B., Whittington M.A., Traub R.D.: Gamma rhythms and beta rhythms have different synchronization properties. Proc. Nat. Acad. Sci. USA. 97(4), 1867-1872 (2000) 17. La Cara GE, Ursino M.: Direction selectivity of simple cells in the primary visual cortex: comparison of two alternative mathematical models. II: Velocity tuning and response to moving bars. Comput. Biol. Med. 37(5), 598-610 (2007) 18. Lu H.D., Chen G., Tanigawa H., Roe A.W.: A motion direction map in macaque V2. Neu-ron. 68(5), 1002-1013 (2010) 19. Priebe N.J., Ferster D.: Direction selectivity of excitation and inhibition in simple cells of the cat primary visual cortex. Neuron. 45(1), 133-45 (2005) 20. Rodrigues S., Chizhov A.V., Marten F., Terry J.R.: Mappings between a macroscopic neu-ral-mass model and a reduced conductance-based model. Biol Cybern. 102(5), 361-71 (2010) 21. Saul A.B., Humphrey A.L.: Spatial and temporal response properties of lagged and nonlagged cells in cat lateral geniculate nucleus. J. Neurophysiol. 64(1), 206-224 (1990) 22. Saul A.B., Humphrey A.L.: Evidence of input from lagged cells in the lateral geniculate nucleus to simple cells in cortical area 17 of the cat. J. Neurophysiol. 68(4), 1190-1208 (1992) 23. Ursino M., La Cara G.E., Ritrovato M.: Direction selectivity of simple cells in the primary visual cortex: comparison of two alternative mathematical models. I: response to drifting gratings. Comput. Biol. Med. 37(3), 398-414 (2007) 24. Vigeland L.E., Contreras D., Palmer L.A.: Synaptic mechanisms of temporal diversity in the lateral geniculate nucleus of the thalamus. J. Neurosci. 33(5), 1887-1896 (2013) 25. Yakimova E. G., Chizhov A. V.: Experimental and Modeling Studies of Orientational Sen-sitivity of Neurons in the Lateral Geniculate Nucleus. Neuroscience and Behavioral Physi-ology. 45(4), 465-475 (2015)

REVIEW This paper proposes a reevaluation of the impact of NMDA conductance recruitment during visual activation of V1 cortical neurons in vivo and questions the validity of previous decomposition methods used by different authors to dissect out the timing and the influence of excitatory and inhibitory synaptic components in visual responses. On the basis of simple conductance-based simulations, they point out the fact that inhibitory contribution may have been under-estimated during phases of strong depolarization, due to the absence of the NMDA component in the decomposition. They conclude that this over-simplification in the decomposition puts in jeopardy the existence of push-pull organization in Simple cells. This report was written with Dr Cyril Monier. The methodological issue raised in this commentary paper is certainly a valid one, but the conclusions of the author remains largely exaggerated if not unsubstantiated. On the one hand, it is a fact that most conductance decompositions of synaptic responses do no account for the contribution of NMDA receptor activation, which may indeed participate to the expression of functional properties of sensory cortical neurons in vivo. On the other hand, the functional impact of this hypothetical recruitment is critically linked to the assumed ratio between NMDA and AMPA excitatory components, and unfortunately the author does not provide new experimental evidence that would justify his claim. In a way the present paper pushes already open doors and it would be unfair to assume that this issue was totally ignored by the pioneers of conductance decomposition of visual responses in vivo. This problem was partially tackled in the PhD thesis of Cyril Monier and the choice of voltage clamp measurements provided a way of minimizing the masked contribution of NMDA currents triggered in current clamp experiments. Cyril Monier, in his thesis, mentioned explicitly the fact that ”it is possible that in the presence of strong excitatory conductances, one underestimates the inhibitory conductance, by not taking into account voltage-dependent and NMDA conductances recruited by the global excitatory process”. We review below in succession a few factors that may weaken significantly the claims of the author: - Voltage clamp vs Current Clamp: It is important to distinguish in the published literature between measurements made in voltage clamp (VC) (Fregnac’s team; Borg-Graham et al, 1996; 1998; Monier et al, 2003; 2008) and current clamp (CC) (the later work of Ferster, Hirsch, Priebe’s teams among others: Anderson et al, 2000; Priebe et al, 2005). Not only do the membrane time constants differ between the two recording modes, but the local depolarization state of the recorded neuron is better controlled in VC mode than in CC mode (in spite of imperfect space clamp in vivo). Consequently in this latter mode (VC), the possible impact of NMDA recruitment can be ignored for holding potentials more negative than -40 mV (where only 10% of NMDA receptors are activated). No. The authors should provide distinct simulations for each of the recording modes (CC and VC). - Linearity (or not) of IV curves: Monier et al (2008) report that IV curves (for voltage values between -90 and -40 mV) are more linear when measured during the evoked state than during the ongoing state. For this range of potentials, NMDA activation will correspond to a weak negative conductance component with a linear impact. If inhibitory conductances are strong (which we reported in vivo), the effect can be neglected. The authors should provide distinct simulations for different level of inhibitory conductances. - Amplitude of NMDA conductances : The choice made by the author of the balance between AMPA and NMDA in the proposed simulations is very high (compared to Stern et al, 1992). This raises the question of the criticality of that ratio, and whether the impact of NMDA is a linear function of this ratio. This may constitute a side issue that could be studied parametrically through additional simulations. The main argument in favor of a strong NMDA component in cat cortex in vivo comes from an abstract of Kevin Fox (1990), where the iontophoretic effect of NMDA blocker (D-APV) is described as mostly developmental (before and at the peak of the critical period) and layer-dependent, being stronger in superficial layers. In granular and deep layers (IV, V, and VI), D-APV affected the visual response in young animals but only spontaneous activity in older animals. Fox study does not comfort the present proposal since the push-pull organization is the hallmark of Simple cells in layer 4. Note also that the fact that anti-correlated excitatory and inhibitory conductances can be observed in Simple Cells for other types of stimuli than gratings, such as dense noise (Baudot et al, 2013), whose kinetics are too fast for NMDA recruitment, remains an additional argument for the existence of push-pull organization in cat V1. Note finally, that the fact that the push-pull organization in Simple cells is observed only in the cat and not in the rodent visual coretx, suggests that this organization results from a structural anlage which is specific of the species under study rather being the trivial consequence of a methodological artifact (the consequence of which should be the same in both studies). - GABA-B involvement : Since the authors are looking for a more exhaustive description of the impact of each possible neurotransmitter, it could have been in the scope of this paper to try to separate also GABA-A and GABA-B influences. This may be particular relevant for the type of stimulation used for the simulations. A drifting grating triggers only a transient activation of the GABA-A, which because of different time constants disappears a few tens of ms after the stimulus onset. Iontophoretic experiments by Albus et al showed a clear temporal dissociation of the effect on visual responses. A more complete study may show that two components can be distinguished in inhibition, one dynamically related to the phase of the stimulus (in an anti-phase fashion in simple cells) and another one acting as a non-specific dc component. But without added experimentation, this discussion will remain purely academic. RESPONSE We thank the editor for the choice of such proper reviewer for our Commentary and reviewer for careful consideration and insightful comments. We agree that the style of the Commentary should not be as if the issue was not considered before, thus we have corrected the text in the light of this aspect. Regarding the main conclusion, it is still the same after our careful consideration of the critics of the reviewer (see Supplementary Material), thus in the Commentary we have kept the same Figure and added some remarks in the text that are supported by figures given in Supplementary Material. SUPPLEMENTARY MATERIAL Voltage clamp vs Current Clamp The voltage-clamp (VC) and current-clamp (CC) modes give similar estimations. We illustrate this point by the Suppl. figures [F1] and [F2]. Linearity (or not) of IV curves Regarding the dependence of the estimation errors on the level of inhibition, we show (compare Suppl. fig. [F3] with [F1]) that the absolute value of the error does not depend on gGABA. Nontrivially, for the range of potentials more negative than -40 mV, where NMDA activation corresponds to a weak negative conductance, the error of inhibitory conductance estimation is large, as seen in Suppl. fig. [F4]. Explanation of this effect follow from the voltage-dependence of NMDA-current, fNMDA(V) (V − VE), shown in Suppl. fig. [F45]. A linear approximation of the current between the hold voltage levels, for instance, Vh¹ = −70 mV and Vh² = −40 mV, corresponds to negative conductance and an effective reversal close to VI. That is why, NMDA-current contaminates more gI than gE. The method is perfect in respect of gI estimation for the case of Vh² = 0, however it gives gE close to gAMPA, ignoring contribution of gNMDA. Usage of the hold voltage equal to the reversal potential of glutamatergic currents can be recommended, however, as known, it might require additional tricks to improve space-clamp. Amplitude of NMDA conductances The underestimation of inhibitory conductance depends on the magnitude of NMDA conductance. The dependence is shown in Suppl. fig. [F5]. In order to justify the ratio of NMDA versus AMPA conductance, set to be 4 in all simulations except those for Suppl. fig. [F5], we have considered experimental data from Stern et al. (1992). The data is well approximated by the profiles of the conductance components shown in Suppl. fig. [F6], which is seen from comparison of reproduced EPSCs shown in Suppl. fig. [F7]B with the experimental figure 12A from Stern et al. (1992) and reproduced in Suppl. fig. [F7]A. The AMPA-component dominates at -100 mV, whereas the NMDA-component is the most pronounced at +20mV, which is observed in both experiment and model. The ratio of the amplitudes of the components is almost the same in the experiment and model. This ratio is equal to 4, as in Fig. 1 of the Commentary and Suppl. fig. [F6]. Further, with this ratio kept, we changed the shape of presynaptic stimulus from pulse to sinusoidal (Suppl. fig. [F8]). In spite of the evidence that in the case of pulse stimulation the maximum of AMPA-conductance is much bigger than that for NMDA, the sinusoidal long-lasting stimulus results in the conductance ratio close to 4. This fact justifies the simulations described in Commentary and in this Response. Similar ratio of NMDA-versus-AMPA conductances has been revealed in our recent studies of interictal discharges in the combined slices of hippocampus and entorhinal cortex of rats (Amakhin et al., 2016). In this study the estimations of gAMPA, gGABA and gNMDA were performed from multi-level recordings in VC-mode. Estimations of gE and gI with conventional method have shown not qualitatively but quantatively different results. Push-pull organization The push-pull organization (Troyer et al., 1998, Fregnac and Bathellier, 2015) may result in antiphase of gE and gI, or may be less pronounced, not resulting in the antiphase. Nevertheless, the qualitative behavior of the conductances is an important feature to reproduce in computational models of the primary vision. However, we think that the antiphase of gE and gI in cats has not been well proved experimentally yet, because of possible contamination of conductance estimations by NMDA effect. Our Commentary pays attention to this problem and calls to improve the method of conductance estimation. We know from Liu et al. (2010) that in rodents EPSC and IPSC recorded at -70 mV and 0 in response to moving gratings are in antiphase, corresponding to in-phase course of gE and gI. According to our study, this derivation is correct, as follows from Suppl. fig. [F4] for Vh² = 0. Unfortunately, we do not know similar recordings in cats. GABA-B involvement GABA-B conductance does not so strongly depend on voltage as NMDA. That is why, it should be approximated as voltage independent. Together with AMPA and GABA, there are three voltage independent components. Unfortunately, without any additional assumptions on different kinetics, it is impossible to evaluate more than 2 such components. This issue is discussed, in particular, in Chizhov et al. // Front.Cell.Neurosc. (2014). Effect of noise In addition, we have considered the influence of noise, important to compare the methods of conductance estimation in VC and CC modes. We show that the noise does not dramatically affect the estimations in VC-mode (Suppl. fig. [F9]), whereas it does in CC-mode (Suppl. fig. [F10]), if noting quite different amplitudes of noise in these simulations. Alternative method of conductance estimation Better method of conductance estimation with explicit consideration of NMDA-component require recordings at 3 levels of the hold voltage in VC-mode. Its quality is seen from Suppl. fig. [F9]. Unfortunately, the method shows much worse precision in the CC-mode (Suppl. fig. [F10]). In VC-mode, the means of estimated gAMPAEst, gGABAEst and gNMDAEst are equal to the exact gAMPA, gGABA and gNMDA. The variation of gGABAEst due to noise in the 3-level method is smaller than not only the variation of gI but also the error of gI relative to gGABA in the 2-level method, as seen from Suppl. fig. [F11]. Conclusion On a base of the above analysis we may recommend for gE and gI estimations to use 2-level method in VC-mode with the hold voltages close to the reversal potentials of GABAergic and glutamatergic currents. For estimation gAMPA, gGABA and gNMDA we recommend the 3-level method in VC mode with the hold voltages -70, -40 and 0 mV. METHODS Neuron model We consider leaky neuron equation model with the AMPA, GABA and NMDA-receptor mediated synaptic currents: C {dV \over dt}&=&-g_L(V -V_L)+g_{AMPA}(t) (V_E-V ) +g_{GABA}(t) (V_I-V )\\ &+&g_{NMDA}(t) f_{NMDA}(V) (V_E-V ) + I + \sigma_V g_L \eta(t)\nonumber where C is the membrane capacity; gL is the membrane time constant; VL, VE, VI are the reversal potentials of the leak, excitatory and inhibitory currents; I is the current through electrode; fNMDA(V)=1/(1 + Mg/3.57exp(−0.062 V)) is the voltage dependence function of NMDA-conductance with the magnesium concentration Mg; η(t) is the gaussian white noise with zero mean and unity dispersion. Estimation of gE and gI in VC-mode Estimations of gE(t) and gI(t) require two recorded currents I₁(t) and I₂(t) at two voltage levels Vh¹ and Vh². The result is g_I(t)&=&(-G_L(V_h^1-V_L)(V_h^2-V_E)+G_L(V_h^2-V_L)(V_h^1-V_E) \nonumber \\ &-&I_2(t)(V_h^1-V_E)+I_1(t)(V_h^2-V_E)) \\ &/&((V_h^1-V_I)(V_h^2-V_E) -(V_h^2-V_I)(V_h^1-V_E)), \nonumber g_E(t)&=&(-G_L(V_h^1-V_L)(V_h^2-V_I)+G_L(V_h^2-V_L)(V_h^1-V_I) \nonumber \\ &-&I_2(t)(V_h^1-V_I)+I_1(t)(V_h^2-V_I)) \\ &/&((V_h^1-V_E)(V_h^2-V_I) -(V_h^2-V_E)(V_h^1-V_I)) \nonumber Estimation of gE and gI in CC-mode In CC-mode, estimations of gE(t) and gI(t) require two recorded voltage traces V₁(t) and V₂(t) calculated for two input currents I¹ and I² according to eq. ([e1]) with known functions of gAMPA(t), gGABA(t) and gNMDA(t). Then, we assume that the voltages approximately satisfy to the equations with unknown gE and gI: C {dV_1 \over dt}&=&-g_L(V_1 -V_L)+g_E(t) (V_E-V_1 ) +g_I(t) (V_I-V_1 ) + I_1 \\ C {dV_2 \over dt}&=&-g_L(V_2 -V_L)+g_E(t) (V_E-V_2 ) +g_I(t) (V_I-V_2 ) + I_2 where C and gL are supposed to be known. The conductances gE and gI are obtained as a solution of the system of algebraic linear equations ([e131], [e132]). Estimation of gAMPA, gGABA and gNMDA in VC-mode Estimations of gAMPAEst, gGABAEst and gNMDAEst require three currents I₁(t), I₂(t) and I₃(t) calculated with eq.([e1]) at three voltage levels Vh¹, Vh² and Vh³. The conductances gAMPAEst, gGABAEst and gNMDAEst are obtained for each time step as a solution of the system of algebraic linear equations: g_L(V_h^1 -V_L)&-&g_{AMPA}^{Est}(t) (V_E-V_h^1 )-g_{GABA}^{Est}(t) (V_I-V_h^1 ) \nonumber \\ &&-g_{NMDA}^{Est}(t) f_{NMDA}(V_h^1) (V_E-V_h^1 ) = I_1(t), \\ g_L(V_h^2 -V_L)&-&g_{AMPA}^{Est}(t) (V_E-V_h^2 )-g_{GABA}^{Est}(t) (V_I-V_h^2 )\nonumber \\ &&-g_{NMDA}^{Est}(t) f_{NMDA}(V_h^2) (V_E-V_h^2 ) = I_2(t), \\ g_L(V_h^3 -V_L)&-&g_{AMPA}^{Est}(t) (V_E-V_h^3 )-g_{GABA}^{Est}(t) (V_I-V_h^3 )\nonumber \\ &&-g_{NMDA}^{Est}(t) f_{NMDA}(V_h^3) (V_E-V_h^3 ) = I_3(t). Estimation of gAMPA, gGABA and gNMDA in CC-mode In CC-mode, estimations of gAMPA, gGABA and gNMDA require three recorded voltage traces V₁(t), V₂(t) and V₃(t) calculated for three input currents I¹, I₂ and I³ according to eq.([e1]) with known functions of gAMPA, gGABA and gNMDA. Then, gAMPAEst, gGABAEst and gNMDAEst are obtained for each time step as a solution of the system of algebraic linear equations C {dV_1 \over dt}&+&g_L(V_1 -V_L)-g_{AMPA}^{Est}(t) (V_E-V_1 )-g_{GABA}^{Est}(t) (V_I-V_1 ) \nonumber \\ &&-g_{NMDA}^{Est}(t) f_{NMDA}(V_1) (V_E-V_1 ) = I_1(t), \\ C {dV_2 \over dt}&+&g_L(V_2 -V_L)-g_{AMPA}^{Est}(t) (V_E-V_2 )-g_{GABA}^{Est}(t) (V_I-V_2 )\nonumber \\ &&-g_{NMDA}^{Est}(t) f_{NMDA}(V_2) (V_E-V_2 ) = I_2(t), \\ C {dV_3 \over dt}&+&g_L(V_3 -V_L)-g_{AMPA}^{Est}(t) (V_E-V_3 )-g_{GABA}^{Est}(t) (V_I-V_3 )\nonumber \\ &&-g_{NMDA}^{Est}(t) f_{NMDA}(V_h^3) (V_E-V_h^3 ) = I_3(t), \\ where C and gL are supposed to be known. Synaptic kinetics The synaptic conductances were described following Chizhov (2013) , i.e. as follows g_{AMPA}(t)&=&_{AMPA} ~m_{AMPA}(t), \\ g_{NMDA}(t)&=&_{NMDA} ~f_{NMDA}(V(t)) ~m_{NMDA}(t), where mAMPA(t) and mNMDA(t) are the non-dimensional synaptic conductances, each being approximated by the second order ordinary differential equation: \tau_1 \tau_2 {dt^2} + (\tau_1+\tau_2) {dt} + m = \tau \varphi(t), \tau&=&{(\tau_1-\tau_2)}/ \bigl( ({\tau_2}/{\tau_1})^{\tau_2/(\tau_1-\tau_2)} - ({\tau_2}/{\tau_1})^{\tau_1/(\tau_1-\tau_2)} \bigr), \\ &&~~~~\tau_1\neq \tau_2, {{\tau_1}{e}}, ~~. \nonumber Here φ(t) is the presynaptic firing rate; $_{AMPA}$ and $_{NMDA}$ are the maximum conductances, τ₁ and τ₂ are the rise and decay time constants. The time scale τ in the form of eq.([e254]) provides independence of the maximum of g(t) on τ₁ and τ₂ in the case of stimulation by a short pulse of φ(t). References Amakhin, D.V., Ergina, J.L., Chizhov, A.V., and Zaitsev, A.V. (2016) Synaptic conductances during interictal discharges in pyramidal neurons of rat entorhinal cortex. Frontiers in Cellular Neuroscience 10:233, DOI: 10.3389/fncel.2016.00233. Chizhov, A.V. (2014) Conductance-based refractory density model of primary visual cortex, J. of Comp. Neuroscience, 36(2): 297-319. Chizhov, A.V., Malinina, E., Druzin, M., Graham, L.J., and Johansson, S. (2014) Firing clamp: a novel method for single-trial estimation of excitatory and inhibitory synaptic neuronal conductances. Front. Cell. Neuroscience 8: 86. Fregnac, Y., and Bathellier, B. (2015) Cortical Correlates of Low-Level Perception: From Neural Circuits to Percepts. Neuron 88, 110-126. Liu, B.H., Li, P., Sun, Y.J., Li, Y.T., Zhang, L.I., and Tao, H.W. (2010) Intervening inhibition underlies simple-cell receptive field structure in visual cortex. Nat. Neurosci. 13, 89–96. Troyer, T.W., Krukowski, A.E., Priebe, N.J., and Miller, K.D. (1998) Contrast invariant orientation tuning in cat visual cortex: thalamocortical input tuning and correlation-based intracortical connectivity. J. Neurosci. 18, 5908–5927.

A number of papers reported data on estimation of excitatory and inhibitory synaptic conductances as input signals into principal neurons of the primary visual cortex (Anderson et al., 2000; Monier et al., 2003, 2008; Priebe and Ferster, 2005). Authors emphasize that they observed the counterphase of excitation and inhibition in simple cells and associate the observations with the push-pull mechanism (Troyer et al. 1998; Fregnac and Bathellier, 2015). However, the procedure of the estimation might occur incorrect and derived conclusions might be controversial, because of using too strong assumption that synaptic conductances are voltage-independent. The role of this assumption has been studied by Monier et al. (2008). It is wrong in presence of NMDA-receptor mediated component, which makes current significantly nonlinear due to magnesium block of the channels. In order to reduce this effect, the recordings are done in the voltage range where the NMDA voltage-current relationship is close to linear. Below we reconsider this issue and demonstrate in simulations how the estimation errors may lead to a doubtful conclusion that excitatory and inhibitory conductances are modulating in counterphase during visual stimulation by moving gratings. Let’s consider leaky neuron equation model with the AMPA-, GABA- and NMDA-receptor mediated synaptic currents: C {dV \over dt}=-g_L(V -V_L)+g_{AMPA} (V_E-V )+g_{GABA} (V_I-V )+g_{NMDA} ~f_{NMDA}(V ,Mg) ~(V_E-V ) + I where C is the membrane capacity; gL is the membrane time constant; VL, VE, VI are the reversal potentials of the leak, excitatory and inhibitory currents, correspondingly; I is the current through electrode; fNMDA = 1/(1 + Mg/3.57 exp(−0.062 V)) is the voltage dependence function of NMDA-conductance with the magnesium concentration Mg. The experimental method of conductance estimations was proposed for both voltage- (Borg-Graham et al., 1996) and current-clamp modes (Anderson et al., 2000). These versions give almost equivalent estimates if the input signals are slow enough in comparison with the membrane time constant (Monier et al., 2008). The synaptic currents are supposed to be linearly dependent on voltage, distinguishing excitatory (gE) and inhibitory (gI) conductances. Thus the equation of current conservation in the voltage-clamp mode is the following: -g_E (V-V_E)-g_I (V-V_I)-G_L (V-V_L)+I = 0 Estimations of gE(t) and gI(t) require two recorded currents I₁(t) and I₂(t) at two voltage levels V₁ and V₂. The result is \nonumber g_I(t)=(-G_L(V_1-V_L)(V_2-V_E)+G_L(V_2-V_L)(V_1-V_E)-I_2(t)(V_1-V_E)+I_1(t)(V_2-V_E)) /((V_1-V_I)(V_2-V_E) -(V_2-V_I)(V_1-V_E)), \nonumber g_E(t)=(-G_L(V_1-V_L)(V_2-V_I)+G_L(V_2-V_L)(V_1-V_I)-I_2(t)(V_1-V_I)+I_1(t)(V_2-V_I)) /((V_1-V_E)(V_2-V_I) -(V_2-V_E)(V_1-V_I)) An example of a comparison of the true AMPA, GABA and NMDA conductances with the estimations is presented in Fig. 1. In simulations we considered sinusoidal gAMPA and gNMDA, mimicking visual stimulation by moving gratings, and constant gGABA. The “measured” currents I₁(t) and I₂(t) were calculated from eq.([e1]) by setting V(t)=V₁ and then V(t)=V₂, correspondingly. As seen, the obtained gE(t) is close to gAMPA(t), whereas gI(t) is very far from gGABA(t). gI(t) oscillates in counterphase with gE(t), as in estimates made in vivo with sinusoidal gratings stimuli (Anderson et al., 2000; Priebe and Ferster, 2005; Monier et al., 2008), whereas the true inhibitory conductance was set to be constant. Monier et al. wrote that they clamped neurons only at membrane potential more hyperpolarized than -40 mV, in order to suppress the NMDA current, but the theoretical estimations show that the effect of NMDA is still strong in that case. Also, the total synaptic conductance gE(t)+gI(t) is underestimated. This contradiction points to the necessity to take NMDA-component into consideration. Estimations of gAMPA, gGABA and gNMDA require measurements of currents at three voltage levels and solving linear algebraic equation composed of the steady-state eq.(1) applied at those three voltages. In our case of absent noise and parameters from Fig. 1, such estimations almost precisely reconstruct the conductances (data not shown). Similarly, in our recent study of interictal discharges (Amakhin et al., 2016) the estimations of gAMPA, gGABA and gNMDA were performed from multi-level recordings in VC-mode. Estimations of gE and gI with conventional method have shown not qualitatively but quantatively different results. Alternative to these multi-level estimations are the conventional two-level estimations of gE(t) and gI(t) with the hold voltages equal to VI and VE, which give gI(t) equal to gGABA(t) and gE(t) close to gAMPA, thus ignoring gNMDA (simulations not shown). The voltage-dependence of NMDA-conductance makes an effect of amplification which exposes as a negative conductance (Monier et al., 2008; Smirnova et al., 2015). The effect might also be necessary to take into account in the techniques of single-trial conductance estimation (Bedard et al., 2012; Chizhov et al., 2014). Thus, we conclude that not accounting of NMDA-component in synaptic conductance estimation compromises evidences of excitation and inhibition counterphase modulation during grating visual stimulation in the striate cortex. ACKNOWLEDGMENTS: The reported study was supported by the RFBR research project 15-04-06234. REFERENCE Amakhin, D.V., Ergina, J.L., Chizhov, A.V., and Zaitsev, A.V. (2016). Synaptic conductances during interictal discharges in pyramidal neurons of rat entorhinal cortex. _Frontiers in Cellular Neuroscience_. 10:233, DOI: 10.3389/fncel.2016.00233. Anderson, J.S., Carandini, M., and Ferster, D. (2000). Orientation tuning of input conductance, excitation, and inhibition in cat primary visual cortex. _J. Neurophysiol._ 84(2), 909-26. Baudot, P., Levy, M., Marre, O., Monier, C., Pananceau, M., and Fregnac. Y. (2013). Animation of natural scene by virtual eye-movements evokes high precision and low noise in V1 neurons. _Front. Neural Circuits_. 7: 206, 1-29. Bedard, C., Behuret, S., Deleuze, C., Bal, T., and Destexhe, A. (2012). Oversampling method to extract excitatory and inhibitory conductances from single-trial membrane potential recordings. _J. Neurosci. Methods_. 210, 3–14. Borg-Graham, L.J., Monier, C., and Fregnac, Y. (1996). Voltage-clamp measurement of visually-evoked conductances with whole-cell patch recordings in primary visual cortex. _J. Physiology_. 90, 185–188. Chizhov, A., Malinina E., Druzin, M., Graham, L.J., and Johansson, S. (2014). Firing clamp: A novel method for single-trial estimation of excitatory and inhibitory synaptic neuronal conductances. _Front. in Cell. Neuroscience_. 8, 86-94. Fregnac, Y., and Bathellier, B. (2015). Cortical Correlates of Low-Level Perception: From Neural Circuits to Percepts. _Neuron_. 88, 110-126. Monier, C., Chavane, F., Baudot, P., Graham, L.J., and Fregnac, Y. (2003) Orientation and direction selectivity of synaptic inputs in visual cortical neurons: a diversity of combinations produces spike tuning. _Neuron._. 37(4), 663-80. Monier, C., Fournier, J., and Fregnac, Y. (2008). In vitro and in vivo measures of evoked excitatory and inhibitory conductance dynamics in sensory cortices. _J. Neuroscience Methods_. 169(2), 323-365. Priebe, N.J., and Ferster, D. (2005). Direction selectivity of excitation and inhibition in simple cells of the cat primary visual cortex. _Neuron_. 45(1), 133-45. Smirnova, E.Y., Zaitsev, A.V., Kim, K.Kh., and Chizhov, A.V. (2015). The domain of neuronal firing on a plane of input current and conductance. _J. Comput. Neurosci._ 39(2), 217-233. Troyer, T.W., Krukowski, A.E., Priebe, N.J., and Miller, K.D. (1998) Contrast invariant orientation tuning in cat visual cortex: thalamocortical input tuning and correlation-based intracortical connectivity. _J. Neurosci._ 18, 5908–5927.