A strategy for mapping biophysical to abstract neuronal network models applied to primary visual cortex

A fundamental challenge for the theoretical study of neuronal networks is to make the link between complex biophysical models based directly on experimental data, to progressively simpler mathematical models that allow the derivation of general operating principles. We present a strategy that successively maps a relatively detailed biophysical population model, comprising conductance-based Hodgkin-Huxley type neuron models with connectivity rules derived from anatomical data, to various representations with fewer parameters, finishing with a firing rate network model that permits analysis. We apply this methodology to primary visual cortex of higher mammals, focusing on the functional property of stimulus orientation selectivity of receptive fields of individual neurons. The mapping produces compact expressions for the parameters of the abstract model that clearly identify the impact of specific electrophysiological and anatomical parameters on the analytical results, in particular as manifested by specific functional signatures of visual cortex, including input-output sharpening, conductance invariance, virtual rotation and the tilt after effect. Importantly, qualitative differences between model behaviours point out consequences of various simplifications. The strategy may be applied to other neuronal systems with appropriate modifications.


Introduction
Theoretical modelling of a neural system can be accomplished with different degrees of biological detail and, conversely, different degrees of mathematical abstraction. Arguably, including more biological detail based on experimental data provides a model which is somehow a more faithful surrogate for the original system. But because detail does not necessarily yield understanding, a more abstract description with fewer parameters and thus fewer degrees of freedom offers the possibility of better insight. It follows that mapping between different categories of models is necessary to relate analytical results obtained from more abstract models, to simulations from detailed models which can hardly be analyzed in their parameter space. Eventually, this sort of back and forth process can facilitate understanding about how specific model assumptions may be linked to detailed model behaviour.
In this paper we address the challenge of parameter reduction and quantitative mapping from original biophysical quantities to simplified model parameters, using the orientation hyper-column, or pinwheel, architecture in the visual cortex of higher mammals as a model system. The orientation hyper-column is a functional unit providing organization of cellular tuning with respect to stimulus orientation, and many theoretical approaches have been used to study this neuronal network. At one extreme are simplified population models that retain only gross features of the modeled system, but allow for systematic analytical and numerical investigations, and thus a formal understanding of underlying mechanisms to the extent that the model is reliable. The canonical model in this sense is the firing-rate (FR) ring model with a single neural population type and a one dimensional ring architecture corresponding to preferred orientation [1,2]. At an increased level of complexity, [3,4] also consider populations based on rate models in a 2D distributed architecture corresponding to the cortical surface geometry. Extensions of population models have taken into account different connection profiles [4] or delays [5], both excitatory and inhibitory neuronal populations [3], or more elaborate rate models [6].
At the other extreme are biophysically-detailed models which attempt to incorporate known physiology and synaptic anatomy of the system, giving a coherent description that bridges established parameters and those which are unknown yet essential. For visual cortex there is a wealth of experimental data including various neuron types and their intrinsic biophysics, the functional architecture of thalamic input to the system, the kinetics of the synaptic conductances, and the schematic of the intra-cortical synaptic connections. Numerical simulations of such models, for example [7][8][9][10], have provided realistic reproductions of experimental data, and can provide explicit predictions for subsequent experimental studies. In general, a complex model can incorporate the available data at will, motivated by the fact that a priori we do not know which details are fundamental for function, e.g. the detailed dynamics of receptive fields.
Here we consider a hierarchy of single layer hyper-column models (Table 1 and Fig 1), at the top with a reasonably detailed model that considers biophysical membrane properties of cortical excitatory and inhibitory neurons with synaptic conductances described by second order kinetics. The synaptic architecture of this model arises from a 2D anatomical distribution of the two neuron types over the cortical surface, with spatial interactions defined by isotropic intracortical connections, and thalamic input according to an assumed pinwheel architecture. This full model has been realized with the help of the conductance-based refractory density approach [11]. Via step by step simplifications of the synaptic architecture and cellular elements, comprising five intermediate models, we arrive at the canonical firing rate ring model. We require a consistency between nearest-neighbors within the hierarchy of models to validate each mapping, in the sense that any qualitative or quantitative change in the behaviour of a given model must be linked to specific assumptions made at each stage of the reduction. Note that this particular hierarchy is not comprehensive, e.g. an arbitrarily elaborate network architecture may be evaluated with a simple cell model, or vica-versa. Rather, the motivation for the models studied here is that several are directly linked to previous published results. Furthermore, the methodology that we develop may be readily adapted to other choices of hyper-column models, as well as network architectures of other cortical areas.

Results
We present the results as follows. We begin by describing the hierarchy of network models and the relevant cellular models from simplest to most complex, thus progressing from the abstract current-based firing rate ring model, until the biophysically detailed 2D cortical model that considers Hodgkin-Huxley membrane channels, synaptic currents and conductances. Next, we describe the mapping steps at the cellular and network levels, now in the opposite direction, from the most detailed network and cell models, until the firing rate ring model. Once the mapping strategy and formulas are described, we present an analysis of the parameters of the single population ring models, taking into account the significance of various biophysical values taken from experiments. Numerical simulations of the different network models are then presented, focusing on the similarities and differences in the response steady-state and dynamics to a stereotyped visual stimulus sequence.

Model hierarchy
We consider two geometries of the orientation hyper-column circuit in primary visual cortex, thus a 2-dimensional geometry that maps the cortical surface to function according to the orientation preference hyper-column map found in higher mammals, and a one dimensional ring geometry, where orientation preference corresponds to the position along the ring. Table 1 and Fig 1 present the model hierarchy, emphasizing the basic characteristics of the models in terms of network architecture and cellular properties. The models are assigned letters according to the presumed hierarchy, ranging from the most elaborate 2-dimensional, biophysical, model A at the top, to the simplest version of the ring model, model G, at the bottom. For clarity, in this section we present the models in order from the simplest to the more complex (model G, model F, model E, etc.).

Current-based firing rate (FR) ring model with instantaneous synapses (Model G).
The simplest network that we consider is the classical FR ring model described by [1] (see also [2]). The network consists of a ring in orientation space, parameterized by θ, of a single population of firing rate neurons whose activity is given by the spike rate ν(t, θ): t FR dnðt; yÞ dt ¼ À nðt; yÞ þ n FR SS ðI ring ðt; yÞÞ where τ FR is a phenomenological time constant of the rate relaxation. The steady-state firing rate of the population, n FR SS , is given by a threshold-linear transfer function: The network architecture is fully defined by the expression for the input I ring (t, θ): where I 0 and I 1 are amplitudes of non-tuned (stimulus background) and tuned thalamic input, respectively, and θ 0 is the orientation of the input; J 0 and J 1 are amplitudes of non-tuned and tuned recurrent intracortical input, respectively. These parameters aggregate the effects of spatial connections, synaptic strengths and dynamics between populations of different types. In general the thalamic input is a function of time, thus I 0 (t), I 1 (t) and θ 0 (t). Note that this model incorporates the simplest synapse model, where input current arising from intracortical activity is directly proportional, and thus instantaneous, with respect to the firing rates of the presynaptic population (which also allows the thalamic firing rates to be implicit in this equation).
In this single population model the effect of inhibitory versus excitatory cortical neurons is implicit. The original interpretation of J 0 and J 1 with respect to excitation and inhibition is given in [1]: "The constant J 0 represents a uniform all-to-all inhibition; J 1 measures the amplitude of the orientation-specific part of the interaction. Neurons with similar preferred orientations are more strongly coupled excitatorally than neurons with dissimilar ones." As originally formulated by [1] (see also [2]), the values of I 0 and I 1 were chosen so that the feedforward input is non-negative, since there is no experimental evidence for inhibitory thalamocortical pathways [12]. In Section 2.2.4 we will present an alternative interpretation of these coefficients that arises from the mapping, specifically that takes into account the implicit inhibitory cortical population of this model.

