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.

#1#2 #1#2 #1#2 #1#2 #1 + #1#2 #1#2#1 / #2 #1#2 #1#1 #1#2 #1#2 #1#1 #1#1 #1 #1 #1 #1 Conductance-based refractory density (CBRD) approach is an efficient approach that describes firing activity of a statistical ensemble of uncoupled Hodgkin-Huxley-like neurons each receiving individual gaussian noise and common time-varying deterministic input. However, the approach requires comparison with experiments and extension to the case of more realistic synaptic connectivity. Here we verify a CBRD model through comparison with experimental data and then generalize the model for the case of lognormal (LN) distribution of input weights within a population. We show that the model with equal weights efficiently reproduces post-spike time histogram and membrane voltage of experimental multiple-trial response of single neurons to step-wise current injection and reveals much more rapid reaction of firing-rate than voltage. Slow adaptive potassium channels strongly affect the shape of the responses. Next, we derive a computationally efficient CBRD model for a population with the LN input weight distribution and compare it to the original model with equal input weights. The analysis shows that the LN distribution does: (i) provide faster response; (ii) eliminate oscillations; (iii) lead to higher sensitivity to weak stimuli; (iv) increase the coefficient of variation of interspike intervals. In addition, a simplified, firing-rate type model is tested, showing better precision in the case of LN distribution of weights. More generally, the CBRD approach can be recommended for complex, biophysically detailed simulations of interacting neuronal populations, and the modified firing-rate type model for computationally reduced simulations. INTRODUCTION Known neuronal population models belong either to the type of very simplified models such as firing-rate (FR) or neural mass models or to those that are based on probability-density approach (PDA). FR-models are usually expressed by algebraic or ordinary differential equations and easy to analyze. However, they are not fully adequate for transient process simulations because of their basic assumption that neurons are always desynchronized. A probability density approach avoids this assumption, however this approach is commonly applied to simplified, one-variable neurons like linear and nonlinear integrate-and-fire, or spike-response model or so , , , , . This limitation has been overcome by a conductance-based refractory density approach (CBRD) that is applicable to regular-spiking, adaptive and fast-spiking neurons described in terms of Hodgkin-Huxley-like approximations . The CBRD approach has been proposed for an infinite set of such uncoupled neurons receiving a common input and an individual gaussian white noise-current. The model was later extended to the case of color noise . Such model is found to be quite efficient for simulation of coupled populations . The approach is worth being further developed and used in analysis of experimental data. In the present paper, the CBRD model is applied to experimental data revealing rapid population response to a weak input. Matching the model to experimental dataset validates the model, which then reveals the influence of neuronal parameters on the response. Further, the model is generalized for the case of lognormal distribution of the weights of the input received by different neurons of a population, which has many physiological implications . In addition, a simplified, firing-rate type model is tested to simulate a population with lognormal weights. METHODS For the purpose of clarifying the ideas of the CBRD-approach and its extension to the lognormal input weight distribution, we present in this section the models for simplified, integrate-and-fire neurons. The details of the model for adaptive Hogkin-Huxley like neurons can be found in or . _Leaky integrate-and-fire neuron with noise_ The LIF neuron is given by the equation C~{dV\over dt}=-(g_L+s(t)) (V-V_{rest}) + I(t) + \sigma ~\xi(t), where ξ(t) is a gaussian white noise characterized by its mean value, <ξ(t)> = 0, and auto-correlation <ξ(t)ξ(t′) > =C/gL δ(t − t′); σ is the noise amplitude. The neuron fires when the potential V crosses the threshold $\VT$. Immediately after the spike the potential V is reset to Vreset. The LIF neuron is characterized by the capacitance C and the leak conductance gL. The input is determined by two signals, the synaptic current I(t) that is measured at the voltage level equal to Vrest, and the total synaptic conductance s(t). The effective membrane time constant is τm = C/(gL + s(t)). We will refer to as a population an infinite number of eq.([e221])-based LIF neurons receiving a deterministic 2-d input (I(t),s(t)) and individual for each neuron noise. The population firing rate is defined as a sum of all spikes, nact, from neurons of the population over a short time window Δt, divided by the number of neurons, N. After taking the limits of N → ∞ and Δt → 0, the firing rate ν is obtained as \nu(t)= {\Delta t} ~(t;~t+\Delta t)}{N}. The steady-state firing rate is given by the solution from \nu^{SS}&=&\Biggl({g_L+s(t)} ~ -U)/\sigma_V \sqrt(2)}^{(\VT-U)/\sigma_V \sqrt(2)} \exp(u^2) (1+\erf(u)) ~du \Biggr)^{-1}, \\ U&=&V_{rest}+I/(g_L+S),\\ \sigma_V&=&{(g_L+s)} \\ where σV is the voltage dispersion; here U is the depolarization potential in the steady-state. _Conductance-based refractory density (CBRD) approach for a population of integrate-and-fire neurons_ As justified in , , the firing rate for a population of LIF neurons is well approximated by the system of equations for the refractory density $\rho(t,\ts)$ and the averaged across noise realizations membrane potential $U(t,\ts)$, where $\ts$ is the time elapsed since the last spike. Thereby, the CBRD method distinguishes neurons only according to their $\ts$-variable. In order words, the state of each neuron is parameterized by this phase variable. For the particular case of a LIF-neuron, its only state variable is V, whereas for more complex neuron models the other state variables, the gating variables of ionic channels, are also parameterized by $\ts$, which is the reduction to an only one phase variable description that rather precisely saves the information about neuronal states because of the same history of input for all neurons . Returning to LIF-neurons, the equations for $\rho(t,\ts)$ and $U(t,\ts)$ are as follows t +{\ts}&=&-\rho~ H(U), \\ C\Biggl(\derc Ut + \derc U{\ts} \Biggr) &=& -(g_L+s(t)) (U-V_{rest}) + I(t), where H is a hazard function which is defined below. The boundary conditions are \nu(t) \equiv \rho(t,0)=\int ^{\infty} \rho H d\ts and U(t, 0)=Vreset, where ν(t) is the population firing rate. When calculating the dynamics of a neuronal population, the integration of eq.([e242]) determines the distribution of not-noisy voltage U across $\ts$. The effect of crossing the threshold and the diffusion due to noise are evaluated by the H-function and affect the equation for ρ. The result of the integration of eq.([e241]) is the distribution of ρ across $\ts$ and the firing rate ν calculated from ([e2425]). The hazard function H is defined as the probability for a single neuron to generate a spike, if known actual neuron state variables. The hazard function H has been approximated in for the case of white noise and in for the case of color noise as a function of U(t) and s(t), and parameters σ, $\VT$ and the ratio of membrane to noise time constants k = τm/τNoise: H&(U)&= A+B,\\ A&=& %A_{0}&=& {\tau_m} e^{0.0061-1.12~T-0.257~T^2-0.072~T^3-0.0117~T^4} %\nonumber \\&\cdot& \biggl(1-(1+k)^{-0.71+0.0825(T+3)}\biggr),\nonumber \\ B&=&- ~ \biggl\lbrack {dT \over dt} \biggr \tilde F(T), ~~~\tilde F(T)= ~{\exp(-T^2) \over 1+\erf(T)},~~~ %\nonumber \\ T={~\sigma_V}, \nonumber where T is the membrane potential relative to the threshold, scaled by noise amplitude; A is the hazard for a neuron to cross the threshold because of noise, derived analytically in and approximated by exponential and polynomial for convenience; B is the hazard for a neuron to fire because of depolarization due to deterministic drive, i.e. the hazard due to drift in the voltage phase space. Note that the H-function is independent of the basic neuron model and does not contain any free parameters or functions for fitting to any particular case. In this aspect, the H-function can be useful not only as a component of a population model, but as well for analysis tasks such as evaluation of neuronal susceptibility , effects of noise etc. _Firing-rate model for integrate-and-fire neurons_ Conventional firing-rate type models often fail to describe well a population activity in transient states. A modified firing-rate model has been proposed in , which improves the transient solutions to step-like or complex-shape inputs . Below such model is given for the LIF-neuron population. The subthreshold voltage U(t) is calculated as C~{dU\over dt}=-(g_L+s(t)) (U-V_{rest}) + I(t), and the firing rate is calculated by the following formula \nu&=&\bigl[\nu^{SS} + \nu^{US} \bigr]_{+}, \\ \nu^{US}&=&{1 \over \sigma_V}{dU \over dt} \exp\Biggl( -{(\VT-U)^2 \over 2 \sigma_V^2} \Biggr), where the function [x]+ is defined as x for x > 0 and 0 otherwise. The term νSS is the steady-state solution ([e223]). The essence of the model modification consists in the term νUS, which evaluates the firing rate in response to fast excitation in assumption of a frozen gaussian distribution of neuronal membrane potentials crossing the threshold. RESULTS Validation of CBRD-model: comparison with analytical and Monte-Carlo solutions To validate the CBRD-model based on eqs.([e221]-[e223]) we compare it with Monte-Carlo (MC) simulation on a test problem of step-wise current stimulation. In MC simulations eq.([e221]) has been integrated 4000 times with different realizations of noise. The firing rate has been calculated from eq.([e2013]). As seen in Fig.[F-1]A, the CBRD and MC solutions are close; they converge with the increase of the number of neurons in MS simulations. For LIF neuron with white noise in steady-state, the analytical solution for the firing rate is known, it is given by eq.([e223]-[e225]). The CBRD-model matches to this solution with a high precision, as demonstrated by Fig.[F-1]B. Important, that the both solutions converge not only for variable input current but also for the input conductance. Comparison of CBRD-model with experimental data According to above introduced definition of a population, the neurons are independent from each other, thus, such population is equivalent to a number of noise realizations for a single neuron. Hence, statistics of the behavior of a neuron receiving both deterministic and stochastic inputs is governed by the same model as a single population. That is why, we compare the CBRD-model with experimental multiple registrations in a single neuron . A weak step-wise current was injected into the neuron. The statistics of spiking response is characterized by the post-stimulus spike-time histogram (PSTH) which corresponds to the firing rate for a population. The model has been taken from . In this model a basic single neuron had a few potassium ionic channels. The adaptation was taken into account in the forms of the M-current and the current of afterhyperpolarization (AHP-current) which effectively approximates the effects of potassium calcium-dependent currents. The spike trains in response to step-wise current injection reveals the similarity of the modeled and registered neuronal responses (compare Figure 1A,D in ) with Fig. [Fig1ABCD]A,C) and the effects of the adaptation. As to population response, the experimental PSTH (Figure 1E in ) is similar to the modeled firing rate (Fig. [Fig1ABCD]D). The response in the model is as fast as in the experiment, the timing of peaks, the amplitudes of the firing rate and the character of relaxation to the steady-state are similar. Also we find out similar rapid reaction to the change of the amplitude of noise (Fig. [Fig3B]). The most affective factors determining the shape of firing response are the level of background activity before the stimulus and the slow potassium channels (Fig. [Fig_SpontActivity_and_WhiteColorNoise]). As seen from comparison of Fig. [Fig_SpontActivity_and_WhiteColorNoise]C and D with Fig. [Fig_SpontActivity_and_WhiteColorNoise]B, the presence of spiking activity before the stimulus is critical for the rapidness of the response. The bigger the firing rate before stimulation the faster the response. This fact is consistent with the findings of the experimental study . Blockage of M- and AHP-currents crucially changes the later phases of the response, especially its steady-state, as seen from comparison of Fig. [Fig_SpontActivity_and_WhiteColorNoise]H with Fig. [Fig_SpontActivity_and_WhiteColorNoise]F. Regarding the time correlations of the noise, their effect is moderate (compare Fig. [Fig_SpontActivity_and_WhiteColorNoise]G with Fig. [Fig_SpontActivity_and_WhiteColorNoise]F). Derivation of CBRD-model for a population receiving lognormally distributed input Our goal is to generalize the CBRD-approach for a population of neurons receiving equal input I(t) and noise, that based on eqs.([e241]-[e243]), to the case of a population consisting of neurons receiving lognormally distributed deterministic inputs Ix(t). The distribution of the current Ix(t) scaled by its mean across the distribution I(t), x = Ix(t)/I(t), is as follows \psi(x)=/{(2~^2)}}\bigr)}{~ ~x} Let us consider neurons parameterized with x, i.e. those neurons that receive the deterministic input current Ix(t) as well as the input conductance s(t) and the noise σ ξ(t) introduced into eq.([e221]). The voltage averaged across realizations of the noise, to be denoted as Ux, is governed by eq.([e242]) written for Ux instead of U and Ix(t) instead of I(t). Such equation is linear in respect to Ix(t), thus the deflection of voltage $U_x(t,\ts)$ from the unperturbed state $U_{free}(\ts)$ (defined for zero input current I(t)=0) is also lognormally distributed in accordance to eq.([e3.1]), with its mean $U(t,\ts)-U_{free}(\ts)$ corresponding to the mean input current I(t). That is why, the voltage $U_x(t,\ts)$ can be found as U_x(t,\ts)=(U(t,\ts)-U_{free}(\ts))~x+U_{free}(\ts). It allows to reduce the consideration by avoiding direct solving of equations in partial derivatives for $U_x(t,\ts)$ on a semiinfinite space for x. Instead, the statistics for $U_x(t,\ts)$ is described by the single equation ([e242]), the distribution ([e3.1]) for x and the expression ([e3.3]). Now we derive equation for the neuronal density. The probability for a neuron with the voltage $U_x(t,\ts)$ to generate a spike is estimated by the hazard function H(Ux, dUx/dt). The distribution of neurons parameterized by x in the phase space $\ts$ is denoted as $\rho_x(t,\ts)$. Note that the distribution of $\rho_x(t,\ts)$ across x is not necessarily lognormal. Calculation of $\rho_x(t,\ts)$ require solving of a continuum of eqs.([e241]) for ρx instead of ρ with H(Ux, dUx/dt). The output firing rate is defined as \nu(t)=^{\infty} \rho_x(t,0) ~\psi(x)~dx The firing rate in the case of lognormal weight distribution of the deterministic input is calculated from eqs.([e241]-[e243], [e3.1]-[e3.2]). Because the extension of the numerical method of the calculations from that for the conventional method from is straightforward, the precision of the CBRD-model in this case is determined only by the parameters of integration. The values of the numerical parameters were chosen such that a residual error was less than 5%. The approximation of the distribution ([e3.1]) required 10-15 points of discretization; hence, the number of numerically integrated 1-d transport equations for ρx was equal to that number of points plus one (for U). Comparison of CBRD models with fixed and lognormally distributed input weights To reveal the difference between a population with equal weights and gaussian noise (P) and a population with lognormal distribution of weights (LNP) we considered again the task of step-wise current stimulation (Fig.[Fig_LNP_P_step]). The response of LNP is more rapid than both the voltage and the firing rate response of P, because of contribution of more sensitive neurons with strong weights from the tail of lognormal distribution. The LNP solution does not oscillate as the response of P, because of overlapping of contributions of the whole lognormal distribution. The steady-state firing rate is greater for LNP, because of nonlinear dependence of the firing rate on the input current for a single neuron (Fig. [F-1]B) and thus bigger contribution of the tail of lognormal distribution. The rapidness of the LNP response increases with the amplitude of stimulation (Fig. [Fig_3_ampl_step_and_sin]A). The sensitivity to weak stimulus is greater for LNP, again because of the contribution of the tail of lognormal distribution. In the steady-state, the lognormal distribution affects not only the firing rate but the interspike interval distribution also. According to the construction of the refractory density model, the distribution of interspike intervals is given by the source term in the eq.([e241]), i.e. ρ H, where ρ is the cumulative density $\rho(\ts)=^{\infty} \rho_x(\ts) ~\psi(x)~dx$. The coefficient of variation (CV) of interspike intervals for the whole population follows from the refractory density $\rho_x(\ts)$ and the hazard function H: CV=^{\infty} ^{\infty} \ts^2~\rho_x(\ts) H(U_x(\ts)) ~d\ts ~\psi(x)~dx %\rho_x(\ts) H(U_x(\ts)) ~\psi(x)~d\ts \\ / (E[\ts])^2 %\biggl( % ^{\infty} ^{\infty} \ts ~\rho_x(\ts) H(U_x(\ts)) ~\psi(x)~d\ts %\biggr)^2 -1 }, where \nonumber E[Z]= ^{\infty} ^{\infty} Z ~\rho_x(\ts) ~H(U_x(\ts)) ~\psi(x)~d\ts \big /~\nu % ^{\infty} ^{\infty} ~\rho_x(\ts) H(U_x(\ts)) ~\psi(x)~d\ts The dependence of the firing rate and CV on the input current I and the lognormal distribution dispersion σLN is plotted in Fig. [Fig_Iind_nu_CV]A. The CV strongly increases with σLN, which corresponds to the widening of the ISI distribution (ρH) shown in Fig. [Fig_Iind_nu_CV]B. Comparison of CBRD and firing-rate models We test the quality of solution obtained by the firing-rate model with lognormal distribution of weights (FR-LNP), based on eqs. ([e231]-[e234],[e223] and [e3.1]). As seen from Fig. [Fig_3_ampl_step_and_sin]A, the initial part of the response as well as the plateau of the steady-state are captured by the model. Good approximation of the rapid response at the initial stage is due to the term νUS in the eq.([e232]). The cumulative across the whole lognormal distribution action of this term provides the sharp profile. The most problematic is the phase of relaxation after the end of stimulation, the FR-LNP model shows too delayed relaxation. Better precision of the FR-LNP model is observed in the case of oscillatory input, as seen from Fig. [Fig_3_ampl_step_and_sin]B. DISCUSSION AND CONCLUSIONS The probability density approach and, in particular, the refractory density approach provides good approximation of firing rate calculation for a statistical ensemble of neurons in transient regimes of activity , , , , . Among these models, only the refractory density approach has been designed for complex, Hodgkin-Huxley like neurons , which was referred to as CBRD-model. Also important to note that the hazard-function derived for this model is more precise, does not contain any free forms or parameters for fitting and universal for different basic neuron models and input signal behavior, because it was strictly derived from the first-time passage problem for linearized voltage fluctuations about the time-variable mean near the threshold. These benefits of the model justifies the necessity of its further validation and development. In this paper some additional validations of the approach have been presented, namely its comparison to the alternative methods of firing rate calculation and experimental data. Further, the approach have been generalized to the case of lognormal distribution of “effective” synaptic weights, the weights of the input into neurons of a population. The CBRD model reveals rapidness and sensitivity of a population response to input changes. These properties are consistent with those found in exponential integrate-and-fire neurons and are physiological important for brain networks , . As pointed in , theoretical analysis identified three basic requirements for fast processing: use of neuronal populations for encoding, background activity, and fast onset dynamics of action potentials in neurons. All these requirements were satisfied in our simulations that have shown fast responses. In particular, the condition of the rapid initiation of action potentials , was fulfilled, because the CBRD-model uses threshold criterion for spike initiation, thus the initiation is instantaneous by the construction . The effect of background activity has been studied and found to be consistent with the experimental observations . Higher sensitivity is revealed in the case of lognormal distribution, which is explained by contribution of neurons from the tail of the distribution that are sparse but powerful. However this effect of the distribution is less pronounced than that for the interspike interval distribution. The coefficient of variation of intervals between spikes generated by the whole population as a measure of the interspike interval distribution strongly depends on the dispersion that characterizes the lognormal input weight distribution. Nevertheless, this coefficient does not relate directly to the dispersion of interspike intervals for a single neuron, which mechanism was studied in , etc. What is also important for simulations of interacting neuronal populations is mathematical complexity of a model and efficiency of its numerical solution. Regarding the CBRD approach, even in the case of realistic neurons and lognormal distribution its mathematical complexity is reduced to tens of one-dimensional transport equations, which are quite solvable numerically. Hence, computational efficiency of the method is comparable to that for population models with one-dimensional neurons and complex synaptic weight distributions, like in . Note, however, the difference between the proposed and the mentioned approach from . The synaptic weights were uncorrelated in time in that approach, whereas the proposed approach considers fixed distribution of weights for the deterministic input and temporally uncorrelated or colored individual noise. These assumptions seem to be more realistic. The above consideration allows one to conclude that the CBRD approach can be recommended as useful for cortical tissue simulations, because it provides significant gain in precision in comparison to conventional mean-field models, computational effectiveness in comparison with direct network simulations and the detail consideration of biophysical mechanisms in comparison with Fokker-Planck-based models; it also admits mathematical analysis of the system of equations. Regarding more simplified, firing-rate type population models, generally they fail to simulate correctly the transient regimes, however some modifications help to improve their quality , , . One of such models that was described in Methods was generalized for the case of lognormal distribution of input weights. Its comparison with the CBRD-model has revealed that the modified FR model provides more or less satisfactory approximation in the case of lognormal distribution. Thus, this model can also be recommended for simplified simulations. Acknowledgement This work was supported by the Russian Science Foundation (project 16-15-10201)

