Firing-rate based network modeling of the dLGN circuit: Effects of cortical feedback on spatiotemporal response properties of relay cells

Visually evoked signals in the retina pass through the dorsal geniculate nucleus (dLGN) on the way to the visual cortex. This is however not a simple feedforward flow of information: there is a significant feedback from cortical cells back to both relay cells and interneurons in the dLGN. Despite four decades of experimental and theoretical studies, the functional role of this feedback is still debated. Here we use a firing-rate model, the extended difference-of-Gaussians (eDOG) model, to explore cortical feedback effects on visual responses of dLGN relay cells. For this model the responses are found by direct evaluation of two- or three-dimensional integrals allowing for fast and comprehensive studies of putative effects of different candidate organizations of the cortical feedback. Our analysis identifies a special mixed configuration of excitatory and inhibitory cortical feedback which seems to best account for available experimental data. This configuration consists of (i) a slow (long-delay) and spatially widespread inhibitory feedback, combined with (ii) a fast (short-delayed) and spatially narrow excitatory feedback, where (iii) the excitatory/inhibitory ON-ON connections are accompanied respectively by inhibitory/excitatory OFF-ON connections, i.e. following a phase-reversed arrangement. The recent development of optogenetic and pharmacogenetic methods has provided new tools for more precise manipulation and investigation of the thalamocortical circuit, in particular for mice. Such data will expectedly allow the eDOG model to be better constrained by data from specific animal model systems than has been possible until now for cat. We have therefore made the Python tool pyLGN which allows for easy adaptation of the eDOG model to new situations.