Conductance-based shunt FR and leaky integrate and fire (LIF) ring models with instantaneous synapses (Models E, F).
We now elaborate the single population currentbased ring by incorporating synaptic conductances as well as currents, and by including a membrane current noise term. We consider firing rate and spiking versions of the ring, based on a shunt firing rate neuron (model F), and a leaky integrate and fire neuron (model E), respectively.
The conductance-based LIF ring model E explicitly considers the effect of spikes on the membrane potential, V. The current equation for a single compartment LIF neuron at location θ along the ring is given by: where C is the cell capacitance, and g L is the leak conductance. Synaptic input is defined as the total synaptic conductance, S(t, θ), and the total synaptic current, I(t, θ), the latter measured with the cell held at the resting potential (note that a current term representing external stimulation may be included without loss of generality). ξ(t) is a Gaussian white current noise process characterized by its mean, hξ(t)i = 0, and auto-correlation, hxðtÞxðt 0 Þi ¼ s 2 I ðC=g L Þ dðt À t 0 Þ, and σ I is the noise amplitude. When V crosses a defined threshold voltage V Th , the neuron "fires" and the membrane voltage is immediately reset to a defined voltage V Reset . For simplicity, in particular to allow a direct mapping to a firing rate model, we neglect an absolute refractory period under the assumption that the evoked firing rates are well below maximal.
We note that a noisy population of LIF neurons with instantaneous synapses can be evaluated either with Monte-Carlo methods, or with a probabilistic approach based on the Komogorov-Fokker-Planck (KFP) equation. We use the latter approach, described in the Methods (Section 4.1), to derive the simulations presented in the Results. Thus, in the presentation below of the input terms for model E we will refer to specific formulae of the KFP approach.
The average steady-state firing rate of the LIF neuron defined by Eq 4 provides a hybrid "shunt FR" network model (model F), situated between the previously described currentbased FR model G, and the conductance-based LIF model E, that takes into account noise as well as the impact of synaptic input on the membrane time constant. For a given input I and S, the average steady-state firing rate of the LIF model, n LIF SS , is given analytically by [13] (see also Section 4.2): where the asymptotic potential, x(I, S), the steady-state voltage dispersion, σ V (S), and the effective membrane time constant of the LIF neuron, τ m (S), all depend on the synaptic conductance S: The shunt FR model F is then defined by using Eq 1 of the current-based FR model G, now with the steady-state solution given by n LIF SS of the LIF model E: t FR dnðt; yÞ dt ¼ À nðt; yÞ þ n LIF SS ðIðt; yÞ; Sðt; yÞÞ ð8Þ The form of the synaptic current and conductance network connections of models E and F are identical; we assume they recapitulate the constant plus cosine form of the current-based FR ring model G (Eq 3) and behave identically as a function of θ: Iðt; yÞ ¼Ĩ 0 ðtÞ þĨ 1 ðtÞ cos ðy À y 0 ðtÞÞ þ 1 2p  B, C, D). We next consider two population ring models that include explicit excitatory (E) and inhibitory (I) neuron populations, based either on the LIF cell model, or more complex Hodgkin-Huxley (HH) type cell models. The two population LIF-based ring models include a single compartment inhibitory cell model, with either a one compartment LIF excitatory cell model (model D), or a two compartment LIF excitatory cell model (model C). Note that by definition, the standard LIF model as used here does not exhibit adaptation. At a next level beyond the LIF description, HH-type models are the standard framework for describing cellular dynamics tied to explicit biophysical mechanisms, including voltage-dependent membrane channels [14]. The HH-based ring model (model B) includes an adaptive two compartment excitatory cell model and a non-adaptive single compartment inhibitory cell model. The intracortical weighting functions for the excitatory and inhibitory pathways in all three of these two population ring models, fitted to anatomical data, follow the constant plus cosine form presented above.
The two population models allow distinct cellular and synaptic properties for excitation and inhibition, which in turn motivates the consideration of more elaborate synaptic kinetics. Each cortical population, therefore, has three types of synapses: those arising from thalamic excitatory (Θ), from intracortical excitatory (E), and from intracortical inhibitory (I) pathways. Accordingly, each connection in the network is denoted by a double index ij where i and j indicate the pre-and post-synaptic population, respectively. The complete synaptic input to a target population of type j is then given by the total synaptic current measured at rest, I ij , and the total synaptic conductance, S ij , summed over the different synapse types arising from presynaptic populations i.
The excitatory and inhibitory synaptic currents measured at rest for each cell type j at location θ, and the associated synaptic conductances, are thus: where V E and V I are the excitatory and inhibitory reversal potentials, and V Rest j is the resting potential of the target population.
For the one compartment LIF and HH cell models (both cell types in model D; the inhibitory cell type in models B and C), we refer to the total input, thus: I j ðt; yÞ ¼ I Ej ðt; yÞ þ I Ij ðt; yÞ S j ðt; yÞ ¼ S Ej ðt; yÞ þ S Ij ðt; yÞ The membrane current equation for the single compartment LIF and HH models of cell type j at location θ is then: where I Act is the current via active, voltage-gated channels for HH neurons, and ignored for LIF cells.
For the two compartment HH and LIF excitatory cells of models B and C, the post-synaptic sites are pathway specific. Thus we refer directly to Eq 11 for the excitatory input exclusively on the dendrite, I EE (t, θ) and S EE (t, θ), and for the inhibitory input exclusively on the soma, I IE (t, θ) and S IE (t, θ). The two compartment models are an approximate solution of two boundary problems for a spatially distributed passive dendrite and a point passive soma, following the description in [15] (see also [9] and [16]). The problems correspond to the two cases of, first, dendritic synaptic current estimation under somatic voltage clamp (the reverse voltageclamp problem) and, second, somatic voltage measurement under current-clamp in response to the previous estimates of the dendritic synaptic current.
The membrane current equations for the soma voltage V(t) and the dendrite voltage V D (t) for this excitatory cell model at location θ are: where t 0 m is the resting time constant (¼ C=g 0 tot , where g 0 tot is the total conductance at rest), ρ D is the ratio of the dendritic to somatic passive membrane conductances, and L is the dendritic length in units of the characteristic length (see Methods, Section 4.3).
Notably, spike generation in the HH cell models used here is governed by an explicit threshold mechanism, reminiscent of the LIF model, that replaces the spike generating current I Na . This modification, described in detail in the Methods (Section 4.2) as the conductance-based refractory density (CBRD) approach, allows for efficient evaluation of neuron populations that still considers the impact of the various HH currents that comprise I Act (Eqs 13 and 14). It is important to note that this approach was taken for its computational advantage alone; using a standard HH cell model in the framework of the network simulations will give very similar results [17] (see also Section 4.2).
The connectivity rules for all the two population models are written in terms of pre-synaptic firing rates, φ ij , to facilitate incorporating synaptic kinetics in the translation of pre-synaptic activity to the post-synaptic conductances g ij (t) that figure in Eq 11. The relation between φ ij (t) and g ij (t) is given by a second order kinetic model that implicitly accounts for mechanisms such as axonal delays, receptor dynamics and dendritic integration: with the scaling "synaptic capacity" c ij described in Section 4.5. The thalamic inputs driving g Θj (t, θ) of the excitatory and inhibitory cortical populations are taken to be identical, and follow the constant plus cosine form of the single population ring models: where φ 0 is the un-tuned stimulation background, φ 1 is the tuned part of the cortical input, and θ 0 is the stimulus orientation angle. As mentioned, the thalamocortical pathway is exclusively excitatory. Intra-cortical connections driving g ij (t, θ) of the two populations are given by: where the weight function, w ij (ϕ), again recapitulates the constant plus cosine form of the single population models, with the parameter q ij to take into account the orientation tuning between different types of neurons: Note that, implicitly, 0 � q ij � 1. This approximation will be developed in Section 2.2.1.

Two population, two dimensional HH model with the second order synapse kinetics (Model A).
As the most elaborate model, we consider an anatomical 2D cortical network of the two compartment excitatory, and the single compartment inhibitory, HH cell models developed in the previous section. As for the two population ring models, the current and conductance inputs are driven by pre-synaptic firing rates.
As opposed to the angular functional index θ of the ring models, in the 2D model each population is parameterized by its position on a cortical surface map of adjacent, radially-symmetric hyper-column pinwheels [3] (Fig 2). The 2D architecture is captured entirely by new expressions for the firing rates φ Θ and φ ij that are now functions of the x and y coordinates of the target cells. Thus, the inputs of the 2D model recapitulate Eqs 11-15, with arguments (t, θ) of the ring models replaced by (t, x, y). An orientation column (Fig 2) corresponds to a hyper-column sector, containing populations with similar preferred orientations at varying distances from the hyper-column center. The functional architecture of each hyper-column pinwheel arises from a thalamic input whose orientation tuning corresponds to the angle of the target location with respect to the center of its hyper-column [8,3]. Hyper-columns are arranged so that those with a clockwise progression of orientation columns are adjacent to those with counterclockwise progression, thus adjacent orientation columns from different pinwheels have the same orientation preference. The centers of the hyper-column pinwheels are distributed on a rectangular grid with pinwheel radius R and indexed by i PW and j PW . Here we consider a square containing 4 pinwheels on the cortical surface. The coordinates of the pinwheel-centers are x PW = (2i PW − 1)R, y PW = (2j PW − 1)R. The orientation angle for the point (x, y) which belongs to the pinwheel (i PW , j PW ) is defined as θ = arctan((y − y PW )/(x − x PW )). The progression is determined by the factor ðÀ 1Þ i PW þj PW .
Thus, the thalamic input in terms of effective firing rate is described with an elaboration of Eq 16: It is important to note that, as with the ring models, this constant plus cosine form of the thalamic input constitutes an implicit model of the retino-geniculo-cortical pathway; this is equivalent to each cortical neuron receiving an oriented receptive field corresponding to a single sub-region of a classical simple cell. The orientation dependence of the input in particular characterizes the degree of symmetry breaking at the input stage, parameterized here by φ 1 and φ 0 . The homogeneous case, where φ 1 = 0, represents thalamic input with no dependence on stimulus orientation. The maximum selectivity of this model for the input stage, where φ 1 = φ 0 , is a pure shifted cosine dependence. Intra-cortical connections between populations are local and isotropic with an exponential decay [12], with the effective pre-synaptic firing rate φ ij (t, x, y) given by: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À x 0 Þ 2 þ ðy À y 0 Þ 2 q =d ij Þ dx 0 dy 0 R R exp ðÀ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where, as above, i and j correspond to the cell type (E or I), ν i (t, x 0 , y 0 ) is the firing rate of presynaptic cortical population i at location (x 0 , y 0 ), and d ij is the characteristic anatomic length of the i to j pathway. Note that in contrast to the explicit tuning of intracortical connections in the ring models (parameterized by θ), these connections have an implicit functional orientation tuning since neighboring orientation columns receive similar oriented thalamic inputs.

Mapping of network and cellular parameters
We now describe the quantitative mappings between the various models at the network and cellular levels, which allow a path from the parameters of the full HH EI 2D model A, to the different LIF networks models, and finally to the current-based FR ring model G. In contrast to the previous section, we present the mappings from the more complex to the less complex models. The essential goal of these cellular and network mapping steps, as diagrammed in Fig  1, is to systematically map the high-dimensional full model, with multiple non-linearities, to a low-dimensional model with a single threshold non-linearity.

