Contrast response in a comprehensive network model of macaque V1

The response to contrast is one of the most important functions of the macaque primary visual cortex, V1, but up to now there has not been an adequate theory for it. To fill this gap in our understanding of cortical function, we built and analyzed a new large-scale, biologically constrained model of the input layer, 4Cα, of macaque V1. We called the new model CSY2. We challenged CSY2 with a three-parameter family of visual stimuli that varied in contrast, orientation, and spatial frequency. CSY2 accurately simulated experimental data and made many new predictions. It accounted for 1) the shapes of firing-rate-versus-contrast functions, 2) orientation and spatial frequency tuning versus contrast, and 3) the approximate contrast-invariance of cortical activity maps. Post-analysis revealed that the mechanisms that were needed to produce the successful simulations of contrast response included strong recurrent excitation and inhibition that find dynamic equilibria across the cortical surface, dynamic feedback between L6 and L4, and synaptic dynamics like inhibitory synaptic depression.

each Magnocellular LGN neuron, as in Zhu et al., 2009). The phase at location x is given by p(t, x) = wt + x, wu (where , is the inner product). Constants are chosen so that background firing rate of an LGN cell is ∼ 20 spikes/sec, and when driven by a grating at full contrast, peak firing rate is ∼ 100 spikes/sec (Hawken et al., 1996). The contrast-response function of the LGN was simply modeled as the amplitude vs contrast from the measured values in macaque M ganglion cells (Kaplan and Shapley, 1986) and LGN Magnocellular neurons (Kaplan et al., 1987).

The basic integrate-and-fire equations for 4Cα neurons
4Cα neurons are modeled as conductance-based, integrate-and-fire point neurons; the membrane potential V of each evolves according to where C is membrane capacitance; g L is the leak conductance; g AMPA and g NMDA are the time-dependent, excitatory conductances corresponding to AMPA and NMDA receptors; g amb is a conductance term representing the modulating effects of multiple substances not modeled; g GABA is the time-dependent inhibitory conductance from GABA A receptors; and V L , V E , and V I are the leak, excitatory, and inhibitory reversal potentials. The equations governing the time evolution of the conductances are given below. For simplicity the voltage is rescaled so that V L = 0 and the spiking threshold is at V = 1. Capacitance is also moved to the RHS and absorbed in the conductances. As a result, instead of using equation (1) directly, we rewrite it as where, in rescaled units, τ L = 20 ms for E cells and 15 ms for I cells (Beierlein et al., 2003), the spiking threshold is V = 1, upon spiking V is immediately reset to 0, V E = 14/3 and V I = −2/3 (Koch, 1999). Following a reset after each spike, V is kept at 0 for a refractory period of 2 ms for E cells and 1 ms for I cells. When the refractory period is over, the evolution of V resumes according to equation (2).
A remark on units: Membrane capacitance was normalized to 1 (McLaughlin et al., 2000). Therefore, conductance was proportional to the inverse of the membrane time constant and has the same units, sec −1 . A normalized conductance of 1000 sec −1 was equivalent to ∼ 20 nS. Membrane current was conductance scaled by the difference between membrane potential and the reversal potential for the particular conductance channel.
The currents also were calculated in normalized units where, for current, sec −1 is equivalent to ∼ 300 pA.
To state the conductance equations, we introduce the following notation. The synaptic conductance terms g AMPA , g NMDA , and g GABA are driven by the spiking activity of presynaptic cells. For neurons of type Q ∈ {E, I} and Q ∈ {E, L6, LGN, amb, I}, we denote by S QQ the synaptic strength of a spike when the presynaptic neuron is of type Q and postsynaptic neuron is of type Q. The ambient conductance g amb is driven by a Poisson-timed train with rates ρ E amb for E cells and ρ I amb for I cells. The time courses for EPSC's and IPSC's are denoted by functions G AMPA , G NMDA , G GABA and G amb , which are given below, following the conductance equations. The model parameters ρ E amb , ρ I amb , and S QQ are given in section 0.5. For a layer 4 neuron of type Q, the equations governing the time-evolution of its conductances in equation (2) can now be stated as follows: Here t E i are the arrival times of E spikes from all neurons in layer 4Cα presynaptic to the neuron in question; t L6 As for the conductance time-courses, for R = AMPA, NMDA, GABA or amb, the time courses are

Synaptic and adaptation dynamics
Here we modify the equations of section 0.2 to include synaptic and adaptation dynamics.
To model synaptic depression of I cells (Galaretta & Hestrin, 1999;Gibson et al., 1999), i.e., the phenomenon where I cells become less effective as they increase in spike rate, we vary S EI as follows. If an I cell spikes at time t 0 and subsequently at time t 1 , the synaptic strength of the second spike is whereS QL6 = S QL6 is a constant and the function g is given in section 0.5.

To model the adaptation of E cells in layer 4 (Yuan et al., 2005), we assume the firing threshold is elevated
after each spike, to be relaxed back in time to its original value. If an E neuron spikes at time t 0 , then its spiking threshold follows where the values ofθ and τ adp are given in section 0.5.

Determination of feedback from L6 to 4Cα
Layer 6 neurons are made to spike at a rate which is modulated by the local 4Cα firing rate activity. We assume there is a response function f (R) so that when R is the local layer 4 firing rate, f (R) is the local layer 6 response. The function f is chosen to be sigmoidal and to fit the data at background and high contrast (Ringach et al., 2002).
We stress that this response function does not depend on cortical location or on the nature of the stimulus, i.e., we assume L6 response at any location under any condition depends solely on L4 activity at that location at that point in time. The function we used in our simulations is given at the end of section 0.5.
In more detail, let r 4 (x,t) denote the instantaneous firing rate of a layer 4 neuron at location x and time t. We compute the local firing rate R 4 (x,t) by averaging over r 4 within a disc of radius 37.5 microns and over time, weighting the influence of the past by a temporal kernel K(t): where |D| is the area of the disc, and K is a sigmoidal function in the form Here t 1/2 is a constant describing approximately how long a time-window we are averaging over. This average layer 4 rate is then used to compute where f (R) is the response function above. Finally, to simulate the effects of longer range connections in layer 6, we take a weighted average of R 6 (x,t) over like orientation domains: the firing rate r 6 (x,t) of a layer 6 cell at location x and time t is computed as where L(x) is a set of locations including x as well as locations in neighboring HC's with the same orientation preference. The final computed rate of the layer 6 cell at location x and time t is r 6 (x,t).
The implementation of this time-dependent firing rate by a L6 neuron is realized by taking a baseline spike train and copying or canceling spikes with a probability consistent with the firing rate r 6 (x,t). A new value r 6 (x,t) is recomputed every 2.5 ms.

Parameters
The parameters in Table 1 (4), as discussed in section 0.3. B) The response function of layer 6 to local layer 4 activity, as introduced at the beginning of section 0.4. Here R 0 is the local firing rate of simple cells. We have elected to index layer 6 response to simple cell firing rate because they are more stable. Mean local firing rate is ≈ 1.6 times R 0 .
The following references are cited in the Appendix: