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.

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.

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 \(\mu\)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].

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 \(i_{PW}\) and \(j_{PW}\). The adjacent columns owing to different pinwheels have the same orientation preferences. The coordinates of the pinwheel-center are \(x_{PW}=(2\;i_{PW}-1)R\), \(\;y_{PW}=(2\;j_{PW}-1)R\). The orientation angle for the point (\(x\), \(y\)) of V1 which belongs to the pinwheel (\(i_{PW}\), \(j_{PW}\)) is defined as \(\theta(x,y)=\arctan((y-y_{PW})/(x-x_{PW}))\). The progression is determined by the factor \((-1)^{i_{PW}+j_{PW}}\). Finally, the input firing rate is

\begin{equation} \varphi_{Th\;E}(x,y,t)=\iint d\tilde{x}d\tilde{y}\,D^{LGN-V1}(x,y,\tilde{x},\tilde{y})\,L_{LGN}(\tilde{x},\tilde{y},t-\delta(x,y,\tilde{x},\tilde{y})),\\ \end{equation}where

\begin{equation} D^{LGN-V1}(x,y,\tilde{x},\tilde{y})=\frac{1}{\pi\;\sigma_{pref}\sigma_{orth}}\exp\left(-\frac{x^{\prime 2}}{\sigma_{pref}^{2}}-\frac{y^{\prime 2}}{\sigma_{orth}^{2}}\right),\\ \end{equation}where \(L_{LGN}\) is the activity of the non-lagged LGN neurons at (\(\tilde{x},\tilde{y}\)), (\(x_{cf},y_{cf}\)) are the coordinates of the centre of the footprint of V1 neuron in LGN, \(x_{cf}=x_{cf}(x,y)\), \(y_{cf}=y_{cf}(x,y)\); \(\delta\) 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:

\begin{aligned}
L_{LGN} & & (\tilde{x},\tilde{y},t-\delta(x,y,\tilde{x},\tilde{y}))\approx \\
& & \left[L_{LGN}(\tilde{x},\tilde{y},t)-\delta(x,y,\tilde{x},\tilde{y})\frac{\partial L_{LGN}(\tilde{x},\tilde{y},t)}{\partial t}\right]_{+}\nonumber \\
\end{aligned}

\begin{equation}
\delta(x,y,\tilde{x},\tilde{y})=30~{}ms,\quad\hbox{if}\quad(-1)^{i_{PW}+j_{PW}}x^{\prime}>0;\quad 0,\quad\hbox{otherwise},\\
\end{equation}

\begin{equation}
x^{\prime}=(\tilde{x}-x_{cf})\cos\theta-(\tilde{y}-y_{cf})\sin\theta,\\
\end{equation}

\begin{equation}
y^{\prime}=(\tilde{x}-x_{cf})\sin\theta+(\tilde{y}-y_{cf})\cos\theta.\\
\end{equation}

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: \(\bar{g}_{ampa\;Th\;E}=0.4\) mS/cm\({}^{2}\), \(\bar{g}_{ampa\;E\;E}=0.4\) mS/cm\({}^{2}\), \(\bar{g}_{ampa\;E\;I}=0.4\) mS/cm\({}^{2}\), \(\bar{g}_{gaba\;I\;E}=0.8\) mS/cm\({}^{2}\), \(\bar{g}_{gaba\;I\;I}=0.2\) mS/cm\({}^{2}\), \(\bar{g}_{nmda\;Th\;E}=0.8\) mS/cm\({}^{2}\), \(\bar{g}_{nmda\;E\;E}=1.6\) mS/cm\({}^{2}\), \(\bar{g}_{nmda\;E\;I}=1.6\) mS/cm\({}^{2}\), \(\tau_{r}^{AMPA}=1.7\)ms, \(\tau_{d}^{AMPA}=8.3\) ms, \(\tau_{r}^{GABA}=0.5\) ms, \(\tau_{d}^{GABA}=20\) ms, \(\tau_{r}^{NMDA}=6.7\) ms, \(\tau_{d}^{NMDA}=100\) ms, \(Mg=2\) mM, \(V_{AMPA}=V_{NMDA}=0\), \(V_{GABA}=-77\) mV. Region geometry: cortex region is 1 mm x 1.5 mm, containing 2x3 regularly distributed pinwheels. Numerical parameters: tim