2.2.1
Network mapping: 2D cortical architecture to the cosine ring models. The first mapping step at the network level is to reduce the 2D geometry of the pinwheel architecture in anatomic space to a 1D ring geometry in orientation space (solid green arrows in Fig 1). We consider the synaptic strength function for the intracortical connection between two pointsr 1 andr 2 on the cortical surface (corresponding to the locations (x, y) and (x 0 , y 0 ) in Eq 19): f w ðjr 1 Àr 2 jÞ ¼ exp ðÀ jr 1 Àr 2 j=dÞ where d is the characteristic length of the connection. The weight of the connections of two elementary areas ΔS 1 = r 1 ΔθΔr 1 and ΔS 2 = r 2 ΔξΔr 2 at the points shown in Fig 3A is: Dx are the areas of the sectors corresponding to two radial vectors emanating from the center of a pinwheel with radius R and having the orientation angles θ and ξ (Fig 3A). Combining the last two expressions, we integrate over the vectors to obtain the total connection weightw: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðr 1 cos y À r 2 cos xÞ 2 þ ðr 1 sin y À r 2 sin xÞ 2 q � � We now findw for the angle difference α = ξ − θ. By rotating the coordinates such that θ = 0 and α = ξ, we obtain: The resulting normalized weight is shown in Fig 3B by  . We also considered an exponential form w ij (ϕ) = 1 + q ij exp(|ϕ|/ϕ ij ) ( Fig  3B, dashed line), which gives a more precise mapping of the anatomically based 2D cortical geometry [18] to the ring. This relationship between the exponential profile in 2D plane and the constant-plus-exponential profile in the ring of orientations is novel.

Cellular mappings: HH to LIF neurons, and 2 compartment to 1 compartment LIF excitatory neuron.
The mappings of the HH neuron models to their respective LIF models were made as follows. For the population of one compartment HH inhibitory neurons, the steady-state transfer function between current input and firing rate (f/I) of a one compartment LIF model, with the same passive parameters of the HH neuron, was fitted to the f/I curve of the HH model by adjusting the LIF reset potential, while retaining the threshold potential V Th of the HH inhibitory model (dashed purple arrows in Fig 1). A similar fitting procedure was done for the two compartment HH excitatory neurons, but importantly using a non-adapting version of this model. Thus, as a first step, the two adaptation currents of the excitatory neuron, I M (intermediate time-scale voltage-dependent current) and I AHP (slow voltage and calcium-dependent current), were blocked, and the f/I was obtained. The f/I curve of a two compartment LIF model was then fitted to the non-adapting HH f/I curve, with the same passive parameters and threshold potential V Th of the HH excitatory model (dashed green arrow in Fig 1). The comparison of the spike timings in response to constant input is shown in Fig 4. To map the two compartment LIF excitatory cell model to single compartment LIF model (dashed cyan arrow in Fig 1), we fixed the passive input conductance of the latter to the somatic input conductance of the former, and decreased the membrane time constant from 14.4 to 10.3 ms, reflecting the fact that the effective time constant of a two compartment model is smaller than intrinsic membrane τ m . Because of different expressions for the input conductance in two and single compartment models, expressed by the ratio of the dendritic and somatic input conductances, ρ D = G in,d /G in (Eqs (63) in Methods, Section 4.3), the capacitance of the one compartment model changed slightly from 0.25 to 0.27nF. The synaptic input to the dendrite of the two compartment model (excitatory inputs from the thalamus and the  excitatory cortical population) were adjusted for the one compartment model by ρ D so that the sizes of the synaptic conductances relative to the post-synaptic input conductance were maintained (ref. Fig 4A for a comparison of the spike times of the two models).

Cellular mapping: 2-dimensional LIF steady-state to 1-dimensional threshold-linear f/I relation.
We now describe the mapping of the steady-state expression for the current (I) and conductance (S) driven LIF model, n LIF SS ðI; SÞ (Eq 5), to a one dimensional threshold-linear model, thus with the current and conductance inputs combined a priori (dashed orange arrows in Fig 1). This transformation is used twice. First, this is used for eliminating an explicit inhibitory population when going from the two population LIF ring model D, to the conductance-based single population models E and F (Section 2.2.4). Second, this step is used for mapping from the conductance-based shunt FR model F, to the classical current-based FR model G (Section 2.2.5).
We first make a non-standard assumption that the dispersion of noisy current, s 2 I , is linearly dependent on the synaptic conductance: with an only free coefficient, which is the noise level at rest σ I0 . This approximation reflects the fact that the synaptic current fluctuations typically increase as more synaptic channels are activated. The ionic channels provide independent additive stochastic currents, thus the variation of the total current fluctuations is the sum of variations of the currents from separate channels, which is roughly proportional to the number of activated channels. The later is characterised by the total conductance. Importantly, for the LIF neuron this leads to constant voltage dispersion due to noise, σ V (S) (Eq 6): Independence of the voltage dispersion on the total conductance corresponds to the case of thermal noise in RC circuits [19]. The scaling Eq 23 and, consequently, the constant σ V were also noted by Chance et al. [20], who considered a balanced excitatory and inhibitory synaptic input with the synaptic conductance dominated over the leak conductance. The proportionality of s 2 I on S holds true for the noise generated by a Poisson process. So either an addition of the free term s I 2 0 to this proportionality in Eq 23 or the assumption of synaptic conductance domination, S >>g L , as in [20], lead to constant s 2 V . Given the form of Eq (5), this in turn allows a recasting of the steady-state solution n LIF SS as a one dimensional function: We note that this reduction of a 2 to 1 dimensional solution is unique for this noise model, and is an important analytical property of the noisy LIF model, obtained due to the novel assumption Eq 23. In the notations from [20], the formula of Eq 24 corresponds to the relationship r = g f(I/g) with the total conductance g, the input current I and the firing rate r. With respect to the modulation of the f/I curve gain, we see that the shunt g provides a divisive effect on the gain f 0 (I/g) only if the second derivative of f is positive. Next, we approximate n LIF SS ðI; SÞ for S = 0 as a one dimensional threshold-linear function, parameterized by the slope k LIF and the rheobase current I LIF Off : The approximation parameters for the LIF models introduced in the previous Section 2.2.2, thus k LIF E and I LIF Off;E for the excitatory LIF neuron, k LIF I and I LIF Off;I for the inhibitory LIF neuron, are obtained by fitting to n LIF SS;E ðI; 0Þ and n LIF SS;I ðI; 0Þ, respectively. The precision of this approximation for the excitatory cell type can be seen from Fig 5, which compares the threshold-linear approximations for noisy LIF and HH cell models plotted for parameters given in Section 2.3.
Combining the last two expressions obtains an approximate one dimensional threshold-linear dependence of the rate on I and S: A similar to Eq 26 approximation was described in [6], in the context of fitting thresholdlinear functions to the numerically obtained steady-state transfer function of a HH neuron with a slow potassium A-current, specifically approximating the impact of the cell conductance as a rightward shift in the transfer function. Brizzi et.al. [21] showed that this approximation was justified for spinal motor neurons in vivo, and Persi et al [22] generalized the approach for HH neuron models with current noise. In contrast, our novel interpretation of Eq 26 arises from the application of the similarity law (Eq 24), which is exact for the particular noise model (Eq 23) and the LIF neuron.
Thus, the effect of the synaptic conductance is to modulate the threshold current of the LIF neuron I LIF Th : Note, in the aspect of the f/I curve gain problem, we see from Eq 26 that the channel-to-channel independent noise results in shunt-independent gain of f/I curve, i.e. the shunt provides a pure subtractive effect. As shown in the Section 2.3, the quality of these approximations allows to obtain consistent mapping of the model parameters. Note that the choice of a threshold-linear transfer function approximation for the steady-state is not essential for the mapping from the two population to