ABSTRACT Minimal level of detail needed to model extracellular potentials originating from activity of interacting neuronal populations is an important but still open question. Whereas in regimes with moderate ionic concentration changes the neuronal contribution into the local field potential in a layered cortical or hippocampal tissue can be evaluated as the difference between the somatic and dendritic membrane potentials of pyramidal neurons, hyper-synchronous states of activity reguires taking into account a contribution of glia. We propose an effective consideration of the glial contribution that is reduced to a low-pass filtered neuronal firing rate. Our formalism is applicable to the study of cortical activity using two-compartment neuronal population models. INTRODUCTION In experimental neuroscience, most electric signals are picked up from the electrolytic solution that constitutes extracellular (interneuronal) space, yet the origin of these potentials is incompletely understood. Our study touches the related open questions: - What level of detail is needed to model extracellular electrical phenomena, and which formalism can be employed? - How valid is the cable equation, or need slow concentration changes of the major ionic components involved be modeled explicitly? - How can LFP calculations be incorporated in simulation software? - Can single-compartmental model neurons generate local field potentials at all? Interpretation of extracellular potentials is unclear citeBuzhaki et al. 2012, citeBedard et al. 2006, citeBedard and Destexhe 2011, citeBedard and Destexhe 2013, citeBedard and Destexhe 2014 even in the case of synchronous neuronal epileptiform events the relationship between simultaneously registered extracellular and intracellular signals is non-trivial citeAlvarado-Rojas et al. 2015. Assuming an homogeneous and isotropic extracellular space (but see citeRanck and BeMent 1965), an extracellular potential is commonly described by the Poisson equation with macroscopically averaged transmembrane currents as distributed sources citeNicholson 1973, citeNicholson and Freeman 1975 or by the Laplace equation with the transmembrane currents distributed along all membrane surfaces citeLogothetis 2007, citeLinden et al. 2011, citeSchomburg et al. 2012, citeGold et al. 2006, citeEinevoll et al. 2013, citeIbarz et al. 2010. If neuronal activity is modeled by multiple compartment neurons, then calculation of the local field potential (LFP) with the Poisson equation is straightforward. However, in most cases biophysically detailed mean-field models also deals with somatic signals of neuronal activity exclusively, and the question of how to estimate the LFP from these somatic signals is unclear citeChizhov 2014, citeBaladron et al. 2012. A simple formula linking intracellular parameters of neuronal activity to extracellular field potential has been derived in citeChizhov et al. 2015, neglecting ionic dynamics and glial contributions. In our reduced consideration of the LFP interpretation problem we assume a layer of homogeneously distributed neurons and glial cell. As shown in , in the first approximation, the relationship between local field potential (LFP) and intracellular signals in layered neural tissue is governed by the balance of postsynaptic, leak, capacity and trans-dendritic currents. This contribution is approximated by a proportionality of the LFP to the dendritic and somatic voltage difference. In the second approximation, the currents during action potentials contributes into the LFP, giving an additional term proportional to the neuronal firing rate. However, this model does not take into consideration the contributions of glia and slow extracellular ionic concentration changes. Taking this into account might be especially important in cases of rather intensive synchronized neuronal activity such as epileptic discharge. In the resent paper, we assert that the next important factor (third-order approximation) is the change of the extracellular potassium concentration, [K+]e. The increase of [K+]e in time scales of tens and hundreds of milliseconds is caused by ionic flows through neuronal membranes, mainly due to voltage gated potassium channels, compared to slow effect of active and passive transporters. The increase is then buffered by glial cells, which effectively filter the impact of the transneuronal potassium flow . The glial membrane potassium currents contribute into the LFP, in addition to the neuronal membrane currents. As a result the LFP model has two terms. Neuronal contribution into LFP According to , the LFP at the somatic level, φN(t), is determined by the neuronal membrane voltage distribution and the firing activity as follows: \varphi_N (t)={2\, \sigma \, r_{i} } \bigl(V_d(t)-V_s(t)\bigr) + k_1 ~\nu(t) where Vs(t) and Vd(t) are the somatic and dendritic neuronal membrane potentials; neurons are homogeneously distributed with the density p per area; ri is the specific intracellular resistivity; σ is the mean conductivity of extracellular medium, assumed to be a constant; ν(t) is the firing rate; k₁ is a coefficient. Glial contribution into LFP The glial contribution into LFP, φG(t), is proportial to the glial membrane currents of potassium ions, thus can be approximated as follows: \varphi_G(t)=k_G \biggl[ {d[K^+]_e \over dt} \biggr]_{Glia}, where $\Bigl[ {d[K^+]_e \over dt} \Bigr]_{Glia}$ is the potassium concentration change due to glia; kG is a coefficient. External potassium was modeled in . Glial uptake was modeled in . To simplify the consideration, we proposed the following phenomenological model. Potassium currents through neuronal membranes are most pronounced during the spiking activity. That is why, approximately these currents are just proportional to the firing rate of principal neurons, ν(t). [K+]e is then buffered by glial cells. Concequently, the equation for the rate of potassium concentration change due to glia should look like the refractory equation, as well as the associated extracellular potential φG: \tau {\varphi_G(t) \over dt } = - \varphi_G(t) + k_2 ~\nu(t), where τ is the characteristic time constant; k₂ is a coefficient. The time constant cumulatively express the kinetics of neuronal and glial voltage-gated channels, as well as the scales of glial intracellular processes evoked by neuronal activity. Cumulative LFP model The contributions of potassium currents through mebranes of both neurons and glial cells are additive. Thus taking into account e.([e1]), we get \varphi(t) = k_1 ~\nu(t) + k_3 ~\bigl(V_d(t)-V_s(t)\bigr) + \varphi_G(t) where k₃ = p/(2 σ ri). A set of equations ([e3]) and ([e4]) determines the LFP model. According to these equations, LFP is governed by the somatic and dendritic membrane voltages Vs(t) and Vd(t) as well as the firing rate of principal neurons ν(t). The model might be matched to experimental data by fitting four coefficients, k₁, k₂, k₃ and τ.

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.