Introduction
Visually evoked signals pass the dorsal geniculate nucleus (dLGN) on the route from retina to primary visual cortex in the early visual pathway. This is however not a simple feedforward flow of information, as there is a significant feedback from primary visual cortex back to dLGN. Cortical cells feed back to both relay cells and interneurons in the dLGN, and also to cells in the thalamic reticular nucleus (TRN) which in turn provide feedback to dLGN cells [1,2]. In the last four decades numerous experimental studies have provided insight into the potential roles of this feedback in modulating the transfer of visual information in the dLGN circuit [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Cortical feedback has been observed to switch relay cells between tonic and burst response modes [20,21], increase the center-surround antagonism of relay cells [16,17,22,23], and synchronize the firing patterns of groups of such cells [10,13]. However, the functional role of cortical feedback is still debated [2,[24][25][26][27][28][29][30].
Several studies have used computational modeling to investigate cortical feedback effects on spatial and/or temporal visual response properties of dLGN cells [31][32][33][34][35][36][37][38]53]. These have typically involved numericallyexpensive dLGN network simulations based on spiking neurons [31-33, 35, 38] or models where each neuron is represented as individual firing-rate unit [36,37]. This is not only computationally cumbersome, but the typically large number of model parameters in these comprehensive network models also makes a systematic exploration of the model behavior very difficult.
In the present study we instead use a firing-rate based model, the extended difference-of-Gaussians (eDOG) model [39], to explore putative cortical feedback effects on visual responses of dLGN relay cells. A main advantage with this model is that visual responses are found from direct evaluation of two-dimensional or three-dimensional integrals in the case of static or dynamic (i.e., movie) stimuli, respectively. This computational simplicity allows for fast and comprehensive study of putative effects of different candidate organizations of the cortical feedback. Taking advantage of the computational efficiency of the eDOG model, we here explore effects of direct excitatory and indirect inhibitory feedback effects (via dLGN interneurons and TRN neurons) on spatiotemporal responses of dLGN relay cells. In particular we investigate effects of (i) different spatial spreads of corticothalamic feedback and (ii) different corticothalamic propagation delays.
Our analysis suggests that a particular mix of excitatory and inhibitory cortical feedback agrees best with available experimental observations. In this configuration an ON-center relay cell receives feedback from ON-center cortical cells (ON-ON feedback), consisting of a slow (long-delay) and spatially widespread inhibitory feedback combined with a fast (short-delay) and spatially narrow excitatory feedback. Here the inhibitory and excitatory ON-ON feedback connections are accompanied by excitatory and inhibitory OFF-ON connections, respectively, following a phase-reversed arrangement [38]. For one this feedback organization accounts for the feedback-induced enhancement of center-surround antagonism of relay cells as observed in experiments [16,17,22,23,38]. Further, it seems well suited to dynamically modulate both the center-surround suppression and spatial resolution, for example, to adapt to changing light conditions [40].
Moreover, a longer thalamocortical loop time of ON-ON inhibitory feedback loop compared to ON-ON excitatory feedback may contribute to temporal decorrelation of natural stimuli [41], an operation that has been observed accomplished at the level of dLGN in the early visual pathway [42]. At the same time, the rapid excitatory feedback may contribute to linking stimulus features by synchronizing firing of neighboring relay cells [10,19].
Previous experimental studies have focused on cat, monkey, and ferret dLGN, and the present model was adapted to neurobiological findings from cat. However, the last years have seen a surge of interest in mouse visual system, where new optogenetic and pharmacogenetic methods provide new tools for precise manipulation of identified neurons in the thalamocortical circuit [43][44][45][46][47][48]. Such data will expectedly allow for a detailed adaptation of the eDOG model to rodent dLGN, likely much better constrained by biological findings than what has been possible until now for cat. To facilitate this we have made the Python tool pyLGN (http://pylgn. rtfd.io) which allows for easy modification and evaluation of the eDOG model to new situations.

Spatiotemporal receptive fields
Spike responses of neurons in the early visual pathway are most commonly described in terms of receptive fields. Mathematically, the spatiotemporal receptive field is defined by an impulseresponse function W(r, t). This function describes the firing-rate response to a tiny (δ-function) spot positioned at r = 0 which is on for a very short time (δ-function) at t = 0. If linearity is assumed, the response to any stimulus S(r, t) can be found by convolving the impulseresponse function with the stimulus [39, [49][50][51][52]: or written more compactly Here S(r, t) is a spatiotemporal stimulus function describing, e.g., the light intensity on a screen as a function of time and position. R(r, t) is the response of a neuron with its receptive-field center at r. The spatial integral goes over the whole visual field, i.e., over all two-dimensional space. For mathematical convenience we have chosen the temporal integral to go from τ = −1 to +1. Since a stimulus input cannot affect the response in the past, it then follows that W(r, τ < 0) = 0.
In Fourier space the convolution in Eq (2) corresponds to a product Rðk; oÞ ¼W ðk; oÞSðk; oÞ; whereR,W , andS are the Fourier transforms of the neural response R, the impulse-response function W, and the stimulus S, respectively. The tilde symbol ( * ) will be used to denote the Fourier transform of any function throughout this paper. The function argument k is the wave vector which is related to the spatial frequency ν via |k| = 2πν. Correspondingly, the angular frequency ω is related to the temporal frequency f via ω = 2πf. WithW andS known, the neural response can thus always be found by an inverse Fourier transform F À 1 fg, which entails an integral over temporal and spatial frequencies Rðr; tÞ ¼ F À 1 fW ðk; oÞSðk; oÞg: ð4Þ The response model in Eq (1) is an example of a descriptive model where the purpose is to summarize experimental data compactly in a mathematical form [51][52][53]. Here the aim is to find an appropriate impulse-response function, i.e., spatiotemporal receptive-field function, that describes the measured neural response to different visual stimuli [50]. With this approach, however, limited insight is gained into how the neurons and neural circuitry in the early visual system provide such a receptive field. To address this question a mechanistic receptive-field model is needed. For a discussion of the difference between descriptive and mechanistic models in visual neuroscience, see [53,54].

Mechanistic receptive-field models
In mechanistic LGN-circuit models the input from retinal ganglion cells have been described by descriptive models, see, e.g., [36,37,39,52,55]. Likewise, in the present eDOG model the input from retinal ganglion cells is represented by the descriptive impulse-response function (Eq (1)). Here a square grid of retinal ganglion cells with identical, spatially-localized receptive fields are considered (see Fig 1). The activity, i.e., firing rates, of the neurons on the retinal ganglion cell layer then serves as input to the dLGN relay cell layer. This is represented by a spatiotemporal coupling-kernel function K RG , which reflects the direct synaptic input from retinal ganglion cells to dLGN relay cells. The coupling kernel, which is analogous to the descriptive impulse-response function in Eq (1), is assumed to only depend on the relative distance between the cells in the visual field [52].
The response of a relay cell located at r is then given by [39,52]: where R R and R G are the firing-rate responses of relay cells and ganglion cells, respectively. The coupling kernel K RG (r − r 0 , τ) denotes the strength with which the response of a ganglion cell, displaced by r − r 0 from the relay cell, at time t − τ influences the response of the latter at time t. Note that, K RG (r, τ < 0) = 0 due to causality. In Fourier space the relationship in Eq (5) can be written as where we have used the general relationship in Eq (3). The key point here is that a descriptive model for the relay-cell impulse-response functionW R now has a mechanistic interpretation. This relation is given as the product of the impulse-response functionW G of the retinal ganglion cells and the coupling kernelK RG from the former cell type to the latter. In the eDOG model this approach is extended to include the various feedforward and feedback connections affecting the relay-cell response. The result is an expression for the relay-cell impulse-response functionW R in terms of the impulse-response functionW G of the retinal input and the coupling kernels connecting the neurons of the circuit. With such a mechanistic expression forW R , the response to any visual stimulus can be computed by means of the inverse Fourier transform in Eq (4).

Extended difference-of-Gaussians (eDOG) model
Here we derive the impulse-response function for dLGN relay cells for the mechanistic eDOG model [39]. The complete circuit is shown in Fig 2. In this figure each cell type correspond to a two-dimensional layer (or population) of identical cells.
We will in the following focus on the dLGN relay cells with ON symmetry, but a similar model can be constructed for OFF-symmetry cells. These neurons receive feedforward excitation and indirect feedforward inhibition (via intrageniculate interneurons) from ON-center ganglion cells in retina. The relay cells further receive cortical feedback from both cortical ON cells and cortical OFF cells.
Feedforward input from retina. With indirect feedforward inhibition included in addition to the direct feedforward excitation, the expression in Eq (5) generalizes to Here K ON RIG is a spatiotemporal coupling-kernel representing the indirect feedforward inhibition from retinal ganglion cells onto relay cells via intrageniculate interneurons.
In Fourier space this gives a simple expression for the relay-cell impulse-response function, i.e.,W where we have used thatRðk; oÞ ¼W ðk; oÞSðk; oÞ, cf. Eq (3). Feedback from cortex. Next we add effects from cortical feedback onto the relay cell. This cortical feedback can be both excitatory and inhibitory. The excitatory feedback corresponds to direct projections from cortical cells onto relay cells. The inhibitory feedback corresponds to indirect inhibitory action on relay cells mediated by cortical projections onto inhibitory TRN and intrageniculate interneurons. Further, unlike the feedforward projection, the feedback is cross-symmetric, i.e., the activity of ON relay cells are affected both by ON and OFF cortical cells.
In the eDOG model cortical ON and OFF cells are assumed to be driven solely by ON and OFF relay cells, respectively. As the corticogeniculate feedback comes from orientation-tuned cells in layer 6 in cortex, we include a set of N mutually uncoupled, orientation-selective cortical populations C n , (n = 1, 2, . . ., N) for both the ON and OFF pathways. Each population C n responds preferably to stimuli (bars, gratings) with orientation θ n . In Fig 2 only a single cortical population is shown for each pathway even though an arbitrary number of N cortical populations can be considered.
Cortical cells are known to exhibit substantial non-linearities when responding to visual stimuli, and here the response is modeled via a static non-linear function acting on a linearly filtered input [51,56]. More specifically we express the response of the ON or OFF cortical population C n by where K ON=OFF C n R is the feedforward kernel between the relay cells and the cortical cells in population C n . Further, the half-wave rectification function H[x] = xθ(x) is used to enforce non-negative firing rates [57], where θ(x) is the Heaviside step function.
We further assume the input to cortical OFF cells to be the negative of the one for the ON cells [53]. That is Finally, the feedback cross-connection (OFF to ON) is assumed to be phase-reversed compared to the same-sign feedback (ON to ON) [39]: where K OFFÀ X RC n is the cross-coupling feedback from cortical OFF cells onto relay ON cells. In other words, we assume the effect of ON-center and OFF-center cortical cells to be the opposite of each other. However, we do not make any specific assumptions on whether, say, the excitatory or inhibitory feedback is driven by ON-center or OFF-center cortical cells [39].
With the three assumptions in Eqs (9)-(11), the total input to the ON dLGN relay cell is found to be [39, 53] where we have used the mathematical identity: In Fourier space we thus havẽ and in analogy with Eq (8) we find after some simple algebrã In this expression the direct feedforward excitation and the indirect feedforward inhibition via interneurons are represented by the first and second terms in the numerator, respectively. The feedback effects are accounted for in the denominator. The general mathematical expression in Eq (14) for the (Fourier transformed) impulseresponse function for the relay cells is the main feature for the eDOG model [39]. The model provides an analytical formula for (linear) impulse-response function for relay cells, despite the non-linearity of the response of the cortical cells providing the feedback. The simulator presented in this paper uses this expression as basis to compute the impulse-response function and the spatiotemporal responses for user-defined kernels and input stimuli. Once the explicit form of the kernels in Eq (14) are defined, the response of the relay cells to arbitrary stimuli can be calculated using Eq (4).
In the next subsections we describe the choices made in this paper for (i) the descriptive spatiotemporal receptive-field function for the retinal input (W ON G in Eq (14)), (ii) the various mechanistic coupling kernels inside the dLGN circuit (K in Eq (14)), and (iii) the visual stimulus (S in Eq (4)). The coupling kernels are assumed to be space-time separable (e.g., K(r, t) * f (r)h(t)), but space-time coupled kernels can equally be used in the eDOG-model. The same applies to the choice of the receptive-field function of the retinal input [58,59]. For presentational simplicity, we will focus on the ON-pathway and skip the ON-superscript on the connectivity kernels and impulse response functions (K 's andW 's in Eq (14)), but an analogous model can be derived for the OFF pathway.
Impulse-response function of input from retinal ganglion cells. The impulse-response function of the retinal input is modeled as a product of a spatial part F(r) and temporal part H (t). The spatial part is described by means of the difference-of-Gaussians (DOG) model [60]: where the first and second term correspond to the center and surround contribution, respectively. Further, A and B (defined to be positive) are the strengths of the center and surround, and a and b are the corresponding width parameters. In the present paper we have used parameters extracted from fitting the function to retinal-input responses to flashing circular spots [55]. The temporal part of the impulse response of the retinal input is modeled as a biphasic temporal function [37, 53]: where B is the weight for the second phase, and τ is the duration of each phase. The same parameter values as in [53] has been used, which correspond to the mean of the range of values reported by [61].
For an illustration of the shapes of the spatial and temporal impulse-response function, see Fig 3A. Coupling kernels inside dLGN circuit. The kernels K(r, t) are considered to have separable space-time parts, i.e., where f and h are normalized spatial and temporal parts, respectively, and w is the connection weight of the kernel. The latter is positive for excitatory synaptic connections and negative for inhibitory connections. The normalization implies that RR f(r)d 2 r = R h(t)dt = 1 where the integrals go over all visual (two-dimensional) space and all times, respectively.
The spatiotemporal coupling-kernels in the circuit, reflecting how the firing in one type of cell affects the firing in another type of cell through their direct synaptic connections, have not been systematically mapped out. However, a key design principle of the early visual pathway is retinotopy, i.e., that neurons representing neighboring positions in the visual field also are neighbors inside the retina, dLGN, and visual cortex. This implies that the coupling kernels are spatially confined. In this paper we describe the shape of spatial kernels using the mathematically convenient Gaussian function: where a is the width parameter. The temporal part of the kernels is modeled as (delayed) exponential decay in accordance to previous modeling studies [52,53]: where τ is the time constant, and Δ corresponds to a combined axonal and synaptic time delay.
For an illustration of the shapes of spatial and temporal part of the coupling kernels, see Fig  3A.
We next describe the kernel parameters used for the circuit coupling. A detailed list of these parameters is given in Table 1.
Feedforward couplings. Relay cells in the cat appear to receive input from a single or a few retinal ganglion cells [62][63][64][65][66][67][68]. Further, the relay cells receive indirect feedforward inhibition via intrageniculate interneurons which in turn receive input from a few retinal ganglion cells [67,69]. Based on these observations and the known'retinotopic' organization of the early visual pathway, we here use narrow Gaussian functions as coupling kernels between the retinal ganglion cells and dLGN relay cells [39,52]. We assume a larger width parameter for the feedforward inhibitory coupling kernel compared to the excitatory kernel [55], reflecting the observed larger receptive field in intrageniculate interneurons compared to both retinal ganglion cells and relay cells [67].
Feedback coupling. The net feedback coupling from cortex to relay cells are determined by two factors: (i) the spatiotemporal response of the cortical cells providing the feedback, and (ii) the spatiotemporal feedback coupling kernels from cortical to LGN cells. The receptive fields of simple cortical cells arises primarily from convergent input from ON and OFF relay cells [70][71][72]. In order to model orientation-selective cortical populations, the thalamocortical kernels K C n R in Eq (9) must have an elongated shape. In [39] these kernels were, for example, modeled as elliptical Gaussians.
As seen in the denominator of Eq (14), the total effect of cortical feedback is a sum over feedback contributions from all n populations, covering all orientation angles. Thus, the net feedback effect is expected to be essentially circularly symmetric [39]. The net effect of the cortical feedback from all can thus be incorporated in the model via a single circularly-symmetric Table 1. List of kernel parameters. W G is the impulse-response function of ganglion cells, K RG and K RIG are the excitatory and inhibitory feedforward kernels, respectively. K ex RCR and K in RCR are ON-ON excitatory and inhibitory thalamo-cortico-thalamic kernels, respectively. K mix RCR denotes the mixed ON-ON feedback kernel, consisting of an excitatory and an inhibitory term. The subscript 'RIG' refers to the indirect inhibitory input from retinal ganglion cells onto relay cells via intrageniculate interneurons (ganglion ! interneuron ! relay), while the subscript 'RCR' refers to the complete thalamo-cortico-thalamic loop (relay ! cortex ! relay). F represents the DOG function, f represents the Gaussian function, H represents the biphasic temporal function, and h represents the delayed decaying exponential function. The width parameters in spatial functions are given in units of degree, while the temporal parameters are in units of ms. In the present example applications we have kept the time constant τ fixed at 5 ms (comparable to what, e.g., was found in [74]), while the temporal delay parameters Δ have been varied in a range of 5-30 ms. † denotes the default values for parameters that have been varied.

Kernel
Weight Spatial Temporal coupling kernelK RCR P nK RC nK C n R . As for the feedforward couplings, we for simplicity model the feedback coupling kernels as product of a Gaussian function of space (Eq (18)) with a delayed exponentially-decaying temporal function (Eq (19)).
The structure of the eDOG model is indifferent to whether the cortical feedback is excitatory, inhibitory, or even a mix of excitatory and inhibitory feedback. For excitatory feedback the weight parameter w in Eq (17) is positive, while for inhibitory feedback it is negative. For mixed feedback the coupling kernelK RCR consists of a sum of excitatory and inhibitory feedback terms. Note that in all cases the ON to ON couplings are accompanied by OFF to ON couplings with the opposite sign, i.e., a phase-reversed arrangement as described in Eqs (10) and (11).
A few experiments give some hints about how the feedback may be organized. In [3] a center-surround feedback configuration was reported in cats where feedback was excitatory when the cortical and relay cell receptive field centers were close to each other and inhibitory when they were further apart. This observation was later supported by [16], where they found in primates a center-surround configuration for feedback, with a facilitatory bias to center and inhibitory surround (but see also [18]). Further, in [73] a particular cross-symmetry organization was observed where a same-symmetry inhibitory feedback was accompanied by an excitatory feedback with opposite symmetry, e.g., ON-ON inhibitory feedback accompanied by OFF-ON excitatory feedback.
In this paper we will study three different spatial organization of the cortical feedback as shown in the list below and illustrated in Fig 3C. In this list ON-ON refers to feedback from ON-center cortical cells to ON-center relay cells, while OFF-ON refers to feedback from OFFcenter cortical cells to ON-center relay cells.
• ON-ON excitatory feedback (K ex RCR ) combined with OFF-ON inhibitory feedback. • ON-ON inhibitory feedback (K in RCR ) combined with OFF-ON excitatory feedback. • Mixed ON-ON excitatory and inhibitory feedback (K mix RCR ). The OFF-ON feedback is also both excitatory and inhibitory.
The superscripts 'ex' (excitatory) and 'in' (inhibitory) refer to the sign of the ON-ON feedback, and the subscript 'RCR' refers to the complete thalamo-cortico-thalamic loop (relay ! cortex ! relay). These three scenarios are illustrated in Fig 3C. The second scenario corresponds to the configuration observed experimentally in [73], while the last configuration is inspired of the center-surround configuration suggested by data from [3,16]. For simplicity we will in the following refer to ON-ON excitatory feedback as just excitatory feedback and ON-ON inhibitory feedback as inhibitory feedback. It is then implicitly assumed that the influence from the OFF-ON feedback has the opposite sign.
The influence of each of these feedback configurations on the relay cell responses is investigated for a range of feedback strengths w, width values a for the Gaussian functions (Eq (18)), as well to temporal delays Δ of the delayed exponential functions (Eq (19)). In Fig 3B the interpretation of these parameters is illustrated.
Visual stimuli. With the general eDOG relay-cell impulse-response function expression from Eq (14), specified by the coupling kernels above, all that is needed to compute the relaycell response by means of Eq (4), is a mathematical expression for the stimulus S(r, t). The two main visual stimuli considered in the present work are (i) circular patch gratings and (ii) fullfield gratings. For a full-field drifting grating, specified by k g and ω g , the relay-cell response is essentially given by Fourier-transformed impulse response in Eq (14) [39].
For a circular patch of drifting grating, the stimulus can be described mathematically as [39,75] Sðr; tÞ ¼ C pg cos ðk pg r À o pg tÞ½1 À yðr À d pg =2Þ; ð20Þ where k pg and ω pg are the wave vector and the angular frequency of the patch-grating, respectively, d pg is the diameter of the patch-grating spot, and C pg is a measure for the contrast of the grating. In all calculations presented in this paper C pg = 1. Note that a static circular patch (spot) is obtained for k pg = ω pg = 0. In the limit d pg ! 1 the Heaviside function in Eq (20) is always zero, and we obtain the simple harmonic function representing a full-field grating.
In addition, natural stimuli (images and movies) are also used. The stimulus is then given as an array of numbers, and the Fourier transform of the stimulus is calculated numerically.

Implementation in pyLGN
In order to allow for easy exploration of the eDOG model and in particular effects of cortical feedback on relay-cell responses, we have developed an efficient, firing-rate based simulator of spatiotemporal responses in the early visual system. The simulator is named pyLGN and is written in Python. The design goals for pyLGN are to provide a software framework for studying the cortical feedback effects that is easy to use, extensible, and open. To facilitate usability, pyLGN has its own documentation page including installation instructions, several usage examples, and technical aspects (http://pylgn.rtfd.io). To achieve extensibility, objectoriented programming is used, making it possible for the user to define new connectivity kernels and input stimuli. Lastly, to support openness pyLGN is both open-source and multiplatform.
All calculations presented in this paper have been tracked using the Python software Sumatra [76], which is an automated tracking tool for computational simulations and analysis. The source code for all presented simulations is available at (https://github.com/miladh/ edog-simulations).

Results
The result section is divided into two distinct parts. In the first part, results for the effects of cortical feedback on the spatial response properties of relay cells are presented. The cortical feedback effects on temporal aspects are presented in the second part.

Effect of cortical feedback on spatial properties
Spatial receptive fields. We start our study of the dLGN network model by characterizing the effects of cortical feedback on the spatial aspects of the relay cell's receptive field structure. From Eq (14) we see that even with separable kernels the impulse-response function, in general, remains non-separable in space and time. However, with a static stimulus, the spatial response properties can be studied in isolation.
Mathematically, the Fourier transformS for a static stimulus is / δ(ω). The convolution in the response integral in Eq (4) is then given by where F À 1 spatial is the spatial inverse Fourier transform, andSðkÞ is the Fourier transform of the spatial part of the stimulus.
Using the kernels shown in Table 1 we then find that the static relay-cell impulse response functionW R ðk; 0Þ is given bỹ where we have used thathð0Þ ¼ 1 and thatH ð0Þ is a constant. For simplicity will we hereafter refer toW spatial R ðkÞ as the spatial impulse-response function of relay cells. Examples of spatial receptive fields, found by an inverse Fourier transform of this function, is shown in Fig 4. As seen in this figure the center-surround receptive field structure of the retinal ganglion cells is qualitatively preserved, in accordance with the notion that cortical feedback has a mainly modulatory effect on response properties of relay cells [77][78][79]. A close inspection of the right panel in the figure reveals that while the response at the receptive-field center (peak value) is increased for excitatory feedback, it is reduced for inhibitory feedback. The cortical feedback effects outside the receptive-field center, on the other hand, are less clear-cut for the examples in the figure.
Next, we investigate the spatial impulse-response function in more detail. In particular we study how the spatial responses depend on the weights (w) and Gaussian width parameters (a) of the connections, see top panel of Fig 3B. We characterize the spatial receptive-field structure  weight w RIG and width a RIG are shown in the top row of Fig 5. The clear tendency is that narrow kernels with high weights most effectively reduce the center excitation and surround inhibition. The largest reduction in the receptive-field size is also observed in this situation. Another observation is that inhibitory kernels widths a RIG similar to the width b G *1.3 deg of the DOG surround of the ganglion-cell input, combined with large weights, give a large surround inhibition.
For the situation with feedback inhibition only (Fig 5, middle row) an overall similar tendency is observed. However, the effects of feedback inhibition are somewhat weaker compared to feedforward inhibition for these example parameter ranges. Finally, we see from bottom row in Fig 5 that strong and narrow excitatory feedback strongly increases the center excitation and surround inhibition. Larger widths, however, reduce the surround inhibition significantly and also results in larger receptive-field sizes.
To see the influence of a mixed cortical feedback on the spatial receptive-field properties, we show in Fig 6 the effects of increasing cortical feedback weights and widths. A configuration consisting of a narrow excitatory and a broader inhibitory feedback both increases the excitation in the center and the inhibition in the surround. A large reduction in the receptive-  Table 1).
https://doi.org/10.1371/journal.pcbi.1006156.g006 field center size is also seen with this configuration, specially for inhibitory width values a in RCR close to one.
The bottom row of Fig 6 shows the effect of increasing feedback weights for a narrow excitatory central core projection and a wider inhibitory projection. Strong excitatory feedback combined with a weak inhibitory feedback increases excitation in the center and inhibition in the surround. In contrast, strong inhibitory feedback combined with weak excitatory feedback, reduces the center excitation. Note however that the effects due to strong excitatory feedback are more significant than the ones due to the inhibitory feedback. This is specially obvious in the surround inhibition which is nearly completely dominated by the excitatory feedback strength. The size of the receptive field decreases with increasing excitatory and inhibitory feedback strength.
In conclusion, these results show that the cortical feedback is well suited to modulate the center-surround organization of relay-cell receptive fields. Excitatory and inhibitory inputs have opposite effects on a relay cell's spatial response: while excitatory feedback can increase the center excitation and center size, inhibitory feedback can do the opposite. Depending on the width of the feedback projection, both excitatory and inhibitory feedback can either increase or decrease surround excitation. A mixed feedback configuration consisting of a combination of narrow excitatory and a broader inhibitory feedback, both increases the excitation in the center and the inhibition in the surround.
In previous experiments the removal of cortical feedback has been observed to give larger receptive-field center sizes [17,19]. The results in Fig 5 demonstrate that inhibitory feedback gives results in accordance with this, while excitatory feedback has the opposite effect. Likewise, this experimentally observed receptive-field center shrinkage effect by cortical feedback is also observed (to various degrees) for the mixed feedback configurations depicted in Fig 6. Area summation curves. A common way to experimentally probe the center-surround organization of cells in the early visual pathway is to measure area-response curves, i.e., the response to circular stimulus spots as a function of spot diameter [16,17,22,55,59,78,80,81]. In the top row of Fig 7 we correspondingly show area-response curves for relay cells responding to static bright-spot stimuli for different feedback configurations. Here the receptive field of the cell is set to be concentric with the spot.
The top left panel in Fig 7 shows that an increasing excitatory feedback enhance the excitatory response to stimuli restricted to be within the receptive field center. It also reduces the suppressive effects of stimuli in the surround area. Inhibitory cortical feedback, on the other hand, reduces the response to optimal patch diameter and enhances the suppressive effects for large patch sizes (top center panel in Fig 7).
In the top right panel of Fig 7, the mixed feedback situation with a combination of narrow excitatory and a broader inhibitory feedback, as suggested by experimental findings [3,16], is considered. Here we observe that an increased feedback strength both (i) enhances the excitatory response to stimuli restricted to be within the receptive field center, and (ii) enhances the suppressive effects of stimuli in the surround area. Stronger feedback also reduces the receptive-field center size, i.e., the spot diameter giving the maximum response.
The area-response curves in Fig 7 are for static spot stimuli, but area-response curves are also commonly recorded for patch-grating stimuli [17,22,38,80]. In our formalism such response curves are readily obtained by use of the circular patch-grating stimulus function S in Eq (20). The resulting area-response curves typically resemble the static-spot curves shown in the top row of Fig 7, and we do not show any example curves here.
However, in the bottom row of Fig 7 we summarize results for area-response curves both for static-spot and patch-grating stimuli. Here, the stimulus size giving the largest response (corresponding to the receptive-field center size for static spot stimuli) and center-surround suppression index are shown as a function of feedback strength. This suppression index, α s , is here defined as: where R max is the maximum response, and R plateau is the response when the large-diameter plateau is reached(see bottom left panel in Fig 7). The figure shows that the suppression index for static spot stimuli (bottom right panel, solid lines) is increased with stronger feedback weights both for inhibitory and mixed feedback. The same qualitative trend is also observed for patch-grating stimuli (dashed line). Here the suppression index without feedback is fairly small (*0.4), but increases more strongly with feedback strength than for static-spot stimuli. This relative difference in suppression index between spot and patch-grating stimuli is qualitatively in agreement with experimental observations [17,22]. With excitatory feedback on the other hand, the suppression index is  reduced with increasing feedback strength. The largest suppression indices are found for mixed feedback, again illustrating that such feedback is particularly suited for modulating center-surround antagonism.
The optimal stimulus size (Fig 7, bottom center panel) is seen to be the same for static-spot and patch-grating stimuli. For both stimulus types this size is seen to decrease with increasing feedback strength both for inhibitory and mixed feedback, while the opposite is true for excitatory feedback.
Note that in experimental measurements of spot area-response curves, 'flashing' spots rather than static spots have been used [22,59,78]. This means that the spots were 'flashed' on and the subsequent response, which contained both a transient and a sustained response, were used to compute the area-response curves [78]. Here we have focused on the spatial response properties of relay cells in isolation by using Eq (22) where the responses are directly obtained without any assumptions about the temporal properties of the circuit, such as the time constants and delays of the synaptic connections. Therefore the static-spot area-response curves in Fig 7 corresponds to the sustained response. However, calculations with flashing spot has also been performed, and the observed feedback effects on the tuning curves are similar to the case with static spot stimulus (see S1 Fig in Supporting information). For the patch-grating experiments in [17,22] drifting patch gratings with a temporal frequency of only ω*6 Hz was used, so that in the 'fast-loop limit' (i.e., assuming sufficiently short propagation times around the thalamocortical loop) the expression in Eq (22)  For the smaller patch the frequency characteristic corresponds to a low-pass filter for all feedback configurations, cf. upper row panels in Fig 8. Increasing feedback strength leads as expected to higher response values for excitatory feedback, but also for the mixed feedback. For inhibitory feedback the opposite is the case.
For the larger patch size, (effectively corresponding to a full-field grating), relay cells have band-pass characteristics in all cases, cf. lower row panels in Fig 8. Excitatory feedback is seen to overall increase the response as well as shift the frequency giving the maximum to smaller frequencies. Inhibitory feedback is seen to have opposite effects. For mixed feedback an interesting combination of these effects are seen, i.e., the maximum-frequency response is shifted towards higher frequencies, but the maximum amplitude is also increased. This shift from low-pass to bandpass characteristics when changing the grating size can be understood when considering the center-surround organization of the receptive field. When the stimulus only covers the center of the relay-cell receptive field, the filtering of the circuit is effectively Gaussian-like, i.e., a low-pass filter. However when the stimulus also covers the surround region, the circuit filter is effectively an antagonistic center-surround filter with bandpass characteristic.
In [11] it was found that in cat, cortical feedback enhanced this band-pass feature observed for full-field gratings, that is, enhanced the suppressive surround effects at low spatial frequencies. This was demonstrated by computing, both with and without cortical feedback, the percentage response reduction that occurs when both the center and surround of the receptive field are simulated compared to when only the center is stimulated [11,Figs. 8 and 9]. For the lowest frequencies they observed about a 65% reduction with cortical feedback for X cells, while only about a 45% reduction was observed when cortical feedback was removed. For the model example results in Fig 8 we find 80% reduction with mixed feedback (w mix RCR ¼ 1:0) and 70% reduction without feedback for |k pg | = 0.25/deg, but the exact values of these reduction factors will depend on the model parameters.
A putative benefit of the shift of the response towards higher frequencies observed for our mixed feedback, can be alluded to in the context of information theory and efficient coding. In natural scenes there are usually extensive spatial correlations since neighboring regions often have similar luminance values [29]. This leads to a power spectrum of the input with large contributions from low spatial frequencies. An antagonistic center-surround organization dampens the low-frequency components and enhances the higher frequency components of the image and reduces the redundancy in the signal conveyed to the cortex. The shift towards higher frequencies sharpens the spatial receptive field of the relay cells and thereby increases the saliency of edges. To illustrate this point, we show in Fig 9 the response map for the relay cells for different circuit configurations with a natural image as stimulus. Here, we see that inhibitory feedback reduces the response magnitude to both high and low spatial frequency components in the image compared to the case without feedback. The opposite is seen for excitatory feedback: both high and low frequencies are enhanced in magnitude compared to  Table 1 have been used for fixed parameters. https://doi.org/10.1371/journal.pcbi.1006156.g008 Effects of cortical feedback on response properties of dLGN relay cells the no-feedback case. With mixed feedback on the other hand, the high frequency components (i.e., edges) are slightly enhanced while the low frequency components (i.e., areas with small luminance changes) are reduced in magnitude. Thus, a representation highlighting contrasts in the image is obtained. Mixed feedback has different effect on low and high frequency components of natural scenes in contrast to pure excitatory or inhibitory feedback. Each subfigure shows activation of a layer of relay cells in response to the input image, shown as a logarithmic color map from blue to red (from reduced to increased response). the responses are normalized with respect to the maximum response in the case without cortical feedback. The red and blue circles mark representative parts of the image with high and low spatial frequency, respectively. Default values from Table 1

Effect of cortical feedback on temporal properties
Temporal receptive fields. We have so far focused on influence of cortical feedback on spatial response properties of relay cells. Next, we investigate the effect of cortical feedback on temporal properties. The relay cell impulse-response function in Fourier space in Eq (14) can be written as Here,F andH represent the assumed spatial and temporal response functions for the retinal input, i.e., Eqs (15) and (16), respectively. Further the coupling kernels have been expanded into products of spatial (f ) and temporal (h) functions, cf. Eq (17).
Illustrations of the real-space version W R (r, t) of this impulse-response function is shown in Fig 10 using the kernel parameters listed as default parameters in Table 1. In this figure the temporal evolution of the spatial structure of the receptive field is shown (top panel), in addition to an x-t plot of the impulse-response function, summarizing how the one-dimensional spatial organization of the receptive field changes with time. These figures illustrate the biphasic nature of the center and surround responses, as has been observed experimentally [61,82]. For t between 0 and 50 ms the impulse-response function exhibits a bright-excitatory center, i.e., an increased firing to a tiny bright test spot placed in the receptive-field center. However, for times later than 60 ms, the polarity of the center response is reversed, and becomes dark-excitatory, i.e, increased firing-rate for dark spots. A similar behavior is observed for the surround, but with opposite polarities.
Experimental studies have reported several distinct temporal receptive-field profiles among relay cells, including monophasic and triphasic responses in addition to the more common biphasic response seen in Fig 10 [ 71,83]. To explore scenarios where these different response profiles may arise, we next study how different model parameters change the shape of the realspace temporal impulse response. In accordance with previous experimental and computational studies we use the biphasic index (I BP ) and peak response latency (t peak ) as measures to characterize the temporal properties of the impulse-response function [53,61,83]. The biphasic index is defined as the ratio between the peak magnitude of the (negative) rebound phase and the peak magnitude of the first (positive) phase, and thus measures how biphasic the response is (Fig 10, bottom right). A biphasic index equal to one means a perfect biphasic response, while zero corresponds to a monophasic response. Fig 11 illustrates the dependency of the temporal part of the impulse-response function, as well as I BP and t peak , on various circuit configurations and model parameters. Each temporal coupling kernel is described by two parameters: the time constant τ of the exponential decay and the parameter Δ accounting for delay in the propagation of the signals between the different neuronal population. In the present examples we keep the τ fixed and instead focus on how different values of the Δ affect the temporal impulse-response function. The parameter range for Δ is chosen based on the findings in [84], where a large variation in axonal conduction times was found for corticothalamic axons, from a few milliseconds to several tens of milliseconds.
In Fig 11 it is seen that the effects of feedforward and feedback inhibition on the temporal impulse are qualitatively similar. In both cases the depth of the second phase is increased for delayed (large Δ) inhibitory input, i.e., the biphasic index I BP is increased. The peak response latency t peak is seen to be substantially reduced with feedforward inhibition, but less so for feedback inhibition. Finally, we see from the middle left panels in Fig 11 that for strong delayed inhibitory inputs, a triphasic impulse-response may arise.
Focusing on the middle and bottom row of panels in Fig 11, we see that the effect of excitatory feedback is essentially the opposite of that for inhibitory feedback. The biphasic index I BP is decreased by the excitatory feedback and increased by the inhibitory feedback, in particular for delayed inputs. Further, the peak response latency t peak is increased with excitatory feedback, in contrast to with inhibitory feedback.
In Fig 12 the effect of temporal kernel parameters on the biphasic index I BP and t peak is shown for the more complex situation with a mixed excitatory and inhibitory feedback. Here the spatial spread as well as the relative weight of the excitatory and inhibitory feedback is kept Effects of cortical feedback on response properties of dLGN relay cells fixed while the delay parameters (D ex RCR , D in RCR ) are varied. With long-delay excitatory feedback and short-delay inhibitory feedback (D ex RCR ¼ 30 ms, D in RCR ¼ 5 ms), the first positive phase of the impulse-response function is only modestly affected by the feedback. In this case the feedback mainly affects the second negative phase, which generally is reduced in depth. Thus the response becomes more monophasic, as reflected in smaller values for the biphasic index ( Fig  12, rightmost panel).
In the case with long-delay inhibitory feedback and short-delay excitatory feedback (D ex RCR ¼ 5ms, D in RCR ¼ 30ms), the rapid excitatory feedback boosts the initial positive peak while the former boosts the following negative peak. Further, this feedback combination gives multiphasic responses, i.e., distinct responses also after the initial biphasic response. parameter dependence of two impulse-response measures t peak and biphasic index I BP . The biphasic index is normalized with respect to the value for the case with feedforward excitation only (I BP = 0.35), while the t peak plots show the difference in peak time in milliseconds compared to the corresponding value for feedforward excitation only (t peak = 29 ms). Default parameters have been used for the fixed parameters (see Table 1). https://doi.org/10.1371/journal.pcbi.1006156.g011 The figures in the right panel of Fig 12 show that the peak response latency t peak and the biphasic index I BP are reduced for long-delayed excitatory feedback. Long-delayed inhibitory feedback combined with short-delayed excitatory input, on the other hand, increases both t peak and I BP .
Temporal frequency tuning curves. To investigate the effect of cortical feedback further, we next explore the temporal frequency tuning of relay cells. This is done by computing the response of the relay cells to a full-field grating stimulus for a range of different temporal frequencies [85][86][87][88]. In this case the response is given by the magnitude of the Fourier-space impulse-response function in Eq (24) evaluated with a varying angular frequency ω g combined with a fixed spatial wave vector k g [75]. In the present example the wavenumber is kept fixed at |k g | % 1 deg −1 .
In Fig 13 (top row) the frequency tuning is shown for inhibitory(left) and excitatory feedback(right). As expected, rapid inhibitory feedback (D in RCR ¼ 5 ms) reduces the overall response. In addition, it is also seen to shift the peak frequency to slightly higher values. Longdelayed inhibitory feedback (D in RCR ¼ 30 ms), on the other hand, gives sharper tuning curves, enhancing the band-pass characteristics of relay cells. In this case increased feedback weights both sharpen the resonance and shift it to higher frequencies.
Excitatory and inhibitory feedback essentially have opposite effects on the tuning properties of relay cells (Fig 13, top row, right panel). Rapid excitatory feedback shifts the tuning curve to higher response values, while long-delayed excitatory feedback blunts the tuning profile. In all cases the peak frequency is shifted to lower frequencies.
In the bottom row of Fig 13 we next investigate the effect of mixed cortical feedback on temporal frequency tuning properties of relay cells. In particular, we consider three different cases: (1) long-delay inhibitory feedback combined with short-delay excitatory feedback, (2) long-delay excitatory feedback combined with short-delay inhibitory feedback, and (3) synchronized feedback where excitatory and inhibitory feedback are received at the same time. Long-delayed inhibition (combined with rapid excitatory feedback) leads to sharper tuning with increasing feedback. This contrasts the case with long-delay excitation (combined with rapid inhibitory feedback), where a more flat spectrum is observed with increasing feedback. Synchronized excitatory and inhibitory feedback does not change the tuning properties by much when comparing with the case without cortical feedback. This reflects that in this case  Table 1). the excitatory and inhibitory effects largely cancel each other. Note that this cancellation is not perfect since the weight and spatial properties of the two feedback types are not identical, cf. Table 1.
In conclusion, the tuning curves from Fig 13 show that both temporal low-pass filtering and band-pass filtering can arise from cortical feedback. The detailed spectral shape will depend on both the relative weight and relative delay of the excitatory and inhibitory feedback contributions. While long-delay inhibitory feedback sharpens the temporal frequency tuning of relay cells giving more band-pass-like characteristics, long-delay excitatory feedback makes the tuning more low-pass like. These tuning behaviors can be related to the temporal impulseresponse functions depicted in Figs 11 and 12: multiple phases in the temporal response leads to band-pass filtering, while monophasic responses have low-pass filter characteristics.
Note, finally, that the results in Figs 11 to 13 are all obtained with temporal kernels with a fixed, relatively short, time constant τ of 5 ms, i.e., a relatively short duration of the feedback, cf. Fig 3B. With a longer duration, i.e., larger value of τ, qualitatively similar results for both the impulse-response function and temporal frequency tuning are found, though the curves were found to be more blunt. . Full-field grating is used as stimulus (|k g | % 1deg −1 ) and default values from Table 1 is used for fixed parameters. https://doi.org/10.1371/journal.pcbi.1006156.g013 Effects of cortical feedback on response properties of dLGN relay cells Decorrelation of naturalistic stimuli. In natural visual scenes there are, in addition to extensive spatial correlations, also large temporal correlations [41,89]. It has previously been shown that the biphasic temporal response, seen in retina and dLGN, decorrelates the incoming signal in time, resulting in a more efficient representation of the information in the naturalscene images [41, [90][91][92]. It has also been suggested that the cortical feedback may control the degree of temporal decorrelation in relay cells depending on the signal to noise ratio [37,41,93].
Here we investigate the putative role of cortical feedback on the decorrelation of visual input by calculating the temporal autocorrelation function for relay cell responses to natural movies for different feedback arrangements. The movie was recorded by a camera mounted on the head of a cat exploring the environment [94,95]. The average stimulus autocorrelation function in time and the averaged response autocorrelation function for the relay cells, is shown in Fig 14. The correlation has been calculated for both cases with and without cortical feedback. For the case with cortical feedback the three mixed-feedback scenarios from Fig 13 (bottom row) are considered: (1) long-delay inhibitory feedback combined with short-delay excitatory feedback, (2) long-delay excitatory feedback combined with short-delay inhibitory feedback, and (3) synchronized feedback where excitatory and inhibitory feedback are received at the same time. Fig 14 shows that even without cortical feedback, the correlations in the relay-cell response are significantly lower than the correlations in the stimulus. We observe that the response correlations are further reduced by synchronized cortical feedback, and even more so when the feedback inhibition is delayed. When the excitation is delayed, the response correlations are instead increased compared to the no-feedback case for short time lags. These results thus show that cortical feedback may influence the temporal decorrelation of naturalistic stimuli,  Table 1 have been used for fixed parameters, and the temporal feedback parameters are the same as in Fig 13 (bottom row). Right: Frames from the complex naturalistic movie used as stimulus. This movie was recorded by a camera mounted on the head of a cat exploring the environment (forest) [94,95]. The red circle marks the receptive-field center size for the relay cell at the center.
https://doi.org/10.1371/journal.pcbi. 1006156.g014 and that the degree of decorrelation depends on the spatiotemporal configuration of the feedback. In particular, our mixed feedback configuration seems to be particularly suited for reducing the temporal information redundancy in the signal.

Discussion
In the present work we have developed a firing-rate based simulation tool (pyLGN) to compute spatiotemporal responses of cells in the early visual system to visual stimuli. The simulation tool is based on the extended difference-of-Gaussians (eDOG) model. This model provides closed-form expressions for (Fourier transformed) responses of both dLGN cells and cortical cells, also when cortical feedback projections to dLGN are explicitly included [39]. A main advantage of pyLGN is its computational and conceptual ease. The computation of visual responses corresponds to direct evaluation of two-dimensional or three-dimensional integrals in the case of static or dynamic (i.e., movie) stimuli, respectively, contrasting numerically expensive LGN network simulation based on spiking neurons [31,35,38,96,97] or models where each neuron is represented as individual firing-rate unit [36,37]. This computational simplicity of pyLGN allows, for example, for fast and comprehensive exploration of a wide range of candidate scenarios for the organization of the cortical feedback.

Spatial effects of feedback
As a first example application we focused on the effect of cortical feedback on the spatial response properties of dLGN cells. A specific focus was on so-called area-response curves, i.e., responses to circular spots and patch-gratings as a function of stimulus size, which has received substantial experimental attention. In particular, studies have reported several effects of cortical feedback including sharpening of the receptive field by enhancing the center-surround antagonism of relay cells, increased receptive field center size with removal of feedback, and increased peak response to a an optimal diameter stimulus [17,22,23]. Other studies have reported a more diverse influence from corticothalamic feedback, including both facilitatory and suppressive effects on dLGN cell responses, and changes in the receptive-field structure [18].
Further, several studies have indicated that cortical feedback plays a role in extraclassical suppression in LGN [8,98,99]. A recent study confirms that extraclassical suppression indeed involves extraretinal mechanisms with a delay that is consistent with polysynaptic inhibition [100]. However, the early phase of the suppression is unlikely to involve cortical feedback, as the onset of the suppression is too fast, even though the cortical feedback may contribute in later phases of suppression.
Our model demonstrated that cortical feedback can, depending on the feedback configuration, both enhance and suppress center-surround antagonism and both increase and decrease the receptive-field center size of relay cells. While the receptive-field center size decreases and the center-surround antagonism (as measured by the suppression index α s ) increases with increased (indirect) cortical inhibitory feedback, the opposite is seen for excitatory feedback (Fig 7).
These results support that a phase-reversed arrangement of the cortical feedback, where the ON-ON feedback is inhibitory while the OFF-ON feedback is excitatory, as suggested by data from [73], is more effective to enhance the center-surround antagonism of relay cells as observed in experiments [16,17,22,23,38]. However, with this arrangement a reduction in response to the optimal diameter stimulus in the size tuning curves was observed in our model (Fig 7), in contrast to some experimental studies where an increase in response was reported [17] Here we also considered the more complicated mixed phase-reversed feedback situation with a spatially broad ON-ON inhibitory feedback (combined with a corresponding OFF-ON excitatory feedback) and a spatially narrow ON-ON excitatory feedback (combined with a corresponding OFF-ON inhibitory feedback). Such a center-surround spatial organization of the feedback with excitatory bias to center and an inhibitory bias to the surround has been seen experimentally [3,16]. In our model studies such mixed feedback was seen to give increased center-surround antagonism compared to the situation with ON-ON inhibitory feedback alone. Further, this configuration could also both reduce the size of the optimal stimulus diameter, as well as increase the magnitude of the response to the optimal stimulus diameter (Fig 7). Correspondingly, with this configuration a sharper band-pass property of the spatial-frequency spectra was observed (Fig 8).

Temporal effects of feedback
As for the spatial response properties, the effects of ON-ON inhibitory and ON-ON excitatory feedback (accompanied by the corresponding phase-reversed OFF-ON feedback) are seen to be quite distinct. While delayed inhibitory feedback makes the impulse response more biphasic, the opposite is the case for delayed excitatory feedback (Fig 11). Likewise, while the temporal frequency tuning becomes sharper with delayed inhibitory feedback, it becomes blunter with delayed excitatory feedback (Fig 13, top row).
These features, i.e., increased biphasic index and sharper temporal frequency tuning, are maintained also for the case of mixed cortical feedback as long as the thalamocortical loop delay for the inhibitory feedback is much larger than for the excitatory feedback (Figs 12 and  13, bottom row). Such a mixed-feedback configuration is also found to be particularly suited to remove temporal correlation in the stimulus and thus reduce the temporal redundancy in the neural signals that are sent from dLGN relay cells to cortex (Fig 14).

Spatiotemporal organization of cortical feedback
Our results concerning the spatial and temporal feedback effects suggests that a situation with a mixed organization of cortical feedback consisting of a slow (long-delay) and spatially widespread ON-ON inhibitory feedback, combined with a fast (short-delay) and spatially narrow ON-ON excitatory feedback may have particular advantages. Here the inhibitory and excitatory ON-ON feedback connections are accompanied by excitatory and inhibitory OFF-ON connections following a phase-reversed arrangement [38]. This specific prediction of a mixed organization of feedback could be tested experimentally, for example, by pharmacological inactivation (or other means) of specific feedback connections in the circuit.
This feedback organization seems well suited to dynamically modulate both the center-surround suppression and spatial resolution, for example, to adapt to changing light conditions where the most efficient neural representation of the stimulus is expected to vary depending on the signal-to-noise ratio [40]. In particular, for high light levels (i.e., high signal-to-noise) a band-pass like spatial spectrum (as obtained with our model for certain parameter choices) is expected to provide the most efficient coding, while for low light levels (low signal-to-noise) a low-pass spatial spectrum (as obtained with our model for some other parameter choices) seems better (see Sec. 3.6.1 in [101]) Further, a longer thalamocortical loop time of ON-ON inhibitory feedback compared to ON-ON excitatory feedback assures that temporal correlations in the natural visual stimuli are reduced in the relay-cell responses (Fig 14). This temporal feedback arrangement gives a large biphasic index (Fig 12) which previously has been shown to provide temporal decorrelation of natural stimuli [41], a feature that has also been seen in experiments [42]. An increased biphasic index in the presence of cortical feedback, has also been seen in other modeling studies [37,102]. In particular, a study using a predictive coding model, where the phase-reversed configuration emerged during training, reported that stronger biphasic response in LGN may result from (predictive) cortical feedback interactions [102].
While the slow inhibitory feedback is key for providing temporal decorrelation, the rapid excitatory feedback may have a role in linking stimulus features by synchronizing firing of neighboring relay cells to provide a strong input to cortical target cells [10,19]. Interestingly, a recent study found a large variation in axonal conduction times for corticothalamic axons, from a few milliseconds to many tens of milliseconds [84]. This suggests that differences in feedback delays indeed may have a functional role.
It should be noted that the eDOG-model [39] on which pyLGN is based, assumes that the cortical feedback has a phase-reversed arrangement where each ON-ON feedback connection is accompanied by a phase-reversed OFF-ON feedback connection (Eq (11)), i.e., a push-pull arrangement as experimentally observed in [73]. An alternative is a phase-matched arrangement where relay cells receive feedback only from cortical cells with the same symmetry, including both the direct excitatory feedback and the indirect inhibitory feedback. However, such an arrangement is not only at odds with the observations in [73], but also fails to explain the experimentally observed cortical-feedback induced increase in center-surround antagonism [38].

Application of pyLGN
The aim of present work has been to introduce the simulation tool pyLGN and to provide a variety of example results to illustrate its use. The simulation tool is obviously not limited to the presented example applications, and here we discuss future applications of the tool, as well as its limitations.
Visual stimuli. A key advantage of pyLGN is its computational efficiency and the possibility to use any type of visual stimuli to drive the network. Specifically, natural images (Fig 9) and movies (Fig 14) can be used as stimuli, and an interesting future study would be to compare pyLGN model predictions with experimental data based on natural and other complex stimuli.
X vs. Y channels. A noteworthy feature of early vision is the existence of distinct channels for encoding the visual information [85,103]. In cats, two parallel channels (X and Y) are distinguished based primarily on the spatial summation properties of the cells. Cells in the X pathway show linear spatial summation properties, while cells in the Y pathway show nonlinear summation [30]. As the eDOG model implemented in pyLGN assumes linearity at the geniculate level (even though half-wave rectification was assumed at the cortical level), the present examples were thus geared toward the linear X pathway.
However, the non-linearity of the Y channels is also present in the retinal input [59,78]. A viable strategy for applying pyLGN to Y channels is therefore to apply a suitable non-linear transformation to the response of the retinal ganglion cells driving the dLGN cells in the model. In monkey, the P pathway is known to be fairly linear, while the M pathway may be more non-linear [30]. The same approach may be thus be used to model the monkey M-type dLGN cells.
Adaptation. The response of dLGN cells are known to adapt to changes in stimulus contrast and correlations [104,105]. Such non-linearities at the retinal level can be also incorporated by suitable adaptation of the response of the retinal ganglion cells, for example by using a series of filters [106]. A similar approach can be used to include contrast gain mechanisms originating in the retina [49,107,108]. The eDOG circuit model itself assumes the coupling kernels to be non-adapting, that is, constant over time. Long timescale adaptation effects can still be used in a quasi-stationary manner where the responses during different epochs are computed for fixed connectivity kernel parameters, but where the kernels parameters themselves 'adapt', that is, differ between the different epochs.
Lagged vs. non-lagged cells. A subclass of dLGN relay cells, the so-called lagged cells, has been found to have delayed response onset to visual stimuli and an initial suppression of response until it reaches a maintained firing level [65]. This response profile is contrasting the fast and strong transient response found for the presently modeled non-lagged cells. An application of pyLGN could be to investigate putative advantages of having lagged cells in the dLGN circuit [41].
Bursting vs. tonic firing. Depending on the level of depolarization of the resting membrane potential, dLGN cells can fire in two distinct firing modes: tonic and bursting [2,20,109]. The eDOG model, on which pyLGN is based, assumes the response properties of dLGN cells to be linear, mimicking the firing behavior in the tonic mode [109]. To investigate network behavior when dLGN cells are in burst mode, spiking neuron models including calciumlike slow conductances [35,110] or biophysically detailed models with Hodgkin-Huxley style conductances [38,97] are likely necessary.
ON/OFF asymmetry. The assumption of ON/OFF symmetry is essential for the derivation of the eDOG model, where it is assumed that the input to cortical OFF cells is the negative of the one for the cortical ON cells (see Eq (10)). The exploration of effects of asymmetry in the ON/OFF pathways, for example due to cross-channel mixing as observed in the retina [111,112], would in general require other more complicated network models. However, asymmetry in the ON/OFF pathway can still be studied with pyLGN in a plain feedforward model where cortical feedback is neglected.
Cortical cell model. The eDOG model assumes the response of the cortical feedback cells to be given by a simple half-wave rectifying functions of the input, an approximation most suitable for cortical simple cells [57]. Studies of feedback from cells with other non-linear characteristics such as complex cells can thus not be done in pyLGN.
The choice of half-wave rectifying cortical cells in eDOG assures that the response of dLGN cells is linear, i.e., that the different frequencies of the dLGN response are independent. Thus cortical-feedback induced coupling of neural activity of dLGN cells at different frequencies, such as the slow oscillations induced by visual response in [35], cannot be modeled directly with pyLGN.
Modeling of cortical responses. pyLGN has been designed primarily for modeling of responses of dLGN cells, and in particular effects of cortical feedback. Although not demonstrated in the present examples, it can also be used to explore cortical responses to visual stimuli. Given the computational efficiency of the simulation tool, cortical responses to natural images and movies can readily be evaluated. For the cases when cortical feedback effects are neglected, the cortical output from pyLGN can also be passed through additional static nonlinearities [51].
Descriptive modeling of visual responses. pyLGN has been designed primarily for mechanistic modeling, that is, the modeling of response properties based on the underlying circuitry. However, the tool can also be used to evaluate visual responses for descriptive models, such as the DOG models used to describe responses of cells in retina and dLGN or the Gaborfunction models used to describe cortical cells [51,Ch. 2]. In such models the stimulus is convolved with a (user-defined) descriptive impulse-response function, and thereafter (possibly) passed through a static nonlinearity, to estimate the response. In such applications, custom-made kernels in pyLGN can be used as impulse-response functions and the software interface can be used to efficiently compute the convolution integrals.

Outlook
Studies over the last decades have shown that neurons in dLGN are not simple relays, and play much more important roles in information processing than previously appreciated [19,30,44,100]. (With this in mind, referring to dLGN neurons as 'relay cells' could be misleading, and a renaming could be in order. If so, a more neutral term to use could for instance be 'thalamocortical cells'.) Compared to the primary visual cortex (V1), i.e., the next station in the early visual pathway, the dLGN has received relatively little attention from computational neuroscientists [113]. From a modeling strategy point of view, this is somewhat unfortunate as progress towards a mechanistic understanding of the function of the dLGN circuit seems more attainable given that (i) the dLGN circuit involves much fewer neuron types and is more comprehensively mapped out [1,30], and that (ii) the dLGN has much fewer neurons making simulations computationally less intensive (18000 neurons in dLGN vs. 360000 neurons in V1 in mouse [114,115]). Further, the strong recurrent interactions characteristic for cortical networks (which make them difficult to understand and analyze) appear absent between the principal cells, that is, the relay cells, in the dLGN, even if key circuit network motifs such as feedforward and feedback interactions are present. Thus a focused and comprehensive effort on mechanistic modeling of the dLGN circuit would not only be of interest in itself, it would also likely be a very useful stepping stone for later attempts to model the visual cortex.
While network simulations based on integrate-and-fire type neuron models (e.g., [96,110]) and biophysically-detailed neuron models (e.g., [38,97]) for entire dLGN nuclei are becoming computationally feasible with modern computers, there will still be a need for conceptually and mathematically simpler network simulation tools such as the present pyLGN tool based on the eDOG model. Such models are important to gain intuition about how the different circuit components may affect the overall circuit behavior, and will also be important for guiding the choice of the numerous, typically unknown, parameters in more comprehensive dLGN network simulations. Thus we envision that a future mathematics-based understanding of the dLGN circuit will be of a 'multiscale' nature and be based on a set of interconnected models at different levels of biophysical detail.
The mouse seems particularly suitable as model animal since construction and testing of the multiscale models can be greatly facilitated by the ever more sophisticated techniques for controlling gene expression in mice as well as the possibility for optogenetic activation [114,116]. There are however differences in the way the mouse dLGN is organized compared to monkey and cat, which the current study has been based on, that should be considered [30]. For instance, although there are regions in the mouse dLGN clustered by morphology and ocular dominance, there are no clear laminations (except for a core and a dorsal shell region), and it has been difficult to identify functionally distinct classes of thalamocortical cells [117]. With respect to the model presented in this work, the main challenge is the lack of experimental data on whether the cortical feedback in rodents is organized in a phase-reversed arrangement (ON-ON feedback connections accompanied by OFF-ON connections) as observed in cats [73]. Experimental recordings to characterize this will be of particular interest since several computational studies have based their model on this organization, and in one study where the biphasic responses in LGN was explained using a predictive coding model, this arrangement of feedback emerged during training of the network [37, 93,102].
We envision that at a close and targeted collaboration between modelers and experimentalist, in particular direct assessment of predictions from mechanistic models in targeted experiments, holds great promise for unraveling mechanisms of visual information processing at the different levels of advancements in visual information processing from the earliest sensory systems to more complex computations of the higher cortical areas.
Supporting information S1 Fig. Area response curves for flashing spot stimulation. Top row: predicted area-response curves of relay cells for different arrangements of cortical feedback. Excitatory (left) and inhibitory feedback (center): solid lines and dashed lines correspond to short-delay (5 ms) and long-delay (30 ms) feedback, respectively. Mixed feedback: solid lines correspond to delayed inhibition (D ex RCR ¼ 5 ms, D in RCR ¼ 30 ms), while dashed lines correspond to delayed excitation (D ex RCR ¼ 30 ms, D in RCR ¼ 5 ms). Bottom row: optimal size and suppression index (α s ) are shown as a function of cortical feedback weight for different feedback configurations with different delays. Line styles correspond to different delays as described above. The values on the x-axis represent factors multiplied with the default values for w in RCR and w ex RCR listed in Table 1. Default values for fixed parameters are also listed in this table. (TIF)