PLOS COMPUTATIONAL BIOLOGY
one-population models, nor between the conductance-based to current-based FR models, and thus a more elaborate one dimensional function may be used. However, the threshold-linear approximation not only simplifies these mappings (e.g. solving for ν I in Section 2.2.4), but it is also necessary for the final mapping to the classical threshold-linear FR model that is especially amenable to analytical results.
Considering the form of the synaptic current and conductance introduced previously, we can simplify the mapping further. Now we explicitly indicate the LIF cell type j. First, we express the threshold voltage of the LIF neuron in terms of its rheobase current I LIF Off;j and input conductance g L,j : We now define DV Th ij as the driving force relative to the target's threshold potential: where V Rev i is the reversal potential of synapse type i. With this definition of DV Th ij , we introduce the network input, recalling the definitions of I j and S j (Eqs 11 and 12), and finally rewrite Eq 26 for n LIF SS;j : 2.2.4 Network mapping: Conductance-based two population ring to single population ring. We now map the current and conductance inputs I j (t, θ) and S j (t, θ) (Eq 12) of the two population ring model D, to the I(t, θ) and S(t, θ) inputs (Eqs 9 and 10) of the single population LIF and FR with shunt ring models (models E and F; ref. solid red arrows in Fig 1). This mapping effectively eliminates the explicit inhibitory population, replacing it by an instantaneous linear function of the excitatory activity. Note that the cell models in models E and F are derived from the single compartment excitatory LIF model of model D.
The first step is to replace the 2 nd -order synaptic kinetics by instantaneous transfer functions, thus synaptic currents and conductances being directly proportional to the effective presynaptic rates, according to the steady-state solution for Eq 15, which is: Now, the time dependence of the input is implicit: Next, we expand Eq 31 for the current input to excitatory neurons, I E (θ), using the pre-synaptic firing rates of the thalamic (φ Θ , Eq 16) and intracortical (φ EE and ϕ IE , Eq 17) populations: We now replace the inhibitory population rate ν I in this equation by the steady-state equation for the inhibitory LIF model expressed as a linear function of the excitatory population rate ν E , based on three assumptions. First, we assume that the connections between inhibitory neurons in the cortex are local, thus: Next, we ignore the intrinsic dynamics of the inhibitory population and consider only its steady-state response, thus n LIF SS;I ðyÞ, furthermore with the assumption that, on average, the inhibitory neurons around the ring are above threshold during stimulus presentation. Using the one dimensional approximation given by Eq 29 for the LIF inhibitory neuron (importantly neglecting the rectification operation because of the above-threshold assumption), Eq 31 for the current term (subsituting DV Th ij according to Eq 29), and the assumption of inhibitoryinhibitory locality (Eq 33), we obtain: Eq (34) may then be solved for ν I (θ): where:k The termk LIF I represents the reduction in the effective gain of the inhibitory population because of self-feedback, as such parameterized by the inhibitory-inhibitory pathway parameter c II . Assuming that the reversal potential of the inhibitory synapses is beneath the spike threshold, we note that DV Th I I (Eq 28) must be negative, and thereforek LIF I � k LIF I . Finally, we need to express the effective pre-synaptic rate of the excitatory input to the inhibitory population, φ EI (θ), in Eq 35 as a linear function of ν E (θ). For this we assume that the excitatory activity over the ring can be reasonably approximated by a shifted cosine centered at the stimulus orientation θ 0 : We can now re-write Eq 32 for the current I E (θ) as a function of only the excitatory population activity ν E (θ), substituting Eqs 16 and 38 for φ Θ and ϕ EI , respectively: Again using the relation of Eq 37, we obtain the coefficientsĨ 0;1 andJ 0;1 by comparing Eq 39 with the expression for I(θ) (ref. Eq 9, i.e. setting I(θ) = I E (θ), with implicit time dependence): These equations provide the current terms in Eq 9, with the corresponding conductance terms in Eq 10 given by: Note that the C IE term in the above expression forĨ 0 is not, in fact, driven by thalamic input (i.e. the thalamic firing rate term φ 0 ), but rather is constant. Nevertheless, the formats of the defining equations for the inputs (Eqs 9 and 10) that distinguish between terms integrated with the network activity, from those that are "external", oblige the inclusion of the C IE in the "thalamic" component. In comparison to the thalamic terms, the input "driven" by intracortical activity is proportional to bothJ 0;1 andL 0;1 , i.e. there is no non-driven, constant component.
To conclude, Eqs 41 and 42 constitute the mapping between the conductance-based two population ring models with 2 nd -order synapses (models B, C and D) to the conductance-based single population ring models with instantaneous synapses (LIF model E, and shunt FR model F).
2.2.5 Network mapping: Conductance-based single population ring to current-based FR ring. The final mapping is between the current and conductance inputs I(t, θ) and S(t, θ) of the conductance-based ring models E and F, to the current input I(t, θ) (Eq 3) of the currentbased FR ring model G (dashed blue arrows in Fig 1).
As before we exploit the mapping of the two-dimensional steady-state expression for the LIF neuron, n LIF SS ðI; SÞ, to the one dimensional threshold-linear approximation (Eq 26). Here, this step allows the elimination of the explicit conductance input, therefore to obtain the stationary solution n FR SS ðIðt; yÞÞ of the current-based FR model (Eq 1 Note that, unlike the previous use of the one-dimensional transfer function approximation for eliminating an explicit inhibitory population, here there is no assumption that the input is suprathreshold. Comparing Eq 3 for I(t, θ) with Eq 43, we obtain: with the I and S terms on the right-hand side defined by Eqs 41 and 42. These relationships complete the mapping to the current-based FR ring model G, and thus constitutes the main analytical result of the work.

Simulations
We now present comparative simulations of the different models described in Section 2.1 and responding to to the same input. These simulations demonstrate the consistency of the mapping expressions, show how different behaviours of the models reveal implications of each step of the reduction, and highlight essential mechanisms for certain signatures of responses in the visual cortex.
We first set the parameters of the most complex model with excitatory and inhibitory populations distributed in 2D cortical space (model A), according to the experimental literature on single neuron properties, synaptic kinetics and connectivity (Section 4.6). The parameters of the remaining models (Section 4.6) are then calculated according to the formulas derived in Section 2.2.

Overall transient dynamics.
We compared responses of the models to a specific visual stimulus, focusing on the excitatory population firing rate. The initial conditions were the steady-state in response to a homogeneous background input of "contrast" or strength φ 0 . The visual stimulus of contrast φ 1 is presented at t = 0ms with a time-dependent orientation (e.g. oriented bar, grating, edge, or more complicated scene) of 0˚at t = 0ms, shifting to 45˚at t = 100ms. For model A, we simulated the response of 4 identical hyper-columns (Fig 2), with an identical input for each. This stimulus provides the fundamental "step response" of the system from the background state, followed by a step response in the orientation domain, for example following an abrupt change in the visual scene during a saccade.
We now consider the overall dynamics of the responses of the 2D model A. The spatial distributions of the excitatory population firing rate at different times are shown in Fig 6. Initially, the modeled cortical domain responds with a firing activity in the orientation columns that correspond to maximal stimulus (frames at 5-20 ms). The maximum response appears at about 20 ms; its location corrsponding to the stimulus orientation. The activity is homogeneously distributed along a radius of a hypercolumn. At later stages, the firing rate decreases (frames at 30-50 ms) with a peak located near the pinwheel center (frames at 50-90 ms). The input changes its orientation at 100 ms from 0˚to 45˚. The activity rapidly rotates toward the column preferring the orientation of the stimulus (frames at 105-110 ms), and even beyond (frames at 120-130 ms; see also Section 2.3.5). Later, the activity settles again in the column preferring the orientation of the stimulus with a peak near the pinwheel center (frames at 150-190 ms), similar to the activity in response to initial stimulus (compare to frames at 50-90 ms). The time-dependent profiles of the rate and voltage for the populations on a ring of radius 2/3 R centered on a pinwheel (marked as a circle in Fig 6), thus as a function of preferred orientation, are shown in  Fig 7D and 7E, respectively. Importantly, the steady-state location of the maximum activity in each hypercolumn corresponds to the orientation of the stimulus, and the firing is localized in 2D and orientation coordinates, i.e. the cortical response is tuned.
All of the constructed hierarchy of models are compared in Fig 8. Common and different features of the models' portraits constitute an essentially novel result of the work. The population rate as function of orientation and time is shown. When considered along a ring centered on a pinwheel, the dynamics of the 2D adaptive, two population model A, is similar to the adaptive ring geometry model B, although the full dynamics of the 2D model depend on the radius of the ring, as will be discussed below. This similarity is achieved due to the mapping performed in Section 2.2.1. Both exponential and cos-shaped profiles of connection spread on a ring (Fig 3B) result in similar solutions, one of which seen in Fig 8 for model B. In contrast, the transient responses change significantly between the adaptive HH models A and B, and the non-adaptive LIF model C, primarily due to an overall increase in the activity in the latter model. Damped oscillations are observed to varying degrees for models A-D, being strongest  for the adaptive models, and quite weak for the 1 compartment LIF EI model D. Thus, steadystate-like behaviour in the adaptive models A and B is seen only after 200ms. The non-adaptive models C and D show steady-state-like behaviour within approximately 100ms, thus observable for both stimulus orientations.
The reduction step from the 2-compartment LIF neurons of model C to the 1-compartment LIF neurons of model D, which reduces the order of the system of equations, has a more striking effect on the transient response. In particular, as noted above, oscillations following stimulus transitions are substantially reduced for the 1-compartment LIF model. A larger propensity to oscillations of the higher order model C in comparison to the lower order model D is consistent with the observation of a faster voltage response of a spatially distributed neuron to changes in synaptic current due to the electrotonic properties of the implicit cable structure. Following the change in the stimulus at 100ms, both models show a continuous shift of the response in orientation space, or a virtual rotation, as discussed below.
The subsequent model E incorporates two more simplifications by reducing the number of populations from two to one, and assuming instantaneous synaptic kinetics. These lead to quite different transient behavior, with a much more rapid establishment of the steady states, and the disappearance of oscillations and virtual rotation. Why is the response of model E much more rapid than the previous in models, even though the essential cellular time constant is similar, i.e. given by the membrane time constant which scales the rate of the voltage

PLOS COMPUTATIONAL BIOLOGY
integration in response to a change in the input current? To explain, we note that the firing rate of model E effectively depends on both the voltage relative to threshold as well as the change in the voltage [23]. The result is a rapid response of the population to an abruptly changing stimulus [24,25] and explains the near-instantaneous response in our simulations.
Further reduction to the FR-model with shunt, model F, with its time constant equal to the membrane time constant, leads to the delays of the responses. The reduction implies the substitution of the statistically precise consideration of neuronal states in each population by a statistically approximate approach. Namely, the Fokker-Planck-based and CBRD-based methods of models A-D (see Methods) are substituted by the FR-model that is strictly valid only at steady-states of a population. For nonstationary regimes the FR-model provides an approximation that filters the firing rate with some characteristic time constant of the filtering. By chance, given the various parameters of the models, the time constant of model F is close to the synaptic time constants, thus the transient solution is again as smooth as the solutions of models A-D. The small residual difference between the responses of the FR-ring models with and without shunt (models F and G) is explained by the approximation of the dependence of the rate on the current and conductance by the piecewise linear function (Eq 26).

Sharpening of output tuning vs. input tuning.
Our analysis continues with the steady-state response of the networks, taken as the activity profiles at t = 300 ms, in terms of the sharpening between input and output, and of contrast-invariance of the output tuning curves.
A fundamental aspect of the orientation selectivity in primary visual cortex is that this functional property is essentially absent in individual inputs from the thalamus. Various theories have been proposed to account for the emergence of tuning in cortical neurons, but a common assumption is that there is an initial bias in the spatial retinotopic footprint of the feedforward pathways of cortical neurons, and this bias is strengthened in terms of a sharper orientation tuning in the cortical receptive fields ( [26,27,28]). In the context of the models described here, this bias is encapsulated by the constant plus cosine form of the thalamic input.
We consider this input-output sharpening by examining the steady-state solutions of the models (Fig 9). All the models show orientation tuning sharper than the input tuning; furthermore, the sharpening increases as the models simplify.  Table 1, obtained as crosssections at t = 300ms of the plots shown in Fig 8 (un-scaled in A, and normalized in B). The profile with the smallest amplitude corresponds to the simulation by the "adaptive" firing-rate ring model G as described in the text and indicated in Fig 10. The half-widths at half-maximums are 38˚for model A, 36˚for model B, 31˚for model C, 28˚for model D, 27˚for models E and F, 24˚for model G, and 22˚for model G with the gain rescaled to account for adaptation.
The canonical FR ring model G has provided a semi-analytic framework for accounting for this sharpening in terms of the strength and tuning of intracortical connectivity. In this model the non-tuned inhibition is effectively stronger as J 0 decreases, whereas the tuned excitation is stronger as J 1 increases. As presented above, the derived values for J 0 and J 1 correspond to the "marginal" domain (Fig 10) which in turn predicts enhanced output sharpening. In contrast, while the analysis described above suggests that the 2D HH model A is operating in the homogeneous domain, cell thresholds and tuned "mexican hat" inhibitory and excitatory populations (Fig 7) still allow for input-output sharpening. The prediction that adaptation underlies a shift from the marginal to the homogeneous regime of model A during the response suggests that sharpening of the phasic response should be greater than the tonic response, as reported by measurements of orientation-tuning dynamics in the macaque [29].
The impact of spike threshold is also seen in the simulations. When comparing subthreshold voltage and spikes, experimental investigations have reported that the tuning of the former is broader than the latter [30,31,32,28,33] et al., which can be explained simply by firing threshold. This behaviour is replicated in the the HH and LIF models which explicitly consider  adapted from [2]). Amplitude instability corresponds to the state where the activity in the ring increases without any possibility to regulate it. Homogeneous phase ("feedforward" hypothesis) corresponds to a state of weak interactions. The activity in the ring follows directly from the thalamic input, apart from a threshold nonlinearity. Marginal phase ("recurrent" hypothesis) corresponds to a state where only a tuned activity profile is stable, partially but not completely determined by the input shape and dynamics. This state occurs for sufficiently strong recurrent tuned inputs (J 1 ) and, to a lesser extent, with sufficiently strong inhibition (J 0 ). The small square marks the parameters found by the mapping expressions for model G ( Note that the 2D and ring HH adaptive models A and B give small amplitude steady-state profiles in orientation space (Fig 9A) that are qualitatively consistent with the solution for the "adaptive" firing-rate ring model described above. However, quantitative comparison in this case, especially for output sharpening, is problematic, mainly because for the adaptive neuron the rate-current-conductance function is no longer scaled by Eq 24, and is not approximated by Eq 29.

Contrast invariance.
Contrast invariance, which is the qualitative maintenance of the output tuning sharpening with respect to the input tuning over a range of input strengths, is another classical property of neurons in primary visual cortex described by experimental and theoretical studies [34]. Attractor dynamics of the neuronal network have been proposed as the underlying mechanism for this functional property; thus the marginal domain of the canonical FR ring-model model G predicts contrast-invariance, and this property is conserved for the LIF-based ring-models E, D and C. In comparison, the fact that the adaptive version of the G model is in the homogeneous domain suggests that the adaptive models A and B will not show contrast-invariance. Fig 11 confirms these predictions for models A and D.
The models predict therefore that firing adaptation of excitatory neurons is sufficient to cause a dynamic qualitative evolution of the visual response, specifically that contrast-

Fig 11. Contrast invariance for the excitatory population firing rate (red) in models D (left) and A (right), compared with the responses (green) of feed-forward versions of each model (� g EE ¼ �
g IE ¼ 0) to distinguish the contribution of intracortical pathways. The "weak" and "high" contrast values were 0.4 and 0.8, respectively. As described in the text, the higher gain of model D as compared to model A, as seen in the different response magnitudes, account for the marginal phase behavior of the former, thus contrast-invariance, and the homogeneous phase behavior of the latter, thus no contrast-invariance. The intracortical connections serve to increase inputoutput sharpening for both models. invariance should be observed essentially during the transient phase of visual reponses, being reduced or absent during adapted (long term) responses. Similarly for input-output sharpening, the experimental prediction is thus that the phasic response provides the basis for measured contrast-invariance; on the other hand considering the tonic response alone will show reduced or absent contrast-invariance.

Virtual rotation.
In their analysis of the canonical FR ring-model, Ben-Yishai and co-workers [1] introduced the concept of virtual rotation in the marginal phase regime, defined as the non-instantaneous movement of the peak activity following a new stimulus. They suggested that virtual rotation may be implicated in the psychophysical phenomena of mental rotation, the perception of three dimensional attributes after discrete views of an object from different angles. From the point of view of analysis, virtual rotation is one example of attractor dynamics, which in turn may provide a fundamental signature of different cortical architectures.
The stimulus sequence that we use here is particularly appropriate for examining this phenomenon, and can be appreciated in Fig 12, which traces the center of mass of the population firing rate over time. As expected given its position in the marginal domain, the canonical FR ring model G here shows virtual rotation, thus a time constant for the shift of orientation preference of approximately 30ms, as do the LIF and FR-w/shunt models D, E and F. This effect is weaker, thus the shift in orientation faster, in the non-adaptive 2 compartment LIF model C, and is essentally absent in the adaptive HH models A and B.

Tilt after effect.
The tilt after effect [35,36,37] refers to a perceptual bias of an oriented stimulus away from the orientation of a preceding stimulus. Electrophysiological studies indicate that this effect arises, at least in part, from primary visual cortex [38,39], and this hypothesis has been further studied with models of V1 (e.g. [40]).
As the case for virtual rotation above, the stimulus sequence examined here is well suited for studying the tilt after effect. In our models this effect is seen to varying degrees in Fig 8, as an overshoot of the maximum activity beyond the orientation of the second stimulus. Specifically, this effect is seen clearly between roughly t = 110 and 150 ms in the 2D and ring adaptive HH models A and B, whereas removing adaptation in the subsequent LIF ring model C suppresses the effect. The physiological studies of Dragoi and colleagues [41] have confirmed a crucial role of adaptation both for suppressing the response to the new stimulus on the side near the first stimulus, and shifting its peak activity. Inspired by these findings, the modeling study of Jin et al [40] focused in part on the impact of adaptation mechanisms on the time scale of minutes. This can be compared to the adaptation mechanisms used in the present study (M and AHP currents) which have time scales on the order of one second or less; indeed, as Jin et al note, fast adaptation (as studied by [41,42,43,44]) gives a weak tilt after effect [45]. Nevertheless, the qualitative demonstration of the effect in the models examined here is consistent with additional and longer lasting adaptation mechanisms, and provide a basis for exploring this phenomena with biophysically-constrained population models.
2.3.6 Spatial distribution effects. The consideration of a 2D cortical geometry can be related to two aspects of the model behaviors shown here. First, in the full 2D model A there is a gradient as a function of distance from the pinwheel centers in the steady-state or quasisteady-state responses to both stimulus orientations, thus at t = 90ms and t = 190ms (Figs 6  and 7). This effect can be explained by the distribution of inhibition along the vector towards a pinwheel center. Thus, neurons far from the pinwheel center receive nearly circularly symmetric inhibition, whereas those closer to the center receive inhibition increasingly limited to the side away from the center. Of course, as stated this effect derives directly from the assumption of homogenous cortical connectivity in the model, and thus may be altered with a different connection scheme. Thus the model predicts higher activity near pinwheel centers if the input connections are homogeneously distributed along the radius on a scale of a single pinwheel. Second, the tilt after-effect is more pronounced for neurons remote from the pinwheel centers, as seen from curved isolines in the frames at 110, 120, 130 and 150ms. By construction, these features cannot be reproduced in any of the ring models.
In summary, the residual differences between the adaptive HH 2D model A at a fixed pinwheel radius, and it's immediate ring version, the adaptive HH ring model B, are most likely due to the qualitative reduction of the 2D to 1D geometry as well as quantitative approximation of the synaptic projections to the ring by the cosine profile.

Parameter analysis of the ring models
We now present an analysis and implications of the mapped parameters of the reduced ring models, especially the canonical ring model G, in terms of the original biophysical parameters of the most detailed 2D model A.

Reformulation of the mapping expressions.
For the sake of analysis, we define several terms to simplify the expressions Eqs 41, 42 and 44, in a manner that collates parameters of the different circuit pathways.
We first develop expressions for theĨ 0;1 ,J 0;1 ,K 0;1 andL 0;1 terms in Eqs 9 and 10 that define the single population conductance-based rings (models E and F, repeated here for convenience): Iðt; yÞ ¼Ĩ 0 ðtÞ þĨ 1 ðtÞ cos ðy À y 0 ðtÞÞ þ 1 2p Referring to Eq 41, for the thalamic input coefficientsĨ 0 andĨ 1 , we define: These terms characterize specific paths within the network. First, A EE describes the explicit feedforward excitatory path from the thalamus. Recalling thatk LIF I (Eq 36) reflects the effective gain of the inhibitory population (less than the actual gain of these neurons because of inhibitory-inhibitory feedback), B IE describes the concatenation of the implicit inhibitory feedforward and recurrent connections, and C IE characterizes the implicit inhibitory recurrent connections.
We define similar terms for the intracortical coefficientsJ 0 andJ 1 : Note the similarities between D EE and A EE , and between E IE and B IE , respectively, in each case differing only by the final c ij term characterizing the source pathway. Thus, D EE describes the explicit recurrent intracortical excitatory path, and E IE describes the concatenation of the implicit inhibitory and excitatory recurrent connections of the conductance-based single population models.
Using these definitions, we now re-arrange Eq 41: Likewise, referring to Eq 42, the corresponding conductance terms are given by: where for convenience, we now introduce the notation "F j DV a " such that all the ΔV terms in some function F are set to α. For the case of Eq 48 (and later equations) we set α = 1 (unitless) to obtain a shorthand expression for the synaptic conductance.
We now complete the mapping by expressing the I 0,1 and J 0,1 terms of the current-based ring (model G) directly in terms of the biophysical parameters (Eq 3, repeated here for convenience): Finally, considering Eqs 44 and 47, we obtain:

Feedforward inhibition.
Conceptually, one can consider distinguishing functional excitatory and inhibitory synaptic pathways in terms of their thalamic, or feedforward, component, versus their intracortical component. In this context, the fact that there are no anatomical inhibitory connections from the thalamus has led to the assumption that a population of cortical inhibitory neurons, driven primarily by thalamic input, act as a surrogate for a true feedforward inhibitory pathway. The mapping to the FR ring model G allows a quantitative interpretation of this assumption.
We recall that in the original formulation of the canonical ring model, the I 0 and I 1 terms (Eq 3) describe the thalamic input. The mapping expressions for these terms in Eq 50, specifically the B Th IE term for the FR ring, provide an explicit quantitative definition of so-called "feedforward inhibition": Recalling the definition of DV Th ij given by Eq 28, we see that B Th IE is negative because of DV Th I E . Notably, this characterizes the component of the functional thalamus-driven pathway, allied with cortically-driven inhibition (parameterized by E Th IE ), that competes against feedforward and intracortical excitation.
Finally, as seen in Eq 50, the relative strength of the un-tuned versus tuned components of the feedforward inhibition depends on the tuning of the intracortical inhibitory-excitatory pathway, parameterized by q I E .

Contrast and tuning of thalamic input.
We now compare the thalamic input across the models, characterizing them in terms of a "contrast" parameter O and a tuning parameter Γ. We define O as the average strength of the input to the network, thus across orientation space. The response of many visual neurons is strongly related to the visual contrast of the stimulus, e.g. the response scales with the amplitude of a grating, after adaptation to its average value. The tuning parameter Γ characterizes the orientation dependence of the retino-geniculo-cortical, thus thalamic, pathway to the cortex. Indeed, as noted before, this pathway is explicitly established by the constant plus cosine form of thalamic input for all the models. For the two population models driven by thalamic firing rates, we have: Given that firing rates must be non-negative, note that the constant plus cosine form (1 + cos (θ)) imposes the constraint 0 � Γ φ � 1. We follow the same form for defining input contrast and tuning of the single population conductance-based ring models, but now introducing separate definitions for the current and conductance inputs: With these definitions, we now establish their relation with the contrast and tuning of the thalamic firing rate input for the two population models. ExpandingĨ 0;1 (ref. Eq 47), and considering that φ 1 = Γ φ O φ : We see that there is not a one-to-one correspondance between the measures for the different architectures. First, for zero "true" thalamic input defined by firing rates, the equivalent contrast of the current and conductance inputs is −C IE and À C IE j DV 1 , respectively, which reflect the subthreshold constant terms that result from eliminating the inhibitory population. Second, these offsets make the tunings of the current and conductance inputs dependent on the contrast of the thalamic firing rates, whereas in princple tuning is defined purely by anatomical connectivity, regardless of contrast. Finally, the coefficients for the input rates φ 1 and φ 0 in the tuning expressions are not equal unless B IE ¼ B IE j DV 1 ¼ 0 (e.g. zero thalamic input to the inhibitory population of the two population models, thus c Θ = 0).
If we assume that inhibition is purely shunting (V I = V Rest and therefore B IE = C IE = 0), then the contrast and tuning for the current input, at least, simplifies: In this case, compared to the two population models A-D, we see that the "contrast" of the conductance-based single population models E and F are the same apart from a scaling term, and more importantly the tunings are identical. We now consider contrast O FR and tuning Γ FR for the canonical FR ring model G, again following the definitions for the other models: To simplify the previous expressions for I 0,1 (ref. Eq 50), we define: Applied to the expressions for I 0 and I 1 , we obtain: We re-express these equations, now in terms of the original definitions of contrast and tuning (Eq 51): Given the definitions of contrast and tuning for the canonical FR ring (Eq 52), we obtain: As before with the input terms of the conductance-based single population models (except for the case of pure "shunting" inhibition), the contrast of the input to the FR ring model is a nonlinear function of the contrast as defined for the original two population models.
We can supply some constraints on A 0 and A 2 by asserting that intracortical connection weights are non-negative, thus for the constant plus cosine approximation 0 � q i,j � 1. This gives A 0 � A 2 , and therefore: Note that Γ FR diverges as O φ approaches B 0 /A 0 , which in turn is equal to the threshold firing input: The distortion of the input contrast and tuning when compared to the thalamic firing rates is due to B 0 , which encapsulates the impact of the original two firing thresholds (of the HH excitatory and inhibitory neurons) on the final mapping. As one exercise, note that since C Th IE / I LIF Off;I , at the limit of ðI LIF Off;E ; I LIF Off;I Þ ¼ 0, then B 0 = 0. Of particular interest are the quantitative values of the input to FR ring model G, specifically that there is a strong negative current away from the preferred orientation. This aspect runs contrary to the general assumption of the canonical FR ring model, thus that the excitatoryonly thalamic pathway to the cortex constrains the (current) input to this model to be nonnegative. This result arises from the effect of the cortical connections on the thalamic input terms as defined by the mapping. This is first seen in the mapping of the two population conductance models to the single population conductance models, specifically as expressed in the equations describing the cortical input, thus forĨ 0;1 (Eq 41) and, by extension,K 0;1 . In particular, these expressions include negative terms that depend on the strength of the intracortical connections between inhibitory and excitatory populations, thus directly proportional to c IE (for bothĨ 0;1 andK 0;1 ) and the weighting parameter q IE for the cosine approximation (forĨ 2 andK 2 ). Note that these same dependencies hold for the intracortical connectivity parameters J 0;1 (ref. Eq 41) and by extensionL 0;1 . This "crosstalk" onto the thalamic input is conserved in the final mapping between the single population conductance models and the current-based FR ring, as seen in the equations for I 0,1 (ref. Eq 44) which depend directly onĨ 0;1 andK 0;1 . As with the single population conductance models, there are similar dependencies for the intracortical parameters J 0,1 (Eq 44), which depend directly onJ 0;1 andL 0;1 .
We propose that the quantitative mapping developed here provides a novel explicit interpretation of the concept of "feedforward-inhibition", as introduced in Section 2.4.2, thus corresponding to the negative input current shown here.

Intra-cortical connections and the dynamical steady-state behaviour of the network.
We have shown that the diagram of the steady-state solutions mapping formulas provide an interpretation of the coefficients J 0 and J 1 of the FR ring model in terms of quantitative biophysical, thus synaptic and neuronal, parameters. This then allows an analysis in terms of the phase-plane of the FR model, as described by [2], shown in Fig 10. From this plot we see that the parameter values of the derived ring model G place it in the marginal phase, and thus allowing for attractor dynamics.
The mapping expressions (Eq 50) show that all the firing-rate ring model parameters, I 0 , I 1 , J 0 and J 1 , are proportional to the gain of the excitatory neuron population, k LIF E . This underlines the fundamental role of k LIF E , a purely local cellular, and not network, property, with respect to the phase space of the canonical FR model: a linear trajectory in this space necessarily begins in the homogeneous regime (which includes the origin). Depending on the relative values of the network components D Th EE and E Th IE (parameterizing the recurrent excitatory and inhibitory pathways, respectively, and weighted by the q ij parameters in the case of J 1 ), the trajectory subsequently passes either into the marginal or unstable regime and remains there, or first into the marginal and then the unstable regimes.
As constrained by the biophysical parameters, this dependence is significant across the model spectrum described here, most crucially on J 1 . Indeed the removal of adaptation between models B and C significantly increases the excitatory cell gain (from 0.069 Hz/pA for the adaptive neuron to 0.23 Hz/pA for non-adaptive neuron, ref. Fig 5), thus k LIF E . Derivation of an "adaptive" version of model G underlines this significance, where the value of k LIF E derived from the adaptive HH excitatory cell model puts the model in the homogeneous state (Fig 10). As presented later, this qualititative difference is manifested in several ways by the simulations. The effect of rate-dependent adaptation can be qualitatively compared to the effect of tuned inhibition which decreases J 1 , both providing negative contributions to the terms with the rate in the right hand part of Eq 1. The consequences of an increase of the adaptation or the tuned inhibition are similar, however the adaptation is more concentrated than the cosine-like tuned inhibition, and hence is more efficient.
We now examine D Th EE and E Th IE in more detail. As the case for constraints on the biophysical terms A Th EE , B Th IE and C Th IE described earlier, reasonable assumptions on the reversal potentials for the excitatory and inhibitory synapses imply: Thus, the expressions for J 0 and J 1 show explicitly the direct competition between the net strengths of excitatory and inhibitory recurrent pathways in determining the behavior in phase space, specifically both the quadrant that the network is restricted to (given by the signs of J 0 and J 1 ) and the sensitivity of the location in phase space as a function of the excitatory gain k LIF E . For J 0 the competition between the pathways is independent of the anatomical spread of the various connections. In contrast, for J 1 the contribution of the excitatory recurrent network is weighted by its characteristic anatomical tuning parameter q EE , while the inhibitory recurrent network term is weighted by the product of q IE and q EI , reflecting the effective convolution of the recurrent excitatory-inhibitory and inhibitory-excitatory pathways.
The sign of J 1 establishes the qualitative relation between the intrinsic stable attractor states and the input. Thus, a positive J 1 causes the attractor steady-state to line up with the stimulus, while a negative J 1 tends to cause the attractor steady-state to be orthogonal to the stimulus. Note that this applies strictly to the limit case of homogenous input, because the actual steadystate of the network is established by both the inherent attractor properties and the quantitative value of the stimulus. For illustration, we consider the limit cases for a network with flat excitatory recurrent connections (i.e. q EE = 0), flat inhibitory recurrent connections (i.e. q EI = 0 or q IE = 0), or both. Noting that J 0 is unaffected by these parameters, for the first condition: and thus the network is constrained to the left half of the phase space. For a network with flat inhibitory recurrent connections: and thus the network is constrained to the right half of the phase space. These limit cases provide a more general interpretation, thus that the effective width of intracortical inhibition must be broader than that for intracortical excitation to allow the attractor properties to complement the input, and vica-versa. Finally, if both intracortical pathways are flat, then J 1 = 0, and the marginal regime is inaccessible. To our knowledge the possible relevance of a negative J 1 to actual biological circuits has not been explored.
In the present case, the biophysical parameters establish that J 0 < 0 and J 1 > 0, thus limiting the network to the fourth quadrant of the phase space as a function of k LIF E . The specific mapping of the canonical FR ring parameters for the intracortical connections predict that the network is operating in the marginal phase, which is of particular interest in the context of attractor networks. More specifically, this has been used to interpret the contrast invariance of the network. However, the distortion of the input parameters of the canonical model described previously, with respect to the original thalamic input parameters φ 1 and φ 1 , complicate the interpretation. As described earlier, a simple change of contrast O φ = φ 0 for the full model while maintaining the same tuning (defined as Γ φ = φ 1 /φ 0 ), maps to an input for the canonical FR ring model with a different tuning.

Mapping assumptions
We can speculate on the implications of the several steps to eliminate specific non-linearities in order to achieve the final linear model (apart from the FR threshold). First, the role of inhibition may be under-estimated for two reasons. The experimental evidence on the impact of synaptic shunting on neuronal transfer functions is mixed, with some studies reporting a pure subtractive effect [46], while others showing a mixture of a threshold shift and gain reduction [47,48]. The linear mapping here necessitates neglecting any gain change, implying that the gain of the current-based FR ring may be over-estimated, particularly in the face of strong recurrent inhibition. Furthermore, the anatomical spread of inhibitory-inhibitory connections in cortex is non-negligible on the scale of the hyper-column. The mapping requires that this spread be ignored to allow the elimination of the explicit inhibitory population, again with the result that the impact of inhibition in the final model may be under-estimated.
Second, the role of excitation within cortex may be overestimated for the following reason. A necessary assumption here is that the inhibitory population is always suprathreshold, allowing the actual threshold of this population to be ignored. However, given the large dynamic range of cortical dynamics, it is likely that different populations of neurons will be sub-threshold at different moments, especially for relatively small, but realistic, stimuli. In the context of the single population models, this implies that with small stimuli these models will transform inhibitory paths to excitatory (e.g. at the tails of the activity profile), effectively equivalent to a broader anatomical tuning of the excitatory recurrent connections.
Finally, the post-synaptic response in biological neurons to an incoming spike train saturates for several reasons, most importantly because of the reduction in driving force as the post-synaptic potential approaches the reversal potential, the finite number of post-synaptic receptor/channel complexes, not to mention mechanisms such as synaptic adaptation [49]. Indeed, only the first mechanism is considered in the two population models (in part in the single population LIF model for the excitatory population). The net result for the final FR ring is that the activity in the ring may be overestimated, particularly the peak at the preferred orientation.

Interpretations from the mapping
The proposed mapping allows a step-by-step fitting of parameters of increasingly abstract models to experimental data. At the most abstract level, the mappings allows the interpretation of the canonical firing-rate ring model in a slightly uncommon but more rigorous way. We see, for example, that although I 0 and I 1 are defined as thalamic input, when mapped from the full model, both parameters depend on the inhibitory-to-excitatory intracortical connections � g I E (Eqs 50, 49, 66 and 28), and therefore can not be distinguished as pure background and tuned inputs. On the other hand, study of information processing in the visual cortex by the means of the detailed model, also requires fitting of the canonical firing-rate ring model to the experimental data, in order to limit the domain of physiologically meaningful parameters' values. This work exceeds the frame of the present paper and will be done in the future.

Transient behaviour of the models
After the parameters of the full model A are established, the mapping procedure provides parameters of the progressively more abstract models that are explicitly constrained by replicating the steady-state regime. At the same time, we note that the dynamics of the full model A show a rich variety of the solutions, including oscillations, waves etc. not seen in the more abstract models. Nevertheless, the systematic mappings provide constrained reference points in the parameter spaces of the models, which give essentially similar behaviors between nearest neighbors among the model hierarchy.
Generally, models with more complex dynamics possess a larger variety of solutions, like the FR ring model with a delay [5] in comparison to the canonical ring model [2]. We note that some predictions of the marginal domain of the canonical FR ring-model are maintained throughout the model hierarchy, specifically sharpening of the output tuning, while others fail for more complex models, specifically contrast invariance and virtual rotation. Summarizing the most significant differences between the models, we underline that the reduced models do not reproduce the ripples of activity provided by the synchronization of spikes and refractory effect after stimulus presentation, nor the effects of firing adaptation. We failed to observe any significant delay in the detailed models due to virtual rotation effect predicted by the reduced models, instead the reaction was fast. For the parameter sets corresponding to more prominent virtual rotation in the canonical FR-ring model we observed oscillations in the detailed models. These observations warn that an oversimplification may lead to misleading conclusions in regards to transient behaviour of the reduced models.

Comparisons of mean-field descriptions, and limitations with respect to discrete neuron network models
Refractority. The consideration of refractority is one the most significant factors that lead to discrepancies between discrete neuron networks and most mean-field approaches [50]. Discrete network models explicitly describe each neuron in the network, and thus all their states over time. At any instant different neurons are generally in different refractory states, mainly because of different timing of their last spikes, and hence a discrete model explicitly describes the distribution of neurons by their refractory states. In contrast, most mean-field models, including our models Fand G do not, instead assuming that the activity is asynchronous. While this assumption is violated for most parameter settings, if the networks are not too synchronous such mean-field models are in good agreement with simulations of discrete networks [51]. In contrast, the more elaborated KFP and CBRD mean-field approaches, including our models A-E describe the refractority correctly (and explicitly) for all activity regimes (Sections 4.1 and 4.2).
Number of synapses. In discrete neuron network simulations, the finite number of neurons and synaptic connections determine the synaptic current fluctuations. Mean-field models generally assume a weak coupling limit [52,53], which can lead to results different from simulations of discrete neurons [54]. The models using the HH or LIF cell model (models A-E) approximately take into account the scaling of synaptic noise with the number of synapses, by scaling the noise term by the mean synaptic conductance (Eq 23), thus avoiding the weak coupling limit.
Number of neurons. The same refractory density framework described here can also take into account a finite number of neurons per se [55], though this extension is beyond the scope of the present paper.
Sparse strong connections. By definition, our mean-field framework cannot take into account the activation for each individual connection. In our experience, as an example, this can be important for networks with a fraction of strong connections, e.g. those with long-tail weight distributions, and only for spontaneous activity, not the one driven by stimulus [56]. In that work, we compared the effects of different connectivity profiles on the population dynamics of discrete neuron networks. For long-tailed synaptic weight distributions, we showed that both strong sparse and very weak connections play crucial roles for spontaneous activity, whereas they have little influence on driven population activity. Moreover, neurons with strong connections form multiple chains which are the core of spontaneous activity, with neurons with weak connections only partially enhancing the population activity. These observations are in agreement with the other works focused on clustering in networks [57]. Therefore, the stimulus evoked activity may be reasonably simulated in neglect of the sparse strong connections, assuming that the majority of connections have approximately normally distributed weights. Because the population equations of the KFP and CBRD models well approximate the mean firing rate of neurons with most common intermediate connection strengths, and because we consider only driven activity regimes, we derive the validity of our approaches for the problems solved in the present paper.
Numerous strong connections. In the case of a sufficiently heavy long-tail synaptic weight distribution, the mean-field approach requires an extension. For the case of a lognormal distribution, this extension for the CBRD approach was proposed [25], therefore the CBRD models described here can in principle work with such connection profiles. For the sake of simplicity, however, heavy long-tail distributions of synaptic weights were not the focus of the present work.

Conclusions
We conclude that the consistency of the activity patterns of the considered models support the validity of the obtained mapping expressions, and thus permits understanding the significance of the various assumptions made during the derivation of the models' equations. The constructed hierarchy of models can therefore serve as a useful instrument for the fitting of mathematical models [58][59][60] to experimental data and their subsequent analysis. In particular, the architecture of the visual cortex ring model has been recapitulated in models of other systems, either generalized [61] or specific, for example head direction cells [62,63] and prefrontal working memory [64] in rodents, as well as navigation circuits in the fly [65,66]. The methods developed here should be amenable towards constructing biophysically constrained abstract models of these systems.

Methods
To build a model of interacting populations of neurons we explore the probability density approach [11,67,68]. It is based on the equivalence of the consideration of a stochastic differential equation for a single neuron to the probabilistic consideration of neuronal density evolution in the phase space of neuronal state variables. In this framewok, simulations of the stochastic model is equivalent to consideration of a statistical ensemble of similar neurons receiving similar input signals.

Kolmogorov-Fokker-Planck (KFP) approach for LIF neurons with instantaneous synapses
Here we review the probabilistic evaluation of an infinite population of 1-compartment LIF neurons with Gaussian membrane current noise, used to obtain the simulations of the single population LIF ring model (model E) [68]. Thus, Eq 4 can be considered as a Langevein equation for a single neuron. This equation is equivalent to the KFP equation written for the probability density function, or neuronal density, in the voltage phase-space, ρ(t, V) (the temporal dependencies of the synaptic current and conductance input, I and S, are implicit): with the effective membrane time constant of the LIF neuron, τ m (S), given earlier by Eq 7, and the average firing rate, ν LIF (t), given by: The set of Eqs 53 and 54 is equivalent to an infinite set of Eq 4. We compare the current step response of different models (Fig 13), in particular against simulations of individual LIF type neurons using a Monte-Carlo framework which provides a "gold-standard".

Conductance-based refractory density (CBRD) approach for multiple state variable neurons
Because of their multiple state variables, evaluating noisy neural populations with non-instantaneous synapses, multiple compartments and/or active HH channels with a KFP-based method, as described above, cannot be evaluated analytically, and thus a KFP (or Monte-Carlo) approach requires significant computational resources [11,69]. To address this problem, we employ a CBRD approach that can evaluate a broad class of models in a one dimensional phase space parameterized by the refractory time, or time since the last spike, t � [17,25,70]. In the present case, this approach is used to evaluate all the two population models, including the HH models A and B, and the LIF models C and D.
In conjunction with a given neuron model, the CBRD approach considers a hazard function H which defines the probability density of firing, and describes the dynamics of a population as a whole by the probability density of neurons distributed according to their refractory times, ρ(t, t � ). In particular, the CBRD model substitutes the phenomenologically derived hazard function developed in the original refractory density approach [68] by its rigorously derived expression [17,70]. The population model for the neuron model is expressed by a set of transport equations, one for each state variable, in partial derivatives with t � and time t as independent variables, and thus the number of equations is directly proportional to the complexity of the model.
In order to obtain an expression for H in a probabilistic framework, e.g. in the presence of noise, the CBRD approach requires that the state variables be well described at any given t � by a linear model around their mean value. This constraint holds for most of state variables of the HH model given their voltage-dependence and a typical dispersion of the membrane voltage of several millivolts. For the LIF model, this constraint is automatically satisfied. In contrast, the strong nonlinear properties of sodium channels of the HH model near threshold do not meet this criterion. To address this problem, and as mentioned in Section 2.1.3, we use a threshold HH cell model without sodium channels, with spikes defined by an explicit dynamic threshold and renewal function. This model is based on three important assumptions: first, that the probability of firing can be well predicted by the dynamics of a neuron without an explicit sodium current, that the influence of strongly non-linear channel dynamics, namely the sodium current, is only significant during the spike, and third, that there is a fixed (i.e. renewal process) impact of each spike on the state variables. The last assumption allows imposed conditions on the state variables following a spike, for example reset values for the voltage and fast potassium channel gating variables (e.g. for I DR ), and increments for slow potassium channel gating variables (for I AHP and I M ), as appropriate. We have shown previously that these approximations are valid in a range of neuron models [17]. The associated threshold model compares quite well to the full conductance model with respect to interspike potentials and spike time moments.
The dynamics of a given population in the CBRD model are thus described by a set of equations of the form dY(t, t � )/dt = F(), where Y() includes the distribution ρ(t, t � ) of neurons with a given value of t � , the average soma voltage U(t, t � ), the average dendrite voltage U D (t, t � ) for two compartment models, and finally the average gating particle states x(t, t � ) for HH cell models. Since the refractory time t � between spikes is proportional to t, the total time derivative dY(t, t � )/dt may be replaced by a sum of partial derivatives in t and t � , i.e. d/dt = @/@t + @/@t � .
The source term in the right-hand side of the equation for ρ(t, t � ) is the fraction of neurons crossing the threshold, given by ρ(t, t � ), multiplied by the hazard function H: The remaining terms consist of the model details. The expressions corresponding to the cell voltages derive from those presented in Section 2.1.3, specifically the sets of Eqs 13 and 14, with the following replacements of the average voltages U, U D and ξ, for V, V D and 0, respectively: The boundary conditions for ρ (which also defines the population firing rate, ν(t)) and the model voltages are: where Dt AP j is the duration of the action potential for cell type j (= 0 for the LIF model; for the HH cell models see Sections 4.3 and 4.4); likewise, the somatic reset voltage V Reset j is defined also for each cell type.
The CBRD expressions for the gating particle states x(t, t � ) of the HH channel models, including their boundary conditions, are presented in the next sections (Sections 4.3 and 4.4), which in turn reference the average somatic voltage U.
The dynamics of the entire neural population are found by the integration of Eqs 13 and 14, with the replacements indicated above for the average voltages U and U D , which then defines the distribution of cell voltages across t � at each time t (for the HH models similar equations are solved for the gating particle state variables). The effect of threshold crossing in response to the input and noise is then taken into account by the H function, with the integration of Eq 55 giving the distribution of ρ across t � as well as the firing rate ν at time t.
As the source term in the density equation, the H function is the solution of the classical first-time threshold crossing problem for arbitrary history of stimulation. The H function has been semi-analytically derived from the KFP equation for voltage fluctuations, based on the following assumptions: (i) Away from threshold, voltage fluctuations due to individual noise realizations can be described by a linear equation given the mean voltage U and the mean membrane conductance, which are in turn given by the associated transport equations and the synaptic input. Furthermore, the evolving probability density distribution of voltage fluctuations about U can be described by a KFP equation.
(ii) The flux across threshold described by H are due to two additive underlying processes: diffusion along the voltage axis due to noise, and transfer in response to the input.
(iii) Diffusion due to noise dominates the firing when U is constant, e.g. when the input is constant. In this case the problem is described by the KFP equation with a leak and a constant threshold, denoted here as A. For white noise the governing equation is reduced to an ordinary differential equation with an analytical solution for the spike generation probability expressed in Kummer functions [17]. In the case of correlated (colored) noise, A is obtained numerically and tabulated for a range of values for U and the time constant of the noise correlations τ [70].
(iv) Transfer based firing dominates when excitation starts abruptly, i.e. U increases with infinite rate. The initial condition is a stationary Gaussian distribution of neurons in terms of their potential V(t) − U(t). This implies that the state of zero distribution immediately following the previous spike is forgotten at the time of the next threshold crossing. As U(t) increases, the moving boundary at the threshold of fluctuations, V Th − U(t), crosses the distribution. The flux through the boundary determines the probability of spike generation, which we denote in this case as B, and is obtained algebraically.
In our previous work we formulated an approximation for H for white Gaussian noise, as a function of the time-varying quantities that characterize the cell, including U (and its derivative with respect to time t at a given t � ), σ V , V Th and the effective membrane time constant τ m = C/g tot , where g tot is the total membrane conductance: Note that the hazard function above has no free parameters, because they reflect the approximation of the first time passage solution. The dynamics of a neural continuum are thus evaluated at each time step by solving the set of one dimensional partial differential equations in terms of the refractory time, including the probability density of neurons, and the population averages of the membrane potential and channel gating variables. The CBRD model well approximates the firing rate of an infinite set of biophysically detailed neurons for an arbitrary stimulus, e.g. oscillatory input, and recurrent connections, when compared with Monte-Carlo simulations of individual neurons (see [17] and [70]). Fig 13 compares the response of a population of noisy LIF neurons, evaluated by Monte-Carlo simulations, a KFP approach (Section 4.1), and the CBRD method, with the response of the FR model with and without a time constant. Fig 14 compares the steady-state f/I characteristics as a function of synaptic conductance for a population of LIF neurons, with and without noise (Eq 5), with the result from the CBRD evaluation.

Excitatory two compartment HH neuron model
Here we complete the description of the excitatory neuron model as described by Eq 14, with the average voltage U substituting for V. As stated earlier, this model incorporates several HH membrane currents, including the voltage-dependent potassium currents responsible for spike repolarization, I DR and I A , the voltage-dependent potassium current that contributes to spike frequency adaptation, I M , the voltage-dependent cation current, I H , and calciumdependent potassium current that also contributes to spike frequency adaptation, I AHP . The • Synaptic parameters: � g YE ¼ 0:011 mS, � g Y ¼ 0:002 mS, � g EE ¼ 0:053 mS, � g EI ¼ 0:01 mS, � g I E ¼ 0:11 mS, � g I I ¼ 0:02 mS, t r • Spatial connections: d EE = 100μm, d EI = 200μm, d IE = 200μm, d II = 200μm.
• Numerical parameters: time step 0.1ms, spatial grid 40 × 40, the t � -space was limited to 100ms and discretized by 100 points. The absolute refractory period was introduced by setting the term A(T) of the hazard function to zero for t � < 6ms.
• Numerical parameters: ring was discretized by 40 points.

LIF Models C (LIF 2E EI Cos Ring), D (LIF 1E EI Cos Ring), E (LIF Cos Ring).
Excitatory model: The parameters of V Reset LIF;E ¼ À 90mV, C E = 0.27 nF are determined in Section 2.2.2, V Reset LIF;I ¼ À 82mV was obtain by fitting as demonstrated by
• Numerical parameters: ring was discretized by 40 points, voltage space ranged from -105 (to be more negative than any expected polarization in neurons) to V Th = −57mV was discretized by 250 points; time step 0.02ms.

Model G (FR Cos Ring).
• The characteristic time constant of the FR model was assumed to be equal to the membrane time constant of the excitatory neurons, t FR ¼ t 0 m;E given in Section 4.3. • Connections: I 0 = −20pA, I 1 = 43pA, J 0 = −0.35pA/Hz, J 1 = 2.7pA/Hz were calculated by Eq 44.
• Numerical parameters: ring was discretized by 40 points; time step 0.02ms.